在石油与天然气工业中，随着常规油气资源的逐渐枯竭及能源需求量的持续增长，页岩气^[1]、煤层气^[2]、致密油^[3]等非常规油气资源的勘探与开发对保障我国能源安全，实现“碳中和”目标具有重要的战略意义^[4]。水力压裂是非常规油气开采的核心技术，有利于提高油气产量及采收率^[5]。压裂过程中，压裂液携带固体支撑剂通过高压管汇输送至井筒，在深部低渗透底层中形成人工微裂纹。高压管汇作为压裂作业的核心设备，不仅处于高达140 MPa的高压环境中，而且受到压裂液中固体颗粒的冲蚀破坏，从而引发管道的破裂刺漏，甚至导致更大的安全事故和经济损失。所以，对高压管汇的失效分析显得尤为重要。地面高压管汇主要由弯头、三通、直管、单流阀和控制阀等组成，而三通结构的失效占高压管汇总失效的比例较高^[6,7]。三通是高压管汇中改变流动方向的重要部件，同时承受高压、固体支撑剂的冲击及压力波动等作用，主要失效形式是固体颗粒的冲蚀。所以，有必要对三通管件的冲蚀进行研究，现有研究方法主要有实验和数值模拟方法两种。Hong等^[4]利用混合模型、离散相模型(DPM)和冲蚀模型对压裂过程中盲三通、90°弯头、直角管的冲蚀规律数值模拟，结果表明：管壁最大冲蚀率随流体流速的增加呈指数增长趋势；相同工况下，盲三通管道的抗冲蚀性能优于90°弯头和直角弯管。Shan等^[8]通过数值模拟、实验和支持向量机模型(SVM)研究了T型三通管道中颗粒的运动轨迹及冲蚀情况，结果表明：管道的冲蚀与含水率、含气量和含砂量成正比。Zhang等^[9]利用DPM模型数值模拟高压管汇T型、Y型和歧型三通的冲蚀规律，结果表明：三通结构分流时的高压区主要集中在接头处和下游，合流时高压区主要集中在下游；Y型三通冲蚀最严重，歧型三通的冲蚀程度最低。Khan等^[10]利用数值模拟方法分析油气生产过程中T型管、45°Y型管、30°Y型管和15°Y型管的冲蚀性能，结果表明：颗粒在Y型管中相交处撞击并反弹一次，而在T型管中观察到2~3次，从而导致更大的冲蚀率。祝效华等^[7]采用DPM模型研究高压管汇Y型三通和歧型三通在液固两相流中的冲蚀规律，结果表明：两种结构中冲蚀最严重的部位均在相交处，Y型三通的最大冲蚀率随空间夹角的增大而增大，歧型三通在空间夹角为60°时冲蚀最为严重。黄宝华等^[11]利用稠密离散相模型(DDPM)模拟高压管汇岐型三通的冲蚀磨损机理，结果表明：利用DDPM模型模拟精度优于DPM模型。

以上基于三通管道的冲蚀研究均在常压状态下进行，事实上，水力压裂过程中高压管汇在高脉动压力作用下受到冲蚀磨损和拉伸疲劳的共同影响，所以实际工况中，高压管汇不可避免地受到波动载荷作用。虽然Wang等^[12]提出一种考虑静载作用的冲蚀方程，并分析了弯管在高压作用下的冲蚀机理。Yang等^[13]提出适用于静载作用下的冲蚀磨损模型并利用该模型数值模拟不同工况下压裂管道弯头的冲蚀行为，但所考虑的高压情况为静载，而非动载。Dai等^[14]虽提出一种动载作用下的冲蚀磨损模型并数值模拟不同流速、内压、粒径和浓度的颗粒对双弯头的冲蚀作用，但所提出的冲蚀模型中的碰撞角函数在拟合高碰撞角时存在峰值且未将粒径对冲蚀的影响考虑在内。以上学者仅研究载荷作用下弯头、双弯头结构的冲蚀行为，未见用于流动复杂的三通管道中。

因此，本文通过动载作用下的冲蚀试验，提出一种考虑粒径且修正碰撞角函数的新的冲蚀模型，并利用用户自定函数(UDF)嵌入数值模拟软件中。采用DPM模型数值模拟不同空间角度的三通结构在不同内压、流速和颗粒质量流量下的冲蚀规律。所得结论可为高压管汇的安全运行提供理论基础。

