LNG气化器管内气泡行为及两相流型转换

徐少杰, 高文学, 严荣松, 王 艳, 杨 林, 张 欢

(1.天津大学 环境科学与工程学院，天津300072;2.中国市政工程华北设计研究总院有限公司城市燃气热力研究院，天津300384)

1 概述

在中小型LNG供气厂站中，空温式气化器是比较常见的一种气化设备，其利用空气作为热源，对翅片管内的LNG进行加热。目前关于LNG空温式气化器的研究，大多数集中于单根星型翅片管的传热问题，内容从管内的流动沸腾到管外的自然对流与湿空气结霜[1-2]。虽然不同流型的转换机理对于流动沸腾传热特性的研究具有重要意义，但是目前对于翅片管内天然气气泡行为以及气液两相流型转换的研究远远不足。竖直管内的流动沸腾沿管程方向依次表现为泡状流动、弹状流动、环状流动以及弥散状流动[3]。李祥东等[4]通过分析研究两相流型常用理论，认为常用理论很难准确地判断低温流体的流型变化。刘珊珊等[5]使用经典Lee模型[6]研究了两相流型变化与局部传热系数之间的联系，表明局部换热系数随气相分数呈现先增后减的趋势。文珏[7]的研究结果发现，与Lee模型相比，S.Hardt模型[8]对气液的蒸发沸腾过程描述更为准确。

虽然近年来采用计算流体力学方法对空温式气化器内外场耦合传热特性的研究已有诸多报道[9]，但是更多是从宏观层面探索空温式气化器的传热规律，然而，气泡在沸腾表面的形成、生长以及脱离对于流动沸腾传热有非常显著的影响。为了研究管内LNG沸腾过程气泡的行为和两相流型的转化，本文借助于Fluent 19.0软件，采用流体体积法(Volume of Fluid，简称VOF)模型并耦合S.Hardt模型，计算空温式气化器翅片管内的气液两相流动沸腾过程，分析不同LNG入口流速、LNG入口温度以及管外表面传热系数等条件对管内流动沸腾换热的影响。

2 问题描述及假定说明

2.1 问题描述

空温式气化器由呈阵列形式布置的铝合金翅片管簇组成，LNG由翅片管底部进入，自下而上流动，在此过程中不断从管壁吸收来自环境空气的热量。当LNG从入口温度上升到饱和温度后，管内流体进入蒸发沸腾阶段，形成气液两相流动。随着管内流体温度升高，天然气的体积分数不断增大，管内气液两相流型将出现显著变化。本文的研究对象为翅片管内流动沸腾过程中气泡行为和两相流型，分析不同LNG入口流速、LNG入口温度以及管外表面传热系数等条件下，两相流型的转换特性。

2.2 假定说明

① 管道导热为各向同性。

② 管道出口天然气体积分数不随时间改变时，认为管内流动沸腾传热过程发展到稳态。

③ 管外空气表面传热系数沿管长方向为定值。

④ 管外空气温度为定值。

3 几何模型描述

本文研究的重点为翅片管内流动沸腾过程中气泡行为和两相流型的变化，管外翅片对于管内流动沸腾过程的影响仅仅在于传热问题。因此，翅片所带来的传热影响将通过等效的方式作用于管道外壁面，即将翅片管外的传热系数折算为竖直光管面积所对应的传热系数，从而将翅片管外壁处理为无翅片的圆光管，圆光管的直径与翅片管外直径相同。为了实现沸腾与固体基质的传热耦合，保留了管壁的厚度。

由于竖直光管具有对称性，因此本文建立的几何模型为二维模型，见图1。LNG从管道的底部进入，受热气化后从顶部流出。几何坐标轴的原点O在管道底部入口的中心处，x轴沿管道的径向方向，y轴与管道中心轴线重合，方向与LNG的流动方向一致。管道的内直径为0.02 m，外直径为0.026 m，管道的长度为5 m。

图1 几何模型

4 网格划分

本文将计算域分为LNG流体域和铝合金管道固体域，计算域网格划分见图2。首先在ICEM 19.0软件中建立几何模型，并创建块文件，通过两次切分，在几何模型与块之间形成映射，最终生成结构化四边形网格。对铝合金管内外壁面附近法线方向上的LNG流体区域进行边界层加密。经过网格无关性检验后确定，在管内侧采用BiGeometric分布规律，节点数为131，第1层的高度为1.0×10-4m，增比Ratio=1.02；铝合金管壁内的网格采用Uniform节点平均分布，节点数为5，每1层的高度为7.5×10-4m。沿铝合金管高度方向，网格高度相同，每层高度为2.5×10-3m。最终所得的网格总数为1 379 655。

