绕管式换热器壳侧沸腾模拟时的时间松弛参数研究

吴志勇， 刘 洋， 高 阳， 蔡伟华， 姜益强

(1.辽宁石油化工大学 石油天然气工程学院，辽宁 抚顺 113001；2.哈尔滨工业大学 能源科学与工程学院，黑龙江 哈尔滨 150001；3.哈尔滨工业大学 市政与环境工程学院，黑龙江 哈尔滨 150090)

绕管式换热器壳侧沸腾模拟时的时间松弛参数研究

吴志勇1， 刘 洋1， 高 阳1， 蔡伟华2， 姜益强3

(1.辽宁石油化工大学 石油天然气工程学院，辽宁 抚顺 113001；2.哈尔滨工业大学 能源科学与工程学院，黑龙江 哈尔滨 150001；3.哈尔滨工业大学 市政与环境工程学院，黑龙江 哈尔滨 150090)

利用ANSYS Fluent模拟软件下的流体体积(VOF)模型，对液化天然气用绕管式换热器壳侧沸腾过程进行数值模拟，针对控制壳侧沸腾强度的时间松弛参数取值问题展开研究，提出采用汽化当量法确定时间松弛参数。结果表明，壳侧制冷剂为乙烷、丙烷时的时间松弛参数为3 s-1，并由沸腾换热模拟偏差得到验证。汽化当量法可用于各种沸腾仿真计算，对准确模拟传质过程具有使用价值。

天然气液化； VOF模型； 沸腾传质； 汽化当量法； 汽相体积分数

1 绕管式换热器介绍

绕管式换热器作为一种特殊类型的换热设备广泛应用于天然气液化工艺，其特点是结构紧凑、耐低温、抗高压以及多种介质可以同时换热。这种换热器是由许多小管径换热管以螺旋形状分层缠绕在换热器芯筒上而制成的，相邻两层换热管的旋向相反，层间采用隔条保持间距[1-2]。绕管式换热器内部结构见图1[3]。

对于天然气液化工艺，壳侧烷烃制冷剂首先进入管侧进行过冷，天然气与需要过冷的制冷剂在不同的换热管内形成管侧多股流，并一同向上流动。过冷后的管侧制冷剂在换热器外部经节流后导入到壳侧并向下流动，对管侧天然气和需要过冷的制冷剂实现制冷。烷烃制冷剂在壳侧自上向下流动时，吸收管侧热量而汽化，先后经历两相流动和过热流动两种状态。

在绕管式换热器流动与换热研究领域，已公开发表的文献比较少。文献[1]介绍了绕管式换热器结构参数间的耦合关系，给出了壳侧流通面积的计算方法。B.O.Neeraas[4]建立了绕管式换热器管侧实验装置，以实验方式获得了烷烃混合工质在管侧冷凝时的压降、换热数据，并对管侧流动、换热影响因素进行了分析。B.Aunan[5]建立了液化天然气用绕管式换热器壳侧实验台，对乙烷、丙烷在壳侧的流动与换热特性进行了实验研究。M.Moawed[6]以实验方式研究了定热流密度下缠绕管曲率比、扭曲比与壳侧换热系数之间的关系，指出当缠绕管曲率比、扭曲比变大时壳侧单相换热系数明显增大。

图1 绕管式换热器内部结构

在绕管式换热器数值模拟研究方面，文献[2]研究了壳侧气相流动时的压降与换热特点，指出了换热器优化设计的方向。贾金才[7]建立了绕管式换热器三维计算模型，研究了换热管管径、缠绕角、径向比、轴向比对壳侧换热系数的影响，但所有模拟工况均为气相流动，与壳侧两相流动为主的实际工况存在差别。X.Lu等[8]针对绕管式换热器壳侧气态换热过程展开了实验与计算机仿真研究工作，建立了具有三层绕管的实验装置以及CFD仿真计算模型，实验证实了所用CFD模型具有较高的计算精度。

