彭文山，马力，侯健，曹学文，韩明一

（1.中国船舶重工集团公司第七二五研究所 海洋腐蚀与防护重点实验室，山东 青岛 266237；2.中国石油大学（华东）储运与建筑工程学院，山东 青岛 266580；3.廊坊中油朗威工程项目管理有限公司，河北 廊坊 065000）

随着我国油气开采由陆地向海洋转移，海底管道投用量持续增加[1]，海洋混输管网得到不断应用。混输管道中流体流动复杂，流体对管壁冲刷严重，油井出砂量及原油含水量随着油气田开采时间的增长而增加。油气开采及运输过程中，管内固体颗粒随油气流动，极大地加剧了管内介质对管壁的冲击，极易造成管道冲蚀破坏，导致危险事故[2-3]。因此，提出准确的多相流冲蚀计算方法对于预测管道可能出现的冲蚀最严重位置，进而保障管输安全、减少经济损失意义重大。

泡状流在垂直上升气液两相流中是常见的流型[4]，例如海底垂直上升管道以及海洋立管中。含固体颗粒泡状流中，颗粒受气液流动影响，空间分布随机性大，而管道冲蚀最严重位置与颗粒的分布直接相关，要准确仿真多相流冲蚀，需考虑气体-液体-颗粒-管道壁面之间的多向耦合，研究具有挑战性。国外部分学者开展了含固体颗粒多相流冲蚀实验，在实验结果基础上提出了一些经验和半经验冲蚀计算公式[5-7]，但这些公式仅可计算最大冲蚀速率，局限性较大。Chen 等人[8]提出采用将气液两相简化成混合流体相来计算管道内颗粒的冲蚀，计算获得的冲蚀结果偏小，且无法获得管道内部气液分布及颗粒运动情况。因此以上方法仅适应于估算管道冲蚀速率，并不能正确解释管道冲蚀机理。因此，提出了一种基于流体体积模型（VOF）模型和离散相模型（DPM）模型的多相流瞬态冲蚀仿真计算方法，研究管内气液分布、颗粒运动对冲蚀的影响，并与经验模型计算结果、简化计算流体动力学（CFD）计算结果及实验结果进行对比，验证该模型准确性，为含固体颗粒管道的分散泡状流冲蚀预测提供参考。

1 分散泡状流冲蚀求解

管道内为气液固三相流，气液两相为连续相，固体颗粒为离散相，采用DPM 模型完成固体颗粒的计算。分散泡状流冲蚀求解分为以下三步：流场计算、颗粒追踪及冲蚀速率求解。采用Eulerian-Lagrangian方法[9]，求解气液两相连续相流场和离散相颗粒受力方程，进行颗粒轨迹计算，冲蚀计算模型提取前两步计算获得的信息，计算得到冲蚀速率。

1.1 连续相流场计算

标准k-ε 模型在计算管道内部不同截面处流体速度的误差较小[10]，选用该模型计算管道内部流场。使用VOF 模型实现管内分散泡状流的求解。另外，考虑表面张力的作用，采用Brackbill 等人[11]提出的连续表面力模型（CSF）。利用这个模型，VOF 模型中的附加表面张力通过源项的方式添加到动量方程中。

1.2 离散相控制方程

颗粒轨迹追踪过程中不考虑颗粒之间的相互作用，单位质量固体颗粒在分散泡状流气液两相流中的受力方程为：

1.3 冲蚀计算模型

目前已有许多冲蚀模型被应用于冲蚀计算[13]，Liu等人[14]对冲蚀计算研究发现，广岛大学Oka 等人[15-16]提出的冲蚀模型能较为准确地预测多相流条件下管道的冲蚀速率。另外，模型中考虑了较多的影响因素，适应范围较广。使用该模型进行冲蚀计算：

式中：ρw是管材密度，kg/m3；α 为颗粒冲击角度，rad；ρw为管材密度，kg/m3；HV 是管材维氏硬度；uC为参考颗粒速度，m/s；dC为参考直径，m；k0、k1、k2、k3、n1、n2、dC、uC为常数，见表1[15-16]。

表1 Oka 等冲蚀模型参数Tab.1 Parameters of erosion model proposed by Oka et al.

2 仿真模型建立

2.1 管道基本参数

Birchenough 等[17]进行冲蚀实验，获得一组竖直管道中空气-水系统分散泡状流冲蚀的实验结果，在进行数值分析时，与该组实验结果进行了对比验证。Birchenough 等人的实验管道管径为49 mm，工作压力为0.2 MPa，管道材料为碳钢，仿真参数值见表2，D 为管道直径，R/D 为弯径比，VSG、VSL分别为气相和液相的表观流速，ρl、ρg、ρp分别为液体、气体、砂粒密度。弯管模型由进口直管段、弯头和出口直管段组成。为使得管内流动充分发展，并考虑模型数值计算消耗，入口直管段长度取为50 倍管径（50D），出口直管段长度取为25 倍管径（25D），如图1 所示。水的表面张力为0.073 N/m。

