绕管式换热器壳侧汽化模拟时的液滴粒径研究

吴志勇 张安昊 国慧敏 建伟伟 蔡伟华

1.辽宁石油化工大学石油天然气工程学院 2.哈尔滨工业大学能源科学与工程学院

绕管式换热器(spiral wound heat exchanger，SWHE)作为一种结构复杂的换热器广泛用于天然气液化领域，具有抗高压、耐低温、多种介质可同时换热的特点。绕管式换热器内部构造见图1[1-2]，直径较小的换热管以螺旋方式分层缠绕在中央芯筒上，相邻两层换热管的缠绕方向不同，层间用金属隔条控制间距[3]。

绕管式换热器用于天然气液化工艺时，天然气以多股流方式在换热管内自下向上地流动；而在换热器壳侧，烷烃制冷剂则自上而下地流过。管、壳两侧流动都存在相变换热，其中烷烃制冷剂以液相进入壳侧吸热汽化，先后经历降膜流、剪切流，并以过热气形式流出换热器。

在绕管式换热器壳侧换热性能研究方面，可借鉴的相关资料较少。Aunan[4]搭建了LNG绕管式换热器壳侧实验装置，对烷烃制冷剂在壳侧的流动与换热性能展开了研究，指出了干度、热流密度、压力、质量流率对换热的影响程度。Wu等[5]对绕管式换热器壳侧降膜流沸腾工况做出了数值模拟，基于Ansys Fluent软件下的VOF模型开展了传质时间松弛系数的研究，给出了乙烷、丙烷冷剂下的传质时间松弛系数的合理取值。李剑锐等[6]利用Fluent模拟出绕管式换热器壳侧降膜流下的沸腾换热特性，指出冷剂干度对换热系数影响大。季鹏等[7]对多个绕管式换热器壳侧换热关联式进行了对比分析，指出修正后的Abadzic关联式在天然气预冷段最为准确。Wu等[8]以数值模拟方法研究了绕管式换热器壳侧降膜沸腾时的换热特性，分析了利用VOF模型模拟剪切沸腾的不适用性。Ren等[9]对滚动工况下的绕管式换热器壳侧降膜流沸腾换热进行了数值模拟研究，认为滚动振幅与周期对换热的影响大于质量流率和压力的作用。

鉴于Wu等[8]指出了使用VOF模型模拟绕管式换热器壳侧剪切流沸腾换热的不适用性，所以本研究建议使用Ansys CFX软件下的欧拉模型对绕管式换热器壳侧剪切汽化流动进行数值模拟研究，并对欧拉模型下影响流动与换热强度的液相平均粒子直径这一关键参数的取值问题做深入研究，以便能够准确模拟换热器壳侧剪切流沸腾换热过程，进而掌握其流动与换热特性。

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

1.1 壳侧几何模型

实际应用的绕管式换热器体积巨大，而内部结构却细致紧凑，在对其进行数值模拟研究时必将导致几何模型网格单元数量巨大，计算耗时长且难以收敛。因此，本研究按照文献[4]中实验装置建立壳侧几何模型，并通过文献[4]中的实验数据对模拟结果予以验证。

文献[4]中的实验装置是对绕管式换热器壳侧的简化，换热管分三层螺旋缠绕，中间层采用断面完整的换热管，而内、外两层使用断面为半圆的换热管，其结构见图2。中间缠绕层下部并行的4根换热管壁面为加热壁面，由管内电加热丝控制热流密度，其他壁面均是绝热壁面。烷烃冷剂由换热器顶部30个分流孔进入壳侧，经过换热器内部多排换热管缓冲后进入到底部4根换热管测量段。测量段的加热壁面温度和流道断面温度由热电偶多点实测获得。该实验装置的几何参数见表1。

表1 壳侧模型几何参数几何参数参数值换热管外径/mm12.00换热管轴向间距/mm13.94换热管径向间距/mm15.91换热管缠绕角/(°)7.938换热管缠绕层数/层3 缠绕层并管数目/根3,4,5缠绕直径/mm96,127.82,159.63 模型高度/mm160 分流孔直径/mm10 分流孔数目/个30

