目前海底管道埋设工程中应用最广泛的是喷冲式海底挖沟，利用高压水射流对海底泥面进行喷冲，实现基础开挖和埋管。这类开沟机不仅可用于管道铺设阶段，还可用于管道维修阶段的开挖作业，且对土壤适用性也强。足够的开挖深度和宽度是管道埋设作业的重要指标，为了获得开挖后的泥沙床面形态，需要对射流冲刷、泥沙运动进行深入分析。

AKASHI等^[1]纪录了不同冲刷时刻的冲刷坑形态；ADERIBIGBE等^[2]、CHAKRAVARTI等^[3]及HOU等^[4]分析了射流速度、喷嘴直径、距泥面高度、射流角度、泥沙粒径等不同参数对于冲刷效果的影响。由于涉及到泥沙和水流2种不同物质形态，数值模拟这类问题难度较大。目前的数值研究方法主要分为只计算水流运动继而建立泥面变形公式的间接法和同时计算水流、泥沙运动的直接法。间接法中利用泥沙输运方程计算出不同水流条件下床面泥沙的位置，并改变网格生成新的计算域。该方法适用于几何模型简单、冲刷不剧烈的情况^[5-6]。直接法还可进一步分为连续介质多相流法^[7-8]和粒子法^[9-10]。连续介质多相流法将泥沙视为连续介质，计算效率较高、能满足工程快速预报。在多相流模型里，泥沙相本构模型有粘性模型^[7]、粘弹性模型^[8]、弹塑性模型^{[9, 11]}和粘塑性模型^[12-13]等。VAN等^[8]和LAGRÉE等^[14]提到，由于泥沙相被视为流体，数值计算中会产生不符合物理现象的情况，如水沙界面波动、泥沙运动不停止的现象。

本文之前采用简易的数值方法对一套配置抽吸装置的喷射式挖沟机进行了计算并指导了工程设计^[15]。本文进一步基于Fluent软件，选用欧拉多相流法，结合泥沙起动、边坡稳定理论，利用用户自定义函数(user-defined function，UDF)修正动量源项，定量计算了冲刷坑形态，实现了合理的二维泥沙冲刷数值计算。

1 欧拉多相流模型及其修正模型 1.1 控制方程

欧拉多相流模型是对泥沙、流体相建立各自的控制方程，求解得到两相的流场特征：

1) 连续方程。

$ \frac{\partial }{{\partial t}}({\alpha _{\rm{q}}}{\rho _{\rm{q}}}) + \nabla ({\alpha _{\rm{q}}}{\rho _{\rm{q}}}{\mathit{\boldsymbol{v}}_q}) = 0 $

(1)

式中：q为f、s，分别代表流体相和泥沙相；α_q为各相的体积分数；ρ_q为各相密度；v_q为各相速度。

2) 动量方程。

水流相：

$ \begin{array}{l} \frac{\partial }{{\partial t}}({\alpha _{\rm{f}}}{\rho _{\rm{f}}}{\mathit{\boldsymbol{v}}_{\rm{f}}}) + \nabla ({\alpha _{\rm{f}}}{\rho _{\rm{f}}}{\mathit{\boldsymbol{v}}_{\rm{f}}}{\mathit{\boldsymbol{v}}_{\rm{f}}}) = - {\alpha _{\rm{f}}}\nabla p + \\ \;\;\;\;\;\nabla \cdot{\mathit{\boldsymbol{\tau }}_{\rm{f}}} + {\alpha _{\rm{f}}}{\rho _{\rm{f}}}\mathit{\boldsymbol{g}} + {K_{{\rm{sf}}}}({\mathit{\boldsymbol{v}}_{\rm{s}}} - {\mathit{\boldsymbol{v}}_{\rm{f}}}) \end{array} $

(2)

泥沙相：

