辽宁科技大学材料与冶金学院 刘 邓
中国科学院大学工程热物理所 武广龙
【摘 要】利用有限容积法，针对纯制冷剂R134a和非共沸混合制冷剂R407c，在水平管和U型管组合的物理模型内沸腾换热过程进行数值模拟，得到了温度、速度、干度及换热系数的分布，并分析了管内工质的流速和热流密度的变化对沸腾换热系数的影响。此外，还对R134a在U型回转式弯头的水平光滑管内的沸腾换热进行非稳态的数值模拟，可以清晰地观察到气泡的产生、跃离、生长和破碎的气液两相流型变化，并且模拟结果与Baker流型图也较吻合。
【关键词】气液两相流 沸腾传热 数值模拟
Abstract:Adopting the method of Finite Volume Method to simulate the saturated boiling heat transfer in a smooth horizontal tube of R134a and R407c, and to get the distribution of temperature, velocity, dryness and heat transfer coefficient of R134a and R407c in the process of boiling heat transfer in tubes under the condition of different velocity and different heat flux density. In addition, the boiling heat transfer in tubes of R134a under the condition of unsteady three-dimensional numerical simulation is conducted separately, the changing process of the gas-liquid two phase flow pattern in tubes is observed, and the simulations give good agreement with the flow regimes expected from the Baker chart.
Key words:Gas liquid two-phase flow, Boiling heat transfer, Numerical simulation
1 引言
在80年代，科学界确认氟氯烃（CFC）是引起臭氧破坏和温室效应的危害物质之一。淘汰了R11、R12、R113、R114等普遍使用的CFC类制冷剂，即使R22等HCFC类对环境的污染破坏影响相对较小一些，最终也将被禁止使用。所以，研究并开发环保的新型制冷剂在制冷领域内是一个重要的研究课题。R134a是第一个被提出来的非臭氧破坏物质（ODP=0），具有较好的热物理性质，其温室效应潜能值（GWP）很低，不易燃，无毒，在冰箱、冷柜、汽车空调等制冷设备中已经开始使用。R407c具有清洁性、低毒性和不燃性等特点，尤其不破坏臭氧层（ODP=0），是目前使用较多的中高温环保制冷剂。所以研究制冷剂R134a和R407c的换热特性对今后的制冷领域有着重要作用。
随着在社会生产力提高的同时，也开始引起人们对研究气液两相流的注意。郭雷[1]等着重研究在竖直矩形的细通道内水的沸腾换热情况；欧阳新萍[2]等在水平管外R404a和R407c的沸腾换热进行了实验研究；李静[3]模拟研究的是在竖直圆管内的沸腾换热；吴晓敏[4]等对CO2/丙烷混合工质在水平管内的沸腾换热进行分析。而本文选用的是将水平管和U型管两者组合起来的物理模型，针对纯制冷剂R134a和非共沸混合制冷剂R407c在带有翅片的光滑管内沸腾换热过程进行数值模拟，并进行对比分析。此外，还对在U型回转式弯头的水平光滑管内的沸腾换热过程单独做非稳态的数值模拟，以此观察气液两相流的流型变化，并与Baker流型图做了对比。
2 数值模拟的模型
2.1 物理模型及参数
图2.1为模拟用的水平弯管的示意图。其中管直径D为25 mm，壁厚为0.5 mm，管长为12.89 m。选取R134a和R407c两种工质。由于本文重点研究的是制冷剂在管内的流动换热，所以创建几何模型时，忽略管壁，只是针对管内区域划分网格，然后进行了相关网格的独立性检验，最后得到了符合条件的计算网格。本文采用固定热流密度（30 kW/m2）的方法，对管壁使用壁面边界条件，并且为0.5 m/s的入口流速，自由出流的出口流速。
图2.1 水平弯管示意图
非稳态模拟采用的边界条件与稳态模拟相同。此外，利用软件Refprop8.0计算分别得到了R134a和R407c的热物性参数，分别如表2.1和表2.2所示。
表2.1 R134a的热物性参数
表2.2 R407c的热物性参数
2.2 数学模型
控制方程可以表示成如下的通用形式[5]：
式中ф为通用变量；гф为广义扩散系数；Sф为广义源项。
制冷剂液体在管内蒸发时，流体被加热达到所在压力下的饱和温度以后，会从液相变为气相，所以两相间必然会有质量传递。质量传递过程采用De Schepper等[6]提出来的方程，不同相之间的质量源项如下：
当Tliq≥Tsat时，制冷剂液体蒸发
液相： 
气相：
气液两相间的能量传递是伴随着质量传递过程进行的，蒸发过程中能量源项的计算公式为：
其中，Tsat为蒸发温度；Tliq为制冷剂液相的温度；αliq为液相容积率；αvap为气相容积率；ρliq为制冷剂液相密度；
为单位时间单位体积内向液相传递的质量，
为单位时间单位体积内向气相传递的质量；ΔH为汽化潜热。以上模型公式针对的是纯制冷剂工质，而对于混合制冷剂工质，拟采用Lee模型[7]定义不同相之间的质量传递：
液相：
气相：
制冷剂在蒸发过程中，气液两相间的能量传递与纯制冷剂工质的能量源项的计算公式类似：
其中，Tb为泡点温度；Td为露点温度；τl为液相的时间松弛系数，τv为气相的时间松弛系数，建议取值为0.1[7]。
3 数值模拟结果及分析
3.1 R134a非稳态模拟结果分析
管内沸腾换热的两相体积分数分布如图3.1所示，其中液相的体积分数用红色表示，气相的体积分数用蓝色表示。气弹状流、细泡状流、波状流和分层流的两相流特定流型能在图中很明显的展现出来。分析：液体在进入刚管内时被加热，所以换热管开始部分并没有气泡，只能观察到液相。首先被加热到饱和温度的液体开始蒸发，在这些部分首先开始形成气泡，受重力影响会逐渐上升到管子顶部，随着细泡状流的形成，一系列其他两相流流型也开始逐渐发展。液相的速度小于气相的速度，当气相的数量逐渐增加时，管内的两相流动也变得愈发混乱。
图3.1 管内沸腾换热的两相体积分数分布
当液体继续向管内流入时，气相的体积分数和速度都会增加，也会观察到气弹状流型。气弹状流区域的混乱程度比细泡状流和气塞状流的混乱程度大，但其两相交界面曲率与后两者相比较会变缓和。第二管道后部，第三和四管道分别出现了出现气弹状、波状流和分层流流型。许多采用水平管件的换热器的设计必须用180°弯头（回转型弯头）。由于回转型弯头的出现，气相在管道弯头处会重新分布，如图3.2（a）所示。在一些绝热弯头处，还可能出现气相的再冷凝现象。两相的温度场如图3.2（b）所示，液体在管内加热，当达到沸点283.15 K时，液体开始蒸发，气相逐渐形成。就如图3.2（b）所示，累积的气相由于过热其温度有所上升，然而由于液体的替换流动和气相在管顶部，所以只是壁温有所上升。气液两相的速度分布具有明显的差别如图3.2（c）所示。而且，气相速度往往大于液相速度，速度的差别导致了气弹状流和波状流的形成和发展。
图3.2 气液两相在管内的沸腾换热过程
在上述对非稳态模拟结果分析的基础之上，本文利用相关的相应的流型图与模拟得到的结果进行比对。图3.3为在管道中心界面上选取的流型不同的6个点，与目前使用较多的Baker流型图较为吻合，如图3.4所示。
图3.3 选取的6个不同流型点
图3.4 Baker流型图[8]
3.2 R134a稳态模拟结果分析
由于换热管较长，故选取具有特征的截面。如图3.5所示可以看出，当换热管管壁为恒定的热流密度时，R134a从换热管入口流动到出口进行沸腾换热，管壁温度高于工质的温度，工质被加热，其温度逐渐升高，近壁面处的工质温度高于工质中心部分的温度，故温度与温度梯度沿管热管管壁面的法线方向逐渐逐渐减小。
图3.5 R134a在截面上的的温度场分布
如图3.6所示，R134a从换热管入口流动到出口进行沸腾换热时，在入口速度一定的情况下，由于受管壁加热的作用，工质由粘性系数相对较高的液相制冷剂逐渐沸腾转化为粘性系数相对较低的气相制冷剂，工质的沿程阻力也相应地逐渐减小，而且在同一压力梯度下，蒸汽的密度较小，速度较大。在沿壁面法线方向，由于壁面处是首先发生沸腾换热的地方，传质阻力影响比较明显，阻力沿法线方向起初较大，但逐渐减小。在边界层内出现很明显的速度梯度变化现象，即中心部分流速大，壁面附近流速相对较小，故出现沿管长方向R134a速度逐渐增大的现象。
图3.6 R134a在截面上的速度场分布
图3.7为R134a在截面上的含气率分布。由于管道的几何的对称性，气体的体积分数的分布因此也具有很好的对称性。R134a从换热管入口流动到出口时，在加热作用下，壁温升高迅速，当获得能够沸腾的过热度时，管内开始产生气泡，同时气泡的数量越来越多，故气相体积分数也逐渐增大。在刚刚开始沸腾时，同一截面上，由于浮力和拽曳力等外力，壁面上产生的气泡从壁面脱跃离，逐渐进入到主流中，故贴近管壁面处含气率比管中心部分较高。随着沸腾换热过程的停止，大部分液相都转化为气相。此外，从图中可以看出，在弯头部位，其截面的体积含气率往往比其前后的截面含气率较低，出现气相再冷凝的现象。
图3.7 R134a在截面上的含气率分布
本文截取典型区域的图像来展示管壁面上的换热系数分布情况如图3.8所示。由于工质在进口时已经处于饱和状态，所以直接就开始饱和沸腾，没有了过冷沸腾阶段。饱和沸腾换热时，由于管道上壁面的液相体积率远小于气相，随着换热系数的降低，换热过程逐渐减弱。
图3.8 R134a的局部换热系数分布
3.3 R407c模拟结果分析
如图3.9所示，当换热管管壁为恒定的热流密度时，R407c从换热管入口流动到出口沸腾换热时，管壁温度高于工质的温度，工质被加热，其温度逐渐升高，近壁面处的工质温度高于工质中心部分的温度，故换热管中心区的温度分布较为均匀，出现温度与温度梯度沿管热管管壁面的法线方向逐渐减小的现象。
图3.9 R407c在截面上的温度分布
如图3.10所示，R407c从换热管入口流动到出口进行沸腾换热时，在入口速度一定的情况下，由于管壁的加热作用，工质受热开始发生沸腾，由液相向气相开始转变，气相体积分数逐渐增大；在管长方向上，速度也随之逐渐增大。在沿管径方向，管壁受加热作用温度迅速上升，壁面处的液体工质首先出现气泡，自然其处的传质阻力的影响比较明显，但是沿管径方向逐渐减小，导致了中心部分流速大，壁面附近流速相对较小，故出现沿管长方向R407c速度逐渐增大的现象。
图3.10 R407c在截面上的速度分布
图3.11为R407c的截面含气率分布。由于管道的几何的对称性，除弯头部分外，气体的体积分数的分布也具有很好的对称性。R407c从换热管入口流动到出口时，在加热作用下，壁温升高迅速，当获得能够沸腾的过热度时，管内开始产生气泡，同时气泡的数量越来越多，故气相体积分数也逐渐增大。在沸腾开始阶段，同一截面上，气泡在壁面上产生，在各种力的作用下，从壁面跃离，变大，逐渐进入到主流中，故贴近管壁面处含气率比管中心部分较高。随着沸腾换热过程的停止，大部分液相都转化为气相。此外，与图3.7相比较，可以看出，由于流速增加，液体的替代作用增强，相似横截面含气率分布的截面出现的位置向后延迟出现。在弯头部分，两相会重新分布并且会有一定的再冷凝现象。
图3.11 R407c在截面上的含气率分布
图3.12是模拟得到的R407c的局部换热系数沿管长方向上的分布和变化情况。由于换热管细长，所以本文只是截取典型区域的图像来展示管壁面上的换热系数分布情况。由于工质在进口时已经处于饱和状态，所以直接就开始饱和沸腾，没有了过冷沸腾阶段。饱和沸腾换热时，由于管道上壁面的液相体积率远小于气相，随着换热系数的降低，换热过程逐渐减弱。
图3.12 R407c的局部换热系数分布
3.4 R134a和R407c沸腾换热特性的对比分析
3.4.1 R134a和R407c在管内沸腾换热的温度分布
图3.13中为制冷剂R134a和R407c温度随干度的变化情况。从（a）可以看出，在发生沸腾后流体的温度并不是保持不变，而是有一定的升高，大约有2 K，这是由于此时传热的的热量大部分用于补充制冷剂发生相变时的所需要的潜热。虽然理论上讲，流体发生相变时温度保持不变，但是此结果与文献[8]的结果有较好的吻合。从图（b）可以看出，随着干度的增加，流体温度也逐渐增加，由于混合制冷剂R407c的非共沸的特性使其与纯制冷剂相比具有非等温特性，所以其温度增加的幅度要远大于同条件下的R134a的温升。
图3.13 制冷剂的温度随干度的变化情况
3.4.2 R134a和R407c的沸腾换热系数对比
如图3.14所示，R134a和R407c随着干度的增加沸腾换热系数先增加后减小。
图3.14 制冷剂的沸腾换热系数随干度的变化情况
对比分析图3.15（a）、（b）得，R407c的平均沸腾换热系数大于R134a的平均沸腾换热系数。
图3.15 制冷剂在不同流速下的平均沸腾换热系数
3.4.3 不同的流速与沸腾换热系数的关系
随着制冷剂流速的增加，由于气泡沿壁面滑移速度逐渐变快，使得气泡从壁面跃离，故平均沸腾换热系数随着制冷剂流速的增加而增加，如图3.16所示。
图3.16 制冷剂在不同流速下的沸腾换热系数
3.4.4 不同的热流密度与沸腾换热系数的关系
在泡状沸腾时，热流密度对换热的影响要大一些，使得沸腾换热系数增大，故出现沸腾换热系数会随着热流密度的增加而增加的现象，如图3.17所示。
图3.17 制冷剂在不同热流密度下的沸腾换热系数
4 结论
本文对R134a和R407c在水平光滑管内的沸腾换热过程进行数值模拟。结果表明：
（1）由于水平光滑管内除弯头部分的横截面的对称性，使得R134a的速度、干度及温度分布也具有对称性。沿流动方向，工质的温度逐渐升高，从非稳态模拟结果中可以看出气泡数量逐渐增加，并且不断跃离壁面、生长、破碎进入主流区。壁面边界层内，由于粘性，速度由0沿管径方向逐渐增加到主流速度，而速度梯度则与之相反。在管弯头部分，由非稳态模拟结果可以看出，气相和液相会重新分布，而且还可能会伴有气相的再冷凝现象。
（2）另外，混合制冷剂R407c的非共沸的特性使其与纯制冷剂相比具有非等温特性。
（3）当制冷剂流速或热流密度增大时，换热管内的沸腾换热系数也会增加。
（4）将R134a在非稳态条件下的模拟结果与Baker流型图进行了验证，两者结果比较吻合，模拟计算满足了一定的合理性，为进一步模型优化提供了理论依据。
参考文献：
[1] 郭雷,张树生,程林. R134a竖直矩形细通道内谁沸腾换热的数值模拟[J].热能动力工程, 2011,(1):31-35.
[2] 欧阳新萍,刘超,林梦. R404A和R407C在水平强化管外的凝结换热实验研究[J].制冷学报, 2015,(4):72-77.
[3] 李静.三维竖直圆管内流动沸腾的数值模拟[J].广东化工, 2011, 38(3): 75-76.
[4] 吴晓敏,魏兆福,莫少嘉,等.CO2/丙烷混合工质管内流动沸腾换热的数值模拟[J].工程热物理学报, 2011, (6):1019-1022.
[5] 陶文栓.数值传热学[M].第2版.西安:西安交通大学出版社,2001:5.
[6] De Schepper S C K, Heynderickx G J, Marin G B. Modeling the evaporation of a hydrocarbon feedstock in the convection section of a steam cracker[J]. Computers & Chemical Engineering, 2009, 33(1): 122~132.
[7] Lee W H. Pressure iteration scheme for two-phase flow modeling[J]. Multiphase Transport: Fundamentals, Reactor Safety, Applications, 1980: 407-432.
[8] De Schepper S C K, Heynderickx G J, Marin G B. CFD modeling of all gas–liquid and vapor–liquid flow regimes predicted by the Baker chart[J]. Chemical Engineering Journal, 2008, 138(1): 349-357.
备注：本文收录于第21届暖通空调制冷学术年会论文集。版权归论文作者所有，任何形式转载请联系作者。