表2 分散泡状流仿真参数Tab.2 Simulation parameters of dispersed bubble flow

图1 竖直弯管几何模型Fig.1 Geometry model of the vertical pipe bend

2.2 网格划分

为提高VOF 模型与DPM 模型耦合计算的稳定性，计算模型采用规则四边形面网格，采用六面体弯管体网格，网格划分如图2 所示，竖直弯管三维模型网格总数为524 000。

图2 竖直弯管网格划分Fig.2 Computational mesh of the vertical pipe bend

2.3 边界条件

模型入口采用速度入口边界，速度大小设置为气液混合速度，出口采用压力出口。入口水力直径为49 mm，湍流强度经计算为3.2%。DPM 模型中进口和出口采用逃逸（Escape）条件，壁面采用反弹（Reflect）条件，采用面射流源喷射固体颗粒，固体颗粒初始速度与流体入口速度相同。管壁壁面粗糙度常数设置为0.5，并设定为“静止壁面”和“无滑移壁面”。壁面设置恢复系数来计算固体颗粒碰撞壁面前后的速度变化，计算过程中使用Grant 和Tabakoff[18]提出的随机颗粒-壁面碰撞反弹模型的恢复系数，形式如下：

式中：en为法向恢复系数；et为切向恢复系数。

2.4 数值求解

使用FLUENT 软件完成冲蚀计算，基于VOF 模型和DPM 模型的CFD 冲蚀计算中，压力-速度耦合采用PISO 方式，压力离散采用PRESTO!方式，动量方程、湍动能方程、湍流扩散率方程均使用QUICK格式，体积分数离散采用Geo-Reconstruct 方式，来提高计算精度。为了使收敛比较平稳，并保证模拟结果的可靠性，时间步长设置为0.001 s。颗粒同样采用瞬态计算，计算开始前打开离散相模型，加入离散相粒子，追踪时间步长与连续相计算步长相同。在瞬态计算过程中，每计算一步，固体颗粒在入口处“喷射”一次，入口处有500 个网格单元，则入口处每一时间步要入射500 个固体颗粒。总的计算时间为8 s，时间步长为0.001 s，则数值计算共喷射进入管道的颗粒数目为4 000 000 个。颗粒轨迹追踪时，假定颗粒不影响气液两相流动，仅考虑气液两相对颗粒运动的影响，即单向耦合计算。为了便于分析管道内部气液固多相的运动情况，分别在上、下游直管段中间位置以及弯头45°截面处设置监控面，对瞬态计算结果进行监控保存如图1 中红色虚线位置。

3 结果分析

3.1 冲蚀结果对比

基于VOF 模型和DPM 模型的CFD 冲蚀计算中，冲蚀量随时间增长是持续增大的，冲蚀速率的大小为冲蚀量随时间变化的快慢，选取计算过程中冲蚀量随时间变化最快的计算结果为最大冲蚀速率。FLUENT 软件数值模拟直接计算获得的冲蚀单位为冲蚀速率（kg/(m2·s)），为了便于同Birchenough 等人的实验结果对比，冲蚀速率的单位转换方式如下：

ER（m·kg-1）=ER（kg·m-2·s-1）/ρw（kg·m-3）/mp（kg·s-1）

式中：mp为固体颗粒质量流量。

将瞬态模型计算结果与 Salama 经验模型[5]、Bourgoyne 经验模型[6]以及简化CFD 模型计算结果做对比。

Salama 经验模型：

式中：Sm为管道形状系数，m2/s2，取2000；Vm为气液混合流速，Vm=VSG+VSL，m/s；ρm为混合流体密度，kg/m3。

Bourgoyne 经验模型：

式中：Fe为砂粒侵蚀系数，取为0.01g/kg；ρp为固体颗粒密度，kg/m3；ρw为管壁材料密度，kg/m3；Wp为固体颗粒体积流量，m3/s；Apipe为管道横截面积，m2；HL为液体体积分数。

简化CFD 冲蚀计算思路为：冲蚀简化计算方法假设固体颗粒均匀地分布在气液两相中，将分散泡状流假定为混合均匀的单相流体，采用混合密度、混合黏度以及混合速度完成稳态的冲蚀计算，见式（7）、（8）。该方法与Chen 等人[8]提出的分散泡状流冲蚀求解思路一致。

式中：μg、μl、μm分别为气体、液体、混合流体的黏度，Pa·s。

Birchenough 等人实验结果为4.6×10-9m/kg，采用冲蚀预测值/实验值来直观地反映不同模型的计算精度，如图3 所示。由图3 可知，基于VOF 模型和DPM 模型的CFD 冲蚀计算结果最接近实验值。经验模型可直观快速地确定冲蚀速率大小，但是不能确定管道内部的流动过程及冲蚀细节。简化CFD 计算模型由于将流体简化，其管道最终冲蚀分布并不能反映实际冲蚀情况。

