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

1 黄河口泥沙异重流基本控制参数的数值试验摘要: 为了定量地评价泥沙异重流基本控制参数对泥沙异重流流动的影响, 本文设计了一系列数值试验来展示各种参数与泥沙异重流流动特征之间的定量关系。采用 ADI 法求解三维方程经垂向积分后的平面二维泥沙异重流方程。所选取的基本控制参数包括水下斜坡坡角, 总阻力系数, 水力挟带系数, 泥沙饱和浓度和泥沙沉速。通过数值试验得出水下斜坡坡角θ、泥沙饱和浓度 Es 越大, 总阻力系数 cd、泥沙沉速ω越小,相应地泥沙异重流流速、流动厚度以及悬沙浓度越大。其中斜坡坡角对泥沙异重流影响最大。当泥沙异重流流动处于超临界流状态,水力挟带系数 Ew 大于零,这时由于有一定水量的卷入,所以流动厚度沿程增加。关键词:黄河口泥沙异重流控制参数数值试验 ADI 法泥沙异重流及形成的地质体在沉积学和石油地质中有 2 重要意义。在黄河口,泥沙异重流是泥沙搬运的重要方式, 是河口区沉积动力的重要过程, 是认识河口三角洲发育演化的重要途径。为使研究深入一步, 我们进行了泥沙异重流的数值模拟。近年来, 国内对异重流数模研究较少, 仅见方春明等(1997) 关于泥沙异重流潜入时的立面二维数值模拟以及邱晨霞(1995) 对盐水异重流的二维两层数值计算。 1 平面二维泥沙异重流控制方程及数值模拟目前,包含三个基本控制方程( 流体质量、动量守恒方程和泥沙质量守恒方程) 的泥沙异重流数值模拟被广泛应用,本文把一维方程推广到平面二维的情形,成为四方程模型,如下所示[1] 3 (1) (2) (3) (4) 4 其中 u,v 分别为 x,y 方向速度,θx,θy 分别为 x,y 方向坡度, Ewx 与 Ewy 分别为 x、y 方向的水力挟带系数, τx、τy为x 方向与 y 方向的剪切摩擦力, kx、 ky 分别为 x 方向与 y 方向的悬沙扩散系数, t 为时间, h 为异重流厚度, s 为深度平均的悬沙含量,ρs,ρw,ρf 分别为泥沙颗粒、水体和泥沙异重流的密度, Es 为泥沙饱和浓度, ω为泥沙颗粒的沉降速度, α为泥沙运动饱和系数, g 为重力加速度, f 为柯氏力参量。在方程(1) 中, 水力挟带系数 Ewx 与 Ewy 是流动 Richardso n 数(Rix 、 Riy) 的函数,可用若干经验式表达[5] 。以 x 方向为例,它们是 Ashida et al.(1975): Ewx= ; Parker et al.(1987): Ewx=/(1+) ; Fukushim et al.(1985):Ewx=/(+Rix) ; 以及俞维升(1991) : Ewx= 等。 5 在方程(2) 、(3) 中,剪切摩擦力τx、τy 是作用在异重流上下界面摩擦力的总和,它们由下式所示。其中 cd=cf(1+ α′), 被称为总阻力系数,范家骅等(1980) 的cd 值约 , Parker et al.(1987) 的 cd 值在 ~ 之间, 俞维升(1991) 的 cd 值在 ~ 之间,可见 cd 值有很大的不确定性。在方程(4) 中, 泥沙饱和浓度 Es 是描述泥沙异重流流动特征参数的函数,对黄河口可用下式表示蕌中国水科院(1997) 公式[2, 3] s≤ 15kg/m3,Es=123(V2/gh)(V/ ω)-(1- Δ h/h) 6 s15kg/m3,Es=(V2/gh)(V/ ω)(1- Δ h/h) V 为合流速,h 为水深,ω为泥沙沉速, Δh 为潮差。张青玉经验式 Es=(V2/h) 定解条件: 边界条件:陆地边界: Vn = 0(n 为岸线的法线方向), s′=0( ′表示陆地) 水边界: h*(t)=h(t),s*=s(* 表示水边界) ;河口边界: u=(m/s) s=50(kg/m3) 初始条件: u=0 v=0 s=0 h=(m) 采用有限差分法中的 ADI 法求解方程组(1) 、(2) 、(3) 、(4) 。事实上, 我们应当同时求解描述泥沙异重流流动的方程, 因为它们是一组相互耦合的方程组。然而同时求解相当复杂,我们便采取一种叫冻结系数法[4] 的近似耦合法来求解。计算黄河口泥沙异重流时,泥沙异重流起点在黄河入 7 海处。数值计算区域如图 1 所示, x、y 方向空间步长均取 1/3 千米,时间步长取 24 秒;计算区海底地形如图 2 所示; 柯氏参数 f=2Ω sin φ,Ω取 × 10-5s-1 ,φ取 ° ;泥沙密度ρs取 /