重沸器内三甘醇-水两相流动与沸腾换热特性研究

2021-07-14沈志恒翟占虎

节能技术 2021年3期

李巍，沈志恒，董超，吴尧，翟占虎

(1.海洋石油工程股份有限公司，天津 300452; 2.哈尔滨工业大学能源科学与工程学院，黑龙江哈尔滨 150001)

0 引言

海洋气田中心平台在运过程需要使用三甘醇(TEG)脱水装置对天然气进行干燥，对于脱水装置中重沸器内三甘醇-水的沸腾过程进行数值模拟计算，对于准确预测混合溶液中的三甘醇浓度、确定合理运行温度区间、防止局部温度过高等具有重要意义[1-2]。传统对于三甘醇脱水过程的数值模拟计算大多采用ASPEN等软件进行工艺分析，获得不同流程上的状态参数，该方法对于如重沸器这样的设备的局部温度特性的预测尚无法满足精度要求。

李天斌[3]采用HYSYS软件对某海洋气田中心平台正在运行的TEG脱水装置进行模拟计算，对贫TEG循环量及其质量分数、再沸器温度、汽提气流量等参数进行模拟优化，获得了推荐的优化参数。张明震等[4]以某油田TEG脱水装置为基础，利用ASPEN HYSYS建立TEG脱水装置模拟流程，对TEG脱水装置运行参数进行了优化对比并提出了TEG脱水装置的最优运行参数。杨延明[5]在ASPEN PLUS软件中将离子液体作为脱水溶剂并对其天然气脱水过程进行模拟计算。研究发现增大TEG贫液再生温度、循环量和吸收塔塔板数可使干气水含量减小，同时使贫液再生热负荷以及TEG的损失量增大，并认为离子液体是TEG的潜在替代者。Ahmad等[6]使用前馈人工神经网络(FANN)预测TEG脱水过程中天然气的平衡水露点，研究表明FANN可以很好地预测TEG脱水过程中天然气的水露点。Tatar等[7]利用智能建模技术根据流体的TEG浓度和温度预测天然气流中的平衡水露点。林志军[8]计算了天然气TEG脱水时TEG对应的定性温度、重沸器工作时对应的负荷热及其传热面积，基于此设计了重沸器火管和壳体尺寸并进行校核。可以看出，目前国内外对于TEG脱水过程的数值模拟侧重于采用ASPEN等软件进行工艺上的计算，对于实际脱水过程中TEG溶液内部的流动、传热和蒸发过程未见详细研究。此外，对于重沸器内部加热管与流体的流动换热特性也未见相关报道。

本文基于VOF方法，考虑加热过程中流体的相变过程对温度的影响，通过数值模拟计算获得了重沸器中的流场和温度场分布，并提出了一种用于评估加热器表面温度分布的方法。

1 计算模型

VOF方法通过两种或多种流体(或相)没有互相穿插这一特点，对于增加到模型里的每一附加相，就引进一个计算单元里的相的容积比率。在每个控制容积内，所有相的体积分数和为1。基于某一项的局部值，可以适当的将属性和变量在一定范围内分配给每一控制容积。

1.1 容积比率方程

跟踪相之间的界面是通过求解一相或多相的容积比率的连续方程来完成的。对第q相，这个方程如下

(1)

对于默认情形，上式右端的源项为零，除非给每一相指定常数或用户定义的质量源项。容积比率方程不是为主相求解的，主相容积比率的计算基于如下的约束

(2)

1.2 动量方程

通过求解整个区域内单一的动量方程，作为结果的速度场是由各相共享的。如下所示，动量方程取决于通过密度ρ和粘度μ的所有相的容积比率

(3)

1.3 能量方程

能量方程表示如下

(4)

VOF模型处理能量E和温度T是作为质量平均变量来处理的

(5)

这里对每一相的Eq是基于该相的比热和共享温度来确定。密度ρ和keff(有效热传导系数)被各项共享。源项Sh包含辐射的贡献，也有其他容积热源。

1.4 湍流模型

湍流模型采用RNGk-ε，该模型通过大尺度运动和修正后的粘度项体现小尺度的影响，从而将小尺度运动系统地从控制方程中去除。此外，RNGk-ε湍流模型通过在ε方程中增加一个反映主流的时均变率项，考虑了平均流动中的旋流进而修正湍流粘度。此模型表示如下

(6)

(7)

式中k——湍动能;

t——时间;

ui——时均速度;

αk，C1ε，C2ε——模型常数;

Gk——由平均速度梯度引起的湍动能k的产生项;

ε——湍动能耗散率。

1.5 界面传热模型

本文基于Nu数来进行流体中对流换热系数的定义，传热系数hfg的值可以与Nu数相关联

(8)

式中κf——流体的热导率;

dg——分散相的直径。为了确定公式(8)中的Nu数的大小，这里引入Ranz-Marshall[9-10]关系式。

(9)

式中Ref——基于分散相直径与两相间的滑移速度而定义的雷诺数;

Pr——主相的普朗特数。

1.6 界面传质模型

依据Hertz-Knudsen公式，可以得到基于分子动力学的界面上的相变流量

(10)

式中p——温度为T时可凝结气相的分压;

