黄河口泥沙异重流基本控制参数的数值试验

摘要：为了定量地评价泥沙异重流基本控制参数对泥沙异重流流动的影响，本文设计了一系列数值试验来展示各种参数与泥沙异重流流动特征之间的定量关系。采用ADI法求解三维方程经垂向积分后的平面二维泥沙异重流方程。所选取的基本控制参数包括水下斜坡坡角，总阻力系数，水力挟带系数，泥沙饱和浓度和泥沙沉速。通过数值试验得出水下斜坡坡角θ、泥沙饱和浓度Es越大，总阻力系数cd、泥沙沉速ω越小，相应地泥沙异重流流速、流动厚度以及悬沙浓度越大。其中斜坡坡角对泥沙异重流影响最大。当泥沙异重流流动处于超临界流状态，水力挟带系数Ew大于零，这时由于有一定水量的卷入，所以流动厚度沿程增加。

关键词：黄河口 泥沙异重流 控制参数 数值试验 ADI法

　　泥沙异重流及形成的地质体在沉积学和石油地质中有重要意义。在黄河口，泥沙异重流是泥沙搬运的重要方式，是河口区沉积动力的重要过程，是认识河口三角洲发育演化的重要途径。为使研究深入一步，我们进行了泥沙异重流的数值模拟。近年来，国内对异重流数模研究较少，仅见方春明等(1997)关于泥沙异重流潜入时的立面二维数值模拟以及邱晨霞(1995)对盐水异重流的二维两层数值。

1 平面二维泥沙异重流控制方程及数值模拟

　　目前，包含三个基本控制方程(流体质量、动量守恒方程和泥沙质量守恒方程)的泥沙异重流数值模拟被广泛应用，本文把一维方程推广到平面二维的情形，成为四方程模型，如下所示^[1]

其中u，v分别为x，y方向速度，θ_x，θ_y分别为x，y方向坡度，E_wx与E_wy分别为x、y方向的水力挟带系数，τ_x、τ_y为x方向与y方向的剪切摩擦力，k_x、k_y分别为x方向与y方向的悬沙扩散系数，t为时间，h为异重流厚度，s为深度平均的悬沙含量，ρ_s，ρ_w，ρ_f分别为泥沙颗粒、水体和泥沙异重流的密度，E_s为泥沙饱和浓度，ω为泥沙颗粒的沉降速度，α为泥沙运动饱和系数，g为重力加速度，f为柯氏力参量。

在方程(1)中，水力挟带系数E_wx与E_wy是流动Richardson数(R_ix、R_iy)的函数，可用若干经验式表达^[5]。以x方向为例，它们是

　　Ashida et al.(1975): E_wx=0.0015/Ri_x；

　　Parker et al.(1987): E_wx=0.075/(1+718Ri_x^2.4)^2.5；

　　Fukushim et al.(1985):E_wx=0.0015/(0.0204+Ri_x)；

　　以及俞维升(1991)：E_wx=0.0034/Ri_x^2.9等。

在方程(2)、(3)中，剪切摩擦力τ_x、τ_y是作用在异重流上下界面摩擦力的总和，它们由下式所示。

其中 c_d=c_f(1+α′),被称为总阻力系数，范家骅等(1980)的c_d值约0.003，Parker et al.(1987)的c_d值在0.001～0.058之间，俞维升(1991)的c_d值在0.003～0.004之间，可见c_d值有很大的不确定性。

　　在方程(4)中，泥沙饱和浓度E_s是描述泥沙异重流流动特征参数的函数，对黄河口可用下式表示

?　水科院(1997)公式^[2，3]

s≤15kg/m³,E_s=123(V²/gh)^0.36(V/ω)^-0.33(1-Δh/h)^0.2

s>15kg/m³,E_s=9.7(V²/gh)^0.01(V/ω)^0.16(1-Δh/h)^0.22

V为合流速,h为水深,ω为泥沙沉速，Δh为潮差。

　　张青玉经验式 E_s=9.83(V²/h)^0.23

定解条件：

　　边界条件：陆地边界：V_n＝0(n为岸线的法线方向), s′=0(′表示陆地)

　　水边界：h^*(t)=h(t),s^*=s(*表示水边界)；河口边界：u=1.5(m/s) s=50(kg/m³)

