基于CFD-DPM耦合方法的多级串联弯管内液固两相流冲蚀特性研究

高凌霄，杨贺琦，李小平，敬登伟，马利静

(1.中海油(天津)管道工程技术有限公司，天津 300452;2.西安交通大学动力工程多相流国家重点实验室，陕西 西安 710049)

随着我国陆地油气产量逐年萎缩与国家对石油的需求量稳步增长，油气资源蕴藏量相对较为丰富的海上油田成为我国石油产量的增长及产量接替的重要组成部分[1]。研究表明，我国绝大多数海上油田分布在疏松砂岩地层中，这使得海上油气开采过程的出砂问题尤为严重[2]。出砂引起的各种腐蚀危害不仅增加了开采难度，也导致尤其生产成本的大幅增加[3]。此外，作为油气田开发过程中的油气输送工具，管道中的砂粒侵蚀会导致石油和天然气开采过程中的严重设备故障，并可能限制最大生产率。因此，预测油气生产管道中砂粒冲蚀引起的腐蚀情况是降低管道失效率的关键因素。

关于单相和两相系统中侵蚀率的预测，已有大量研究，但关于气-液-砂多相流在不同管道结构中的流动侵蚀情况研究还相对较少[4]。Liu等人提出了一个简单的基于CFD的模型，用于预测环形流中的侵蚀率。在他们的工作中，假设弯管部位液膜厚度不变，模拟中间接考虑了薄膜液体的缓冲效应，从而观察到环状流弯管中的侵蚀明显减少[5]。Parsi等人建立CFD模型研究了标准弯管处两种不同粒径颗粒的侵蚀率，发现液膜厚度在侵蚀程度中起着重要作用[6]。Zahedi等人分别使用VOF和多流体VOF方法模拟了76.2 mm标准弯头中低液体流速和高气体流速的空气-水流动，作者同时模拟了颗粒轨迹和流动，并预测了几种表观气速和液速下的材料腐蚀情况[7]。结果表明，与VOF模型相比，多流体VOF模型存在一些收敛问题，需要更多的计算时长。相比之下，VOF模型更稳定，需要更少的计算机时。Ogunsan等人将欧拉多流体VOF模型与拉格朗日方法耦合，证明模拟结果与实验数据吻合良好，并发现最大侵蚀率出现在标准弯头的45°处[8]。Sedrez等人用超声波传感器测量了携带300 μm砂粒的空气-水多相流弯管中的壁厚损失。通过将欧拉-拉格朗日方法对不同的流动条件进行的CFD模拟结果与实验数据进行比对，证明了CFD模拟的有效性[9]。

上述已有的数值模拟研究主要针对简单的单弯管工况，而未考虑实际管道输运过程中更为复杂但可能有效降低含砂流体冲蚀的结构，例如多级串联弯管。本研究的主要目标即是开发一个计算模型，以便更准确地预测复杂结构弯管中含砂颗粒的侵蚀速率。通过计算油气连续相流场、管壁侵蚀率和颗粒轨迹，预测侵蚀最严重区域，重点分析颗粒运动轨迹以及颗粒在管壁处的碰撞规律，总结出管道侵蚀与砂粒运动间的关联关系。通过对比单个弯管与多级串联弯管内液固两相流冲蚀特性的差异，为未来含砂多相流体输运管道的设计及腐蚀防护提供理论指导。

本文利用Ansys Fluent软件对90°弯管冲蚀进行数值模拟研究，具体包括连续相(液相)流场模拟、离散相固体颗粒轨迹跟踪和弯管冲蚀速率计算三个部分。本文首先将模拟带颗粒的湍流液体在管道中的流动形态，然后根据拉格朗日离散相模型，分析流场中离散相颗粒所受力情况，从而求出颗粒的运动轨迹。最后，在颗粒与管壁发生碰撞时，利用冲蚀预测模型预测颗粒对管壁的冲蚀范围。

1.1 连续相模型