试验在实验室自建的液固两相射流冲蚀试验台上进行，该装置可对试样施加静载或循环动载作用，装置示意图如图1所示，主要由动力部分、液固混合及循环输送部分、射流冲蚀部分、伺服加载部分和测控部分组成。泥浆泵为整个系统提供动力源，混砂罐通过带叶片的搅拌电机将液相水和固相砂以一定的比例混合成均匀的水砂两相浆体。试验过程中，泥浆泵将混砂罐中的浆体通过循环管道输送至冲蚀腔并以射流形式对试样进行冲击。冲蚀腔主要由喷嘴、伺服缸、拉力传感器、试样夹头、支撑筒等组成，伺服缸与试样左夹头连接，拉力传感器与试样右夹头连接，可通过调节试样夹头的方向调整试样和喷嘴的喷射角度。伺服液压系统主要由高低压溢流阀、电机、齿轮泵、蓄能器、伺服阀、力传感器等组成。测控部分可对试验过程中的数据进行监控，并通过软件控制加载应力的幅值和频率。冲击试样后的浆体重新返回混砂罐中，实现循环流动。为保证试验准确度，一次冲蚀试验之后需更换砂水。

试验中液相为水，固相为陶粒支撑剂。试样材料为35CrMo，其基本参数见参考文献[15]。试样按照GB/T 3075-2008设计加工，其结构如图2所示。喷嘴的直径为8 mm，与试样的距离为20 mm。

试验过程中，伺服液压系统可对试样施加动载拉应力，实现动载与冲蚀的耦合。应力公式为：

式中，σ为拉伸应力，MPa；σ_m为平均拉应力(0~500 MPa)；σ_a为应力幅值(50 MPa)；t是加载时间，s；φ是波动应力的初相位。

试验过程中，施加动载的频率为5 Hz，冲蚀角度为15°、30°、45°、60°和90°，冲蚀速度为7.5、10、12.5、15、17.5和20 m/s，固相浓度为1%、2%、5%、10%和15%，在以上条件下对试样施加 式(1)的拉应力作用，进行冲蚀试验研究。

基于上述试验，在E/CRC冲蚀模型的框架上重新拟合^[16]。该模型已用于多种试验的冲蚀磨损中^[17]，并且广泛用于弯管、三通等管件的冲蚀预测^[18]。模型方程如下：

式中，ER表示冲蚀率；BH为靶材的布氏硬度，HB；F_s为颗粒形状系数，球形颗粒为0.2，半球形颗粒为0.53，不规则颗粒为1；μ_p为颗粒的速度；n为速度指数(n = 2.41)；θ为碰撞角，rad；f(θ)为碰撞角函数；C为常数(C =2.17 × 10^-7)。

无应力情况下，在 式(2)的基础上结合试验数据构建以下冲蚀公式：

式中，d_p为颗粒粒径，c_p为颗粒浓度，f(d_p)为粒径影响系数，f(c_p)为颗粒浓度影响系数，n′为修正后的速度指数；假定颗粒为均匀球形，故F_s = 0.2；ER₀为无应力下的冲蚀率；f′(θ)为修正后的碰撞角函数，修正为与实际碰撞角函数变化趋势相符合的四次函数 式(4)的形式。

f(d_p)和f(c_p)表示为：

将流速、碰撞角、粒径和浓度进行多元非线性回归分析得到无应力状态下的冲蚀方程。其中，C =1.28 × 10^-7，n' = 2.043，f’(θ)、f (d_p)和f (c_p)中的常数项数值分别为R₁ = 3.5、R₂ = -6.28、R₃ = 4.05、R₄ = -0.833、a₁ = 15.868、a₂ = 2.99、a₃ = 3.4、b₁ = 10.76、b₂ = -0.745和b₃ = 0.025。

由于高压管汇工作过程中不可避免的产生应力，所以在 式(3)的基础上引入应力加速项F_σ 来考虑应力对冲蚀的影响，应力条件下的冲蚀率ER _σ 表示为：

其中，f_σ 为应力加速因子；f_σ (θ)、f_σ (d_p)和f_σ (c_p)为应力情况下的角度、粒径和浓度修正系数，表示为：

通过多元非线性回归分析可得常数项的具体数值分别为g₁ = 1.17、g₂ = 7.13 × 10^-6、g₃ = 15.865、h₁ = 0.033、h₂ = 0.887、h₃ = 3.184、w = 0.007、d₁ = 1.401、d₂ = 0.339、d₃ = 3.239、c₁ = -0.803、c₂ = 0.017、c₃ = 1.427和c₄ = 0.007。

图3为不同流速及角度下试验结果(点值)与模型值对比。由图可知，新建立的冲蚀模型与试验值的误差较小，可用于冲蚀模拟的计算。

本次工作中，采用基于Euler-Lagrange方法的DPM模型，并考虑颗粒相和连续相的双向耦合作用，模型方程如下。

液相视为连续相，采用雷诺平均Navier-Stokes方程求解。连续相方程和动量方程如下^[19,20]：