$ \begin{array}{l} \frac{\partial }{{\partial t}}({\alpha _{\rm{s}}}{\rho _{\rm{s}}}{\mathit{\boldsymbol{v}}_{\rm{s}}}) + \nabla ({\alpha _{\rm{s}}}{\rho _{\rm{s}}}{\mathit{\boldsymbol{v}}_{\rm{s}}}{\mathit{\boldsymbol{v}}_{\rm{s}}}) = - {\alpha _{\rm{s}}}\nabla p - \Delta {p_{\rm{s}}} + \\ \;\;\;\;\;\nabla \cdot{\mathit{\boldsymbol{\tau }}_{\rm{s}}} + {\alpha _{\rm{s}}}{\rho _{\rm{s}}}\mathit{\boldsymbol{g}} + {K_{{\rm{fs}}}}({\mathit{\boldsymbol{v}}_{\rm{f}}} - {\mathit{\boldsymbol{v}}_{\rm{s}}}) + {\mathit{\boldsymbol{S}}_{{v_{\rm{s}}}}} \end{array} $

(3)

式中：p为水压力；p_s为泥沙相的压力，当泥沙相的体积分数α_s小于其最大允许值时计算；K_sf=K_fs为相间动量交换系数；S_vs为源项；τ_q为切应力张量，表达式为：

$ \begin{array}{l} {\mathit{\boldsymbol{\tau }}_{\rm{q}}} = {\alpha _{\rm{q}}}{\mu _{\rm{q}}}(\nabla \cdot{\mathit{\boldsymbol{v}}_{\rm{q}}} + \nabla \cdot\mathit{\boldsymbol{v}}_{\rm{q}}^T) + \\ \;\;\;\;\;\;\;{\alpha _{\rm{q}}}\left( {{\lambda _{\rm{q}}} - \frac{2}{3}{\mu _{\rm{q}}}} \right)\nabla \cdot{\mathit{\boldsymbol{v}}_{\rm{q}}}\mathit{\boldsymbol{I}} \end{array} $

(4)

式中：μ_q为各相剪切粘度；λ_q为各相体积粘度。对于泥沙相，其颗粒的剪切黏度表达式为:

$ {\mu _{\rm{s}}} = {\mu _{{\rm{s, col}}}} + {\mu _{{\rm{s, kin}}}} + {\mu _{{\rm{s, fr}}}} $

(5)

式中λ_s=0。式(5)分别考虑了碰撞、运动及摩擦的影响。

紊流模型采用Standard，近壁面用标准壁面函数处理，参数设置为软件默认k-ε数值。压力速度耦合迭代采用SIMPLEC算法求解。动量、湍流粘度、湍流动能和湍流耗散率采用二阶迎风格式进行离散，体积分数用一阶迎风格式离散。

数值模型的边界条件为：喷嘴处为速度出口，水沙界面设为质量动量交互的交界面，距离喷嘴很远的上部水面为压力出口，其压强为标准大气压，其余为壁面边界条件。

1.2 泥沙模型 1.2.1 起动模型

目前判别泥沙起动条件的理论主要有以下3种：起动拖曳力、起动流速和起动功率。其中泥沙开始运动时床面水流切应力称为起动拖曳力τ_c，是最本质的描述泥沙起动条件指标。Shields根据4种不同密度的泥沙试验结果得到了泥沙雷诺数Re_*和τ_c/[(γ_s-γ)D]之间的关系，根据这些数据以及其他研究资料，可以得到一条平均曲线，称为Shields曲线^[16]。已知泥沙颗粒的物性，可根据Shields曲线，计算得到泥沙起动时的起动拖曳力。

数值计算里提取流体内部的速度值比较方便。近泥沙床面，切应力和临界剪切流速U_*c之间的关系可表示为τ_c=ρU_*c²，可使用U_*c作为判断泥沙运动状态的指标，具体算法为

1) 泥沙堆积区：当泥沙的体积分数超过摩擦堆积限值α_sf且水流速度U_w小于临界剪切流速U_*c时，泥沙没有起动或者由于颗粒间摩擦作用力增大而停止运动，泥沙相呈固态特征，无流动；

