绕管式换热器壳侧降膜流动和相变传热的数值模拟

李剑锐，陈杰，浦晖，李恩道，丁国良，庄大伟，胡海涛

（1上海交通大学制冷与低温工程研究所，上海200240；2中海石油气电集团技术研发中心，北京100028）

引 言

绕管式换热器作为天然气液化过程主要使用的低温换热器，受到大型陆上天然气液化工厂和海上浮式天然气液化平台 （FLNG）的广泛应用［1－2］。绕管式换热器具有结构紧凑、适用地域广［3－4］、工况范围大［5］的特点，除去运输困难的缺点，可以组合为超大型换热机组［6］。作为主要的低温换热器，换热器成本占FLNG总投资的20%以上，同时绕管式换热器占到总热量损失的20%［2］，其换热性能及运行效率对天然气开发平台的实际产液率有极大的影响［7］。而绕管式换热器壳侧的传热和压降损失都较大［8－9］，因此针对绕管式换热器壳侧的设计优化和性能改善是十分必要的。

在绕管式换热器的壳侧，两相制冷剂由顶部经过分配器均匀流下，落在管外形成液膜，与管内相对高温的天然气进行换热，发生传质，造成制冷剂侧发生干度变化［7］。而随着干度沿流程上升，气液相间相互作用发生变化，壳侧流体在不同的流型间发生转变。绕管式换热器在不同的流型阶段的换热性能、压力损失都不相同［9］，因此为了能够预测绕管式换热器壳侧性能，开发反映绕管式换热器壳侧流动、传热和传质的数值模型十分必要。

但鉴于绕管式换热器生产成本十分高昂，传统的整机实验优化方法较难实现，因此绕管式换热器的实验研究多采用小型模型。实验研究也多针对管侧流动和换热，而绕管式换热器壳侧相关研究和相关实验研究较少。Fredheim［7］建立了以丙烷为工质的绕管式换热器壳侧实验方法，建立了传热和压降模型，通过单相流和两相流实验分别得到特定流型下的传热和压降关联式；Aunan在Fredheim的基础上使用液氮、甲烷、乙烷和甲烷乙烷混合工质为壳侧工质进行实验，对比了单相流、剪切流、降膜流动下的传热和压降特性，讨论了已有的混合物工质两元传热模型的局限性［9］，并在实验基础上验证了已有模型的精度［10］，同时对原有降膜流动的传热模型进行改进，适用于5000＜Re＜8000范围内的几种特定烷烃类工质［11］。现有的实验研究局限于特定烷烃类工质的传热压降研究，研究过程耗费巨大，并且难以进行换热器的工质的替代、结构上的性能优化和海上恶劣工况的性能预测的研究。如果采用数值模拟的方法对绕管式换热器进行性能模拟，则能够大幅节约设计开发成本。

由于绕管式换热器内部壳侧包含气液两相，为管外降膜流动，而针对壳侧的数值模拟研究也较少，仅Lu等［12］对卧式绕管式换热器进行单相空气流模拟，难以满足壳侧多相降膜流动的要求，因此需要运用针对降膜流动的数值模型。目前已有关于管外降膜流动的数值模拟研究如下：针对模拟模型的选择，Killion等［13］对比了管外的溴化锂液滴流动过程的模拟结果与实验结果，证明模拟管外液膜流动三维模型比二维轴对称模型结果更精确；针对液膜流动的偏移，雷贤良等［14］也采用VOF方法建立降膜流动模型，分析横向气流对液膜流动偏移的影响，并建立了一个量纲为1的方程计算液膜的偏移量；针对管外液膜厚度的分布，王小飞等［15］采用VOF模型研究了不同结构的降膜蒸发器对液膜厚度的影响；邱庆刚等［16］通过降膜流动模拟分析了喷洒密度、管的尺寸对管壁上薄膜厚度的影响因素，并与相关实验数据进行了比较，结果表明液膜同时在管的顶部和底部区域最厚，同时最薄出现在圆形角约120°处；针对不同管型的降膜流动，罗林聪等用二维数值模拟方法对几种管型 （包括光管和滴形管、椭圆形管）进行研究，结果表明圆管、滴形管和椭圆形管的液膜厚度依次减小，最小值出现在周向125°、170°、160°处，而异形管的热边界层更薄，传热性能更优［17－18］，随后通过实验验证得到圆管和椭圆管液膜最薄处角度分别为105°和135°［19］；针对降膜流动流型的转变，Pise等［20］通过数值模拟研究了Reynolds数、Galileo数、几何参数对流型转变过程的影响，得到了低流速下流型转变判别方程；Chen等［21］以流量为判定流型转变的依据，通过模拟得到了特定结构下层状流动、柱状流动、滴状流动的临界流量。尽管降膜流动模拟研究已有很多，但是大都以降膜蒸发器为背景，集中在研究液膜的流动和液膜分布方面，流型研究也在低Re的范围内，几乎没有关于多相流传质的模拟研究，甚至少有不含传质过程的降膜传热研究。

