彭 涛， 张菲菲， 方含之

(长江大学石油工程学院， 武汉 430100)

地下储气库工艺作为季候调峰和能源储备基础设施，其完整性研究对天然气的高效、安全、经济生产有着显著的意义[1]。储气库工艺有其特殊性：井筒多处于交变载荷工况；井筒中流体系持续性周期流动、流量大、作业压力大；储气遇水易酸化；砂岩地层(弱胶结岩层)改造时，持续拉伸作用等易导致砂砾弹塑性变化和储层完整性下降，进而加剧井筒内出砂，严重出砂更会导致砂埋气层，钻柱磨损失效，降低生产压力及产量[2-3]。因此，储气库储层出砂问题，亟须各研究者加以重视。

关于研究固体颗粒对筒壁和设备的冲蚀作用，中外已有大量实验分析和报告，主要依托Magneé[4]结合材料性质、微观结构、冲击角或硬度等建立的广义侵蚀磨损定理或其延伸。而对一般弯管下，影响因素主要包括弯管曲率、弯管直径、气相流速、颗粒密度直径等，其中最关键的是颗粒直径，文献[5-8]总结了一般性结论：在同一携砂流速下，颗粒直径越大，磨蚀作用越大。而考虑井筒冲蚀过程中，出砂基岩直径较小且近似，但不同砂砾形状对冲蚀影响研究较少。

因此，现基于实际采气井工况，结合计算流体力学(computational fluid dynamics, CFD)思想及磨蚀工况设定方法，选取适宜湍流模型并建立冲蚀模型，提出不同储层出砂形状特性，及其对井筒磨蚀率的影响。

1 组分运移模型

在流体动力学计算中，应先分析流场内各物质分布、流动形态等。主要考虑储气由储层进入生产管并运移至井筒造斜点的流动特性，其间储气流动机制主要是由压力梯度所引起的黏性流动。因此，针对造斜点间多组分运移，只考虑较为简单的气相携运模型，其固相与气相间并无能量传递，仅存在气相湍流作用携砂和砂砾冲击壁面，并且携砂模型下，也不需考虑各气相间相互吸附作用。对于较小颗粒与壁面间的相互作用，分子力较惯性力占主导，故存在粒子由于凝聚力被壁面吸收而没有反弹或滑移的力学现象。

2 数值模型建立

由于冲蚀问题作为一个长周期性行为，依托较高模拟精度的CFD方法也是分析该类问题的重要途径[9-10]。在考虑数值方法时，不得忽略物质连续变形特性，而针对刚体物质，如砂砾，一般不考虑其变形。将流场内流体做连续相介质(流体各连续相控制方程均依托于连续性方程和N-S方程)。将砂砾作为离散项介质，即以欧拉法形式定义固相微团，通过对各个微团研究以概括整个流场工况，并对各流固耦合力学方程进行数值求解，模型中假设包括：

(1)设置侵蚀模型下粒径函数为常数，取常规设定值1.8×10-9。

(2)不考虑颗粒间相互作用力。

(3)不考虑冲蚀壁面时产生的几何形状变化。

取短曲率水平井造斜段方法：设计造斜率为10(°)/m(曲率半径5.73 m)。旨在考虑携砂气流对管段各位置磨蚀情况的机理性分析，并且若取实际长径比，对于计算机运行能力要求过高，因此依据动力学相似原理，将井筒原型定义为如图1所示的几何模型。

图1 造斜几何模型及网格划分Fig.1 Deflecting geometric model and meshing

前期利用前处理软件绘制网格模型，再建立冲蚀模型进行模拟。计算时考虑稳态、等温及湍流条件。对各方向动量、连续性方程等采用离散形式计算时，考虑到储气可压及建立结构化网格方法，使用二阶迎风格式。选取求解器时，考虑模型运动边界条件设定方法及两相设定形式，取压力分离形式求解器。该求解器下压力-速率耦合形式有以下算法：SIMPLE、SIMPLEC、PISO及COUPLED。韦鲁滨等[11]研究表明，PISO算法对于各变量均进行一次预测及两次修正，其对计算值更为准确。并且该算法对于非耦合形式的动量方程及标量方程(如热边界条件为常数时)，较SIMPLEC及SIMPLE算法效率更高，故采用PISO算法。计算时对描述流场运动状态各控制方程进行迭代计算，得出解析解，进而定性描述流固耦合问题。迭代过程中取平均磨蚀量和各方程残差设定值为收敛依据。

2.1 湍流模型定义

定义湍流模型时，考虑到流场内携砂尺寸较小，并被壁面捕获。故流场内不易产生大涡街，因此不考虑大涡模拟(large eddy simulation, LES)方法，而选取一般雷诺时均法(reynolds-averaged Navier-Stokes, RNAS)方法。故选取Realizek-e模型、RNGk-e及SSTk-ω模型。