式中，ρ是流体密度，kg/m³；u是速度分量，m/s；p是压力，Pa；μ是液相粘度，Pa⋅s；μ_t是湍流粘度，Pa⋅s；ρg_i 是重力，N；Su_ip是由于固相而增加的动量源项。因Reynolds平均Navier-Stokes方程中引入了Reynolds应力项，所以需要使用湍流模型来封闭上述方程^[21]。因RNG k-ε模型对于强旋流流动时精度较高，故采用RNG k-ε模型，其方程如下：

式中，C_ε1，C_ε2，σ_k，σ_ε 和β均为常数。

利用牛顿第二定律对颗粒进行追踪，其控制方程如下^[22]：

式中，F_D、F_B、F_V、F_g和F_P分别为颗粒所受的曳力、升力、虚拟质量力、重力和压力梯度力。

压裂作业中常用的三通主要有T型、歧型和Y型三通3种，本文主要对T型三通和不同角度的歧型三通模拟分析。利用Solidworks三维建模软件建立如图4所示的T型三通及不同空间角度α下的歧型三通结构，α的角度分别为30°、45°、60°和75°，其中T型三通是α = 90°时特殊的歧型三通。三通结构垂直放置，重力方向在支流方向上。三通的管径D = 47.8 mm，为使流动充分发展，三通出入口的长度设为5D = 239 mm，按照分流和汇流并结合实际工况要求可将三通在两入一出的工况下设定两个入口I，一个出口O，箭头方向代表流体流动方向。该结构为高压管汇常用材料35CrMo高强度合金钢，密度为7850 kg/m³，Poisson比为0.279，布氏硬度为229 HB^[23]。

以T型三通和α = 45°的歧型三通为例，网格划分情况如图5所示。采用多面体网格，并添加边界层对流动的影响，网格skness小于0.5，质量较好。为提高计算精确性，减小网格数量对模拟结果的影响，以无应力条件下的最大冲蚀率为指标对T型三通结构的不同网格进行表1所示的网格独立性验证。结果表明：当网格数量大于26万时，误差较小，最大冲蚀率趋于稳定，故选择网格数量为268716的网格。

数值模拟过程中，采用基于压力求解器的半隐式算法，动量和湍动能采用二阶迎风模式。模拟过程中连续相为水，密度为998.2 kg/m³；离散相为陶粒支撑剂，密度为1800 kg/m³，粒径为400 μm，从入口表面垂直均匀射入。连续相边界条件设置中，入口为速度入口，出口为压力出口，壁面为标准壁面函数；颗粒相出入口边界条件为逃逸(escape)，壁面为反射(reflect)，反射角函数采用Grant和Tabakoff^[24]提出的壁面反弹模型。假设颗粒相和连续相的入口速度相同，并添加随机游走模型对颗粒运动的影响。冲蚀模型采用UDF中DEFINE_DPM_EROSION宏嵌入数值计算。

进行数值模拟之前，对模型的有效性验证是必要的。以试验中喷嘴冲蚀平板试样为数值模型，将试验过程中的冲蚀率转化为数值模拟软件中单位面积在单位时间内被冲蚀掉的质量，为能提前预知冲蚀，数值模拟过程中选用试样在不同碰撞角下的最大冲蚀率，如图6所示。由图可知，试验值和模拟值在不同碰撞角下的冲蚀率变化趋势相同，且数值模拟的最大冲蚀值始终大于试验过程中的冲蚀值，可预防因预测冲蚀率小于实际冲蚀率而带来的危险，两者的误差为10.09%，可用于冲蚀模拟的预测。

以流速7.5 m/s、固相颗粒质量流量0.1 kg/s、大气压出口为固定条件，对不同三通结构的流场进行分析。图7为该条件下不同结构的压力场分布。由图可知，不同三通结构的两个入口处压力较大，分布较为均匀，随着两分流流体汇入主流，流动方向的改变使压力场突变，从主流区流出之后，压力逐渐平稳。空间角度30°、45°、60°、75°和T型三通的压降分别为0.323、0.333、0.356、0.418和0.464 MPa，可以看出，随着三通空间角度的增加，管内流体的压降逐渐增大，直至成为T型管时，压差达到最大。这是因为角度越大，流体流动方向变化越大，需要更大的能量完成流体携砂方向的转变，所消耗的能量越多，压降越大。

图8为该条件下不同结构的速度场分布。由图可知，当流体从两分流入口流入主流时，主流区域流量增大，而通过相同管径的管道时需要更大的速度，从而流速增加。在流动惯性力、离心力和压力梯度力等的共同作用下，主管上部形成一个小的低速区，出现小范围的二次流动^[12]。空间角度30°、45°、60°、75°和T型三通的最大流速分别为22.00、22.02、22.33、23.21和23.70 m/s，可以看出，随着三通空间角度的增加，流速逐渐增加。结合图9所示的不同截面的速度迹线图可知，随着三通角度的增加，二次流作用逐渐增强，迪恩涡流线越来越密集，尺寸越来越大，湍流作用逐渐增强，速度逐渐增大。从横向看，沿着流体流动方向，迪恩涡逐渐加大，二次流作用加强，故需一段较长距离的管道使流动稳定。

