中间介质式气化器LNG跨临界传热特性模拟与分析

姚宝龙，唐建峰，王国瑞，王剑琨，于 笑，许 涛

（1. 中国石油大学（华东） 储运与建筑工程学院，山东 青岛 266580；2. 中国石化天津液化天然气有限责任公司，天津 300457；3. 中国石化青岛液化天然气有限责任公司，山东 青岛 266580）

液化天然气接收站作为LNG产业链中的重要环节，近年来在我国得到大力发展。截止2020年，LNG接收站数量已达22座[1]。接收站储罐内的LNG使用时经过气化器气化后外输。中间介质式气化器（Intermediate Fluid Vaporizer，IFV）由于其没有结冰问题，对海水质量要求低、占地面积小等优点，在沿海LNG接收站及海上浮式气化场景中有广阔的应用前景。

IFV结构简图如图1所示，根据换热过程可将IFV分为蒸发区（E1）、凝结区（E2）和调温区（E3）。IFV热源介质为海水，中间介质为丙烷。LNG在E2段吸收管外丙烷的热量后气化成为NG，最后经E3升温后外输，丙烷在E2段壳侧释放热量凝结成液，在重力作用下进入E1段与管内海水换热蒸发后返回至E1段，形成丙烷的“蒸发-冷凝”内部循环过程。近些年随着气化器国产化进程的加速，国内已有部分厂家掌握了IFV的设计及制造，但现场应用经验不足。国内学者对IFV的研究更多集中在材料、结构及运行操作等方面，对其传热机理研究较少。刘丰等[2]开展了IFV在海上浮动工况下的适用性研究，从结构、材料等方面进行了论证。陈双双等[3]对比了不同IFV的布置方式，采用MATLAB软件分析了整体式IFV与分体式IFV的换热性能。王博杰等[4]利用CFD数值模拟对LNG在IFV内的流动传热特性进行研究，得到了IFV内部流动传热与温度分布特性。IFV涉及多流体的不同换热过程，其内部复杂的换热机理仍未明晰。其整个传热过程有两个难点，一是中间介质工作状态的确定，二是凝结区管内LNG的跨临界流动传热过程。IFV在使用过程中工作压力一般在6~10 MPa范围内，LNG在与管外丙烷换热后由亚临界态转变至超临界态。目前，针对流体跨临界流动传热的研究主要集中在H2O和CO2等常规工业流体中[5-8]，针对特定组分LNG的研究较少，同时文献中研究大多针对水平圆管和竖直管内流体的跨临界流动传热过程[9-11]，无法对IFV内U型管工况提供参考。

图1 IFV结构Fig. 1 Structure of IFV

本文对LNG在IFV中的跨临界流动换热过程进行数值模拟，对超临界状态下LNG热物性进行计算，探究管内LNG在跨临界升温过程中的传热特性与流动特性，分析操作参数对IFV运行情况的影响。

1 超临界态下LNG热物性计算

在进行数值计算时需要准确获得实际LNG组分下各物性随温度的变化规律。LNG通常含有含量（物质的量分数，下同）为90%以上的甲烷，此外还有少量的C2H4、C3H6、C5+及N2等。表1为国内沿海某接收站贫LNG组分。根据表1中LNG组分，利用工质物性计算软件Refprop 9.1对所需工况下LNG的热物性进行计算，以9 MPa为例计算结果如图2所示。由图2(a)可知，LNG在跨临界过程中各项热物性随温度连续变化，气液分界面消失，不存在明显的相变现象。由定压比热容曲线可以看出，LNG在跨临界升温过程中定压比热容存在明显的极值点，该点对应的温度称为拟临界点。由图2(b)可知，压力越高，定压比热容曲线越平缓且峰值越低。本文采用变物性的湍流模型对此过程进行描述分析[12]，利用上述计算得到的LNG热物性随温度变化的规律，通过FLUENT软件自带的线性差值（Piecewise-liner）方式输入到计算过程中。

图2 9 MPa时LNG热物性随温度的变化(a)和压力对LNG定压比热容的影响(b)Fig. 2 Variation of LNG thermophysical properties with temperature at 9 MPa (a) and effect of pressure on specific heat capacityof LNG (b)

表1 贫LNG组分Table 1 Lean LNG components

2 数值计算模型

2.1 物理模型

国内沿海某LNG接收站使用的国产化IFV，具体尺寸参数及性能指标见表2。其凝结区内有多根LNG传热管，各传热管结构相同，选取单根传热管作为研究的计算域。IFV冷凝区内传热管规格如表3所示。利用ICEM对计算域进行建模及结构化网格划分。对弯管部分进行网格加密，以更好地捕捉弯管区域流体流动特性。网格划分如图3所示，其中Z轴起始点为管段入口处。