图2为3种湍流模型不同流速下沿管壁位置最大磨蚀率验证结果，其中气相速率较大时三种方法对应冲蚀模型分析结果近似，而在流速较小时RNGk-e方法沿程振荡幅度较大，并且进入竖直管段时，幅度最大，因此有结论：RNGk-e及SSTk-ω方法效果均较好，其中SSTk-ω模型综合了k-e模型在近壁区的计算特点，并加入横向耗散项，且在定义湍流黏度时加入剪切应力影响，故更适于黏性流动机制模拟。因此使用SSTk-ω湍流模型。

图2 流速6 m/s和3 m/s时不同湍流模型下的 沿程磨蚀分布Fig.2 The wear distribution of different turbulence models along the path under 6 m/s and 3 m/s

2.2 冲蚀模型定义

基于广义侵蚀磨损定理不同参数设定方式，提出Oka模型[12]、Peng模型[13]和Vieria模型[14]等较好的预测模型,其中对于弯管形式冲蚀模型，Vieria等[14]和刘琦等[15]建立了与实验值更为接近的模型，其基本表达式为

(1)

式(1)中：vp为颗粒相对速度，m/s；b(v)为颗粒相对速度函数(或速度指数函数)；α为颗粒与壁面间冲击角,rad；C(dp)为颗粒粒径函数，设为常数，取1.8×10-9；Aface为壁面面积，m2；mp为颗粒质量流量，kg/s；f(α)为冲击角函数；其中参数C、b以及f均为定义壁面的边界条件，以常数、多项式、分段多项式或分段线性方法等定义。

可通过取适宜的碰撞模型和磨蚀参数，再定义携砂气流初始条件及边界条件。

2.3 壁面函数定义

在提出对于各边界条件设定时，主要定义冲击角函数、反弹系数设定、粒径函数设定(上文已讨论)、速度指数函数等方法设定。考虑冲蚀问题时，颗粒务必贯穿壁面边界层，其边界流动可视为层流，因此对于整个流场设定湍流模型时，需考虑壁面区域。利用无量纲u+及以定义y+边界层中三个子层(黏性底层、过渡层、对数律层)，划分图示如图3所示。

(2)

式(2)中：u为流场时均流速,m/s；Δy为流体至壁面距离,m；ρ为流场混合密度，kg/m3;uτ为壁面摩擦速度,m/s。

黄诗嵬[16]研究表明，对近壁面区域物理量计算方法，一种是采用半经验的壁面函数来过渡壁面与核心湍流区域；另一种是对湍流模型进行修改，对边界网格加密，使其适用于近壁面区域。

未加密之前，壁面y+为100～500。随后，采用自适应形式对壁面网格加密，如图4所示。图示各点为沿管各网格节点y+值最大最小范围，将y+定义到50～400范围值，将网格第一个内节点位置修改至对数律层成立的范围，使其配置到旺盛湍流区域，进而在使用SSTk-ω模型模拟壁面区域运动时更为精确。

图3 壁面三子层划分及对应速度分布Fig.3 Division of three sublayers of wall surface and corresponding velocity distribution

图4 自适应前壁面和后壁面部分y+Fig.4 y+ at adaptive front and rear wall segment

2.4 壁面反弹系数定义

设定壁面反弹系数时，可取葛令行等[17]以依托粒子跟踪测速技术(particle tracking velocimetry，PTV)对碰撞前后速度分解，能够获得切向和法向反弹系数。或采用经典反弹模型(已知入射角和入射速率、反弹速率，利用数学模型计算)定义反弹系数，由于需配套设备要求高，且操作复杂。鉴于实验平台有限，考虑采取指数多项式定义各反弹系数。由于小粒径砂砾法向反弹系数与入射角正相关，且变化明显，故颗粒与延展性壁面碰撞时，各反弹系数较大。取李洁琼[18]对280～400 μm粒子与壁面发生弹性碰撞实验所得各弹性系数多项式，并考虑砂砾入射角度和通过防砂筛管逃逸砂砾直径，进行插值计算，得出反弹系数多项式如下。

切向反弹系数：

et=0.997 6+0.942 1α-1.642α2+0.706 7α3

(3)

法向反弹系数：

en=0.874 1+0.964 1α-1.752 4α2+0.961 8α3

(4)

2.5 速度指数定义

对于一般冲蚀模型，其速度指数函数的设定，主要与冲蚀壁面材质、壁面温度梯度有关，由于一般套管钢质壁面导热系数较大，各处温差幅度小，故壁面温度梯度影响可忽略。对于延展性材料，参考文献[19-20]，冲蚀模型中速率指数为2；Stack等[21]认为速度指数应介于2～2.5，甚至后期已开发模型的速度指数超过3。依托冲蚀模型，并考虑与Stack方法下的模型存在几何相似，取速度指数为2.8。