纵观以上文献发现，绕管式换热器流动与换热性能研究可以通过实验和数值模拟两种方法来实现，但数值模拟均是在单相流动条件下进行的，其原因是在数值模拟条件下相变时的传质强度难以准确把握，导致相变换热模拟难度加大。本文对控制绕管式换热器壳侧沸腾强度的时间松弛参数取值问题做深入研究，期望对相变模拟问题具有指导作用。

2 绕管式换热器壳侧模型的建立

2.1 壳侧几何模型

液化天然气用绕管式换热器整体尺寸巨大，其外部直径为3～5 m，高度为10～50 m，加之内部结构复杂，极不利于开展数值模拟研究。本文在研究过程中以文献[5]中的实验数据作为模拟计算检验标准，因此所建壳侧几何模型与文献[5]中的实验设备一致。

文献[5]中的实验设备为LNG绕管式换热器的简化模型，由3层换热管以交替的缠绕方向缠绕而成，每层的换热管并管缠绕数目由内向外依次为3、4、5。由完整的换热管缠绕而成的中间层及两侧流道为实验研究的核心部位；内、外两侧缠绕层均由半剖换热管缠绕而成，并作为内、外边壁使用。换热器简化模型如图2所示。文献[5]对管-壳两侧换热实行解耦，以给定热流方式研究壳侧换热特性，其中内、外缠绕层及中间缠绕层上部为绝热壁面，中间缠绕层下部4排换热管为受热壁面，通过管内电加热丝实现定热流加热。烷烃制冷剂由顶部30个直径为10 mm的分流孔进入壳侧，经过上方多排缠绕管束缓冲后进入到压降与换热测试区。

图2 换热器简化模型

按照文献[5]中的换热器几何参数，本文利用UG NX 6.0建模软件建立了壳侧模型，其几何参数见表1。

表1 壳侧几何模型参数

由于模拟相变换热需要精细的边界层网格，因此在保证计算精度的前提下，壳侧模型网格数量达到1 000万以上，导致模拟计算耗时巨增。考虑到该模型具有轴对称性，本文将模型沿轴向切割36°作为模拟研究对象，使网格数量降至100万左右，此时网格数量与计算精度能够得以兼顾。切割下来的壳侧36°几何模型如图3所示。

(a) 正视图 (b) 侧视图

2.2 控制方程

烷烃制冷剂流入壳侧时，换热管管束被液膜覆盖，壳侧形成降膜流动，适宜采用两相流VOF模型进行计算。ANSYS Fluent下的VOF模型控制方程见式(1)—(4)[9-10]。

连续方程：

(1)

(2)

动量方程：

(3)

式中，p为压力，Pa；g为重力加速度，m/s2；μ为动力黏度，Pa·s。

能量方程：

(4)

式中，μg、μl为汽、液相动力黏度，kg/(m·s)；E为比内能，J/kg；Γe为能量传输源项，W/m3。

控制方程组中的连续方程和能量方程源项通过下列模型进行计算：

(1)相变传质模型。根据W.H.Lee[11]的传质模型，当流体温度高于汽化温度即T≥Tsat时，单位控制单元内由液相转变为汽相的传质量(汽化量)即为Γm，并有：

(5)

反之，当T

(6)

式中，βm为相变传质过程中的时间松弛参数，s-1，其为控制汽化强度的决定性参数，其取值问题是本文研究的重点。

(2)能量传输模型。当汽、液相间传质量Γm确定后，相应的能量传输源项Γe可用式(7)来表达：

Γe=-hLHΓm

(7)

式中，hLH为流体汽化潜热，J/kg，式中的负号表示需要外界向控制单元输入热量。

另外，VOF模型用于湍流时，本文使用RNGk-ε湍流应力模型作为封闭方程，其形式见式(8)—(9)[12]。

湍动能方程：

(8)

耗散率方程：

(9)

式中，k为湍动能，m2/s2；ε为耗散率，m2/s3；Rk为平均速度梯度引起的湍动能产生项，kg/(m·s3)；Eij为湍流时均应变率，s-1；cμ、ck、cε、c1ε、φ0、φ0均为常数，且cμ=0.084 5，ck=1.39，cε=1.39，c1ε=1.42，c2ε=1.68，φ0=4.377。

2.3 边界条件