鉴于该实验装置具有轴对称性，因此本研究在几何模型建立之后，将模型沿轴向切割36°作为仿真计算研究对象，其目的是降低网格数量，减少计算耗时。切割下来的壳侧36°几何模型见图3。

1.2 控制方程

在绕管式换热器壳侧，当冷剂入口干度x>0.2时，由于气、液相密度比很小，所以汽相占据壳侧绝大部分空间，液相以离散相液滴形式存在于汽相之中以及出现在换热管壁面上。在汽相吹扫作用下，换热管壁面会出现暂时的干涸现象。此时，汽、液间作用力明显，相间滑速比增大，在这种状态下适宜采用欧拉模型进行模拟计算。

欧拉模型控制方程组如下[10]：

连续方程：

(1)

(2)

动量方程：

(3)

(4)

能量方程：

(5)

(6)

(1) 相间传热模型

单位控制单元内液相向汽相的传热量可通过粒子模型进行求解，将液相定为离散相粒子公式如下[10]：

Qlg=αlgAlg(Tl-Tg)

(7)

(8)

(9)

(10)

(11)

(12)

式中：αlg是汽、液相间换热系数，W/(m2·K)；Alg是汽、液相界面面积浓度，m2/m3；dl是液相粒子直径，m。

(2) 相间传质模型

(13)

Qsg=αgAlg(Ts-Tg)

(14)

Qsl=αlAlg(Ts-Tl)

(15)

(16)

(17)

式中：Qsg、Qsl分别为单位控制单元内相界面向汽相和液相的传热量，W/m3；αg、αl分别为汽相和液相的换热系数，W/(m2·K)；Ts为相界面温度，K。当液相粒子采用零热阻换热模型时，有Ts=Tl。

(3) 相间作用力模型

(18)

CD=max{24(1+0.15Re0.687)/Re,0.44}

(19)

此外，欧拉模型在湍流计算时，汽相作为连续相使用标准k-ε湍流模型，液相以离散粒子形式使用零方程模型进行简化处理，此时对应的湍流方程分别为[10]：

汽相湍动能方程：

(20)

汽相耗散率方程：

(21)

液相湍流黏度：

(22)

式中:σk=1.0，σε=1.3，c1ε=1.44，c2ε=1.92，σ=0.9。

1.3 边界条件

对于绕管式换热器壳侧36°几何模型，冷剂从上部3个分流孔进入，在底部端面流出，两个轴向端面设为对称边界，中间缠绕层下部4根换热管的外表面是加热壁面，其余壁面属于绝热壁面，边界条件设置见表2和图3。

表2 边界条件模型部位边界条件类型顶部3个分流孔质量流量入口底部端面敞开流动轴向端面对称边界1、2中间层底部4根管壁加热壁面给定热流其余壁面绝热壁面

2 离散液相平均粒子直径推导

绕管式换热器壳侧剪切汽化流动模拟时，汽、液相分别设置成连续相和离散粒子，以便描述汽、液相间动量、能量交换，其中式(8)和式(18)的计算需要引入液相平均粒子直径，但该参数尚无计算方法，给剪切流汽化模拟带来困难，因此，需对液滴平均粒子直径作理论研究。

壳侧剪切汽化流动时，汽、液相温度并不相同，汽相温度高于液相温度，流场内弥散状液滴通过相界面从汽相获取热量实现汽化。液滴平均粒子直径可以通过微元控制体能量平衡推导出来。流场内部微元控制体见图4，图中蓝色部分表示液相，其余部分表示汽相。

假设沿x、y、z方向流入微元控制体的单位体积能量为hV,x、hV,y、hV,z，相应的流速为u、v、w。则x方向流入能量为：

Hx=hV,xudydz

(23)

同理，y、z方向流入能量为：

Hy=hV,yvdxdz

(24)

Hz=hV,zwdxdy

(25)

x方向流出能量为：

(26)

同理，y、z方向流出能量为：

(27)

(28)

