1 前言

加氢反应流出物空冷器(Hydrogenation Reactor Effluent Air Coolers, REAC)系统输送的介质为气液两相流，气相中富含的NH₃、H₂S和HCl在冷却过程中会发生可逆反应^[1]，生成强腐蚀性的铵盐(NH₄Cl和NH₄HS)。工艺上采用在REAC系统上游设置注水点，溶解并冲洗沉积在管壁的铵盐。然而，若冲洗水和输送介质未经充分混合，或多相流动过程中水相分布不均匀，则易引发铵盐的沉积堵塞和垢下腐蚀。因此，通常需要在REAC系统的进口管道增设静态混合器。

在REAC系统入口静态混合器的设计选型过程中，应重点考虑以下两个因素：1)提高工艺注水的混合效果；2)在混合器出口处，壁面剪切应力低于其冲蚀破坏的临界值，避免管道冲蚀。由于静态混合器内部结构各异，流场较为复杂，难以进行理论求解，采用数值计算进行混合器内部流场和混合性能分析已得到广泛认同。Hirschberg等^[2]计算获得了结构优化后静态混合器SMX内部的内部压力场、流体混合效果和停留时间分布；Rahmani等^[3]对SMX和螺旋静态混合器的混合效果进行对比分析；郭瓦力等^[4]通过建立SK型静态混合器的甲苯-水两相混合物理模型，获得了浓度场云图和不均匀系数；张春梅等^{[5, 6]}等通过CFD软件分析静态混合器流场特性，并分析了结构特性和流场边界条件对混合效果的影响。在多相流模拟研究方面，王嘉骏等^[7]在VOF多相流模型的基础上，分别采用欧拉-示踪剂法和拉格朗日-示踪粒子法研究弯曲通道内液体混合过程；毛羽^[8]等采用欧拉多相流模型，对加氢裂化反应器内三相流动、传热、传质开展数值模拟计算；翟之平等^[9]基于Mixture多相流模型对叶片式抛送装置内的气固两相流场进行数值模拟。

本文在加氢REAC系统流动腐蚀机理研究的基础上，针对空冷器入口内置叶片式混合器的结构和介质特性，对静态混合器内多相流动特性进行数值分析，获得流速、湍流强度、水相分率、壁面剪切应力等流体动力学参数，研究混合机理。并且，通过水相分率和不均匀系数对混合效果进行表征。研究成果有望为静态混合器的设计和选型提供参考。

2 REAC系统流动腐蚀机理

为了避免REAC系统中铵盐的结晶沉积堵塞，在空冷器前设置工艺注水以洗涤气相冷凝析出的NH₄Cl和NH₄HS晶体。铵盐颗粒在微量液态水存在的工况下极易吸湿潮解，对碳钢空冷器管束造成碱性环境下的垢下腐蚀，如(1)、(2)所示：

$ \text{Fe}+2\text{N}{{\text{H}}_{4}}^{+}+2\text{C}{{\text{l}}^{-}}\to \text{F}{{\text{e}}^{2+}}+2\text{N}{{\text{H}}_{3}}(\text{l})+{{\text{H}}_{2}}\uparrow +2\text{C}{{\text{l}}^{-}} $

(1)

$ x\text{Fe}+y\text{N}{{\text{H}}_{4}}^{+}+y\text{H}{{\text{S}}^{-}}\to \text{F}{{\text{e}}_{x}}{{\text{S}}_{y}}+y\text{N}{{\text{H}}_{3}}+y{{\text{H}}_{2}}\uparrow $

(2)

在反应(2)中，铵盐溶液与碳钢基体反应生成的腐蚀产物膜Fe_xS_y，可有效阻止腐蚀介质对碳钢的进一步腐蚀。由于多相介质流动时在近壁面存在沿管束径向的速度梯度，当流动产生的剪切应力足够大时，腐蚀产物膜局部破损流失，碳钢基体裸露于腐蚀介质中继续腐蚀，形成破损-腐蚀-破损直至穿孔失效的自催化加速腐蚀体系。因此，对于REAC系统而言，其入口管道静态混合器的设计既要考虑混合后的水相分率分布均匀，又要确保局部最大剪切应力小于冲蚀临界剪切应力，避免冲蚀。

3 叶片式混合器计算模型 3.1 物理模型