对于绕管式换热器的传质过程而言，由于绕管式换热器管侧和壳侧温差不超过10K［4－5］，管侧与壁面间的温差更小，因此不会出现过度沸腾和膜态沸腾的情况。目前大多数蒸发、沸腾的模拟研究集中在两方面：一是池沸腾过程，包括平面核态沸腾的气泡形成、生长、脱离过程［22－24］，并加以实验验证［25－26］，着重于气泡脱离直径和局部热流变化，也有少数圆管外膜态沸腾的研究［27］；二是管内流动沸腾， 包 括 竖 直 管［28］、 水 平 管［29－30］管 内 相 变 引起的流型转变，重力作用下蛇形管［31］、螺旋管［32］管内流动沸腾的气液分布不均问题，微通道内部的气泡成核研究［33］。而关于管外流动沸腾和降膜蒸发的数值模拟则没有相应的模型可以使用。因此，针对壳侧的降膜蒸发过程需要适用于壳侧的传热传质模型。

1 数值模型

1.1 模型对象描述

绕管式换热器内部结构如图1所示。管侧下进上出，沿途发生天然气液化过程；管外壳侧上进下出，沿途为制冷剂蒸发过程［7，9］。壳侧的制冷剂在流动过程中分为如下几步：①入口处的制冷剂处于低干度饱和状态，在压力作用下喷入绕管式换热器上方的缓冲区，然后经过一个分配器，使制冷剂液体均匀漏下，覆盖整个排布面；②由分配器垂直流下的制冷剂液体附着在管外壁面，沿管间空隙进行降膜流动，此时制冷剂与管内温度相对较高的天然气发生换热，传入的热量使饱和制冷剂发生汽化相变，制冷剂干度逐渐上升；③在干度达到一定程度后，气相流速加快，而液体的降膜流动在气体的推动下进行，液膜内流动速度加快，形成剪切降膜流；④在干度更高时，制冷剂的液体逐渐随汽化减少，无法在管间形成完整液膜，进而变为柱状降膜流以及滴状降膜流，此时可能出现液膜无法覆盖的干区，换热性能下降；⑤制冷剂完全汽化，气态的制冷剂从下方出口流出。

图1 绕管式换热器壳侧流型示意图Fig.1 Falling film flow patterns at shell side of CWHE

因此，要描述绕管式换热器壳侧制冷剂降膜流动，必须描述不同干度下的制冷剂降膜流动过程、制冷剂汽化相变过程及制冷剂汽液两相相互作用。

制冷剂壳侧的降膜流动是在表面张力、制冷剂黏性力和重力共同作用下的结果［9］。在制冷剂黏性力的作用下，液相趋于附着在管壁表面流动；在表面张力的作用下，液相趋于平摊，覆盖整个管外壁面，形成液膜，以及在不同干度下表现出不同的流型；在重力的作用下，液相沿管壁间垂直向下流动，形成管间液膜。

制冷剂壳侧的汽化相变过程属于流动沸腾过程［8］。流入的制冷剂处于两相饱和状态，与管壁接触的液相和气相部分发生不同的换热过程。液相获得热量的过程为潜热传热，使得液相发生汽化相变，在管外壁面上产生汽化核心，进而形成气泡；气相获得热量的过程为显热传热，使得气相温度上升，而气相温度高于周围液相温度，造成气液相间发生换热，间接使液相汽化，表现为气泡增大。随着气泡增大至脱离直径，气泡会随液膜流动离开，同时聚集变大，突破液膜进入气相。液相和气相的两种不同换热过程最后都导致制冷剂液相汽化，但换热性能上却差异很大，这也是导致换热器内部不同干度区域换热性能下降的原因。