本文采用湍流模型模拟连续相(液相)流场。对于不同的问题，需要考虑流体的可压缩性、计算精度的要求、计算机设备的性能和时间限制，选择最合适的湍流模型。Realizablek-ε湍流模型可以增强在强逆压梯度情况下的边界层性能，且适用于平面射流、圆形射流、旋转流等复杂流动工况[10-12]。本文选用Realizablek-ε湍流模型对管道弯头冲蚀磨损进行数值模拟，其中，湍流动能(k)及其能量耗散率(ε)输运计算方程如下

(1)

其中

式中xi、xj——空间坐标；

ui——速度矢量；

ρ——液体密度；

κ——湍流动能；

μ——分子黏度；

ut——湍流粘度；

σκ、σε——方程k和ε的湍流普朗特数，σκ=1.0，σε=1.2；

GΚ、Gb——湍流动能的产生;

GΚ——由平均速度梯度导致的湍动能项;

Gb——由浮力因素引起的湍动能项；

YΜ——可压缩湍流脉动对总耗散率的贡献；

Ε——用户定义源项；

υ——方向与重力方向平行的速度分量；

Sij——应变变化率张量；

C2、C1ε、C3ε——常数，C2=1.9,C1ε=1.44,对于流动方向与重力方向相同C3ε=1，对于流动方向与重力方向垂直C3ε=0。

1.2 离散相模型

(3)

式中ρ——流体流体密度/kg·m-3；

u——流体相对速度/m·s-1；

up——固体颗粒速度/m·s-1；

μ——流体的分子粘度/Pa·s；

ρp——固体颗粒的密度/kg·m-3；

dp——固体颗粒直径/m；

Re——相对雷诺数;

CD——曳力系数;

“对，去找穿着制服的交通警察。如果你和妈妈分开了，妈妈会很难过很难过的，所以我们一起出门时，你要跟紧妈妈。”

gy——y方向重力加速度/m·s-2,gy=9.81 m/s2。

1.3 冲蚀模型

McLaury模型是ANSYS Fluent软件中一种应用于预测水中固体颗粒冲蚀速率的侵蚀模型，该模型主要应用于泥浆流，其表达式如下[16-17]

E=AVnf(θ)

(4)

A=F(Bh)k

(5)

(6)

式中E——冲蚀率，代表单位时间内单位面积管壁壁面的质量损失；

F——经验常数；

V——固体颗粒撞击壁面速度；

θ——颗粒撞击壁面的角度；

Bh——壁面的布氏硬度值，其大小为156；

k——指数项，由于本模拟选区的是碳钢材料，此处k值选取为-0.59；
其他均为常数，具体数值见下表1。

表1 McLaury模型常数表

固体颗粒在弯管内运动期间，粒子可能会与管壁碰撞，然后反弹回流体域。增强壁面函数(Enhanced Wall Function)可用于处理近壁问题，本文采用Grant和Tabakoff提出的壁面反弹恢复系数公式来处理近壁区域的流动工况[18]

en=0.993-0.037θ+4.75×10-4θ2-2.61×10-6θ3

(7)

et=0.988-0.029θ+6.43×10-4θ2-3.56×10-6θ3

(8)

式中en和et——法向和切向的反弹系数；

θ——固体颗粒的入射速度和表面切线之间的角度。

2.1 弯管几何模型与网格划分

本文选用内径D=50 mm的90°弯管作为研究对象，单弯管由入口直管段L1、弯头部分、出口直管段L2三部分组成，其几何结构如图1(a)所示；
串联双弯管由入口直管段L1、上(下)游弯头、连接直管段L、出口直管段L2五部分组成，其几何结构如图1(b)所示。为确保流体和颗粒得到充分发展，定义进口直管段和出口直管段长度为20D，即L1=L2=1 000 mm，弯头曲率半径R=50 mm。

