基于两相流EWF模型蒸发式空冷器数值模拟张庆;王学生;阮伟程;刘子龙;陈琴珠【摘要】建立了8排蒸发式空冷器管束三维模型,在恒壁温边界条件下,在管束表面施加质量源项模拟喷淋液膜,采用Eulerian Wall Film (EWF)模型与Mixture Species Transport模型耦合来研究蒸发式空冷液膜与空气间热质传递.数值模拟得到的空气出口温度、含湿量的数据与试验数据的误差分别为-0.67％～-0.98％、-4.95％～2.29％.比较了不同喷淋流量下管束表面液膜质量分布,小喷淋流量下液膜主要分布在管壁下半部分,随着喷淋水流量增加,管壁液膜分布趋于均匀,管排水膜温度由上至下先增加后减小.由于空气在管束背风面的流速较低,形成较高含湿量与温度三角区域.数值模拟得到的水膜与空气间的传质系数比试验值小,误差为8.00％～-9.30％.揭示了蒸发式空冷器热质传递机理,为蒸发式空冷器设计改造提供了理论基础.【期刊名称】《华东理工大学学报（自然科学版）》【年(卷),期】2018(044)006【总页数】7页(P928-934)【关键词】蒸发式空冷;EWF模型;热质传递;数值模拟【作者】张庆;王学生;阮伟程;刘子龙;陈琴珠【作者单位】华东理工大学承压系统与安全教育部重点实验室,上海200237;华东理工大学承压系统与安全教育部重点实验室,上海200237;华东理工大学承压系统与安全教育部重点实验室,上海200237;华东理工大学承压系统与安全教育部重点实验室,上海200237;华东理工大学承压系统与安全教育部重点实验室,上海200237【正文语种】中文【中图分类】TQ051.5蒸发式空冷器通过管束底部强制流动的空气和管束外部喷淋水与管内介质进行热交换，同时空气与水进行热质交换。
相较于空冷式热交换器仅利用掠过管束的空气带走管内热量[1]，蒸发式空冷器利用管外水膜的非饱和蒸发潜热，使空气侧传热传质显著增加。
目前，鲜有关于对蒸发式空冷器热质传递数值模拟的研究。
国内有研究者[2-3]进行了蒸发式空冷器单管的流体动力学模拟，但缺乏试验数据比照，结果仅作为参考。
文献[4-5]对蒸发式空冷器数值计算模型进行简化，特别是水膜被简化，因此没有涉及到空气含湿量的变化。
文献[6]研究了不同条件、不同液膜厚度时，采用数值模拟来计算冷态情况下水平管降膜蒸发器管外液膜流动的影响因素，但未涉及到管束热态传热及传质的研究。
而国外研究者中，Gan等 [7-8]先后采用管壁固定与管壁线性热流密度作为边界条件进行蒸发式空冷器的数值模拟，两者均对评估换热性能影响较大。
Facao等[9]将传质关联式作为数值模拟边界条件，结果存在较大误差。
Kachhwaha等[10]对蒸发式空冷器进行了试验和二维数值模拟研究，指出数值模拟方法可用于指导工业蒸发系统设计。
文献[11-12]只是对蒸发式冷却塔中的外部空气流动进行了数值模拟研究。
本文建立了三维湿式蒸发式空冷器管束物理模型，通过在管束壁面施加质量源项和恒壁温边界条件下，研究喷淋液膜成型及其与空气间热质传递过程，将模拟结果与试验数据进行比较，验证本文所用数值模拟模型及边界条件的合理性，本研究对蒸发式空冷器液膜空气热质传递机理研究有重要指导意义。
1 物理模型与边界条件1.1 试验装置蒸发式空冷器中空气自下而上与喷淋水逆流接触，带走管束内热量，在管束间完成热质传递，试验装置如图1(a)所示。
数值模拟中换热管为8排铝-碳钢复合管管束，换热管规格26.6 mm×2.5 mm，管长 2.24 m，换热管呈正三角形排布，管心距50 mm，奇数排13根管，偶数排12根管。
换热管示意图示于图1(b)。
图1 试验装置及管束参数Fig.1 Experimental apparatus and tube parameters 1.2 物理模型与边界条件考虑物理模型对称性，Design Modeler中建立三维模型如图2所示。
使用ANSYS 17.0 workbench meshing 对物理模型进行六面体网格划分及边界设定。
在操作软件Fluent中设定换热管外壁为定壁温条件，在计算区域中，除管壁及进出口以外边界设定为对称边界，进口为空气质量进口，出口为压力出口，模拟中通过对管壁施加质量源项的方式来模拟管束上层喷淋水喷淋在管壁上形成的液膜，用Eulerian Wall Film (EWF)模型来模拟液膜热质传递。
本文中选取湍流模型Standard k-ε模型与非平衡壁面函数(Non-Equilibrium Wall Functions)来模拟光管蒸发式空冷器热力性能，采用SIMPLE算法求解压力-速度耦合、二阶迎风离散，水膜对空气传质过程，采用不含化学反应的组分扩散(Species Transport without Reaction)模型进行模拟，通过能量方程模拟热量传递。
图2 物理模型Fig.2 Physical model2 控制方程及理论基础在液膜与空气的热质传递过程中，空气作为连续相，包含质量、动量、能量与浓度等各种平衡方程。
作为稳态不可压缩流体处理时，其平衡方程如下：(1)式中：为梯度算子，φ为流体变量，具体包括空气平均速度、平均焓值及平均水蒸汽浓度；V为流体矢量速度，m/s；ρ为空气密度，kg/m3；Γφ为扩散系数，N·s/m2；Sφ为连续相源项；为空气与液膜间质量传递附加源项。
薄液膜假设通常用在欧拉液膜模型中，特别是液膜的厚度比壁面曲率半径小很多、沿液膜厚度方向的物性没有差异，此时液膜可以看作是与壁面平行的。
在三维场中对于液膜的质量平衡为(2)其中为表面梯度算子；ρl是液体密度，kg/m3；h是液膜高度，m；Vl是平均液膜速度，是每一单元液滴收集、液膜分离、液膜剥落和相变带来的质量源项。
液膜的动量平衡：(3)其中：PL=Pgas+Ph+Pσ(4)Ph=-ρh(n·g)(5)(6)其中，式(3)等式左边分别代表液膜瞬态变化及对流对动量影响，公式等式右侧第1项包括了气流压力、垂直壁面重力和表面张力影响动量变化；第2项代表平行壁面重力影响；第3项是气液接触处黏性剪切力影响，第4项是液膜黏性力影响，最后一项代表液膜积聚或分离带来的动量改变。
三维场中液膜的能量平衡方程如下：(7)式中：Ts为液膜空气界面处温度，℃；Tf为平均液膜温度，℃；Tw为壁面温度，℃；为壁面液膜源项；为蒸发或冷凝质量；L(Ts)为和相变有关的潜热。
本文采用EWF模型与Mixture Species Transport模型耦合来研究液膜与空气处相变和汽相浓度变化。
相变控制方程如下：(8)式中：D为汽相质量扩散系数；δ为控制单元中心距离壁面距离；Cphase为相变常数；yi为控制单元中心蒸汽相比率；饱和相比率ysat计算如下：(9)式中：P为汽相绝对压力；Mi和分别为汽相及混合相摩尔质量；饱和压力Psat仅为温度的函数；水的饱和汽相压力计算如下[13]：(10)式(8)中相变常数Cphase相对于冷凝和蒸发分别计算如下：(11)式中: Ccon为冷凝常数；Cvap为蒸发常数；h为壁面液膜高度，以保证在仅有液膜的情况下才会产生蒸发。
式(8)中汽相质量扩散系数D计算如下:(12)式中：T为热力学温度，K；p为总压强，Pa；μA、μB为气体A、B的分子量；VA、VB为气体A、B在正常沸点时液态克摩尔容积，cm3/gmol。
以管束壁温作为边界条件，管束壁温由实验数据计算得出，壁温计算如下：(13)式中：Tm为管内介质定性温度，℃；qo为基管的热流密度，W/m2；hi为管内介质换热系数，W/(m2·℃)；do为管子外径，m；di为管子内径，m；dm为管子平均直径，m；b为管壁厚度，m；λw为管材导热系数，W/(m·℃)；ri为管内污垢热阻，(m2·℃)/W。
管内传热系数计算方法采用Dittus-Boelter公式[14]：Nu=0.023Re0.8Pr0.3(14)式中：Nu为努赛尔数；Re为雷诺数；Pr为普朗特数。
水膜到空气的传质系数(hd)见式(5)，同时对传质系数的试验值与模拟值进行了对比。
(15)式中：hd为传质系数，kg/(m2·s)；ma为空气质量流量，kg/s；A为管束外表面积，m2；imasw为水膜温度对应的饱和空气焓值，kJ/kg；imai，imao分别为进出口空气焓值，kJ/kg。
3 网格无关性检验根据试验条件，本文控制风机频率为38～50 Hz，对应空气质量流量为 3.04～3.90 kg/(m2·s)，喷淋水流量分别为3、5、7、9、11 m3/h，每个喷淋水流量下进行3组不同进口风速试验。
选取喷淋水流量11 m3/h，风机频率50 Hz实验工况下进行网格无关性检验，将网格数划分为 160 225、290 810、341 458、432 032、590 445 和 793 152 六组网格。
对于网格数为 793 152 的模型，约计算 7 500 步后，液膜质量保持稳定，此时计算收敛，在6组网格计算中保持迭代步足够大以使计算收敛，当网格数量达到341 458 时，空气出口温度、水蒸气浓度、液膜质量计算结果基本稳定，综合考虑计算成本与计算精度将物理模型划分网格数 341 458。
4 数值模拟结果分析讨论4.1 出口温度与含湿量对比保持喷淋水流量分别为3、5、7、9、11 m3/h，每个喷淋流量下保持风机频率为38～50 Hz内的3个档位，表1给出了15个工况下，空气出口温度与含湿量的模拟值与试验值的对比，见表1。
由表1可知，数值模拟得到的空气出口温度与出口含湿量与试验数据吻合较好，出口温度误差为 -0.67%～-0.98%，出口含湿量误差为 -4.95%～2.29%。
验证了采用EWF模型与Mixture Species Transport模型耦合来研究蒸发式空冷器中液膜与空气间热质传递的可行性与可靠性。
表1 出口温度及含湿量对比Table 1 Comparison of outlet air temperature and humidity ratioCaseConditionsV a/(m·s-1)Ta/KRHin(a)/%Air temperatureTout(ex)/KTout(f)/KError/%Humidity ratiodout(ex)/(g·kg-1)dout(f)/(g·kg-1)Error/%12.803296.9088.8309.59306.93-0.8636.3235.27-2.8922.756297.3088.3309.77307.70-0.6736.6236.920.8232.683296.8088.6310.25307.90-0.7638.1738.290.3142.614296.7088.7310.39308.17-0.7238.5239.402.2852.593298.1084.4311.55309.13-0.7841.5441.04-1.2063.043297.0088.7309.51306.64-0.9336.4034.87-4.2073.043297.1089.4309.85306.96-0.9336.6635.73-2.5482.987298.3083.7310.73307.75-0.9638.7736.85-4.9592.986298.1584.2310.93308.08-0.9238.9038.30-1.54102.830297.9086.4311.02308.32-0.8740.1839.48-1.74113.308296.5589.6307.96305.49-0.8032.6632.24-1.29123.291296.5589.5309.03305.99-0.9835.0434.12-2.63133.213297.1087.8309.51306.70-0.9135.7935.45-0.96143.144297.1087.2309.67306.96-0.8836.0635.73-0.92153.097297.4087.3309.79306.90-0.9335.7736.592.29a—Air；V—Velocity；T—Temperature；RH—Relative humidity；in—Inlet；out—Outlet；d—Moisture content of air；ex—Experimental value；f—Simulated value of fluent4.2 管束流场及液膜分布图3所示为风机频率50 Hz，喷淋水流量为3 m3/h工况下管束间速度与压力云图，可知在空气流经管束表面时，会在管后部形成大尺度的涡，当空气扰动剧烈时，管束后部环流区的大尺寸涡的尺度逐渐变小。