制冷剂壳侧的气液相间存在的相互作用如图2所示，分为如下几种：①气泡形成，并在液相流动作用下的脱离过程；②气泡聚集变大，突破液膜进入气相的过程；③气相膨胀流速增加，对液相的剪切和推动过程。针对气泡和液膜的受力分析如图3所示，除去气液相都存在的表面张力、黏性力和重力外，需要额外考虑液相沿流动方向对气泡的推动力、气泡受到的垂直向上的浮力以及气流对液膜沿气液交界面的剪切力对降膜流动过程的影响。

图2 壳侧降膜流动气液相互作用Fig.2 Interaction between vapor phase and liquid phase at shell side

对于制冷剂降膜流动的过程，最重要的是能够反映出气液两相的不同流动过程以及各个力的实际作用；对于气泡形成和变大的过程，最重要的是计算出汽化传质速率和相变造成的潜热传热的影响；对于气液相间的相互作用，最重要的是对气液的表面进行正确的追踪。

绕管式换热器降膜蒸发模拟的技术路线如图4所示。需要在模型的连续性方程中加入气液两相质量传递源项，在动量方程中加入表面张力源项和剪切力源项，在能量方程中加入潜热传热源项。

1.2 控制方程

降膜蒸发的传质和传热过程中，根据控制单元气液含量可以划分为如下3种控制单元：控制单元1，全为气相的情况；控制单元2，全为液相的情况；控制单元3，为气液相均存在的情况。如图5所示。

根据N－S方程，对于图5所示的3种类型的控制单元，针对降膜蒸发模型可建立如下的基本控制方程［30］。

连续性方程为：

气相

液相

其中气相和液相分别使用VOF模型进行连续性方程的求解，αv和αl分别代表控制单元内部气相和液相的体积分数；右侧为传质质量源项Sm，表示传质过程中气态和液态间的质量变化。

动量方程为：

图4 技术路线Fig.4 Technology roadmap

图5 不同干度下的控制单元Fig.5 Elements of different vapor quality

可以证明，在相界面另一侧表面上的压降取决于表面张力系数、曲率和表面的两个正交方向，其计算公式如下

式中，＾n为相界面函数，θ为接触角。

能量方程为

由于考虑了传热和传质的情形，能量方程存在显热源项Δ·（kΔT）和潜热源项Q。

1.3 传质子模型

关于传质源项Sm，分为蒸发和冷凝两部分进行考虑。制冷剂为共沸制冷剂的情况下，其中一项的值应为0；而非共沸制冷剂的情况下，则需要对两个源项分别进行考虑，以满足非共沸传质过程。

而传质过程仅发生在两个区域：①有温差的相界面；②管外壁面的覆盖液相且达到蒸发温度的区域。

1.3.1 相界面传质模型 在工质有固定蒸发温度的情况下使用Lee的气液传质模型［29］，该模型通过温度判断蒸发和冷凝两种传质过程，如图6所示。传质质量计算如下。

若局部温度高于蒸发温度，发生蒸发 （Tl＞Tsat）

图6 相界面传质模型Fig.6 Model of mass transfer at phase interface

若局部温度低于蒸发温度，发生冷凝（Tv＜Tsat）

1.3.2 管壁附近传质模型 若管壁面流体为液相，而且局部温度高于蒸发温度，液相获得的热量全部用于潜热传热，发生汽化传质，如图7所示。根据壁面的当地的热通量，蒸发传质质量计算如下 （αl＞0，Tl＞Tsat）

图7 管壁面传质模型Fig.7 Model of mass transfer near wall surface

由于管内温度始终比蒸发温度高，壁面附近不发生冷凝现象，没有冷凝传质质量

1.3.3 潜热传热模型 潜热传热源项Q根据质量传质源项Sm进行计算

