高压空冷器入口管束内流动参数分布特性的数值模拟

偶国富， 徐晓峰， 吕文超， 金浩哲， 张 林

(1.浙江理工大学 流动腐蚀研究所， 浙江 杭州 310018； 2.中国石化 武汉分公司 设备工程处，湖北 武汉 430082)

流动腐蚀是石油化工企业在役压力管道及换热设备失效的重要原因，由于服役工况及输送多相流介质的差异，不同石油化工设备的流动腐蚀机理也各不相同[1-3]。在高硫、高酸等劣质原油的加工过程中，高压空冷器(压力≥10 MPa)管束中的流动腐蚀失效形式主要为铵盐沉积及多相流冲蚀[4-6]，由此引起的堵塞爆管、穿孔泄漏等事故不断增加，经济损失很大[7-8]。因此，研究高压空冷器内部的多相流动特性，建立以流动参数为主的预测方法显得尤为重要。

笔者在计算空冷管束内部流场分布的基础上，研究以水相分率、传质系数与壁面剪切应力为核心表征参数的流动腐蚀高风险区域预测方法，从而为石油化工企业高压空冷器的耐流动腐蚀优化设计及风险检验提供了技术支撑。

1 数学模型

1.1 物理模型

图1为高压空冷器的几何模型。该空冷器采用对称布置方式，管箱入口前段为三通管道，使多相流在进入管箱前能够充分混合。在前处理过程中，对空冷器全流场进行拓扑，采用六面体结构化网格进行划分，并对结构突变区域进行局部加密。为保证管束近壁面流动参数计算结果的准确性，设置边界层网格，根据湍流模型选择第一层网格质心到壁面的无量纲距离30 ≤y+≤ 60[15]，计算得到壁面到第一层网格单元中心点的距离Δy=0.0002 m时符合要求。对空冷器内部流场进行了网格无关性验证，发现在4.04×106、6.02×106、8.01×106网格数量下，空冷器管束流动参数十分接近，管束内最高流速分别为6.11、6.15、6.18 m/s，相对误差≤0.57%，认为已达到网格无关性要求，确定采用网格总数4.04×106进行数值计算。

图1 空冷器结构及网格划分Fig.1 Structure of air-cooler and grid generation

1.2 控制方程

高压空冷器管束中的介质为气、油和水三相流，在空冷器入口管段温度较高，多相流相分率基本不变，因此不考虑介质在空冷器入口管段流动过程中的传热及相变问题。由于多相流的相间耦合作用较强，且存在相间相对运动，故选用Mixture模型模拟管束内多相流动过程。

根据质量守恒定律，建立多相流混合相的连续方程：

(1)

根据动量守恒定律，建立动量方程：

(2)

为了研究剪切应力的分布规律，需要对近壁面的流动进行更精确的计算，SSTk-ω湍流模型有效考虑了流线曲率及逆压梯度效应的影响，可以求解近壁区的低雷诺数流动，因此被用来对流场内部物理量进行求解。其中k与ω由相应输运微分方程确定[16-18]。k由式(3)确定：

Gk-Yk+Sk

(3)

ω由式(4)确定：

Gω-Yω+Dω+Sω

(4)

式(3)、(4)中，k为湍动能，J；ω为湍流比耗散率；x为坐标矢量(i,j=1,2,3,分别表示x、y、z3个空间坐标)；u为速度矢量，m/s；Gk为由平均速度梯度产生的湍动能项；Gω为对应的湍动能比耗散率项；Yk、Yω分别为由湍流引起的k和ω的耗散项；Sk、Sω为自定义源项；Dω为交错扩散项；Гk、Гω分别为k和ω的有效扩散系数项，其计算公式为：

(5)

(6)

(7)

式(5)～(7)中，S为应变速率；常数α*、α1分别取1.0、0.31；混合函数F计算公式如下：

F=tanhΦ2

(8)

(9)

式(8)、(9)中，tanh()为双曲正切函数；Ф为近壁面处的流量变量，m3/h。

在空冷器入口管段，铵盐会溶于液态水形成强腐蚀性溶液，溶液中游离的H+在对流传质作用下，会穿过腐蚀产物膜与金属基体发生氧化还原反应，造成电化学腐蚀。为了反映H+在腐蚀产物膜表面的对流传质速率，采用H+的传质系数作为流动腐蚀的表征参数[19]。根据柯尔邦类比，可得：

(10)

解得

Sh=0.023Re0.8Sc-1/3

(11)

(12)

式中，DAB为溶液中H+的扩散系数，取DAB=9.31×10-9m2/s；L为管束内径，m。

1.3 边界条件