　　初始条件：u=0 v=0 s=0 h=0.5(m)

　　采用有限差分法中的ADI法求解方程组(1)、(2)、(3)、(4)。事实上，我们应当同时求解描述泥沙异重流流动的方程，因为它们是一组相互耦合的方程组。然而同时求解相当复杂，我们便采取一种叫冻结系数法^[4]的近似耦合法来求解。

　　计算黄河口泥沙异重流时，泥沙异重流起点在黄河入海处。数值计算区域如图1所示，x、y方向空间步长均取1/3千米，时间步长取24秒；计算区海底地形如图2所示；柯氏参数f＝2Ωsinφ，Ω取7.2722×10^-5s^-1，φ取36.7°；泥沙密度ρ_s取2.65g/cm³，水体密度ρw取1.015g/cm³，泥沙异重流密度ρ_f取1.040g/cm³。

2 平面二维泥沙异重流的参数试验

　　为了定量地评价基本控制参数对泥沙异重流流动的影响，本文设计了一系列数值试验(或称敏感性试验)来展示各种参数与泥沙异重流流动特征之间的定量关系。所选取的基本控制参数包括水下斜坡坡角θ_x、θ_y，总阻力系数c_d，水力挟带系数E_w，泥沙饱和浓度E_s和泥沙沉速ω，进行数值试验时它们各自的输入值见表1。

以试验13为标准进行对比，试验13结果如图3。

表1 数值试验中基本控制参数的输入值

Input values of different basic controlling parameters

2.1 水下斜坡坡角的影响

　　水下斜坡坡角越大，有效重力作用越强，泥沙异重流加速越快，可以运行更远的距离，而且可以一直保持较高的含沙量。通过数值试验13、1、2、3可揭示坡角大小对泥沙异重流流动产生的显著影响。这些数值试验的控制参数除坡角不同外，其余参数均相同,见表1。流速变化图4显示坡角越大,流动速度越大。在河口以东2km处(图4)，试验13的流速为0.80m/s，试验1、2、3的流速分别为1.41m/s, 1.82m/s和2.67m/s。数值试验13(如图3(a))，由于其输入真实海底地形,泥沙异重流从黄河口流出后，便向周围扩散，这是由于实际上口门处地形最高，无论向北向南，还是向东均有一定的坡度。流动厚度变化图5同样显示，坡角越大，其流动厚度越大。试验13，1，2，3的0.8m等厚度线所括范围占整个计算区域的百分比分别为25.0, 66.6, 75.0, 83.3。悬沙含量变化图6显示，坡角越大，愈易保持高的悬沙含量，在河口以东2km处，试验13，1，2，3的悬浮含量分别为18.5，22.3，31.6，37.9g/l。

2.2 总阻力系数c_d的影响

　　总阻力系数c_d的大小主要与边界条件有关，对于不同边界不同性质的泥沙异重流，其总阻力系数c_d的值很不相同。c_d值大小对泥沙异重流影响很大，c_d越大，泥沙异重流在运动过程中动量损失越大，越易发生沿程沉降而导致快速消亡。

　　本文设计了数值试验4、5、13、7来定量分析c_d对泥沙异重流流动的影响。c_d值分别是0.001、0.005、0.01和0.05，依次增大。图4显示随着c_d值的增大，流速逐渐减小，在河口以东2km处，试验4的流速为1.36m/s，试验5，13的流速分别为1.05m/s，0.80m/s；试验7的c_d值高达0.05，几乎是所量测到的c_d值中的最大值，在如此大的阻力下，流速快速衰减，在2km处流速已降至0.32m/s。流动厚度变化图5表明，流动厚度的衰减速度随c?d值的增大而迅速增大，以0.8m等深线为例，在试验4中，0.8m等深线所括范围很大，占据整个计算区的3/4，而在试验13中，它的范围已大大缩减，仅占整个计算区域的1/4左右，试验7的范围更小，占整个区域不到1/8。 流速和流动厚度的减小，必然导致悬沙含量随之减小，如图6所示，试验4、5、13、7在河口以东2km处的悬沙含量分别是 32.6，25.3，18.5和16.4kg/m³。

2.3 水力挟带系数E_w的影响