式中，hfg为制冷剂的汽化潜热。

基于对以上各个模型进行综合建立降膜流动的综合模型，从而实现对降膜流动传热传质过程的完整数学描述。

2 求解方法

本研究数值模拟使用的模型如图8所示，基于商用软件ANSYS中的FLUENT进行流动模拟。采用VOF模型作为模拟两相流模型，连续表面张力模型 （CSF）作为模拟表面张力模型，从而实现对降膜流动液膜分布的模拟；同时通过FLUENT的用户自定义方程 （UDF）分别建立针对壁面及非壁面的传质模型，以实现通过对网格类型的判断采用不同方式计算传质质量，从而实现气泡生成、长大的模拟；采用VOF－CSF模型作为管壁表面的接触角模型；在高干度情况下气相流速较大，此时需采用标准k－ω作为湍流模型。

图8 几何模型Fig.8 Geometry models

模拟中使用ICEM网格划分软件划分的六面体结构化网格，并且对于壁面边界层和气液交界面进行了网格加密。针对网格独立性进行了验证，在最大网格尺寸不超过0.05mm且壁面附近的网格单元尺寸不超过0.0025mm时传热系数及传质质量计算误差低于2%，因此采用0.05mm作为基准网格尺寸及0.0025mm作为壁面附近的网格尺寸，以获得更高精度的模拟结果。同时网格壁面附近的最大网格尺寸需要小于气泡的脱离直径，脱离直径根据模拟采用的热通量和流量等工况采用管道内流动沸腾估算方法进行估算［28，31］，假定值最小值为0.1mm。而最大网格尺寸为0.05mm已经满足要求。

如图7所示，边界条件设定如下：顶部入口为速度进口边界，底部出口为外流边界，管壁面为不可滑移壁面，定热通量边界，左右两侧设定为周期性边界，与管道接触的壁面设为对称边界。

由于使用瞬态模型进行模拟，无法直接通过残差判断收敛。而计算足够时间后，计算区域内流型基本呈现稳定，传热和传质会逐渐趋于定值，因此判断计算收敛条件为在100个时间步长内平均传热系数和平均传质系数稳定。

3 结果分析

3.1 实验验证

为了验证建立的降膜流动模型的正确性，与文献 ［9］中的丙烷降膜流动和换热实验数据进行对比，工况条件见表1，针对以上工况分别进行了数值模拟，流体物性采用0.4MPa饱和温度下的丙烷物性，并假定计算域内局部物性不变。

不同流量下的传热系数数据及实验数据的对比如图9所示。分析实验数据和模拟结果中降膜过程传热系数可知，采用本研究模型模拟降膜流动的过程，89%的传热系数的模拟数据与实验数据偏差不超过25%，模型的计算结果与实验结果的吻合度较好。

3.2 不同干度下流型模拟

基于绕管式换热器流动过程中伴随干度变化，在不同干度下的降膜流动呈现出不同的形态，基于流动模型，可针对不同干度下的降膜流动形态进行模拟。模拟结果如图10所示。

在低干度工况 （0.02），降膜流动呈现为层状流动。此时液膜流速较小，同时气态流动速度也较低，能够形成连续且较为规整的液膜。

干度上升至0.4时，降膜流动流型呈现为柱状流动。液膜的流量降低，无法充满管间间隙形成完整的片状液膜，纵向管间液膜发生破裂，以分散的柱状液流存在。在管道较长的情况下，管间形成多个液柱，液柱及空隙分布位置并不固定，由于管束为非水平方向，液柱在形成后呈现沿管倾斜方向移动的趋势。柱状流动在二维模型中无法实现。

图9 不同入口流量下的传热系数实验验证Fig.9 Experimental verification of heat transfer under different inlet flux

图10 模拟结果Fig.10 Simulation result of flow patterns

而在干度达到0.8的工况下，降膜流动流型呈现为滴状流动，液膜在管间的纵向流量降低，并呈非连续的间断性变化，此时液膜的运动主要表现为受重力作用下沿管倾斜方向的流动。