图10为利用 式(3)数值模拟得到的无应力作用下最大冲蚀率随三通空间角度的变化。由图可知，随着空间角度的增加，最大冲蚀率增大。但空间角度小于60°时，最大冲蚀增长率较小，60°~75°之间存在一拐点值使增长率变大，从而使大空间角度下的最大冲蚀率增大。图11为不同空间角度三通管道的冲蚀率分布，可直观反映图10所示变化规律。由图可知，当空间角度较小时，冲蚀分散分布在三通交汇处，冲蚀率较小，随着三通空间角度的增加，冲蚀率增加，冲蚀位置逐渐向三通相贯线上移动，直至α = 75°和T型管道时，冲蚀几乎集中于主管和支管交汇的相贯线上，造成冲蚀率的较大增加。结合图12所示的颗粒轨迹图可知，随着空间角度的增大，颗粒流速增大，碰撞能量增强，故冲蚀率增大。

以流速7.5 m/s、固相颗粒质量流量0.1 kg/s为固定条件，数值模拟不同三通结构在不同内压下(25、50、75、100和125 MPa)的冲蚀变化，如图13a所示。由图可知，相同内压下，不同结构的最大冲蚀率变化遵循图10所示规律；相同三通结构下，最大冲蚀率随内压的增加逐渐增加，当三通空间角度较小时，随内压增加的最大冲蚀率增速较缓；空间角度较大时，增速较快，T型管时增速最大。以最大冲蚀率变化较大的T型三通为例，不同内压下的冲蚀率云图如图13b所示，可直观反映图13a中的变化规律。由图可知，随着内压的增加，冲蚀形状及位置无明显变化，但最大冲蚀率值有所增加。

以固相颗粒质量流量0.1 kg/s，出口压力100 MPa为固定条件，数值模拟不同三通结构在不同流速下(10、15、20和25 m/s)的冲蚀变化，如图14a所示。由图可知，随着流速的增加，动能增大，流体携带的颗粒速度增大，对壁面的冲击力增强，冲蚀率呈指数增加。另外，流速增加导致的湍动能增加，加剧了颗粒的碰撞，也会导致冲蚀率的增加。当三通空间角度较小时，随流速增加的最大冲蚀率增速较缓；大于75°时，增速较大。图14b为不同三通结构在不同流速下的冲蚀率分布，直观反映图14a所示变化规律。从横向看，不同三通结构下的最大冲蚀率均随流速的增加而增大，冲蚀区域逐渐变广；从纵向看，不同结构在相同流速下的最大冲蚀率变化遵循图10所示规律。

以流速7.5 m/s，出口压力100 MPa为固定条件，数值模拟不同三通结构在不同颗粒质量流量下(0.2、0.3、0.4和0.5 kg/s)的冲蚀变化，如图15a所示。由图可知，随着颗粒质量流量的增加，颗粒数量增多，颗粒与壁面之间的碰撞加剧，最大冲蚀率线性增加。不同结构下最大冲蚀率的增长速率不同，三通空间角度越大，增长速率越快。图15b为不同三通结构在不同颗粒质量流量下的冲蚀率分布，直观反映图15a所示变化。从横向看，随着颗粒质量流量的增加，三通结构的最大冲蚀率增大，冲蚀面积增加；从纵向看，不同结构的最大冲蚀率变化遵循图10。

(1) 根据动载作用下的液固两相射流冲蚀试验，提出一种考虑粒径且修正碰撞角函数的冲蚀模型，该模型与试验值拟合优度较好，可用于数值模拟冲蚀的计算。

(2) 不同三通结构随空间角度的增大，压降、流速及冲蚀率均增加。当空间角度较小时，冲蚀分散分布在三通交汇处且增长速率较低；空间角度大于75°时，冲蚀集中于主管和支管交汇的相贯线上，造成冲蚀率的骤增。因此，三通管汇服役过程中应尤其关注三通相交处的失效破坏，预防危险的发生。

(3) 不同三通结构的最大冲蚀率随内压的增加而增加；随流速的增加影响较大，呈指数形式增加；随颗粒质量流量的增加线性增加，且冲蚀形状及位置无明显变化。三通空间角度越大，增长速率越快。三通结构最大冲蚀率随空间夹角的增大而增大这一规律不随内压、流速及颗粒质量流量的变化而变化。因此，设计三通结构时，在满足基本管道要求的前提下，应尽可能减小三通的空间角度。