2.6 冲击角函数定义

对于不同冲蚀模型下冲击角函数定义方法，可通过调研各有关文献进行赋值。例如Tusla Angle Dependent模型[22]提出经典砂砾冲蚀钢制壁面各冲蚀率改写公式和相应冲击角函数。对于侵蚀性磨损如海底管道磨蚀情况，也提出利用Tabak-off和Grant模型[23]。对地层砂砾和井筒生产套管冲蚀模型各参数，取经典模型，以多个线性段方法定义砂砾冲击角函数：

f(θ)=Aθ+Bθ2+Cθ3+Dθ4+Eθ5+Fθ6+Gθ7+Hθ8

(5)

表1 砂砾冲击角取值

3 模型分析

考虑用圆球度表征对不同储气库出砂形状变化，如基岩地层，其圆球度较大，而针对页岩屑，其圆球度较小。圆球度(SF)为描述研究对象边界复杂程度：

(6)

式(6)中：Ap为研究对象表面积，m2；Vp为研究对象体积，m3。

SF值为0～1，越趋近于1，边界过度越光顺，SF=1时，边界区域属圆形；若边界起伏程度越大，其形状系数越小。计算圆球度表面积时，可取樊骐铖[24]提出的像素投影法。对预测模型影响主要依据颗粒边界与壁面碰撞前后运动方程，将各撞击点前后轨迹整合到对应磨损模型，再进行磨蚀预测。

取砂砾密度2 500 kg/m3，防砂筛管网孔尺寸0.23 mm(忽略持续过流冲蚀对网孔尺寸变化率影响)，直径小于0.2 mm的砂砾从进口(inlet)被携入，其质量流量取1 kg/s，进口流速6 m/s；在设置砂砾圆球度时，考虑基岩砂砾、页岩屑设定等，其他设定参数不变，参考实验数据如表2所示。

表2 颗粒形状系数表

在相同进口条件和边界条件下，对于井筒内出砂，流场对砂砾携砂作用最大位置为图5位置1、2处。位置1主要由于砂砾进入造斜管段，受力方向改变，砂砾沿管内部偏移；对于位置2处，曳力完全沿竖直方向且由于涡街作用，砂砾进一步沿内部偏移，导致在进入竖直管段后对壁面磨蚀作用减小。图6和图7为对于不同携砂形状系数设定(即砂砾边界光顺度)冲蚀行为。

图5 砂砾携运模型图示Fig.5 Gravel transport model diagram

图6 不同圆球度冲蚀及砂砾轨迹模型Fig.6 Model of erosion and gravel track with different sphericity

图7 管壁各处最大磨蚀分布Fig.7 Maximum wear distribution on pipe wall

由图6、图7可知，由于进口直管段沙砾平行于壁面流动，进而对于进口段，其冲蚀率较小；并且储气紊流波动对于颗粒路径也有影响，且在弯管处，多存在紊流增强及二次涡街等，因此加大了流场运动复杂性，进而弯管外侧中心处磨蚀最严重；考虑不同形状系数对磨蚀影响，由于边界越粗糙，其流线轨迹分布越复杂，即对造斜段外侧冲蚀分布更广，受曳力较小，在贯穿黏性底层时能量耗损较大或被底层吸收，即颗粒与壁面间相互作用较弱，壁面磨蚀程度较轻。

颗粒圆球度越小，不规则颗粒在流场中受曳力越大，并且最大磨蚀位置也会向上偏移。形状越不规则的颗粒(如形状系数为0.32、0.56时)，其储气携砂能力越强，颗粒在进入竖直管段时，仍存在一定磨蚀率。

对圆球度与磨蚀率进行线性回归分析时，如图8所示，可决系数R2=0.903 5，因此在其他参数变化时，如气相流速，可考虑只取两个参考数据做出相关拟合曲线方程，对其他不同形状系数粒子求解时只需代入单一值进行求解。

图8 圆球度与最大磨蚀率关系Fig.8 The relationship between SF and maximum erosion rate

4 结论

通过以上研究，得出结论及建议如下。

(1)对于一般气固两相流，NRAS湍流模型均适宜。若不过分考虑计算能力，其中SSTk-ω模型更多考虑近壁面计算特性，且结合自适应网格形式时，该模型更适宜作为冲蚀问题湍流计算方法。

(2)形状系数越小的颗粒，由于边界越粗糙，受曳力较小，其流线轨迹分布越复杂，在贯穿黏性底层时能量耗损较大或被底层吸收，即颗粒与壁面间相互作用较弱，壁面磨蚀程度越轻，且该情况下其储气能力越强。

(3)商业软件对于各边界调节设定只能通过设定形式，故对模型修改、二次开发等具有局限性，后期研究可考虑利用如OpenFoam等开源软件进行设定。