对不同尺寸参数模型的模拟结果进行对比，发现其流型变化规律不同。在管间距较小的结构中，更容易呈现出层状的降膜流动，而柱状流动和滴状流动所需的干度更高。同时滴状的流动也因为窄管间距而更难形成完整的液滴，反而为间歇性的柱状流动形态。

可见结构参数的影响不能够忽略，因此针对不同干度下的流型变化的研究只能基于固定几何参数的结构。

3.3 相变传热传质过程模拟

在引入传热传质模型后，绕管换热器壳侧传热传质过程的模拟结果如图11所示。

图11 降膜流动传热传质模拟结果Fig.11 Simulation result of phase transition

在液膜受到管道加热发生汽化传质后，气泡首先在内部接近壁面处成核并长大，然后跟随液膜流动离开。一部分较大的气泡直接穿过液膜进入气相，同时，由于气泡的产生，汽化后物性改变，造成膨胀增大，液膜表面发生较为剧烈的波动，使得沿管径方向上不同位置出现流量差异，对流动流型产生影响；部分气泡由于气液密度差造成的浮力作用聚集集中到管壁面下方间隙处，并且不断增大，沿管径方向缓慢移动。

由此可见，在降膜流动过程中，因为气液相间的传质对液膜表面带来了更大扰动，加快了降膜流动的流型的转化。

而同时由于存在传质，对传热造成了两方面的影响：一方面潜热传热的存在使得流体温度维持恒定，能够保持温差不变，传热性能上升；另一方面传质造成的波动产生气相与管壁直接换热的情况，汽化后工质物性改变，造成膨胀增大，热导率减小，热容减小，温度上升更快，与壁面的温差减小，弱化了传热。因此传热性能与管壁面液膜包覆面积比密切相关。

在质流密度为70kg·m－2·s－1的工况下对不同干度下的降膜流动进行模拟，得到的换热模拟结果如图12所示，工况条件取自文献 ［9］。

图12 不同入口干度下的传热系数Fig.12 Heat transfer coefficient under different inlet vapor quality

初期，随干度的增加，传热系数明显增大。从降膜流动的原理上分析可知，随干度增大液膜呈现不同流型，流动过程中的扰动增大，同时气相的流速增大，对液膜产生剪切作用，使得液膜在高干度下厚度较低，减小了壁面外部的热阻，在管壁表面仍覆盖液膜的情况下换热能力增强。在干度达到0.7以上时，局部出现液膜非常薄的情形，此时局部的传热系数达到近3万。但厚度的减小造成液膜在较高干度非常容易被撕裂，使得管壁表面出现未被液膜覆盖的干燥区域，壁面与气相直接换热，气相温度快速升高，导致温差减小，换热能力下降。

为了验证不同热通量工况对传热系数的影响，选取干度为0.5的工况，分别在不同热通量工况下进行模拟，结果如图13所示。

由于流体处于两相状态，局部的变化较小，可以认为蒸发温度近似恒定，换热过程中液相温度维持不变，因此在增大热通量的情况下会造成壁面换热温度的上升，这与Aunan使用电加热的测量方法得到的实验数据一致。

而在高热通量的情况下，传热系数则有所下降。这可能是由于更高热通量带来了更多的汽化传质，使得气泡覆盖了更多壁面，同时造成的更大流态的波动使得液膜更容易破裂，造成了更多干区直接与气相换热，换热能力下降。

图13 不同热通量下的传热系数和传热温差Fig.13 Influence of heat flux on heat transfer

4 结 论

（1）基于VOF模型、连续表面张力模型、接触角模型建立了绕管式换热器壳侧降膜流动的数值模拟模型。

（2）针对气液两相间的传质，对相间质量传递方程引入传质质量源项以预测壳侧的管壁壁面降膜蒸发的气泡形成和生长过程，对能量方程加入潜热传热源项以预测壳侧不同工况下的传热系数。

（3）通过与以往壳侧实验的传热系数数据进行对比，对于降膜流动，89%的模型计算结果与实验的平均误差在25%以内，模型能用于预测壳侧降膜流动的传热系数。

（4）对不同干度下的层状降膜流动、柱状降膜流动、滴状降膜流动进行了模拟，结果显示降膜流动的流型不仅与入口工况有关，还受到管束结构参数的影响。