加氢REAC系统入口叶片式混合器的结构组装如图 1(a)所示，包括筒壁、叶片和整流挡板，具体结构尺寸见图 1(b)。混合器长度l为160 mm，管壁内径为123 mm，叶片间夹角为120°，叶片与管壁垂直，沿管内壁一条母线展开，叶片与XZ平面的夹角为13.12°；三个整流挡板竖直放置，板间夹角为120°，与叶片交错分布。计算过程中，在距混合器上游2500 mm处(与实际注水点位置一致)设置注水点，并在混合器出口处增加长度为20倍(3200 mm)的直管段，使混合器进出口的多相流动达到充分发展的状态。

3.2 基本控制方程

加氢REAC系统介质为气-油-水三相流。在混合器上游的注水点附近，液态水由高压喷嘴注入多相流管道，此时的多相流流型为喷雾流。在混合器内，由于叶片和挡板的搅拌和剪切作用，液滴将逐渐均匀地分布在气相中，逐渐实现空间上的均匀分布。在混合器出口处附近，气相流速较高，气相会以气核的形式在管道中心流动，液相则成为贴壁液膜向下流动，此处的多相流流型为环状流和喷雾流共存。针对加氢反应流出物多相流介质，考虑到气液相间的耦合作用以及计算的稳定性，采用Mixture多相流模型进行计算求解^[10]。

Mixture多相流模型是一种简化的双流体模型，该模型假定短空间尺度上局部的平衡求解连续性方程、动量方程和能量方程，适用于模拟有强烈耦合且各相具有不同速度的多相流。在本文中，Mixture多相流模型通过求解混合相的连续性方程、动量方程、第二相的体积分数方程以及相对速度实现气-油-水三相的数值模拟。其中，主相为气相，次相为油相和水相。Mixture模型的控制方程如下：

连续性方程：

$ \frac{\partial }{\partial t}({{\rho }_{\text{m}}})+\nabla \bullet ({{\rho }_{\text{m}}}\boldsymbol{{{v}_{\text{m}}}})=0 $

(3)

多相混合动量方程：

$ \begin{array}{l} \frac{\partial }{{\partial t}}\left( {{\rho _{\rm{m}}}\boldsymbol{{v_{\rm{m}}}}} \right) + \nabla \cdot \left( {{\rho _{\rm{m}}}\boldsymbol{{v_{\rm{m}}}}\boldsymbol{{v_{\rm{m}}}}} \right) = \\ - \nabla p + \nabla \cdot \left[ {{\mu _{\rm{m}}}\left( {\nabla \boldsymbol{{v_{\rm{m}}}} + \nabla \boldsymbol{{v_{\rm{m}}}}^{\rm{T}}} \right)} \right] + {\rho _{\rm{m}}}\boldsymbol{g} + \\ \boldsymbol{F} + \nabla \cdot \left( {\sum\limits_{k = 1}^n {{\lambda _k}{\rho _k}\boldsymbol{{v_{dr,k}}}\boldsymbol{{v_{dr,k}}}} + {\sigma ^p}} \right) \end{array} $

(4)

第二相的体积分数方程：

从第二相p的连续性方程，可以得到第二相p的体积分数方程为：

$ \frac{\partial }{{\partial t}}({\alpha _{\rm{p}}}{\rho _{\rm{p}}}) + \nabla \bullet ({\alpha _{\rm{p}}}{\rho _{\rm{p}}}\boldsymbol{{v_{\rm{m}}}}) = - \nabla \bullet ({\alpha _{\rm{p}}}{\rho _{\rm{p}}}\boldsymbol{{v_{{\rm{dr,k}}}}}) $

(5)

湍流模型和近壁面处理方式：

多相流中气相的体积分数占三相总体积的97%以上，计算气相的雷诺数为1.32×10⁵，据此可以判断叶片式混合器内多相流的流动状态为湍流。湍流模型选择Realizablek-ε湍流模型。由于Realizablek-ε模型是高Re数的湍流模型，而近壁区内的流动，Re数较低，湍流发展并不充分，湍流的脉动影响不如分子黏性的影响大，因此将近壁面(Near-wall treatment)处理方式与Realizablek-ε湍流模型相结合。在沿法线的不同距离上，将流动划分为壁面区和核心区。在后续计算过程中，依据y_p⁺≈ 1的原则，对近壁面的网格进行加密。y_p⁺定义为：

$ y_{\rm{p}}^ + = \frac{{\rho \sqrt {{\tau _{\rm{w}}}/\rho \bullet {x_2}} }}{\mu } $

(6)

3.3 计算模型及边界条件

计算模型根据叶片式混合器的实际几何模型建立，求解器中的主要参数设置为：体积相分率采用一阶迎风格式进行离散，动量、湍动能和湍流耗散率均采用二阶迎风格式进行离散。选择加氢REAC典型工况和操作参数，运用工艺仿真软件ASPEN Plus进行计算，得到内置叶片式混合器入口的多相流介质组成及物性参数，其结果见表 1。