空冷器计算域的进、出口边界条件分别采用速度进口和自由流动出口。采用有限体积法实现计算区域和控制方程的离散，体积相分率、湍动能、动量和比耗散率均采用一阶迎风格式进行离散，压力项采用Standard格式，压力-速度耦合方程的求解采用SIMPLE方法，壁面按固壁无滑移壁面处理。通过Aspen计算得到的空冷器入口多相流物性参数，如表1所示。

表1 空冷入口配管物性参数Table 1 Physical characteristics of fluids in the air-cooler

Q—Volume flow rate;ρ—Density;μ—Viscosity;φ—Volume fraction;v—Velocity

2 结果与讨论

2.1 空冷器偏流分析

图2为空冷器管束位号划分图。由图2可知，介质经入口配管进入空冷器，由三通管道分为两股物料经交界面k1、k2进入管箱，再由管箱流入空冷管束。空冷管束规格为φ25 mm×3 mm，材料为碳钢，其中管箱出口管段设置长度为300 mm的衬管。空冷器第一管程共有两排管束，按从上至下的顺序命名为a管排和b管排，第一管程共有93根管束，其中a排47根，b排46，从左至右分别命名为a1-47、b1-46。其中空冷管箱入口位于a10、a38管束正对位置。

图2 空冷管束位号划分示意图Fig.2 Order division of air-cooler pipes

图3为空冷器内部流场的速度矢量图。从图3可知，由于入口配管相对于空冷器整体为非对称分布，因此介质在流经三通管道时，会发生流体偏流现象，k1、k2处的流动参数见表2。图4为不同时刻管束出口质量流量分布图，结合图4可以看出，a管排的整体质量流量大于b管排，其中同一管排中k1一侧空冷管束中的质量流量大于k2一侧，质量流量峰值均出现在管箱中间及两侧对应位置。图5 为空冷入口配管及管箱水相分率分布图。由图5 可见，由于进口多相流流速较小，各相分层明显，密度较大的水相主要集中在管道底部，而气相则积聚在顶部。介质在流经三通管时，水相偏向于左侧k1面方向，因此k1面水相分率较大，导致流量减小，整体流速小于k2面。

图3 空冷器内部速度矢量图Fig.3 Diagram of velocity vector in the air-cooler

图4 不同时刻管束出口质量流量(Qm)分布Fig.4 Mass flow rate (Qm) distribution of pipe outlets(a) Pipeline a; (b) Pipeline b

图5 空冷入口配管及管箱水相分率(φw)分布Fig.5 Distribution of water volume fraction (φw) in the air-cooler inlet pipes and pipe box

2.2 空冷管束流动腐蚀高风险区域预测

图6为不同时刻空冷入口管束内湍流强度分布图，结合图4可以看出，流动时间为16 s和22 s时的出口质量流量最大误差小于0.67%，可以认为当流动时间大于16 s之后，空冷管束内多相流特性参数保持不变，此时计算结果趋于稳定。其中b管排管束内部流体的湍流强度明显大于a管排管束，整体流动更为复杂，流体与壁面之间的剪切应力也相应较大，得到的管束剪切应力分布如图7所示。

由图7可知，a、b两管排管束的剪切应力分布趋势基本相反，a管排的峰值出现在入口管两侧及管箱中线，对应位号为a8、a13、a22、a35、a40管束位置；b管排的对应位置则均出现谷值，b管排的峰值出现在b4、b10、b19、b24、b31、b37、b46位置，对应a管排的相应管束基本处于谷值。在壁面剪切应力的峰值区域，腐蚀产物保护膜更容易被冲刷破坏，导致内部金属基体重新裸露在腐蚀性流体中，形成自催化体系，进一步加速腐蚀[20]。

表2 不同交界面的流动特性参数Table 2 Flow characteristic parameters of different interfaces

Qm—Mass flow rate;p—Static pressure;ρ—Density;v—Velocity;φw—Water volume fraction

图6 不同时刻空冷管束湍流强度(I)分布Fig.6 Distribution of turbulent intensity (I) in the air-cooler pipes(a) Pipeline a; (b) Pipeline b

从表2中的数据可以看出，在三通管处存在明显的偏流，导致了空冷管束内部流场也存在非对称的现象。图8为空冷管束水相分率的分布图，结合表2可以看出，随着管束编号的增大，水相分率整体呈减小趋势，a管排管束峰值出现在a2、a13、a23、a35、a45区域，b管排管束峰值出现在b4、b20、b30、b46区域，且a管排管束内水相分率最大值大于b管排，这是由于介质进入管箱时，流体冲击管箱底部后流向两侧，在入口两侧区域形成涡流，将大部分水相带至管箱顶部所导致的，而在2个入口对应偏管箱中心区域，为涡流交汇区域，该处a、b管排管束内水相分率分布基本呈完全相反的趋势，在中心漩涡交汇位置，两排管束的水相分率差值达到最大，其中a29及b23处管束内水相分率达到最小值。