（5）对不同干度下的壳侧丙烷降膜流动换热过程进行模拟，结果显示干度上升，传热系数呈增大趋势。

（6）对不同热通量下的壳侧丙烷降膜流动换热过程进行模拟，结果显示换热量的增大造成换热温差的增大，而传热系数有所下降。

符 号 说 明

A——传热面积，m2

coeff——松弛时间 （传质频率系数），s－1

db——气泡直径，m

E——内能，J

Fσ——表面张力，N

g——重力加速度，m·s－2

hfg——汽化潜热，kJ·kg－1

M——摩尔质量，g·mol－1

˙mlv——汽化传质速率，kg·s－1

˙mvl——液化传质速率，kg·s－1

＾n——相界面函数

R——气体常数

Sm——传质源项，kg·s－1

T——温度，K

t——时间，s

α——体积分数

β——气体的适应系数

ρ——密度，kg·m－3

v——速度，m·s－1

μ——黏度，Pa·s

θ——接触角，（°）

下角标

l——液相

sat——饱和态

v——气相

［1］ Shukri T.LNG technology selection ［J］.Hydrocarbon Engineering，2004，9 （2）：71－76.

［2］ Neeraas B O.Condensation of hydrocarbon mixtures in coilwound LNG heat exchangers，tube－side heat transfer and pressure drop ［D］.Trondheim：Norwegian Institute of Technology，1993.

［3］ Pacio J C，Dorao C A.A review on heat exchanger thermal hydraulic models for cryogenic applications ［J］.Cryogenics，2011，51 （7）：366－379.

［4］ Linde A G.Aluminium plate－fin heat exchangers ［R］.

［5］ Crawford D B，Eschenbrenner G P.Heat transfer equipment for LNG projects ［J］.Chem.Eng.Process，1972，68 （9）：62－70.

［6］ McKeever J，Pillarella M，Bower R.An ever evolving technology ［J］.LNGIndustry，2008：44－49.

［7］ Fredheim A O.Thermal design of coil－wound LNG heat exchangers shell－side heat transfer and pressure drop ［D］.Trondheim：University of Trondheim，1994.

［8］ Shokouhmand H，Salimpour M R，Akhavan－Behabadi M A.Experimental investigation of shell and coiled tube heat exchangers using Wilson plots ［J］.InternationalCommunicationsinHeatandMassTransfer，2008，35 （1）：84－92.

［9］ Aunan B.Shell－side heat transfer and pressure drop in coilwound LNG heat exchangers，laboratory measurements and modeling ［D］.Trondheim：Norwegian University of Science and Technology，2000.

［10］ Neeraas B O，Fredheim A O，Aunan B.Experimental shellside heat transfer and pressure drop in gas flow for spiralwound LNG heat exchanger ［J］.InternationalJournalof HeatandMassTransfer，2004，47 （2）：353－361.

［11］ Neeraas B O，Fredheim A O，Aunan B.Experimental data and model for heat transfer，in liquid falling film flow on shell－side，for spiral－wound LNG heat exchanger ［J］.InternationalJournalofHeatandMassTransfer，2004，47 （14）：3565－3572.

［12］ Lu Xing，Du Xueping，Zeng Min，Zhang Sen，Wang Qiuwang.Shell－side thermal－hydraulic performances of multi－layer spiral－wound heat exchangers under different wall thermal boundary conditions ［J］.AppliedThermal Engineering，2014，70 （2）：1216－1227.

［13］ Killion J D，Garimella S.Simulation of pendant droplets and falling films in horizontal tube absorbers ［J］.Journalof HeatTransfer，2004，126 （6）：1003－1013.

［14］ Lei Xianliang（雷贤良），Li Huixiong（李会雄），Yan Libo（颜 利 波 ），etal.Numerical investigation on the offset characteristics of falling liquid columns over horizontal tubes under the interaction of inter－tube gas flows ［J］.JournalofEngineeringThermophysics（工程热物理学报），2010，31 （11）：1875－1878.

［15］ Wang Xiaofei（王小飞），He Maogang （何茂刚），Zhang Ying （张 颖）.Numerical simulation of the liquid flowing outside the tube of the horizontal－tube falling film evaporator ［J］.JournalofEngineeringThermophysics（工程热物理学报），2008，29 （8）：1347－1350.