2) 泥沙起动区：当泥沙的体积分数超过摩擦堆积限值α_sf且水流速度超过临界剪切流速U_*c时，泥沙起动，呈流体特征，产生流动；

3) 泥沙流动区：当泥沙的体积分数小于摩擦堆积限值α_sf时，此时泥沙已被水流卷起，悬浮在水流中，泥沙相间作用力以碰撞为主，为流动状态。

图 1为泥沙的运动区域示意图。通过上述模型的定义，能描述泥沙相的固相特性，泥沙床面不会产生不合理的波动和回填。

1.2.2 水下休止角模型

水下休止角φ是指静水中堆积的泥沙由于摩擦力作用，形成稳定而不坍塌斜坡面的坡角，是泥沙颗粒的重要参数。当堆积角度大于休止角时，泥沙颗粒会失去平衡下滑，直到堆积角度小于或等于休止角。

冲刷坑边界一侧为剧烈变化的射流水流和反射水流区域，受其影响短时间内冲刷得到的坡面角度可能会大于休止角；另一侧主要为泥沙堆积区，泥沙呈固态特性。由于受起动模型作用，冲刷坑边界局部泥沙的速度被限制住，随着计算时间的演变，冲刷坑局部边坡会变陡，如图 2(a)所示。即使在射流停止后，这种陡坡仍然存在。这种不合理的局部现象是由所施加的起动模型造成的，因此应予以修正。

与动网格的沙滑模型不同，欧拉模型的网格固定，故无法通过控制局部网格节点运动来实现坡度的调节。本文通过对界面陡坡处的泥沙施加一定的扰动力，使这部分泥沙滑落直至坡面角度不超过休止角。由于冲刷坑边界两侧的泥沙浓度变化很大，并参照动量方程中重力项的表达式，本文定义的扰动力大小与界面泥沙体积分数的差值、泥沙密度和重力加速度成正比：

$ \left\{ \begin{array}{l} {S_{{\rm{vs}}, x}} = {k_1}\Delta {\alpha _{{\rm{s}}, x}}{\rho _{\rm{s}}}g\\ {S_{{\rm{vs}}, y}} = {k_2}\Delta {\alpha _{{\rm{s}}, y}}{\rho _{\rm{s}}}g \end{array} \right. $

(6)

通过调整最终选定参数k₁=1.15，k₂=0。经过休止角模型的修正，冲刷坑边坡的形态得到改善，不再有不符合物理现象的局部陡坡，如图 2(b)所示。

1.2.3 数值实现

起动模型里，堆积区的泥沙为固体状态，其速度设定是通过UDF修改泥沙相动量方程源项来实现的。为了数值求解的稳定性，源项S_vs表达式一般写成线性形式，即S_vs=A+Bv_s，A为显式部分，Bv_s为隐式部分。通过UDF定义隐式项dS/dv_s为无穷大的数值，由于在方程离散的时候这项归类在dv_s/dt相关的项中，受其影响dv_s/dt≈0，使得堆积区泥沙保持静止。源项S_vs数值设为零，没有引入其他力的作用。

实现步骤为：

1) 在每个时间步计算前，判别当前网格的泥沙体积分数、水流的速度。若为泥沙堆积区的网格，则设置泥沙相的速度为零;

2) 继而进行控制方程的计算，在泥沙相的动量方程中区分出泥沙堆积区的网格，通过定义dS/dv_s为无穷大数值，保持堆积区泥沙静止。

通过UDF遍历水沙交界面处的网格单元，计算周围网格泥沙所形成的坡度。当坡度超过休止角时，相应网格单元内的泥沙相受力下滑直至满足平衡条件停止运动。

为了简化计算，做如下假设：①水沙交界面处，与水平轴夹角等于泥沙休止角的直线上所有点的泥沙体积分数相同，图 3(a)为其示意图；②泥沙体积分数为非零的节点，在其半径小于单个网格长度范围内，各方向的泥沙体积分数为线性变化。