　　当泥沙异重流流动处于超临界流状态，水力挟带系数E_w大于零，这时有一定水量的卷入；当泥沙异重流处于亚临界流状态，水力挟带系数E_w等于零，甚至小于零出现负挟带。黄河口的泥沙异重流一般情况下，其流速小于1.5m/s,悬沙浓度小于60g/l,属低浓度泥沙异重流，通常处于亚临界状态，因此在本文的数值试验中，水力挟带系数E_w大部分取零值，只有试验6和试验8中的E_w分别取0.0015／Ri和0.00034／Ri^2.9。三个试验的流速和悬沙含量基本类似，说明有无水力挟带项以及水力挟带系数的大小对流速和悬沙含量的分布影响不大；然而三个试验的流动厚度大相径庭(图7)，试验13的流动厚度沿程减小，而试验6，试验8的流动厚度沿程增大，E_w值越大，卷入的水量便越多，试验8的E_w略大于试验6的E_w，因此，试验8流动厚度的增加略快于试验6。

2.4 泥沙饱和浓度E_s的影响

　　泥沙饱和浓度Es是影响泥沙异重流数值模拟的关键因素，E_s越大，泥沙异重流越易保持而不消亡。本文设计了试验9与试验13来定量分析Es不同对泥沙异重流所造成的影响。试验9的E_s用水科院关于黄河口输沙力的公式来表达，试验13的E_s用张青玉公式表达，除E_s不同外，两个试验的其余参数均相同。流速变化图4清楚地显示，当E_s值增大时，流速增大，口门以东2km处的流速从0.56m/s增到0.80m/s。在流动厚度变化图5中，试验9的0.8m等值线范围仅占整个区域的1/10，而试验13为1/4左右，即流动厚度增加了2.5倍左右；如悬沙含量图6所示，试验9的悬沙含量在河口以东不到2km处已降至15kg/m³，而试验13的悬沙含量在河口以东2km处仍为18.5kg/m³，说明泥沙饱和浓度越小，泥沙异重流越易发生沿程沉降而导致其快速消亡。

2.5 泥沙沉速的影响

　　本文进行了数值试验10、11、13和12来定量分析泥沙沉速的影响，四次试验的沉速取值分别是0.475cm/s、0.086cm/s、0.036cm/s和0.020 cm/s。泥沙颗粒沉速的大小对泥沙异重流流速分布几乎无影响，图4显示三个试验11，12，13的流速相差不多，只是试验10的流速与其它3个相比偏小，是由于试验10的沉速值是其它3个的10～20倍相差很大造成的。而流动厚度和悬沙浓度对沉速值大小却相当敏感，沉速越大，流动厚度越小，相应空间点的悬沙浓度亦越小。尤其对于试验10，由于沉速过大(只是一假设值)，悬沙含量在河口以东2km处已降至零，说明不到2km悬浮泥沙已全部沉降。

3 结论

　　通过数值试验及分析得出水下斜坡坡角对泥沙异重流影响最大，即泥沙异重流强度对斜坡坡角改变相当敏感，且呈非线性关系，类似于抛物线型。总阻力系数c_d值大小对泥沙异重流影响次之。泥沙颗粒沉速大婿对泥沙异重流悬沙浓度有较大影响，对泥沙异重流的流速几乎无影响。当泥沙异重流流动处于超临界流状态，水力挟带系数E_w大于零，这时由于有一定水量的卷入，所以流动厚度沿程增加；而E_w对流速及悬沙浓度几乎无影响。泥沙饱和浓度E_s是影响泥沙异重流数值模拟的关键因素，E_s越大，泥沙异重流越易保持而不消亡。

参 考 文 献

[1] 庞重光，杨作升。黄河口泥沙异重流的数值模拟，海洋与湖沼，1999，待刊。

[2] “八五”攻关85-926-02-03报告,黄河口演变及整治研究。水科院,1995，35-38.

[3] UNDP支持黄河三角洲可持续CPR/91/144项目专题报告之十一,黄河口泥沙运动数值模拟研究。中国水利水电研究院,1997，12-15.

[4] 杨国录。 河流数学模型。海洋出版社,1993，296-302.

[5] Zeng J，Lowe D R. 1997. Numerical simulation of turbidity current flow and sedimentation: Theory. Sedimentology，44(1): 67-84.

[6] Wright L D, Yang Z S, 1990. Sediment transport and deposition at Huanghe mouth. Continental shelf Research，10(2): 36-44.