［16］ Qiu Qinggang （邱庆刚），Chen Jinbo （陈金波）.Nume－rical simulation of film formation on horizontal－tube falling film evaporators ［J］.JournalofChineseSocietyofPower Engineering（动力工程学报），2011，31 （5）：357－361.

［17］ Luo L，Zhang G，Pan J，etal.Flow and heat transfer characteristics of falling water film on horizontal circular and non－circular cylinders ［J］.JournalofHydrodynamics，2013，25 （3）：404－414.

［18］ Luo Lincong （罗 林 聪），Zhang Guanmin （张 冠 敏），Tian Maocheng （田茂诚），etal.Numerical simulation of liquid film flow on horizontal circular and shaped cylinders for falling film evaporation ［J］.JournalofChineseSocietyof PowerEngineering（动 力 工 程 学 报），2013，34 （4）：710－714.

［19］ Luo L，Zhang G，Pan J，etal.Influence of oval－shaped tube on falling film flow characteristics on horizontal tube bundle［J］.DesalinationandWaterTreatment，2015，54 （11）：2939.

［20］ Pise A T，Korde N U，Salunkhe R H.Numerical simulation of falling droplets／jets over horizontal tubes ［J］.InternationalJournalofFluidMechanicsResearch，2013，40（3）：266.

［21］ Chen J，Zhang R，Niu R.Numerical simulation of horizontal tube bundle falling film flow pattern transformation ［J］.RenewableEnergy，2015，73：62－68.

［22］ Inaoka H，Ito N.Numerical simulation of pool boiling of a Lennard－Jones liquid ［J］.PhysicaA：StatisticalMechanics andItsApplications，2013，392 （18）：3863－3868.

［23］ Gong S，Cheng P.Numerical simulation of pool boiling heat transfer on smooth surfaces with mixed wettability by lattice Boltzmann method ［J］.InternationalJournalofHeatand MassTransfer，2015，80：206－216.

［24］ Lal S，Sato Y，Niceno B.Direct numerical simulation of bubble dynamics in subcooled and near－saturated convective nucleate boiling ［J］.InternationalJournalofHeatand FluidFlow，2015，51：16.

［25］ Sun T，Li W，Yang S.Numerical simulation of bubble growth and departure during flow boiling period by lattice Boltzmann method ［J］.InternationalJournalofHeatand FluidFlow，2013，44：120－129.

［26］ Chen Z，Utaka Y.On heat transfer and evaporation characteristics in the growth process of a bubble with microlayer structure during nucleate boiling ［J］.InternationalJournal ofHeatandMassTransfer，2015，81：750－759.

［27］ Tsui Y Y，Lin S W，Lai Y N，etal.Phase change calculations for film boiling flows ［J］.InternationalJournalof HeatandMassTransfer，2014，70：745－757.

［28］ Kuang Y W，Wang W，Zhuan R，etal.Simulation of boiling flow in evaporator of separate type heat pipe with low heat flux ［J］.AnnalsofNuclearEnergy，2015，75：158－167.

［29］ Lee W，Son G，Yoon H Y.Direct numerical simulation of flow boiling in a finned microchannel ［J］.International CommunicationsinHeatandMassTransfer，2012，39（9）：1460－1466.

［30］ Zhuan R，Wang W.Flow pattern of boiling in micro－channel by numerical simulation ［J］.InternationalJournalofHeat andMassTransfer，2012，55 （5）：1741－1753.

［31］ Wu H L，Peng X F，Ye P，etal.Simulation of refrigerant flow boiling in serpentine tubes ［J］.InternationalJournal ofHeatandMassTransfer，2007，50 （5）：1186－1195.

［32］ Xie L，Xie Y，Yu J.Phase distributions of boiling flow in helical coils in high gravity ［J］.InternationalJournalof HeatandMassTransfer，2015，80：7－15.

［33］ Zhuan R，Wang W.Simulation on nucleate boiling in microchannel ［J］.InternationalJournalofHeatandMass Transfer，2010，53 （1）：502－512.