图2 计算域网格划分

5 各项设置

采用ANSYS Fluent 19.0软件进行计算，求解器设定为：Double precision、Parallel processing(12 processes)、Pressure-based、Absolute、Unsteady、2D。

模型设置：求解质量、动量以及能量控制方程，多相流模型选择VOF模型。湍流模型为标准k-ε湍流模型、近壁处为标准壁面函数。沸腾相变模型选择S. Hardt源项模型，通过用户自定义(UDF)的方式将S. Hardt源项模型加入控制方程中。选用分段线性构造界面(Piecewise Linear Interface Construction)的几何重构法捕捉气液两相界面。

材料设置：LNG和天然气的物性数据采用文献[10]中的方法，分别对气液两相的密度、比热容、动力黏度和热导率进行计算，同样进行多项式拟合，在Fluent中创建LNG和天然气两种物质，其中LNG中各烷烃组分的摩尔分数分别为，甲烷88%，乙烷8%，丙烷4%。固体域铝合金的物性参数采用默认设置，其密度为2 623.3 kg/m3, 比热容为983 J/(kg·K),热导率为227.95 W/(m·K)。

操作条件设置：打开重力场，沿管道轴线方向，即y轴方向，设置重力加速度为-9.81 m/s2。

边界条件设置：对于铝合金管内，入口处设置为速度入口(velocity inlet)，LNG的速度在0.05～0.20 m/s范围内，质量分数为1，入口温度为111～126 K。翅片管出口采用压力出口边界(pressure outlet)，出口处压力设置为0.5 MPa，液相的回流参数设置为0。铝合金管外壁设置为wall边界条件，其表面传热系数变化范围为100～300 W·m-2·K-1，来流环境温度设定为300 K。管内壁面会被Fluent软件自动识别流固耦合传热，在读入网格时通常为流固界面生成一个对应的shadow面，并将这两个面在传热问题上确定为耦合(coupled)关系。

计算工况的设计见表1。

表1 计算工况

续表1

6 求解

各方程的离散格式：Pressure-Velocity Coupling Scheme采用PISO，Gradient选择Least Squares Cell Based，Pressure选择PRESTO，其余方程均采用Second Order Upwind格式进行离散。

欠松弛因子设定：为保证计算过程的稳定性，在计算前期先调低所有欠松弛因子，待计算稳定后再将其设定default，即Pressure、Density、Body Forces、Momentum、Vaporization Mass、Slip Velocity、Volume Fraction、Turbulent Kinetic Energy、Turbulent Dissipation Rate、Turbulent Viscosity、Energy的取值分别为0.3、1、1、0.7、1、0.1、0.5、0.8、0.8、1、0.91。

残差设定：Energy残差为1×10-6，其余为1×10-3。

初始化：先选用Standard Initialization，Relative to Cell Zone方式，Gauge Pressure、X velocity、Y velocity、Z velocity均为0，Turbulent Kinetic Energy为8.437 5×10-5，Turbulent Dissipation Rate为1.657 8×10-3。由于本文采用瞬态模拟计算，需要确保初始化的结果尽可能贴近实际物理现象，因此，通过Patch功能分别设置管内计算域的温度为300 K，压力为5×105Pa，天然气的体积分数为1；固体铝合金管的温度为300 K。

本文研究的内容是从LNG进入铝合金竖直翅片管开始气化到管内的气化状态达到稳定的过程，因此采用瞬态模拟计算，时间步长采用自适应的方式，将步长控制在1×10-5～1×10-4s。每个时间步长内最大迭代次数为200。Fluent实时监测管道出口天然气的体积分数变化，当体积分数变化达到稳定时，即可认为此时管内流动沸腾达到稳定状态。

7 计算结果及分析

7.1 管内流场特性

本文定义工况0：LNG的入口流速为0.15 m/s，入口温度为126 K，管外表面传热系数为300 W/(m2·K)。图3～6展示了工况0条件下，竖直翅片管管内LNG气化过程达到稳定状态后，不同高度区域内气相体积分数云图和速度矢量分布(软件截图)，其中，色标右侧的标值为气相体积分数，速度矢量用黑色箭头线表示。沿管道竖直向上(y轴正方向)，管内不同位置所表现出的流场特性有明显差异。根据天然气气泡的大小、位置、形状以及密集程度，可以划分不同的两相流型，例如泡状流、弹状流以及搅拌状流等。从图3～6可以看出，当LNG刚进入管道后，管内主流区域的速度场分布较为均匀，但是由于受到热壁面的加热，管壁附近有零星的气泡产生，气泡处的流动出现轻微紊乱。随着LNG向上流动，受到管壁持续加热，管内气泡数量增加，气泡直径增大。在气泡存在的区域，气相的速度比液相的速度大，气液两相之间有明显的速度差，液相受到气相的扰动更强。在图5和图6中，一些两相流动区域出现了局部回流的现象。