当汽、液相间换热采用两热阻模型且液相侧使用零热阻时，有αlg=αg，Ts=Tl。稳态时，微元控制体内汽相向液相传递的热量即是外界应向微元控制体输入的汽化热量，也即是微元体的净导入热量，该热量为：

Qlg=αgAlg(Tg-Tl)dxdydz

(29)

稳定流动下根据微元控制体能量平衡可得到：

Hx+Hy+Hz+Qlg=Hx+dx+Hy+dy+Hz+dz

(30)

因dx、dy、dz均趋于0，所以有：

hV,x=hV,y=hV,z=HV

(31)

将式(23)～式(29)及式(31)代入至式(30)，简化后得到：

(32)

将式(8)及式(12)代入至式(32)中，得到：

(33)

在流场范围Ω内对式(33)进行三重积分，有：

(34)

式(34)左侧积分利用高斯定理可转化为：

(35)

式中:n为流场边界曲面数量，ϑx、ϑy、ϑz为边界曲面外法线与坐标轴夹角。将壳侧模型对称面视为壁面并且以无滑移处理壁面后，式(35)右侧可简化为：

(36)

式中:Aheat为换热管加热壁面的面积，m2；q为换热管壁面热流密度，W/m2。

式(34)右侧积分时利用中值定理可变为：

(37)

(38)

将式(36)和式(38)代入至式(34)中，变形后得到：

(39)

3 热相变比的应用

(40)

式中:Δx为壳侧沸腾时的干度增加量。从定义式可以看出，热相变比ζ表征了向壳侧输入热量中，用于液滴汽化的热量占总热量的比例。

将式(40)变形后代入至式(39)，得到：

(41)

表3 乙烷液相平均粒子直径及其工况参数压力/MPa入口干度质量流率/(kg·(m2·s)-1)热流密度/(W·m-2)热相变比平均粒子直径/m换热偏差/%0.20.2144.947 839.130.967 60.000 082-15.540.3152.587 839.130.971 60.000 04010.760.4243.817 839.130.974 40.000 033-19.980.5361.067 839.130.978 10.000 022-16.470.6251.736 315.290.983 10.000 014-18.780.7251.164 822.870.991 60.000 007-15.180.8251.737 839.130.994 70.000 00414.390.9142.977 839.130.997 30.000 002-14.10.30.2161.347 839.130.947 30.000 12117.250.3260.774 744.320.958 40.000 06219.70.4260.773 157.650.959 50.000 044-2.980.5243.536 315.290.969 40.000 0335.030.6244.107 839.130.980 30.000 0202.270.7343.537 839.130.987 70.000 011-16.190.8742.687 839.130.990 30.000 006-17.890.40.2154.274 822.870.941 50.000 13713.050.3261.342 364.310.949 80.000 07614.750.4260.492 364.310.956 20.000 0565.630.5343.534 744.320.961 10.000 0408.310.6259.647 839.130.967 10.000 0245.480.7260.217 839.130.972 70.000 016-9.460.8157.957 839.130.984 20.000 009-13.970.9159.367 839.130.992 10.000 0041.66

本研究依据文献[4]的实验数据，研究了乙烷在变工况条件下做剪切汽化流动时热相变比的取值问题。通过设定不同的ζ值，将模拟工况下的换热结果与实验数据做比对来确定合理的热相变比。模拟工况共选取23个，涵盖了变压力、变干度、变热流、变流量的情况，经过百余次计算，确定出来的ζ值见表3。

通过表3可以看出，壳侧冷剂入口干度(x)对热相变比的影响最为敏感，两者间呈同向变化，即干度增大时要求热相变比相应变大，从而减小液相平均粒子直径，以满足提升换热系数的要求。而当质量流率(G)和热流密度(g)变化时，两因素对热相变比的影响不明显，即无法改变冷剂干度对热相变比的影响走势，因此，热相变比的取值可以忽略质量流率和热流密度的影响。

图5是乙烷23个模拟工况按定压条件绘制出的热相变比变化图。从图5可以看出，热相变比受压力的影响较为明显。在壳侧入口干度相同的条件下，压力越小则要求热相变比越大。其原因在于，高压向低压转变时,相对于其他物性参数而言，汽相密度变化最为明显，汽相密度变小将导致汽相流速增大，随之而来的相间动量、能量传输增强，必然要求液相平均粒子直径进一步变小，即热相变比需要相应增大。