图3 不同冲蚀模型计算结果对比Fig.3 Comparison between calculation results by different erosion models

3.2 颗粒运动轨迹与气液分布关系

利用截面含液率来表示某时刻管道截面气液分布情况，不同时刻管内固体颗粒瞬态分布与气液两相分布之间的关系如图4 所示。由图4 可知，管道内部颗粒稀疏的区域（箭头指向区域），液相含量也相对较少，气相含量高，说明颗粒与液相分布具有密切相关性。为深入研究颗粒分布与气液两相之间的关系，进一步分析了管道三个监控面上的截面含液率与相对颗粒含量之间的关系，如图5 所示。截面含液率的峰值变化规律与相对颗粒含量的峰值变化趋势比较接近，截面处含液率较高时（通过液体较多），截面的颗粒含量也相对较高，说明管内的固体颗粒大部分分散在液相之中。截面含液率和相对颗粒含量在不同截面的波动范围是不同的，如图6、图7 所示。

图4 弯管不同时刻颗粒运动轨迹与管内气液分布Fig.4 Particle trajectories vs gas and liquid distributions in pipe bend at different time

由图6 可知，上游直管段截面处的截面含液率波动幅度小于弯头处和下游直管段截面处。这主要是由于相对于弯头以及下游直管段，上游直管段中气液两相流动比较均匀，固体颗粒分布也较均匀。当颗粒在弯头处发生碰撞后，颗粒的运动轨迹发生改变，并且由于气液两相流经弯头后，分散泡状流流型发生明显变化，在弯头顶部及下游直管段顶部液相含量明显增大，而在弯头及下游直管段底部，气相含量明显增大，气泡与液体的掺混更加剧烈，导致截面颗粒含量的波动较大。截面颗粒含量平均水平同样也是下游直管段和弯头处大于上游直管段，如图7 所示。这主要是由于部分通过弯头的颗粒在进入下游直管段过程中会穿越弯头顶部和下游直管段顶部的流体，而该区域的液相比例较高，固体颗粒进入该区域后，与大量液体进行交互计算，液相密度及黏度比气相大得多。由公式（1）可知，固体颗粒进入该区域后速度减小，形成滞止区，导致颗粒运动减缓，颗粒含量变大。

3.3 冲蚀分布与气液分布关系

基于VOF 模型和DPM 模型的瞬态冲蚀计算获得的冲蚀分布如图8 所示。CFD 计算8 s 时，冲蚀量为2.24× 10-4kg/m2。瞬态仿真计算过程中，管道冲蚀量随时间变化最快时的结果为最大冲蚀速率，求得最大冲蚀速率为4.41×10-9m/kg，冲蚀最严重位置出现在弯头出口位置处。Birchenough 等人的实验结果为4.6× 10-9m/kg，差别很小。气液两相在弯头处分布不均，在弯头弯头出口处（90°处）顶部（A 区域）液相含量较多，拥有较大速度的颗粒在碰撞管壁之前与液体相互作用，使得颗粒碰撞管壁速率减小，液体的存在对弯头起到一定的保护作用，在弯头小于90°角度处，部分位置气相含量较高（B 区域，如60°、45°位置），颗粒与壁面碰撞时能量依然较大。虽然90°附近（A 区域）的颗粒碰撞能较大，但是损失也多，而小于90°区域（B 区域，如60°、45°位置）碰撞能略小，但是损失也较小，综合以上因素，导致冲蚀严重区域较大。

图5 不同截面处含液率与颗粒含量Fig.5 Liquid hold-up vs particle mass concentration at different sections: a) Section of upstream straight pipe in bend; b) 45° section of elbow; c) Section of downstream straight pipe in bend

4 结论

1）基于VOF 模型和DPM 模型的CFD 冲蚀模型计算结果与实验值接近，且与经验模型及简化CFD模型相比，可直观反映管道内部流场及固体颗粒运动与壁面冲蚀之间的关系。

图6 不同截面处截面含液率Fig.6 Liquid hold-up at different sections

图7 不同截面处截面颗粒含量Fig.7 Particle mass concentration at different pipe sections

图8 基于VOF 模型和DPM 模型耦合计算弯管冲蚀分布及气液分布Fig.8 Erosion profile and gas-liquid distribution of the pipe bend based on the VOF and DPM coupled models

2）分散泡状流中，管内固体颗粒大部分分散在液相中，管道不同截面处的含液率与颗粒含量相关性较大，下游直管段和弯头处，截面含液率较大，受到弯头处流速滞止的影响，固体颗粒的含量也较大。

3）分散泡状流中，固体颗粒在弯头部位受到液相的缓冲作用以及气液分布等影响，固体颗粒对管道冲蚀区域较大，冲蚀最严重位置出现在弯头出口处附近。