计算域的入口采用速度进口，出口采用自由流出口，压力-速度耦合采用PRESTO!格式，梯度采用Least Squares Cell Based，壁面按固壁无滑移条件处理。

3.4 网格无关性验证

利用Gambit对叶片式混合器进行网格划分，中间叶片及整流挡板采用非结构化网格，其余部分采用结构化网格，网格最大边长比为2:1，歪斜度均不大于0.1，网格总数为80万。采用表 1中的多相流物性参数和边界条件，分别划分80、100、120、140和160万五种网格密度，计算得到各网格密度下混合器的进出口压力差，如表 2所示。由此可知，相对于80万的网格数量，随着网格密度的逐渐增加，进出口压力差的相对变化值均小于3%，认为已经达到网格无关性要求。

4 计算结果分析 4.1 模拟结果检验

为检验数值计算的可靠性，本文采用环道式多相流实验装置开展叶片式混合器的压降试验^[11]。压降试验选用的多相流介质分别为氮气、白油和水，在实验过程中分别通过隔膜计量泵、磁力传动泵和往复压缩机调节水、白油、氮气的体积流量比率，使其接近实际比率。在混合器进出口设置压力传感器，测试进出口压力差。在不同流量下进行三组比对实验，实验结果如表 3所示。由表 3可知，实验值和模拟值吻合较好。

4.2 混合机理分析

采用多相流物性参数，开展数值模拟。图 2为X=0和Z=0截面内，叶片式混合器内部及出口区域的流速分布图。由图 2可知，叶片和挡板的存在使流道发生明显的改变，流线随之改变，多相流流速迅速增加。混合器内部的最高流速达到20 m·s^-1。并且，由于叶片上下两侧压力差的影响，叶片一、叶片三的下方形成了与流动方向相反的漩涡，可进一步增强三相间的混合作用。由图 3可知，流体在流经混合器时，流体的湍流强度显著增强，促进了三相间的相互渗透和扰动。随着流体逐渐远离混合管段，由于缺少混合叶片的持续作用，湍流强度逐渐减弱。在单相气体条件下，模拟混合器内部的速度和湍流强度分布，并与多相流条件下的计算结果进行对比分析，如图 4、5所示。对比图 2、4可知，上述两种介质下，混合器内流速的整体分布规律基本一致。然而，单相气体条件下的流速较低，最高流速为14 m·s^-1。同样地，由图 3、5可知，湍流强度的整体分布规律也较为相似，但单相气体条件下的湍流强度值较小。这是由于在多相流动条件下，气相和油相、水相间存在滑移速度，相间剪切作用导致湍流脉动更加强烈，导致了湍流强度的增加。

4.3 混合效果分析

图 6(a)~(h)为叶片式混合器及其进出口管道不同截面的水相分率云图。其中，图 6(a)、(b)为进口管道区域，图 6(c)~(f)为混合器所在区域(Y =0~-160 mm)，图 6(g)~(i)为出口管道区域。采用无量纲参数λ定义混合器出口沿程距离L与混合器长度l的比值。由图 6可知，当多相流进入混合器前，水相分布主要集中在管道的一侧，分布很不均匀；进入混合器后，在叶片及挡板的搅拌作用下，多相流流速及湍流强度迅速提高，水相进行了重新分布，完成了从高浓度聚集区向管道中心的扩散过程，均匀性得到了明显改善；离开混合器后，由于混合器出口的漩涡和高湍动能区域的存在，混合作用会持续一段距离；在距混合器出口400 mm后(λ=3)后，管道中的水相分率不再随流动距离的增加发生明显改变，表明其已经达到充分发展状态。

本文进一步采用不均匀系数φ定量计算分布混合程度^[12]。φ的定义如下：

$ \varphi = \frac{\sigma }{{\bar \phi }} $

(7)

其中，σ代表截面上水相分率的分布方差，可通过截面上所有节点的浓度计算得到：

$ \sigma = \sqrt {\frac{1}{N}\sum\limits_{i = 1}^N {{{({\phi _i}-\bar \phi )}^2}} } $

(8)

式中，N代表截面上的节点数，Φ_i代表每个节点上水相的体积分数，Φ是截面上水相分率的算术平方根，由式(9)计算：

$ \bar \phi = \frac{1}{N}\sum\limits_{i = 1}^N {{\phi _i}} $

(9)