鉴于变压力时，汽相密度对热相变比有较大的影响，因此本研究对热相变比进行公式拟合时自变量包含了干度x和汽相密度ρg两个物理量，根据模拟工况下的计算数据拟合出的热相变比公式见式(42)：

(42)

a1=-0.155，a2=0.088，a3=7.606，a4=-0.207，a5=0.034，a6=7.564；

b1=8.122，b2=-14.247，b3=7.879，b4=-13.859

式中：ρ/7.36为变压力时汽相密度发生改变而做出的修正，以乙烷0.4 MPa下汽相密度7.36 kg/m3作为基准进行热相变比修正。上式拟合程度较好，拟合相关系数R=0.993 9。

图6和图7分别为乙烷于换热器壳侧做剪切汽化流动时在轴向剖面、中间缠绕层壁面处的汽相体积分数模拟云图，其模拟工况参数p=0.4 MPa、G=59.36 kg/(m2·s)、x=0.91、q=7 839.13 W/m2。

通过图6和图7可以看出，壳侧冷剂高干度条件下的两相流动趋于汽态流动，液相是以极细的液滴粒子形式弥散于汽相当中，因此，流道内每个网格单元的汽相体积分数接近于1。在中间缠绕层壁面上，由于液相润湿效应而使壁面局部位置处的液相体积分数略有增加。在近似于汽态流动、换热的条件下，只有将热相变比增至0.992 1，液滴平均粒子直径减小到0.000 004 m，才能准确地模拟出换热系数高达12 153 W/(m2·K)的汽化工况，并使模拟偏差低至1.66%。

4 结论

(1) 两相流欧拉模型适用于绕管式换热器壳侧剪

切汽化流动时的数值模拟研究，通过理论推导以及引入热相变比参数，可以确定出离散液相的平均粒子直径，进而实现对壳侧剪切流汽化换热的准确模拟，其换热模拟偏差在±20%范围内。

(2) 以乙烷在绕管式换热器壳侧剪切汽化流动为例，通过换热模拟结果与实验数据相对比的方式，确定出变工况条件下的热相变比取值。经分析发现，壳侧冷剂入口干度对热相变比的影响最为强烈，冷剂工作压力对热相变比的影响次之，而壳侧质量流率、热流密度对热相变比的影响可以忽略不计。

(3) 在分析热相变比影响因素之后，选取壳侧冷剂入口干度和汽相密度作为自变量对热相变比进行公式拟合，给出了相应的计算式，其拟合程度较好，相关系数R=0.993 9。

(4) 数值模拟研究结果表明，通过热相变比确定离散相平均粒子直径的方法，可以使欧拉模型准确地模拟出两相流相变换热过程，解决了相变模拟失真这一突出矛盾。该方法对两相流仿真计算具有实用价值。

符号表符号物理量名称单位Alg汽、液相界面面积浓度m2/m3Aheat换热管受热面积m2Cp定压比热J/(kg·K)d离散相粒子直径mF相间作用力N/m3g重力加速度m/s2h比焓J/kg

续符号表符号物理量名称单位hV单位体积焓J/m3H单位时间总焓Wk湍动能m2/s2m壳侧质量流量kg/m单位控制容积相间传质量kg/(m3·s)Nu努塞尔准则p压力PaPr普朗特准则q热流密度W/ m2Q传热量W/ m3Re雷诺准则T温度Ku、v、w直角坐标系下速度分量m/sv速度矢量m/sV体积m3x干度α换热系数W/(m2·K)β体积分数ƍx、ƍy、ƍz曲面外法线与坐标轴夹角(°)ρ密度kg/m3λ导热系数W/(m·k)ε湍流耗散率m2/s3τ时间sζ热相变比μ动力黏度N·s/m2下角标符号含义符号含义g汽相gl汽相作用于液相l液相lg液相作用于汽相model壳侧几何模型s相界面sat 饱和t湍流