黄河河口海岸二维非恒定水流泥沙数学模型

黄河河口海岸二维非恒定水流泥沙数学模型

曹文洪，何少苓，方春明

(中国水利水电科学研究院泥沙研究所)

摘要：针对黄河河口海岸岸线变化剧烈和含沙量变幅大的特点，开发和建立了适合黄河河口海岸应用的平面二维动边界非恒定水流泥沙数学模型。验证表明，本模型可以较好地模拟黄河河口海岸泥沙输移和冲淤变化，为研究和解决多沙河口海岸的泥沙问题提供技术手段。

关键词：黄河口；挟沙能力；窄缝法；非恒定流；数学模型

收稿日期：2000-01-06

基金项目：国家重点基础研究发展规划项目(G1*******)资助

作者简介：曹文洪(1963-)，男，(满族)，黑龙江省人，中国水利水电科学研究院教授级高工，博士。

自本世纪七十年代以来，由于计算机技术的迅猛发展，国内外相继出现了众多的河口海岸泥沙数学模型[1-7]，有力地促进了河口海岸的泥沙研究的发展。然而，已有的河口海岸数学模型大多是模拟含沙量较低的河口海岸的泥沙运动，而能够模拟多沙和岸线延伸剧烈的河口海岸泥沙数学模型还极为少见。近年来，已有个别学者尝试用泥沙数学模型模拟黄河河口海岸的泥沙运动，如张世奇开发了一套黄河口平面二维泥沙冲淤数学模型，得到了较好的效果[18～20]。为了全面系统地反映黄河三角洲海陆动态交互影响机理和泥沙运动与湿地演替关系，本文开发和建立了径流、潮流和波浪作用下的黄河河口海岸平面二维动边界非恒定流非均匀沙不平衡输沙数学模型。

1 模型结构

1.1 水流运动基本方程

(1)

(2)

(3)

为谢才式中：U、V分别为潮流速在x及y方向的垂线平均值分量；Z为潮位；C

f

系数；F为柯氏系数，F=2ωsinφ，式中ω为自转角速度，φ为地理纬度；h为水深；Z

为海底起始高程。

b

1.2 潮流和波浪共同作用下泥沙运动方程在河口海岸地区，潮流和波浪是泥沙运动的最主要动力。“波浪掀沙和潮流输沙”使河口海岸地区的泥沙运动极为活跃，但也更为复杂。因此，在计算和预报河口海岸地区的泥沙运动和冲淤变化时，仅仅考虑潮流的作用是不全面的，还应考虑波浪的作用。

二维非恒定流不平衡情况下悬移质含沙量随空间和时间的变化规律的偏微分方程式(考虑泥沙粒径不均匀性时，以下各方程式中含沙量和沉速是相应于某一粒径组的)

(4)

式中：l表示非均匀沙中的第l组泥沙，h为水深；U、V分别x、y方向垂线平

均流速；S

l 为第l组泥沙垂线平均含沙量；S

*l

为第l组泥沙水流挟沙力；ω

l

为

第l组泥沙沉速；ε

s

为泥沙扩散系数；α为恢复饱和系数。

为方便起见，在以下分析中省略下角标l.设U

c 、U

w

和V

c

、V

w

分别表示潮流速

和波浪质点速度在x和y方向上的分量，则在流场任意点x和y方向上瞬时流速为：

U=U

c +U

w

(5)

V=V

c +V

w

(6)

将上述二式代入式(4)即得到潮流和波浪共同作用下的悬移质泥沙不平衡输沙方程式如下：

(7)

波浪对泥沙的作用在一个周期内为往返输移，其单向输沙作用一般只表现在半个周期内。由于波浪的周期比较小，一般只有几秒至十几秒；而潮流的周期是半天或一天。因此，当观察的时间或计算的步长大于波浪周期时，就可以对上式进行波浪周期的平均。由于潮流速与波浪质点速度基本无关，而波浪质点速度在一个波周期内的平均值基本为零，对式(5)在波浪周期内取时均后简化为：

(8)

式(8)在形式上与式(4)相同，但差别在于式(8)中的挟沙能力S

*

为潮流和波浪共同作用下的挟沙能力。

由于悬移质的不平衡输沙运动引起海底地形的冲淤变化，其河床变形方程为：

(9)

式中：γ′为海底淤积泥沙的干容重，ε为海底冲淤厚度。

通过上述分析可见，在计算河口海岸大范围泥沙冲淤变化时，只需计算潮流流速场，无需计算波浪流速场，只是在计算挟沙能力时考虑波高因子的影响，这与窦国仁[8]的分析结果是一致的。

1.3 动边界处理[21～22] 在河口海岸的计算区域内，由于潮流的涨落，实际水域边界是变动的。如何模拟变化的水域一直是比较困难的问题，本文第二作者采用窄缝法解决了这一问题，这里不详述。窄缝法是设想在岸滩上各空间步长内存在一条很窄的缝隙，缝内的水和岸滩前的水相连。这相当于把岸滩前的水域延伸到岸滩内。这样就可以把计算边界设在岸滩的窄缝内，成为具有一定水深的固定边界。 1.4 波高计算 波浪衰减过程以破碎点为界分为两部分，破碎点以外以底摩擦为主，破碎点以内主要是波浪破碎。波浪在破碎前波高可以表示为：

H/H 0=K s K r K f

(10)

式中：H 为浅水波高；H 0为深水波高；K s 、K r 、K f 分别为浅水因子、折射因子和底摩擦因子。

浅水波传播速度C 和波长L 服从如下规律：

L=gT 2/2πtanh(2πh/L) (11) C=gT/2πtanh(2πh/L)

(12)

根据微小振幅波理论：

K s =[tanh2πh/L(1+4πh/L/sinh4πh/L)-1/2 (13) K 2r =cos α0/(1-sin 2α0tanh 22πh/L)1/2 (14) 式中：α0为波峰与海岸线的夹角。

根据1984～1994年黄河口的实测资料得到Kf 计算式如下

K f =exp[-0.022H 1/10Δ

/

]

(15) 式中：Δx 为依据波所在位置到计算点的距离；为该距离内的平均水深；Ｈ1/10、L 0分别为依据波的波高和波长。

破波时破波波高Ｈb 与破波水深h b 的关系如下：