根据不均匀系数的定义，φ越小，表示混合效果越好。图 7绘制了不同截面的不均匀系数，其中X轴为管道的沿程距离长度，X轴的负方向为流体的流动方向，X =0 mm为混合器的进口。由图 7可知，多相流在流经混合区域(X=0~160 mm)时，不均匀系数迅速降低，从0.44降低至0.19。随着流动沿程距离的增加，不均匀系数的下降趋势变得平缓，距混合器出口640 mm处(X=-800 mm，λ=3时)，不均匀系数降低至0.13。此后，随着流动充分发展，不均匀系数不再发生明显改变。与图 6对比可知，不均匀系数变化规律与管道横截面的水相分率分布基本一致。

空冷器入口管道的流速范围通常为2~6 m·s^-1。图 8为不同流速下，水相分率的不均匀系数随沿程距离的变化规律。由图 8可知，流速的增加未对混合效果造成明显的影响，证明该叶片式静态混合器对于REAC系统入口管道的多相流混合具有很好的适用性。

4.4 剪切应力分布

剪切应力是衡量碳钢管束是否存在冲蚀风险的重要指标，若局部剪切应力过大，则管壁的腐蚀产物保护膜容易变形脱落，形成自催化加速体系，易导致管壁冲刷减薄。因此，在考虑满足混合效果的前提下，应控制混合器出口剪切应力的数值低于临界值。图 9和图 10分别是当介质为单相气体和气-油-水三相流条件下，混合器各叶片后部管壁剪切应力分布云图。由图 9、10可知，叶片的存在使流场发生了突变，在叶片1和叶片2的下部管壁产生了较大的剪切应力，其中剪切应力较大的部位主要集中在整流挡板所在区域的壁面上。由于液相黏度高于气相，多相流动条件下会在壁面上产生更大的剪切应力，最高达23 Pa，高于冲蚀实验获得的冲蚀临界剪切应力15.2~15.8 Pa^[13]，具有很高的冲蚀风险。并且，由图 2、3可知，混合器出口附近的流速和湍动能较高，流态不稳定，高速气流易将液段击穿，从而造成液膜的失稳。当无稳定的液膜存在时，气相中夹带的液滴将增大对腐蚀产物保护膜的冲刷作用。因此，需要通过在混合器出口加装不锈钢衬套等方式，避免因冲刷腐蚀导致管道泄漏和爆管。

5 结论

本文在加氢REAC系统流动腐蚀机理研究的基础上，运用CFD模拟，计算叶片式混合器内的多相流动特性，并对模拟所得到的混合器进出口压降与实验值进行比对，验证了计算结果的可靠性。通过分析混合器内部的流速、湍动能、水相分率、剪切应力等流动参数，得到以下结论：

(1)多相流流经混合器时，流道变化剧烈，混合区域内的湍流强度明显增强，促进了油、气、水三相间的相互渗透和扰动。最下方叶片的底部会形成较大尺度的漩涡，对混合器出口的流体起到了持续的混合作用；

(2)通过分析管道各截面的水相分率和不均匀系数φ可知，多相流在流经混合器后，水相存在从局部聚集区向管道中心的扩散过程，其不均匀系数φ从0.44降低至0.19，不均匀性得到了明显改善；当流速范围为2~6 m·s^-1时，流速的改变对混合效果没有明显的影响；

(3)采用无量纲参数λ定义混合器出口沿程距离L与混合器长度l的比值，计算发现，当λ=3时，多相流动达到充分发展状态，水相分率的不均匀系数不再发生明显改变。

(4)由于混合区域内流速急剧增加，会在出口管道壁面形成剪切应力较大的区域，接近碳钢材料冲蚀破坏的临界值，需考虑在混合器出口增加衬套，防止壁面冲刷减薄。

符号说明：

F        -体积力，N

ρ_m        -混合密度，kg·m^-3

L        -混合器出口沿程距离，mm

σ        -水相分率的分布方差

l        -混合器长度，mm

ρ        -流体密度，kg·m^-3

n        -相数

ρ_p        -第二相的密度，kg·m^-3

V        -流速，m·s^-1

α_p        -第二相的相分率

V_{dr, k}        -第二相的漂移速度，m·s^-1

λ        -混合器出口沿程距离与混合器长度的比值

V_m        -质量平均速度，m·s^-1

τ_w        -壁面剪切应力，Pa

x₂        -到壁面的距离，mm

μ        -流体的运动黏度，m²·s^-1

Y_M        -可压缩湍流脉动膨胀对总的耗散率的影响

φ        -不均匀系数

φ_i        -每个节点上水相的体积分数

φ        -截面上水相分率的算术平方根

μ_m        -混合黏性，m²·s^-1