如图2所示，运用ICEM CFD软件采用O型划分方法对弯管进行结构化网格划分。为提高模拟计算精度、节约计算成本，依据谢振强等人[19]提出的网格划分建议，本文模型中弯头轴向和周向网格尺寸范围为(10～20)dp(dp为颗粒直径)。此外，为避免由于不同网格类型和大小的划分而导致模拟结果出现偏差，本文进行了网格无关性验证，验证结果如图3所示，以网格数量为自变量，弯管最大冲蚀率为因变量，探究不同网格数量下最大冲蚀率变化情况。由图3可以看出，当选取的网格数量小于13×105时，计算得到的最大冲蚀率随网格数量增加有明显的波动性增加，而当网格数量大于该值时，计算的最大冲蚀率基本不再随网格数量的增加而变化。故为最大限度节省计算资源并降低计算误差，本文选用网格数量为1 377 000的网格划分方法对弯管进行结构化网格划分。经检验，其网格质量均在0.75以上。

2.2 边界条件设置

本文研究中边界条件设置如下：连续相密度定义为998.2 kg/m3，离散相固体颗粒密度定义为2 650 kg/m3。对于连续相计算区域，给出管道的进口速度值，而管道出口的压力进行统计。此处，压力设置为常压而湍流的强度范围设定为5%，并选取SIMPLE为压力-速度耦合方案，空间微分的梯度为Green-Gauss节点，压力设定为PRESTO！且选择离散格式，其他计算项均选用二阶迎风差分格式。壁面为静止壁面且无滑移，采用增强壁面函数(EWF)对近壁区域进行处理；
对于离散相固体颗粒，速度与连续相液体保持一致，管道入口和出口均设置为逃逸(Escape)条件，管壁设置为反弹(Reflect)条件，壁面反弹恢复系数公式见式(7)式(8)。假设管道中入射颗粒均为相互独立且形状均匀的球形粒子，同时忽略粒子与粒子之间的相互作用带来的影响。

3.1 弯管冲蚀分布情况

以管道内径为50 mm、弯管曲率半径为50 mm、流体流速为10 m/s、颗粒质量流量为0.20 kg/s、颗粒直径为28 μm的工况为例，本文首先对弯管内的流场分布进行了计算，获得了弯管壁面冲蚀分布、弯管壁面最大冲蚀率、弯管流场的压力分布、速度分布和颗粒运动轨迹等关键参数信息。

从图4可以看出，对于单弯管的情形，冲蚀集中在弯管外侧壁中心区域，冲蚀区域呈椭圆形分布，最大冲蚀位置在弯管外弧60°～90°范围内且冲蚀区域呈上宽下窄的梭形分布，最大冲蚀率为2.06×10-3kg/(m2·s)，而弯管内壁几乎未发生冲蚀，这与李睿等人数值模拟所得结论一致且在颗粒轨迹图中得到了证实，也从侧面证明了本研究中采用的数理模型的可靠性[11]。从图5可以看出，粒子在进入弯管前其运动轨迹与上游直管轴线近乎平行，后进入弯管区域强烈地撞击弯管外弧管壁并发生反弹，运动轨迹发生显著变化。又由弯头部分局部放大图可知，弯头外侧中心线两侧的粒子发生交错流动，即弯头外侧中心线右边的粒子流向中心线左边，中心线左边的粒子流向中心线右边，这与弯管外侧壁冲蚀分布图案完美契合。在流体作用下，反弹后的粒子并未流向弯管内侧，而是顺应流线流出弯管区域，弯管内侧出现一块无碰撞区域，部分颗粒进入下游直管道后与壁面发生二次碰撞，这便进一步解释了为什么弯管内壁几乎未发生冲蚀。

图6显示了弯管流场的压力分布，图7为速度分布。从两图中可知，弯管外侧压力明显上升，形成高压区，而速度下降；
弯管内侧的压力值大小沿弯管径向逐渐降低，并指向弯管曲率中心形成顺压梯度，在弯管内侧，压力达到最小值，流速最高。同时从弯管0°～90°各截面压力和速度分布可知，流体进入弯头后弯头外侧压力先增加后下降，内侧趋势则与之相反；
而弯头外侧速度逐渐增加，内侧速度呈先增加后下降趋势。针对上述现象，作者推测在离心力作用下弯管内侧流体流向弯管外侧并对弯管外壁进行挤压，造成弯头外侧压力高而内侧压力低，弯管内侧动能由弯管内侧比压能转化而增加，导致弯管内侧流体的流速升高。