表2 IFV设计参数Table 2 IFV design parameters

表3 凝结区传热管规格Table 3 Specification of heat transfer tube in condensation area

图3 弯管及其横截面网格划分Fig. 3 Bend and its cross section meshing

LNG在传热管内的流动受到加热壁面的约束，在近壁区会形成边界层，FLUENT为湍流模型提供了不同的壁面函数来处理边界层内的流动与传热。本文采用标准壁面函数处理近壁面区的流动，标准壁面函数需要将第一层网格节点置于对数律层，引入无量纲化的壁面距离y+定义第一层网格高度，对数律层要求y+的范围在30~300之间[13]，y+的定义式如式(1)。通过式(1)可估算第一层网格高度y。

式中，y为第一层网格高度，m；τw为壁面的切应力，MPa；ν为流体的平均动力黏度，Pa·s；ρ为流体密度，kg/m3。

在符合上述网格划分要求下，划分了4组不同网格数量的换热管模型并进行网格独立性检验，结果如图4所示。由图4可知，随着网格数量的提升，计算结果逐渐趋于稳定，考虑到计算时间成本，选择网格数量为1264784的网格进行后续计算。

图4 网格独立性验证Fig. 4 Grid independence verification

2.2 数学模型

为使得数值计算结果更精确，现对LNG在IFV换热管中的实际气化传热过程进行如下假设：（1）将LNG平均分配到每根传热管，且LNG经过管线到达气化器入口时处于充分发展段；（2）当气化器壳侧工况达到稳定后，管外壁温度与丙烷气相主体温度相差不大，忽略该部分温差，设置模型中传热管的外壁温度为恒定值；（3）管程流体的压降（0.04 MPa）在数量级上远远小于工作压力，可忽略不计。基于上述假设，LNG在传热管内的流动传热过程满足连续性方程、动量守恒方程和能量守恒方程，如式(2)~式(5)。

连续性方程：

式中，为Hamilton算子；为速度向量，m/s。

对于不可压缩牛顿流体，N-S方程如下：

式中，E为热力学能，J；keff为流体的有效导热系数，W/(m·K)；方程右边前两项分别是由导热引起的能量传递和粘性耗散引起的能量传递；τeff为粘性应力张量；Sh为体积热源项，W/m3。

整个传热过程还应包括换热管壁中的热传导：

式中，xi代表i的方向，i= 1, 2, 3；λ为固体壁面的热导率，W/(m·K)。

IFV内部LNG传热管是U型传热管，而RNGk-ε模型在模拟强旋流、弯曲壁面或弯曲流线流动时表现出较好的吻合性。本研究选取RNGk-ε模型，壁面函数选取Standard Wall Functions，该模型的输运方程如式(6)~式(7)。

湍流动能方程：

式中，k为流体的湍动能，m2/s2；αk为k方程的湍流普朗特数；ui为速度，m/s；μeff为有效动力黏度，Pa·s；Gk为由平均速度梯度引起的湍流动能，m2/s2；Gb为由浮力引起的湍流动能，m2/s2；ε为湍流耗散率，m2/s3；YM为可压缩湍流流动脉动膨胀对总的耗散率的影响；Sk为k方程中的源项。

湍流耗散率方程：

式中，αε为ε方程的湍流普朗特数；G1ε和G2ε均为模型常量；G3ε为决定ε方程受浮力影响程度的量；Sε为ε方程中的源项。

本文采用商业CFD软件FLUENT1 9.0对数值过程进行求解。求解方法采用基于压力-速度耦合的SIMPLEC算法，压力项采用Body Force Weighted，动量方程和能量方程采用QUICK格式进行离散，时间步长设置为0.001 s，设置迭代步数为2000步。在计算过程监测出口截面的平均温度，当出口截面平均温度不再变化后，认为计算收敛。

2.3 模型准确性验证

以某LNG接收站夏季某日现场运行数据为基础，对数值模型的准确性进行验证，该日实际运行数据如表4所示，验证结果如表5所示。通过对比模拟数值与实际监测数值可知，模拟计算得到的NG出口温度与实际运行数值误差为0.1%，出口压力与实际数值误差为4.9%，考虑到模型中海水温度、LNG入口温度等与实际运行存在一定的偏差，因此误差在合理范围之内，可利用此数值模拟方法进行进一步的研究。

表4 某LNG接收站气化单元IFV某日实际运行数据Table 4 Actual operation data of IFV of gasification unit in a LNG terminal on a certain day

表5 IFV模拟数值与实际数值对比Table 5 Comparison between IFV simulation value and actual value

3 结果与讨论

3.1 管内温度场及速度场分析