图3 工况0稳定状态时管内y=[0.015 m, 0.047 m] 区域内气相体积分数云图和速度矢量分布(软件截图)

图4 工况0稳定状态时管内y=[0.145 m, 0.177 m] 区域内气相体积分数云图和速度矢量分布(软件截图)

图5 工况0稳定状态时管内y=[0.385 m, 0.417 m] 区域内气相体积分数云图和速度矢量分布(软件截图)

图6 工况0稳定状态时管内y=[0.725 m, 0.757 m] 区域内气相体积分数云图和速度矢量分布(软件截图)

7.2 气泡行为

气泡的生长过程受到多个力的共同作用，主要有体积力、表面张力、来自气液界面的黏性应力、气泡内的蒸气压力以及热壁面对气泡的反作用力。起初，气化核心的半径较小，尽管热壁面的液膜温度较高，但是气泡内的压力大于液相压力，其对应的饱和温度比较高。此时液膜与气泡之间的温差较小，将其视为等温生长过程，在这个阶段气泡生长的推动力主要是气泡内部的蒸气压力与液相作用在气液界面上的应力之和减去气液之间表面张力。当气泡直径逐渐增大时，其受到的体积力迅速增大，气泡内的蒸气压力逐渐减小，温度随之降低，液膜与气泡之间的温差逐渐变大，此后气泡的生长过程便受液膜与气泡之间的传热过程主导。

图7～9为工况0条件下，竖直翅片管管内高度为y=[0.075 m, 0.210 m]区域内热壁面上气泡间的聚合现象。其中，色标右侧的标值为气相体积分数。为了便于观察，使用白色方框框选图中选定的气泡，使用粉红色箭头指示气泡行为的变化历程。

图7 工况0条件下热壁面上气泡间的聚合行为 (软件截图)

随着气泡本身的不断生长，当气泡直径增大到一定范围时，气泡之间就会出现聚合现象，将这种沿管长方向的聚合称为纵向聚合。从图7中可以看出，气泡的聚合行为最先出现在热壁面附近，上游气化核心处产生的气泡多数并未直接脱离，而是受到浮力和液相的拖曳力，沿着壁面向上滑移，同时由于持续受热，其直径不断增大，相邻两个气泡之间的距离减小，从而逐渐发生融合。

从图8可以看到，在时间为0.65 s时，两个气泡间的液膜已经开始逐渐融合，最终在时间为0.67 s时，两个气泡在靠近管道中心轴的主流区实现了完全融合，成为体积更大的气泡。在气泡脱离热壁面进入主流区以后，由于浮力和拖曳力的作用，气泡上升速度不断加快，但是对于处在上下相邻位置的两个气泡而言，下方的气泡会受到上方气泡的尾流作用，其上升加速度大于上方气泡的加速度，使得两个气泡之间的距离不断减小。

图8 工况0条件下主流区气泡纵向聚合(软件截图)

在气泡重点聚合区，其数量逐渐增多导致间距变小，在主流区出现了气泡间的横向聚合，见图9。在时间为1.23 s时，图9中标出的两个气泡基本处于同一水平高度，且两者之间尚存在较厚的液膜，但是在时间为1.26 s时，两个气泡之间的距离缩短，逐渐发生聚合行为。值得注意的是，与热壁面处气泡的上升速度相比，图9所示的主流区的气泡上升速度较小，原因是热壁面附近气相体积分数较大，导致混合相流速高于主流区。

图9 工况0条件下主流区气泡横向聚合(软件截图)

在主流区，气泡失去热壁面的横向束缚，常以球状或者椭球状向上运动。上升气泡的形态和大小与流体本身的物性有较大的关系，主要受到表面张力和上升浮力的影响。奥托斯数[11]正是表征气泡上升过程中表面张力和浮力的相对影响，其计算公式见式(1)。

(1)

式中Eo——奥托斯数

de——与气泡体积相等的圆球的直径，m

g——重力加速度，m/s2

ρl——液相密度，kg/m3

ρg——气相密度，kg/m3

σ——表面张力系数，N/m