具体实现步骤为：

1) 搜寻水沙交界面网格，定义网格P为当前计算网格，其相邻的四个网格分别为W、E、S与N，如图 3(b)所示；

2) 以泥沙堆积坡面斜率为负值时为例，网格P、E、N中心连线的夹角为θ。计算网格P与North中心连线上的P_n点，使其满足与网格E、P中心形成的角度为泥沙休止角。由假设①可知，此时P_n的泥沙体积分数为α_{s, pn}=α_{s, east}；

3) 由假设②，根据网格N中心、P_n点的泥沙体积分数，计算满足休止角度模型时网格P的泥沙体积分数α_{s, 2}；

4) 比较网格P此时的泥沙体积分数α_{s, p}与α_{s, 2}的值：当α_{s, p}＞α_{s, 2}时，说明该处泥沙堆积过多，会使堆积坡度大于休止角，故对网格P的泥沙相动量方程施加扰动力，使泥沙沿着坡面下滑直至坡度小于休止角；当α_{s, p}≤α_{s, 2}时，说明该处泥沙的堆积量正常，不需要施加扰动；

5) 当泥沙堆积斜率为正值时，则计算对象为网格P、W与N，分析过程类似。

通过对水域、沙床的网格进行搜索判别，将泥沙起动、休止角模型载入数值程序的动量方程，对符合条件的泥沙相约束速度或施加扰动力。单个网格的判别流程如图 4所示。

2 轴对称、平面射流冲刷的计算分析 2.1 轴对称冲刷

Aderibigbe等^[2]对水下淹没单喷嘴定点射流冲刷进行过系列实验，本文选择其中的一例作为参照实验。实验为三维模型，根据几何、物性的特点，可以简化为二维轴对称问题进行数值计算，物理模型示意图如图 5所示。

实验初始泥沙厚度为150 mm，喷嘴距泥沙床面高度h₁为103 mm。喷嘴直径d₁为4 mm，射流速度U₁为2.74 m/s。采用疏松非粘性均匀泥沙，其平均粒径D₁为0.88 mm，密度为2 650 kg/m³，孔隙率为0.25，水下休止角为36°。

计算模型为半模，网格数量为2万，计算了历时为60 s的冲刷过程，此时冲刷坑的坡面形状变化已不明显。图 6(a)为数值和实验的冲刷坑形态比对图, 由于Aderibigbe等的实验数据主要集中于坑内(z＜0), 所以主要进行坑内的形态对比。由图 6可知，数值实验获得的冲刷深度为16.5 mm，冲刷坑最大宽度为66 mm，均与实验结果非常接近。

图 6(b)、(c)为冲刷时间t=60 s时喷嘴-冲刷坑附近水流相的速度场。冲刷坑形态以泥沙相体积分数云图来表示。从图 6(b)中可见，水流从喷口垂直射向沙床, 在泥沙床面处受到撞击作用，水流被反射, 沿着水沙交界面向两侧运动。在入射流与反射流混合区域, 形成平行于沙坑界面的小涡旋, 涡旋中心高于沙丘。从图 6(c)可以看出，水流流速在射流过程中迅速衰减，水沙交界面的速度约为0.05 m/s，泥沙床面内的孔隙水速度低于0.05 m/s，无法让泥沙起动。

图 6(d)、(e)为冲刷达到平衡时冲刷坑附近泥沙相的速度场。图 6(d)显示了水沙交界面附近的泥沙运动，可以看出，冲刷坑内的泥沙在反射水流作用力、重力的影响下，主要沿着坡面运动；越过堆丘后水流流速降低，泥沙在自重影响下作沉降运动。图 6(e)发现交界面悬浮着运动的泥沙颗粒，速度值不超过0.047 m/s，而远离交界面处的泥沙相速度为零，即泥沙未起动保持着静止的固体特性。

2.2 平面射流冲刷