对于壳侧36°几何模型，流体从顶部3个分流孔流入，在底部端面流出，两侧剖面均近似处理为对称边界，中间缠绕层底部4排换热管为加热壁面，其余壁面均为绝热壁面，边界条件见表2。

表2 边界条件

2.4 网格划分及无关性验证

使用ANSYS ICEM CFD网格划分软件包对壳侧 36°几何模型进行网格划分，采用非结构化四面体及混合网格，先生成面网格，进而自适应生成体网格。在模型中，缠绕管轴向缝隙 1.94 mm为最小几何尺寸，故以此确定面网格尺寸大小。在确保轴向管缝处存在 3 行网格且网格边长比近似于 2 的前提下，做出3套网格，网格参数见表 3。此外，为了准确计算换热系数，在换热管受热壁面设置精细边界层，边界层初始厚度 0.005 mm，增长因子1.2，共计10 层边界层。

表3 壳侧模型网格方案的网格参数

为了考察网格疏密程度对数值求解的影响并确定网格划分方案，本文以文献[5]中的丙烷降膜沸腾流动实验数据作为依据，进行模拟结果对比。实验工况参数为：p=0.2 MPa；质量流率G=48.6 kg/(m2·s)；入口干度x=0；q=3 951 W/m2。数值模拟的换热系数及其偏差即网格无关性验证数据见表 4。

表4 网格无关性验证数据

由表 4 可知，随着网格数量的增多，换热系数模拟偏差略有下降，但变化趋势趋缓。网格方案 1 的网格数量明显高于网格方案 2，导致网格方案 1 的计算时间远大于网格方案 2 所需时间。所以，在综合考量计算精度与耗时的前提下，本文将网格方案 2 作为最终的网格划分方案。

3 沸腾模拟下时间松弛参数的确定

VOF模型中的连续方程含有质量源项，由式(5)—(6)可知，传质量的多少主要取决于时间松弛参数。数值模拟时，若选取的时间松弛参数过小，则沸腾汽化传质量少；反之，若选取的时间松弛参数过大，则必将导致汽化量增大。因此，合理确定烷烃类的时间松弛参数是准确模拟绕管式换热器壳侧沸腾这一实际过程的关键。

纵观国内外文献，目前尚无沸腾传质时间松弛参数相关研究内容。鉴于此，本文提出采用汽化当量法确定时间松弛参数。汽化当量法是指相对于任一沸腾工况，存在一个汽化当量工况(Equivalent-Evaporation Simulation， EES)与之对应，该工况是以沸腾工况下的理论汽化量及两相总流量为依据而构造出来的，属于绝热工况，其汽、液相入口质量流量分别为：

(10)

(11)

汽化当量法的涵义是：在沸腾工况下，工质汽化后在壳侧出口的汽化量应等于由热平衡计算出来的理论汽化量；并且，在两相总流量不变的前提下，沸腾工况下的壳侧出口汽体量应与汽化当量工况下的汽体量相同。

在VOF模型下，壳侧任一截面处关于分相质量流量及干度的计算数据无法读取，仅能得到汽、液相体积分数计算结果，所以汽化当量法只能通过比对出口端流道截面汽相体积分数的方式来使用。具体而言，在同一沸腾工况下，选取不同的时间松弛参数进行沸腾模拟计算，当汽化量不同时流道出口处汽相体积分数必定不同，而能反映真实汽化程度的时间松弛参数值可通过汽化当量工况下的流道出口汽相体积分数来决定。

本文参照文献[5]的实验工况，对乙烷、丙烷的时间松弛参数进行研究。对于同一沸腾工况，时间松弛参数分别取0.1、1.0、2.0、3.0、4.0、5.0 s-1，经过动态模拟并对收敛结果求均值后在流道截面2处(见图3)得到不同的汽相体积分数，然后与汽化当量工况下的对应值做比较从而确定时间松弛参数值。

模拟工况分为两类，一类用于确定时间松弛参数值；另一类覆盖多压力、多热流、多流率情况，用于检验时间松弛参数值的通用性。在两类工况下，壳侧入口干度均为0，其余参数见表5—6。