3.2 颗粒流量对弯管局部冲蚀的影响

为了研究不同颗粒质量流量对弯管局部冲蚀的影响程度，本文选择直径为28 μm的颗粒，假定液体流速均为10 m/s，研究了0.05～0.25 kg/s流量范围内颗粒质量流量对弯管的局部冲蚀情况，基于冲蚀云图的分析，进一步得到了不同砂粒质量流量下的平均及最大冲蚀率。由图9的云图可以看出，随着颗粒质量流量的不断增加，弯管内部砂粒碰撞产生的冲蚀区域位置及面积大小基本无变化。但在图8中，最大冲蚀率和平均冲蚀率呈显著的线性上升趋势，这与Hong等人的研究结果一致，说明管道中颗粒含量的增加会增强管壁的冲蚀磨损[13]。推测随着颗粒质量流量的增加，管道中的颗粒数量也会随之增加，从而导致单位面积上颗粒对管壁的冲击频率增加，最终使得管壁冲蚀速率显著增加。

3.3 液体流速对弯管冲蚀的影响

同样，在其他影响因素一定的工况下(颗粒直径为28 μm、颗粒质量流量为0.20 kg/s)，本文取不同液体流速进行研究，得到不同液速下弯管冲蚀云图、最大冲蚀率和平均冲蚀率。如图10所示，随着液速的不断增加，弯管冲蚀区域未发生明显变化，主要集中在弯管内壁外侧和下游直管道连接处，而弯管最大冲蚀区域随液速增大而不断扩大。不同流速与弯管最大冲蚀率和平均冲蚀率关系曲线如图11所示。结果与徐上在探究弯管入口速度与冲蚀速率关系时得到的规律基本一致，即弯管的最大冲蚀率随流速的增加而显著增大，且增长速度加快[21]。作者推测上述现象的产生可能受以下两方面因素影响：一方面，流体中的固体颗粒随着流体的流速增加而具有更大的动能，从而以更大的速度对弯管产生碰撞，并导致冲蚀加剧。同时，当流速达到某一程度时，固体颗粒的随流性得到改善，连续相和离散相双向耦合作用提高，颗粒的湍流动能增加，这进一步加剧了颗粒对管壁的冲蚀[21]；
另一方面，随着流速不断增加，固体颗粒与管壁碰撞频率随之增高，弯管内壁在两方面因素的共同作用下冲蚀程度将愈发严重。

3.4 串联弯管连接长度对弯管冲蚀的影响

如前述，目前文献中已有的数值研究主要针对简单的单弯管工况，而未考虑实际管道输运过程中更为复杂但可能有效降低含砂流体冲蚀的结构，例如多级串联弯管。本文通过对比单个弯管与多级串联弯管内液固两相流冲蚀特性的差异，为未来含砂多相流体输运管道的设计及腐蚀防护提供理论指导。为探究串联弯管不同连接长度对弯管冲蚀的影响，本文选择颗粒直径为28 μm、液体流速为10 m/s、颗粒质量流量为0.20 kg/s的工况，取串联弯管连接长度3～18D进行分析，得到了不同连接长度下的最大冲蚀率、平均冲蚀率和串联弯管冲蚀云图。图12依次为不同连接长度(L)的上游弯头的冲蚀云图。可以发现，两个弯头之间的距离从3D增加到18D，但是上游弯头的冲蚀模式几乎没有变化。此外，上游弯管最大冲蚀区始终位于同一位置。推测可能是因为上游直管中的流动条件(包括速度和压降)对于所有连接长度都是相同的，导致上游弯头中的冲击速度和冲击位置相同。进一步观察可发现，串联弯管结构中，最大侵蚀区均出现在上游弯头外弧中心线附近。推测可能原因为：上游直管核心区域的流速更高，更大的曳力使颗粒加速，此外，管壁和颗粒之间的摩擦阻力降低了颗粒的速度，故而核心区域的颗粒比靠近管壁的颗粒运动得更快。当连续相从上游管道进入弯管后，在弯管处发生90°方向变化，而核心区域的颗粒撞击弯管壁面的速度高于靠近管壁的颗粒，因此具有较大动能，从而在该区域产生较为严重的冲蚀。