当离散相为空气，连续相为水时，Eo临界值为40，Eo低于临界值时气泡主要呈现球状或者椭球状，高于临界值时则由帽状转变为不规则的形状[11]。由于Eo临界值是混合相物性参数的函数，因此，本文根据空气-水混合相的Eo临界值，以及LNG-天然气混合相的密度和表面张力系数，推算出LNG-天然气混合相的Eo临界值为7 443。天然气气泡的体积在上升过程中不断增大，其等效直径增大后，从式(1)中可以看出，等式右侧分子部分所代表的浮力将逐渐增大，并开始占据主导地位。当Eo高于7 443时，气泡形状随之发生相应变化，不再呈现规则的球状或椭球状。从图9可以看出，气泡的体积越大，其形状越扭曲。此外，从图4可以发现，气泡在靠近热壁面侧和主流侧的速度相差较大，气泡紧贴热壁面侧的速度大于主流区的速度，因此，形成了拖曳的效果，导致图4和图9中热壁面上的大体积气泡多呈现“半个箭头”形状。

7.3 气液两相流型

本文模拟结果表明，LNG在管内的流动沸腾过程中出现的流型有泡状流、弹状流、搅拌流、雾状流，并未出现环状流。在对流型的划分中，本文部分借鉴了文献[5]的划分方法，具体的划分方法见表2。与文献[5]不同的是，将气相体积分数为0～0.18阶段分为了单液相流和泡状流两种流型。对于单液相流，主流温度尽管未达到饱和温度，但由于存在过冷沸腾的原因，依然会有气泡产生，极少量的小气泡脱离壁面后，受到过冷液体的冷却而快速湮灭，因此，将其作为两相流区的特殊阶段。

表2 LNG竖直翅片管管内气液两相流型划分

从LNG进入管道一直到管内完成气化，即气相天然气的体积分数从0达到1，这个区域统称为两相流区域。各流型在管内分布长度的统计依据是流型的特征和对应的体积分数。各流型在管内分布长度与两相流区域总长度的比值，称为两相流中各流型占比。在上述众多流型中，搅拌流所对应的流动沸腾传热系数较大[5]，故而对不同LNG入口流速、LNG入口温度以及管外表面传热系数等条件下的管内流动沸腾进行模拟研究，目的在于扩大搅拌流区域，增大流动沸腾传热系数，提高空温式气化器的利用效率。

7.3.1 LNG入口温度对两相流型分布的影响

LNG一般经泵加压后进入空温式气化器，进入空温式气化器的LNG并不会即刻蒸发沸腾。若进口LNG的温度太低，将使得气化器管长增加，导致设备体积过大。图10为竖直翅片管管内流动沸腾达到稳定状态后LNG入口温度对两相流型分布的影响，其计算条件是，LNG的入口流速为0.15 m/s，管外表面传热系数为300 W/(m2·K)，入口温度为111～126 K。

图10 LNG入口温度对两相流型分布的影响

从图10可以看出，在不同LNG入口温度的条件下，管内两相流中各流型占比不同。当空温式气化器的入口处温度为111 K时，管内单液相流区域较长，占气液两相流区域的31.4%；其次是泡状流区域，其占比为18.1%。单液相流和泡状流的区域均随着LNG入口温度的升高而减小，其中单液相流区域的变化幅度最大，当温度升高到126 K时，单液相流区域的占比下降到8.2%。这是因为当入口LNG温度较低时，低温的LNG需要在管内流动更长的距离、吸收更多的热量，才能进入泡状流区域。对于泡状流区域而言，在管壁和管道中轴线之间存在较大的温度梯度，靠近热壁面处气泡密度较大，中轴线附近的气泡密度小，由此使得泡状流区域受到一定程度的拉长。随着LNG入口温度的增高，搅拌流区域的占比获得了较大的提高。当入口温度为126 K时，搅拌流和雾状流区域占比的总和达到了70.6%，表明LNG的入口温度越高，管内的流动沸腾换热系数越大，换热效果越好。

7.3.2 LNG入口流速对两相流型分布的影响

随着空温式气化器负荷的变化，入口处的LNG流速会随之发生变化。当LNG的入口流速为0.05 m/s，对应的天然气出口流量(折合成标准状态)为34 m3/h。图11为竖直翅片管管内流动沸腾达到稳定状态后LNG入口流速对两相流型分布的影响，其计算条件是，LNG的入口流速为0.05～0.20 m/s，入口温度为126 K，管外表面传热系数为300 W/(m2·K)。

图11 LNG入口流速对两相流型分布的影响