psat——温度为T时的饱和压力;

R——通用气体常数。考虑到Clapeyron-Clausius方程，在饱和状态附近压力可以与温度关联起来

(11)

式中L——工质的潜热;

vv和vl——气相和液相的比容;

γ1——单位体积内界面蒸发调节系数，表征界面蒸发的强度大小

(12)

利用相似的思路，同样也能得到对于冷凝的表达形式，因此界面上的传质情况可以写成如下的形式

(13)

2 计算模型

2.1 几何模型

重沸器半径为550 mm，长度约为3 700 mm，入口直径约为400 mm，重沸器内液面高度为730 mm。重沸器加热棒直径16 mm，浸入长度为2 137 mm，共30根加热棒，分为5排布置，第一排有5根，最后一排仅有1根加热棒，其中有3根作为备用管，加热棒中心距为24 mm，不考虑用于固定加热棒的挡板，模型简化后如图1(a)所示。

图1 重沸器几何模型及网格划分

考虑到重沸器内部结构复杂，且筒体径向特征尺寸与加热棒径向特征尺寸相差2个数量级，需要使用较多数量的网格才能获取足够精度的数值结果。为了降低计算资源的消耗，本文采用多面体网格，进行网格无关性检验后最终确定的数量为127万，质量在0.25以上，最终的网格如图1(b)所示。

2.2 计算条件

对于重沸器，操作压力为10 kPa，操作温度为377.16 K。采用基于压力的simple算法实施瞬态计算，过程中考察加热棒表面温度分布，当其分布不随时间变化时，认为其温度分布已经能代表真实情况下的膜温度分布情况。

相间作用中质量传输选用蒸发-冷凝模型，相关参数按照冷凝器的冷凝模型计算获得，具体参数见表1。

表1 相关物性

标准工况中，重沸器入口设置质量入口，温度170 ℃，TEG的质量流量为0.889 kg/s，蒸汽的质量流量为0 kg/s；出口设置为压力出口，回流温度204 ℃；重沸器筒壁设为绝热壁面；加热棒壁面设为定热流密度条件，热流密度取20 000 W/m2。

3 计算结果及讨论

3.1 时间无关性检验

再沸器内流体相变和流动是一个动态过程。因此计算初始的一段时间内，不合理的初始条件可能对数据造成一定的影响。因此需要进行时间无关性检验。

从计算的第10 s开始监测壁面温度的分布，并将不同时刻壁面温度在207～215 ℃范围内的壁面占比提取出来，绘制在图2中，可见计算开始10～11 s时，207～215 ℃的各个壁面温度区间的概率有明显上升，说明此时的沸腾还受到较多的初始情况的影响。而13～15 s时刻，壁面温度的概率分布趋于一致，单调递增的趋势已经消失，说明此时加热器周围的沸腾已经进入稳定的周期过程，其温度分布已经具有代表性。图3是壁面最高温度随时间上的散点图，其在16 s后变化不大；为保证分析的准确性,本文中后续的分析均是基于15 s之后的数值结果展开的。

图2 不同时刻温度温度分布曲线

图3 加热棒最大温度分布散点图

3.1 重沸器内流场与加热面温度分布

图4是加热棒壁面温度分布稳定后的重沸器内流体速度分布，可以看到速度最大的位置集中在加热棒束的正上方，这是因为相变产生的气泡在重力作用下上浮携带流体运动导致的。入口处呈现深蓝色，这是由于循环流量变化范围仅为2～4 m3/h，结合重沸器液相容积换算时均流速仅为1～2 m3/h，即使是入口速度也仅为0.01 m/s的数量级；远远小于由于气泡上升携带液体和径向热对流引起的速度(0.3～0.5 m/s)。

图4 重沸器内速度分布

图5所示为重沸器内部径向流体迹线。可以观察到，在加热棒束段，流动以径向对流为主；而在重沸器远离加热棒束的一端，流动以沿轴向的迁移运动为主。结合图2和图3，管束核心区及其上部以外的区域的流体运动速度要比棒束区低1～2个数量级，即重沸器内的流动以径向对流为主导。

图5 重沸器内径向不同位置流体迹线图

图6是监测温度为208 ℃、212 ℃和215 ℃时的加热器表面温度分布，可见热流密度为20 000 W/m2时：(1)加热器表面能够观察到大量区域温度高于208 ℃；(2)温度高于212 ℃时的区域也可以被明显观察到；(3)而监测温度提高到215 ℃时，已经基本观察不到超温区域。结合图2中的温度分布曲线，可以认为超过相应温度的壁面面积占比低于3%时，即观察不到明显的超温区域。

图6 不同计算时刻加热棒表面温度分布

3.2 重沸器壁面温度与热流密度的关系

设置热流密度为8 000 W/m2、12 000 W/m2、16 000 W/m2、18 000 W/m2和20 000 W/m2，其它设置均和标准算例相同，提取各个算例207～215 ℃对应的温度分布概率，并将概率以对数坐标绘制在图7中。可以观察到随着热流密度逐渐增高，相同温度的概率也逐渐增加；例如热流密度从8 000增加到12 000 W/m2时，其概率分布增加了约4倍。

使用最小二乘法对图7中的数据进行关联，得