表5 时间松弛参数确定工况

表6 时间松弛参数检验工况

丙烷和乙烷在流道截面2处的汽相体积分数随时间松弛参数(βm)的变化曲线如图4所示。丙烷与乙烷相应的EES模拟结果也示于图4中。由图4可知，出口汽相体积分数随着βm的增大而变大，即流道内汽化强度增强；βm=3 s-1对丙烷比较适宜，而βm=4 s-1对乙烷较适宜。

图4 出口汽相体积分数随βm变化曲线

当时间松弛参数的取值不同时，加热壁面汽相体积分数及壁温也会不同，由此可能会导致模拟沸腾换热系数发生相应的变化，因此可通过考核沸腾换热系数的方法辅助确定时间松弛参数。βm分别为2、3、4 s-1时换热系数偏差如图5所示。换热系数偏差是根据模拟结果和实验数据而得到的，换热系数计算公式见式(12)。

α=q/[Tw,heat-0.5(T1+T2)]

(12)

式中，Tw,heat为加热管的壁温，K；T1、T2分别为流道截面1、2的平均温度，K。

图5 沸腾换热系数偏差随βm变化

由图5可知，当βm=3 s-1时，丙烷与乙烷均取得了最小换热系数偏差，鉴于乙烷在βm=3 s-1时出口汽相体积分数与EES对应值偏差不大，故将丙烷和乙烷的时间松弛参数统一定为βm=3 s-1，其值将在检验工况下被进一步考察。

βm=3 s-1时，丙烷、乙烷在各检验工况及其相应EES的出口汽相体积分数如图6所示；丙烷、乙烷在各检验工况下出口汽相体积分数与EES对应值之间的偏差如图7所示。相对偏差以EES出口汽相体积分数为基准。由图7可知，丙烷的相对偏差为-15%～20%，乙烷的相对偏差为-20%～25%，表明βm=3 s-1比较理想。

(a) 丙烷

(b) 乙烷

(a) 丙烷

(b) 乙烷

丙烷和乙烷在各检验工况下的沸腾换热系数偏差如图8所示。由图8可知，两种工质的换热系数偏差均控制在-20%～10%，说明从换热角度来看βm=3 s-1也是合适的。

(a) 丙烷

(b) 乙烷

经过汽化当量法以及换热偏差两方面检验后，本文认为βm=3 s-1对以乙烷和丙烷为代表的壳侧冷剂是适宜的，并且在变压力、变热流密度以及变流率情况下具有通用性，在VOF模型下使用此值能够比较真实地反映壳侧沸腾情况。

4 沸腾模拟下的壳侧可视化图像

对各沸腾工况设置时间松弛参数并模拟求解后，即可得到壳侧不同部位的可视化图像。以P1和E1工况为例，在时间松弛参数取值不同时，壳侧不同部位的汽相体积分数分布见表7。可视化部位包括换热管中间缠绕层、壳侧轴向剖面以及流道周向剖面。

从表7可以看出，随着时间松弛参数变大，壳侧汽化量增加，即时间松弛参数控制壳侧沸腾时的汽化强度；同时，通过对流道出口汽相体积分数以及沸腾换热系数偏差的分析，可以确定βm=3 s-1的取值比较适宜。因此，可以认为βm=3 s-1时的可视化图像能够反映实际沸腾过程。

表7 βm变化时壳侧不同部位汽相体积分数分布

5 结 论

(1)汽化当量法能够用于确定时间松弛参数，通过比对模拟工况与汽化当量工况下的汽化量可以明确时间松弛参数的取值。

(2)数值模拟乙烷、丙烷低温制冷剂在绕管式换热器壳侧沸腾时，经过比较模拟工况与汽化当量工况的汽相体积分数，可以确定时间松弛参数βm=3 s-1比较适宜，并且该值适用于变压力、变热流、变流率的不同工况。

(3)时间松弛参数取值变化时，由于汽化强度不同，加热壁面的壁温也会不同，由此导致模拟的沸腾换热系数发生相应变化。因此，可以通过考核沸腾换热系数的方法来辅助确定时间松弛参数。在时间松弛参数βm=3 s-1的条件下，不同工况下的沸腾换热系数模拟偏差均为-20%～10%，进一步证实了乙烷、丙烷时间松弛参数βm=3 s-1的合理性。