图7 空冷管束剪切应力(τ)分布Fig.7 Distribution of wall shear stress (τ)in the air-cooler pipes

图8 空冷管束水相分率(φw)分布Fig.8 Distribution of water volume fraction (φw) in the air-cooler pipes

图9为沿流动方向从衬管出口到450 mm处4个不同截面上的传质系数分布图。对图9分析可知，随着流出衬管出口距离的增加，空冷入口管束内的传质系数整体变化趋势未发生明显改变，但数值逐渐减小。a管排管束内传质系数均大于相邻位置的b管排管束，且同排管束中传质较强区域集中在空冷管箱入口两侧及管箱中线正对位置，对应位号为a12～a13、a24、a35～a36管束，与水相分率较大位置基本相同。在此区域的管道内因液相水较多，易溶解管束中结晶沉积的铵盐从而形成强腐蚀性介质[21]，且该区域内流体与管壁间的对流传质作用最强，介质与管壁间易发生电化学反应，破坏金属基体，在管道内壁形成腐蚀产物保护膜。此外，a管排管束较b管排管束整体水相分率及传质系数更高，发生流动腐蚀的风险更大。

图9 空冷管束不同截面传质系数(kc)分布Fig.9 Distribution of mass transfer coefficient (kc) in the air-cooler pipesZ—Distance along axis; Z/mm: (a) 0; (b) 150; (c) 300; (d) 450

综上所述，a管排管束内部传质系数较大，对流传质作用更强；而b管排管束内壁剪切应力较大，易对腐蚀产物膜造成冲刷破坏。结合传质系数与剪切应力分析发现，流动腐蚀的高风险区域主要位于a管排中位号为a13～a14、a34～a35对应管束，而由于空冷器内存在偏流现象，根据水相分率分布规律，a13～a14管束发生流动腐蚀失效的风险最大。而b管排管束由于传质系数与剪切应力分布区域不一致，因此发生流动腐蚀的风险相对较小。

为进一步分析空冷管束沿流动方向的冲蚀高风险区域，取a13管道单独分析其内部流场，其流动参数分布如图10所示。从图中可以看出，a13管束内的传质系数及剪切应力分布趋势较为相似，均呈现先增大后减小的规律。原因在于：流体在经过衬管后突扩段时，受收缩-扩张结构的影响，在突扩段后局部区域会出现一个传质系数及剪切应力的峰值。随着流动距离的增加，湍流完全发展，两者逐渐减小最后趋于稳定。如图所示，两条曲线峰值之间的R1区域(即Z为11.5～26.4 mm)管束内部传质系数及剪切应力均处于较大值，推断该区域为管束中轴向的流动腐蚀高风险位置。

图10 位号a13管束内部流动参数分布Fig.10 Distribution of flow characteristic parameters in pipe a13kc—Mass transfer coefficient; Z—Distance along axis; τ—Wall shear stress

3 失效案例分析及验证

为验证预测结果的可靠性，对某炼油企业泄漏失效的高压空冷管束进行统计，结果发现，腐蚀失效最严重的区域对应管束位号为a13～a15、a33～a34，与图7和图9中的传质系数及剪切应力最大区域基本一致。其中失效管束的穿孔泄漏位置位于底部，对应仰视图见图11。可以看出，减薄穿孔的轴向位置位于衬管后方R2区域(即Z为10.2～25.5 mm)，与图10中预测的流动腐蚀高风险R1区域基本一致。

图11 高压空冷管束失效位置仰视图Fig.11 Bottom view of failure position in a high pressure air-cooler

4 结 论

(1)在铵盐结晶沉积预测的基础上，采用Mixture多相流模型及SSTk-ω湍流模型对空冷器入口管束内的多相流场进行数值分析，提出以传质系数及剪切应力作为表征流动腐蚀失效的关键参数；

(2)通过数值分析，获得了传质系数和剪切应力的分布特性。其中传质系数和剪切应力最大值的重合位置位于a管排管束位号为a13～a14、a34～a35的R1区域(即Z为11.5～26.4 mm)，是流动腐蚀的高风险区域；

(3)失效案例表明：空冷入口管束失效区域对应管束位号为a13～a15、a33～a34中的R2区域(即Z为10.2～25.5 mm)，与数值预测的流动腐蚀高风险区域一致，验证了分析和预测方法的正确性；

(4)流体经过入口配管的三通管段时存在偏流现象，k1一侧对应空冷管束与k2一侧对应空冷管束中实际流动参数存在差异，两侧管束发生流动腐蚀失效的风险也不相同。