以表4中的工况为例，对LNG在管内跨临界流动换热的瞬态数值模拟结果进行温度场分析和速度场分析，结果分别如图5和图6。由图5可知，管程的前半部分LNG处于液态，与管壁之间温差较大，此时流体升温较快，在温度跨过拟临界点后，流体的热导率迅速降低，且流体温度与壁面温度温差减小，升温变缓。由图6可知，在入口速度为0.48 m/s的工况下，沿管程流体速度不断增大，出口处截面平均速度达到10.35 m/s，出现了明显的流动加速效应。管内LNG沿程温度及密度变化见图7。由图7可知，LNG在管内吸热后密度迅速降低，由424 kg/m3下降至77.93 kg/m3，过程中LNG不断膨胀产生轴向推动力使得流体不断加速。

图5 管内LNG沿程温度场分布Fig. 5 Distribution of LNG temperature field along pipe

图6 管内LNG沿程速度场分布Fig. 6 Distribution of LNG velocity field along pipe

图7 管内LNG温度及密度沿程变化Fig. 7 Variation of LNG temperature and density along pipe

3.2 拟临界点传热特性分析

由图2(b)可知，LNG在不同压力下拟临界温度不同，查询制冷剂热物性计算软件可知，该组分LNG在入口压力为9.26 MPa下的拟临界温度为218.400 K。图8给出了管内LNG沿轴向2~8 m各处截面热物性参数的分布云图。从管内LNG沿程的温度分布与密度分布可以看出，流体在进入换热管后物性变化较为剧烈，在此区域内LNG经历了跨临界的状态转变过程。在Z= 2 m处，由图8(a)可知，圆管顶部小部分的LNG较早达到拟临界温度，由于超临界压力下LNG的热物性对温度有较强的依赖性，由图8(b)可知，在Z= 2 m处，靠近圆管上顶部LNG密度下降较快，截面上密度梯度较大，由于受到重力的作用，圆管各截面的密度并不是呈现“同心圆”式的分布，而是出现了明显的分层现象。随着LNG向前流动，当管程超过6 m后，LNG主体温度均超过拟临界温度，截面上密度变化开始趋缓。图8(c)为不同位置截面上定压比热容变化云图，由图8(c)可知，定压比热容的峰值首先出现在横截面上部，随着LNG主体温度达到拟临界温度，截面上定压比热容的峰值不断向截面中心及下部移动。LNG进入传热管后靠近壁面的部分被加热而密度降低，受到重力的作用，LNG在截面上开始重新分布。

图8 沿轴截面管内LNG热物性分布云图Fig. 8 LNG thermophysical property distribution cloud chart on pipe section

图9为距离LNG入口2 m处截面流线轮廓。从图9可知，靠近管内壁，较轻的流体受到下部较重流体的挤压后沿管壁被推至上部，在圆管截面上形成径向对称的“二次环流”，使得管内流体沿轴向加速的过程中，截面上物性也在重新分布。

图9 距入口2 m处截面流线轮廓Fig. 9 Section streamline at 2 m from entrance

3.3 弯管处流动换热分析

LNG在弯管处受到壁面的阻碍，流动方向发生突变，在浮升力和离心力的共同作用下，弯管处的流动换热情况与直管段有着较大的区别，建立如图10所示的3个观测面（1-1'、2-2'和3-3'），分析LNG在弯管处的流动情况。

图10 弯管处观测面Fig. 10 Observation surface at bend

弯管处3个观测面LNG温度分布如图11所示。由图11(a)可知，相比于弯管外侧区域，靠近弯管内侧区域流体温度更高。由图11(b)可知，进入弯管前，由于物性的重新分布，较重的LNG在管内下部流动，进入弯管后，受到管壁的阻碍以及离心力的作用，较重的LNG沿弯管外侧继续流动。在截面内，Y方向的分速度逐渐增大，外侧较重的LNG受到管壁和上升部分LNG的挤压沿管壁向内侧流动，在截面上形成了两个对称的“迪恩涡”[14]。当LNG运动至2-2'截面时，Z方向的分速度逐渐趋于0，X方向和Y方向的分速度达到最大值，此时在横截面上出现了4个对称的“迪恩涡”，截面流场被完全重构，流体边界层厚度被进一步削弱，显著加强了弯管区域的局部对流换热系数。

图11 弯管观测面处LNG温度分布及流线轮廓Fig. 11 LNG temperature distribution and streamline at observation surface of bend

3.4 入口流量对传热的影响

实际生产运行中，由于受到下游用户季节性需求不同的影响，单台IFV流量波动范围较大。选定单台IFV，以外输量50~200 t/h作为研究区间，根据IFV凝结区LNG换热管尺寸及数量，换算为单管质量流量，如表5所示。针对LNG不同入口质量流量进行模拟计算，分析不同入口流量对LNG流动传热的影响。