图13显示了下游弯头的冲蚀分布情况。该冲蚀云图不同于上游弯头变化模式。随着两个弯管之间的连接长度从3D增加到18D，在下游弯管中可以观察到冲蚀的显著变化。当连接长度L=3D时，下游弯头冲蚀区呈“W”型分布，且下游弯管与连接直管相交位置出现冲蚀斑块群；
当连接长度L=6D时，下游弯头冲蚀区两个最大冲蚀位点关于弯管外侧中心线基本对称；
当连接长度L=9D过渡到L=18D时，下游弯头冲蚀区呈“V”型分布，且最大冲蚀区域向管道外侧中心线集中。

通过与图4单弯管冲蚀分布对比可知，多级串联弯管的上游弯管与单级弯管在相同流动条件下的冲蚀分布图案类似，均为椭圆状。这体现了单个弯管与多级串联弯管内液固两相流冲蚀特性的共性，然而单弯管与多级弯管冲蚀特性的差异主要体现在下游流动中。通过进一步分析图12和图13上下游弯管冲蚀云图，不难发现随着连接长度的增加，串联弯管的最大冲蚀位点发生变化，逐渐由上游弯管转移到下游弯管，这对含砂多相流输运管道设计具有重大指导意义，并且下文对最佳连接长度的确定进行了详细分析阐释。

图14为不同连接长度与最大冲蚀率和平均冲蚀率曲线。可以看出，平均冲蚀率随连接长度增加整体呈下降趋势，而最大冲蚀率呈现先增大后减小而后又逐步增大趋势。连接长度从3D增加到6D最大冲蚀率上升可能因为混流颗粒交换动量空间不足，未得到充分发展，多处于湍流状态，撞击壁面频率增加，导致壁面冲蚀严重。当连接长度从9D增加到18D时，混合颗粒流变成完全发展的流动，类似于上游管道中的流动，故而最大冲蚀率不断增加，这一趋势也在文献的研究工作中得到了佐证[22]。基于上述研究结果，未来在含砂多相流体输运管道的设计中，串联弯管的最佳连接长度建议可选在9～15D范围内，以降低液体中固体颗粒对壁面的冲蚀作用。

本文基于CFD-DPM耦合方法，实现复杂结构弯管中含砂颗粒的侵蚀速率的准确预报。通过模拟连续相(液相)流场、跟踪离散相固体颗粒轨迹和计算弯管冲蚀速率，精准预测更有可能发生严重侵蚀的位点，通过对比单个弯管与多级串联弯管内液固两相流冲蚀特性的差异，为未来含砂多相流体输运管道的设计及腐蚀防护提供理论指导。本文取得的主要结论如下：

(1)研究弯管内的流场分布，弯管外侧压力高速度低，极易发生冲蚀破坏，而弯管内侧存在无碰撞区域，发生冲蚀可能性小，这与已有文献结论基本一致。

(2)管道中颗粒含量的增加会增强管壁的冲蚀磨损，颗粒质量流量与弯管最大冲蚀率和平均冲蚀率呈正相关；
随着流体流速增加，冲蚀加剧，并且流速越大，最大冲蚀率和平均冲蚀率增长斜率越大。

(3)随着串联弯头之间的连接长度增加弯管平均冲蚀率整体呈下降趋势，而最大冲蚀率先增大后减小又逐步增大。未来在含砂多相流体输运管道的设计中，为降低壁面冲蚀磨损，串联弯管连接长度可选在(9～15)D范围内。