[1] 吴志勇，陈杰，浦晖，等.LNG绕管式换热器结构与流通参数计算[J].煤气与热力，2014，34(3)：42-47.

[2] 吴志勇，陈杰，浦晖，等.LNG绕管式换热器壳侧过热态流动的数值模拟[J].煤气与热力，2014，34(8)：6-11.

[3] Steffen H.Near-optimaloperation of LNG liquefaction processes by means of regulation[D].Berlin：Berlin Institute of Technology，2011.

[4] Neeraas B O.Condensation of hydrocarbon mixtures in ciol-wound LNG heat exchangers tube-side heat transfer and pressure drop[D].Trodheim:The Norwegian Institute of Technology，1993.

[5] Aunan B.Shell-sideheat transfer and pressure drop in coil-wound LNG heat exchanger[D].Trodheim：Norwegian University of Science and Technology，2000.

[6] Moawed M.Experimental study of forced convection from helical coiled tubes with different parameters[J].Energy Conversion and Management，2011，52(2)：1150-1156.

[7] 贾金才.几何参数对绕管式换热器传热特性影响的数值研究[J].流体机械，2011，39(8)：33-37.

[8] Lu X，Du X P，Zeng M，et al. Experimental and numerical investigation on shell-side performance of multilayer spiral-wound heat exchangers[J].Chemical Engineering Transactions，2013，35：445-450.

[9] Yang Z，Peng X F，Ye P.Numerical and experimental investigation of two phase flow during boiling in a coiled tube[J].International Journal of Heat and Mass Transfer，2008，51(5)：1003-1016.

[10] 李霞，马贵阳，杜明俊，等. LNG球罐分层翻滚汽化过程的数值分析[J].辽宁石油化工大学学报，2010，30(4)：42-44.

[11] Lee W H.A pressure iteration scheme for two-phase flow modeling[M]. Washington,DC：Hemisphere Publishing，1980：407-431.

[12] Orszag S A，Yakhot V，Flannery W S，et al.Renormalization group modeling and turbulence simulations[C]. Proceedings of the International Conference on Near-wall Turbulent Flows.Arizona：Elsevier Science Publishers B.V.，1993.

An Investigation on the Time Relaxation Parameter Used in Boiling Simulation of Helical Coil Heat Exchanger Shell-Side

Wu Zhiyong1, Liu Yang1， Gao Yang1， Cai Weihua2, Jiang Yiqiang3

(1.CollegeofPetroleumEngineering，LiaoningShihuaUniversity，FushunLiaoning113001，China；2.SchoolofEnergyScienceandEngineering，HarbinInstituteofTechnology，HarbinHeilongjiang150001，China；3.SchoolofMunicipalandEnvironmentalEngineering，HarbinInstituteofTechnology，HarbinHeilongjiang150090，China)

The volume of fluid (VOF) model which was a multi-phase flow model in ANSYS Fluent software was used to simulate the shell-side boiling of helical coil heat exchanger (HCHE), the time relaxation parameter which dominated mass transfer was studied by equivalent-evaporating method, and the value of time relaxation parameter was recommended as 3 s-1for ethane and propane. The value of time relaxation parameter was also verified by heat transfer deviation. The equivalent-evaporating method which was proposed could be used in all kinds of boiling simulations, and had good application value for mass transfer calculation.

Nature gas liquefaction； VOF model； Boiling mass transfer； Equivalent-evaporating method； Vapor volume fraction

1672-6952(2017)05-0053-08

投稿网址：http://journal.lnpu.edu.cn

2016-10-27

2016-11-28

工信部高科技项目“大型LNG绕管式换热器研制”(2013418)。

吴志勇(1977-)，男，博士，讲师，从事天然气液化技术与设备研究；E-mail：1123589558@qq.com。

TE821；TK172

A

10.3969/j.issn.1672-6952.2017.05.011

(编辑 宋锦玉)