表5 单台IFV气化量对应的入口质量流量Table 5 Inlet mass flow corresponding to gasification volume of single IFV

LNG入口压力为9.26 MPa，丙烷饱和温度为288.15 K时，入口质量流量对管内LNG对流换热的影响如图12所示。由图12(a)可知，同一截面LNG温度随着入口质量流量的增加而减小，这是因为随着质量流量的增加，相同质量的LNG经过截面的时间减少，使得LNG与管壁的热量交换减少。由图12(b)可知，同一质量流量下，局部对流传热系数呈现先增大后减小的趋势，在弯管处达到最大值。这是因为在弯管前，LNG主体升至拟临界点温度，定压比热容达到最大值，使得局部传热系数出现峰值。质量流量越大，局部传热系数的峰值越靠后，这是因为随着质量流量的增加，LNG流速增大，相同质量的LNG吸热达到拟临界温度点的换热时间增加，导致峰值位置出现延后。从图12(b)还可以看出，LNG入口质量流量越大，局部对流换热系数也越大，这是因为质量流量变大，增加了管内流体的湍流强度，有效抑制了管壁处流体边界层增厚，降低了对流热阻，从而增强了流体的对流传热能力。

图12 质量流量对管内LNG沿程传热的影响Fig. 12 Effect of mass flow rate on heat transfer of LNG along pipe

3.5 丙烷饱和温度对传热的影响

在凝结区中，气态丙烷在管外凝结成液，其饱和温度的大小直接关系到LNG在传热管内的流动换热情况，进而影响到IFV的气化能力。在LNG入口温度为125 K、入口质量流量为0.042 kg/s下，选 择268 K、278 K、288 K和298 K 4种丙烷饱和温度工况进行模拟计算，结果如图13所示。

图13 丙烷饱和温度对管内LNG沿程传热的影响Fig. 13 Effect of propane saturation temperature on heat transfer of LNG along pipe

由图13(a)可知，丙烷饱和温度越高，相同截面处流体平均温度越高，这在管程的后半段表现更为明显。这是由于LNG在拟临界温度前主要是液态，物性对温度的变化并不敏感，流体升温较快，到达拟临界温度点之后LNG定压比热容随着温度的升高迅速降低，同时较高的丙烷气相温度会增大流体与管壁间的温差，同一截面处的LNG平均温度也随之增加。由图13(b)可知，丙烷的饱和状态对管内LNG的局部传热系数有较大的影响。丙烷的饱和温度一定时，局部对流传热系数先升高，达到峰值后迅速降低，原因可能是当LNG跨过拟临界温度后，发生了传热恶化的现象。结合图 8(a)可知，顶部母线内壁附近流体会比底部母线附近流体更早达到拟临界温度，由于LNG热物性对温度有强烈的依赖性，顶部母线附近的LNG密度迅速下降，体积比底部母线附近的流体膨胀的更快。此时由于密度差及传质的影响，LNG向底部区域流动，顶部沿轴向的质量通量开始减小，同时由图8(c)可以看出，随着温度的升高，顶部母线附近LNG热导率迅速下降，导致传热开始恶化，且丙烷饱和温度越高，这种现象越明显。

4 结论

采用数值模拟研究了LNG在IFV U型管中跨临界流动传热的过程，探究了不同入口质量流量、不同丙烷饱和温度对LNG流动传热的影响，分析了LNG在拟临界点附近的传热机理以及在弯管处的流动传热特性。主要结论如下：

（1）升温过程中，LNG在超临界压力下沿轴向物性变化剧烈，出现了明显的“流动加速”效应。在截面上，由于温度分布不均，密度会出现明显的分层现象，从而形成对称的径向“二次环流”，加强了流体的换热能力。

（2）在弯管处，受到离心力的作用及管壁阻碍，在截面上形成了对称的“迪恩涡”，不同观测面上流线分布差异较大，当LNG运动至弯管的中心截面时，径向流动进一步加强，在截面上形成了4个对称的“迪恩涡”，内部流场被完全重构，流体边界层被进一步削弱，显著加强了弯管处的局部换热系数。

（3）局部对流传热系数沿管程呈现先增大后减小的趋势，在拟临界点附近达到最大值。随着入口质量流量的增加，局部对流传热系数的峰值明显增大，峰值出现的位置向后移动。同一截面处LNG温度呈现减小的趋势，但对出口NG温度影响不大。

（4）随着丙烷饱和温度的升高，相同截面处流体平均温度也越高，但对局部传热系数的影响相反，局部传热系数的峰值随着丙烷饱和温度的升高而降低。随着丙烷饱和温度的升高，U型管内外介质温差增大，容易发生传热恶化的现象。