Akashi等^[1]对二维平面泥沙冲刷的历时变化进行了实验研究。实验在水槽中进行，喷嘴宽度为20 mm，流速为0.74 m/s，距离初始床面距离为100 mm。试验采用中值粒径为0.84 mm的非粘性泥沙，其密度为2 650 kg/m³，孔隙率为0.4，水下休止角为31°。

实际埋管工程中，射流是移动的，以便在海床上开挖出管道埋设需要的沟槽，因此射流不会在某点处进行长时间的定点喷冲。此外，开沟埋管用的射流流速较高，故数值计算时间变化的冲刷现象时，也主要关注冲刷前期，即冲刷坑变化最为剧烈的阶段。本文中对平面射流冲刷模型数值计算了历时8 s的冲刷过程，取t为2、3.5和8 s这3个时刻来进行分析。图 7(a)为各时刻的冲刷坑形态。起动泥沙在水流作用下向冲刷坑两侧运动，随着流速降低泥沙沉降，在冲刷坑两侧形成堆丘。冲刷过程中，冲刷坑深度和宽度逐渐增加，最大冲刷深度分别为5、13、22 mm；堆丘高度和宽度也逐渐增加，至t=8 s时，堆丘高度达到约15 mm。

数值计算中，由于计算模型对称性，在原点处的水流条件无切向分量，故原点附近的冲刷深度略低于最大冲刷深度。实验中存在许多其他因素影响，所以该现象不太明显。

图 7为无因次化的冲刷坑形态，其中横坐标x和纵坐标z分别用最大冲刷坑宽度L_m和最大冲刷坑深度H_m作为无因次化的特征尺度参数。和文献[1]实验比较发现：在冲刷坑内(z≤0 m)，数值计算值和实验值吻合度较好，最大冲刷深度非常接近，能很好描述冲刷坑内形态；在堆丘部分(z＞0 m)，t=2 s时数值模拟得到的堆丘高度和宽度比实验值低；而t为3.5、8 s时，数值模拟得到的堆丘宽度小于实验值，堆丘高度相近。

数值计算中，冲刷坑内形态模拟吻合较好，说明本文模型能合理描述泥沙冲刷起动这一过程，但数值模拟显示，起动后的泥沙大量悬浮在流体中，导致模拟的堆丘形态偏小。随着计算时间增加，越来越多的泥沙沉降下来，堆丘形态愈趋合理，如图 7(c)、(d)所示。

图 8(a)、(b)显示了t为8 s时的水流相的速度场，其中图 8(a)为射流域的流速分布，水流到达坑底沿着冲刷坑向两侧流动，在坑底的流速约为0.05 m/s；图 8(b)为泥沙中的孔隙水流速分布，水沙交界面附近速度值不超过0.04 m/s，在远离水-沙交界面处，孔隙水速度衰减，不足以让泥沙起动。

图 8(c)为t=8 s时的泥沙相的速度场。对应图 8(b)的水流相速度场，远离交界面的泥沙相静止无流动呈现固态特性。在两侧堆丘上方可以看到有部分泥沙以较高的速度随水流相悬浮运动，这部分悬沙随着水流作用减弱将会沉积于床面。

3 结论

1) 本文改进的模型解决了泥沙相在数值计算中产生的非物理的流动问题，能正确将其固相特征反映出来。改进模型不仅仅在模拟冲刷坑形态上较为准确，并且对泥沙起动历程的模拟也有一定的准确性。

2) 本文在数值模拟中，为防止出现陡峭的堆丘边坡，提出对水沙交界面处泥沙施加一定的扰动力，从而形成稳定的堆丘水下休止角。而扰动力的形式具有一定的经验性，其合理的表达形式，对于其他问题的普适性还有待商榷。

埋管工程应用中遇到更多的是移动射流冲刷问题，后续研究将对二维移动射流冲刷现象进行数值模拟，观察分析冲刷坑演变及泥沙运动规律。