由图11可见，随着LNG入口流速的增大，单液相流、搅拌流、雾状流在两相流中的占比逐渐增大；泡状流与弹状流区域则随着入口流速的增大而减小，且弹状流区域减小的幅度最大。当入口流速为0.05 m/s时，单液相流段较小，在两相流中的占比为3.8%。这是由于较小的入口流速，意味着管内的质量流量较小，主流区的LNG经过较短的管程距离便达到了泡状流阶段。当LNG的入口流速从0.05 m/s增大到0.20 m/s时，泡状流区域和弹状流区域的占比(绝对值)分别减少了9.5%和22.9%。这是因为随着入口流速的增加，气液两相之间的相对速度不断增大，气液界面之间的黏性应力增加，使得气泡尺寸被拉长，且管内气液两相的流速增大后，两相流的湍流强度增大，气泡之间更容易挤压碰撞，由此导致管内的泡状流和弹状流区域不断缩小。与雾状流相比，搅拌流的占比仅从28.8%增加到33.7%，增加幅度较小。原因主要为随着气相体积分数的增大，气相的流速剧烈升高，较高的气流速度使得雾状流区域中的液相体积分数升高，增大了雾状流区域在管内的分布长度。

7.3.3 管外表面传热系数对两相流型分布的影响

由于环境空气变化以及管外结霜的原因，管外的表面传热系数会有明显的变化。图12为竖直翅片管管内流动沸腾达到稳定状态后管外表面传热系数对两相流型分布的影响，其计算条件是，LNG的入口流速为0.15 m/s，入口温度为126 K，管外表面传热系数为100～300 W/(m2·K)。

图12 管外表面传热系数对两相流型分布的影响

由图12可见，随着管外表面传热系数的提高，单液相流和泡状流在两相流中的占比快速下降，搅拌流和雾状流的占比则是迅速上升，而弹状流占比的变化较小。当管外的表面传热系数为100 W/(m2·K)时，单液相流在气液两相流中的占比为30.2%，泡状流的占比为33.4%，远高于弹状流和搅拌流的占比。这是因为，管外表面传热系数较低时，在管内壁处的热流密度较小，使得液相中气泡的直径和密度均比较小，气泡之间融合的概率就比较低，气泡之间基本呈现彼此孤立的状态。当管外表面传热系数从100 W/(m2·K)增大到300 W/(m2·K)时，热壁面上气泡受到较高的热流密度加热，无论是气泡直径还是脱离热壁面的频率都会增大，使得搅拌流的占比从8.6%增加到32.3%。

鉴于搅拌流所对应的流动沸腾传热系数较大，增大搅拌流区域，有利于提高管内的流动沸腾传热系数，进而提高空温式气化器的利用效率。对比图10～12中的数据可知，当入口温度从111 K增大到126 K时，搅拌流区域占比增大了14.3%；入口速度从0.05 m/s增大到0.20 m/s时，搅拌流区域占比仅增大了4.8%；而改善管外的表面传热系数使其从100 W/(m2·K)增大到300 W/(m2·K)时，则可使搅拌流区域占比增大23.7%。因此，与增大LNG入口温度、LNG入口流速相比，通过改善管外的表面传热强度更能够提高空温式气化器的换热效率。

8 结论

为了揭示LNG空温式气化器管内流动沸腾特征，对沸腾过程中气泡行为和两相流型的分布和转换进行了研究。研究对象为竖直铝合金翅片管，并将其等效简化为竖直光管，其高度为5 m，管内直径为0.02 m，管外直径为0.026 m。管内LNG的入口流速为0.05～0.2 m/s，入口温度为111～126 K，管外表面传热系数变化范围为100～300 W/(m2·K)，环境温度为300 K，环境压力为101.325 kPa。采用ICEM 19.0软件建立竖直光管的二维模型并划分四边形结构化网格，利用Fluent 19.0软件对管内液化天然气或天然气计算域和管壁固体计算域进行传热模拟，研究不同LNG入口流速、LNG入口温度以及管外表面传热系数等条件对沸腾过程中气泡行为和两相流型分布和转换的影响。采用流体体积法对气液两相界面追踪捕捉，利用用户自定义的方式将表征沸腾过程的S. Hardt模型引入计算流程中。研究结果表明：

① 气泡的聚合及沿热壁面的滑移是在表面张力、浮力以及气液之间剪切应力作用下的重要运动形式。

② LNG入口温度增加可以减小单液相流在两相流中的比例，同时增大搅拌流和雾状流区域；LNG入口流速增大时，两相流中泡状流和弹状流区域快速减少，搅拌流区域小幅增加；管外表面传热系数增大，可显著提高搅拌流的占比。

③ 与增大LNG入口温度、LNG入口流速相比，通过改善管外的表面传热强度更能够提高空温式气化器的换热效率。