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

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

曹文洪，何少苓，方春明

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

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

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

收稿日期：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 的关系如下：

h b =1.28H

b

(16)

破波带内波高衰减过程可由下式计算[23]：

H5/2=0.455h

b 5/2-1.138/β+1h

b

1.5-β(h

b

β+1-hβ+1) (17)

式中：β为表征破波能耗散特征的经验常数，对淤泥质海滩β≈1.6.

1.5 絮凝在河口地区，细颗粒泥沙在盐水中易于发生絮凝，因此，必须考虑絮凝的作用。如以ω

s

和ω分别代表絮凝团粒及单颗泥沙在水中的平均沉速，两者之比值F称为絮凝因子，采用黄建维根据我国河口淤泥沉降试验建立的关系式[24]：

F=ω

s /ω=7×10-4d

50

-1.9(18)

式中：d

50为泥沙中值粒径，以mm计。该式的适用范围为d

50

<0.02mm，大于0.02mm

的泥沙基本不发生絮凝现象。

对于细颗粒泥沙在清水中的沉速一般采用Stockes公式计算：

ω=1/18γs-γ/γgd2/ν(19)

1.6 挟沙能力黄河口地区的泥沙主要来源于径流，并在潮流和波浪的共同作用下向深海输移，因此，单向流(包括径流和潮流)与波浪共同作用下的挟沙能力公式是否符合实际，是泥沙数学模型能否准确模拟河口海岸泥沙运动规律和冲淤变化的关键所在。本数学模型采用笔者[25]建立的挟沙能力公式：

S *=ρ

s

δ/νT+

B

dS

Vm

u4

*

/ω

s

u2

*

k[(1+/C

f

k)J

1

-

/C

f

kJ

2

]/(1-δ

b

/δ

b

)Z

(20)

式中：δ为单位床面面积内的平均猝发面积；T+

B

为近壁区低流速带无因次时间

间隔，大量的水槽试验表明：T+

B ≈100，δ≈0.016；d为泥沙粒径；S

vm

为单位床

面层的极限体积含沙量；v为水的运动粘滞系数；K为卡尔曼常数；b为近底某处距床面距离，一般取δ

b

=b/h=0.01～0.05；

;Z=ω

s /kω

s

. u

*

为波浪和单向流共同

作用下的摩阻流速，采用Bijker[26]的研究成果:

u *=u

*c

[1+1/2(δu

m

/U)2]1/2 (21)

其中，u

*c

为单向流的摩阻流速；U为垂线平均流速。

δ=akC f/(22) Bijker给出a=0.45，本文根据1984～1994年黄河口实测资料得到a=0.18.

u

m

为床面波浪质点运动最大水平分速：

u

m

=πH/Tsinh2πh/L (23)

u

*k

为临界摩阻流速，根据谢尔兹(Shields)曲线确定。

目前，人们对水流分组挟沙能力或者挟沙力级配的机理尚不十分明了，各家采用的方法也各不相同。本次计算采用韩其为的处理方法[27]，即:

S *l =p

l

S

*

(24)

(25)

式中:p

l

为悬沙级配。

2 验证计算

如前所述，本数学模型中计算公式中的

待定参数是采用黄河口1984～1994年实测

资料和水槽试验资料率定的。下面采用黄河

口1996年实测资料进行验证计算。本次计

算域为汊河流路口门附近海区南北长

20km(北京坐标4180.0～4200.0)，东西长

22.5km(北京坐标20692.5～20715.0).计算

网格的空间步长为Δx=250m，Δy=500m，

时间步长为Δt=60s.潮流、潮位和含沙量

验证计算时间为1996年9月5～6日实测海

流过程，与此相应河流入海的流量为1560～

1570m3/s，含沙量为21.7～25.0kg/m3；冲

淤验证计算时间为1996年6月至11月；开

边界含沙量由式(20)计算求得。计算域内各

测点和验证断面的位置如图1所示。计算潮

位和潮流速与实测值的比较如图2和图3所

示，表明计算与实测符合良好。含沙量计算

结果与实测资料的比较如图4所示，表明计

算与实测基本符合。根据1996年汛前6月

和汛后10月的施测的两次地形资料，对数

学模型进行冲淤验证计算。从图5可见，计

算的河口淤积形态、淤积部位和淤积量与实

测值基本吻合。

图1 计算网格划分及测点位置

图2 潮位验证

图3 潮流验证

图4 含沙量验证

图5 冲淤验证

3 河口泥沙输移特征及新生湿地演变预测计算

3.1 计算条件计算的初始地形为1996年10月实测的汊河流路河口及附近海域的地形，选择利津站1990～1994年实测水文系列作为黄河入海的水沙条件，将悬移质泥沙分五组，即<0.007mm，0.007～0.01mm，0.01～0.025mm，0.025～0.05mm，>0.05mm.在此基础上采用泥沙数学模型对未来5年湿地的演进进行了预报计算。

3.2 入海泥沙输移特性进入黄河口的泥沙既有小于0.01mm的细沙又有大于0.025mm的粗沙，究竟多大粒径的泥沙被输入到深海，多大粒径的泥沙淤积在三角洲地区，多大粒径的泥沙既有一部分输入深海又有一部分沉积下来，其各自的数量是多少，这些问题都是黄河三角洲重大和关键的泥沙问题，同时又是迄今为

止没有得到很好解决的问题。此次采用泥沙数学模型进行了计算，结果如表1所示。由表1可见，小于0.007mm粒径组的泥沙有86.6%输入到深海，0.007～0.01mm粒径组的泥沙有77.0%输入到深海，0.01～0.025mm粒径组的泥沙有10.2%输入到深海，0.025～0.05mm粒径组的泥沙有4.1%输入到深海，大于0.05mm粒径组的泥沙全部都沉积在三角洲地区。也就是说小于0.01的细颗粒泥沙大部分(不是全部)输入到深海，少部分沉积在三角洲地区；大于0.025mm的粗沙绝大部分落淤。

表1 不同粒径组泥沙输入深海的比例的计算结果

粒径组/mm 来沙量/(108t) 淤积量/(108t) 排入深海沙量

/(108t)

排沙比(%)

<0.007 5.68 0.76 4.92 86.6

0.007～0.01 1.74 0.40 1.34 77.0

0.01～0.025 6.08 5.46 0.62 10.2

0.025～0.05 5.34 5.12 0.22 4.1

>0.05 3.51 3.51 0.00 0.0

全沙22.35 15.25 7.10 31.8

3.3 新生湿地演变预测黄河三角洲是世界上新生湿地增长最快的地区，湿地演替与泥沙运动密切相关，湿地内动植物分布及生态演替过程与湿地形成的大小、高程、部位的关系亦息息相关。因此，对新生湿地演变的预测是十分重要的。

表2 口门外5年累积淤积情况

淤积厚度/m >10 5～10 1～5

淤积面积/km210.5 60.4 167.8

本次计算得出了未来5年汊河流路河口海岸湿地大小、高程和部位。图6

给了5年后河口岸线的变化，5年共造陆(即陆地湿地)76.0km2，在原有黄河三角洲的基础上形成一个新的小三角洲(沙嘴).图7和表2所示则给出了新生湿地部位、高程和淤积面积，表明大于10m淤积厚度的面积为10.5km2，5～10m淤积厚度的面积为60.4km2，1～5m淤积厚度的面积为167.8km2.

图6 黄河口岸线演变图7 口门外5年累积淤积厚度分布

4 结语

针对黄河河口海岸岸线变化剧烈和含沙量变幅大的特点，开发和建立了适合黄河河口海岸应用的平面二维动边界非恒定水流泥沙数学模型。模型中采用笔者建立的潮流与波浪共同作用下的悬移质挟沙能力公式，适应黄河河口海岸含沙量变幅大的特点，反映径流、潮流和波浪等动力因子对泥沙输移的影响；采用窄缝法解决黄河河口海岸岸线变化剧烈的问题；通过分粒径组计算以反映海流对不同粒径组泥沙的输移特性；河口地区泥沙颗粒较细，模型考虑了絮凝影响。经潮位、潮流、含沙量和泥沙冲淤的验证，计算与实测符合良好。计算结果表明，小于0.007mm粒径组的泥沙有86.6%输入到深海，0.007～0.01mm粒径组的泥沙有77.0%输入到深海，0.01～0.025mm粒径组的泥沙有10.2%输入到深海，0.025～0.05mm粒径组的泥沙有4.1%输入到深海，大于0.05mm粒径组的泥沙全部都沉积在三角洲地区。说明海流对不同粒径组泥沙的输移是极为敏感的，细颗粒泥沙大部分输送到深海，粗颗粒泥沙大部分基本都沉积在河口三角洲地区。对汊河流路未来5年新生湿地形成的大小、高程和部位进行了预报，具有实用价值。

参考文献：

[1] Cole P, Miles G V. Two Dimensional of Mud Transport [J]。ASCE, 1983，109(1).

[2] Holly F M. Dispersion Simulation in 2-D Tidal Flow[J]。ASCE,1984，110(7).

[3] Teision Ch, Latteux B. A Depth-Integrated Bidimensional Model of Suspended Sediment Transport[C]。Proc. of the 3rd Int. Symp. on River sedimentation, The Univ. of Mississippi,1986.

[4] Christopher G. Koutitas. Mathematical Models in Coastal Engineering[M]。Pentech Press,1988.

[5] Davies A M, Jones J E, and Xing J. Review of Recent Developments in Tidal Hydrodynamic Modeling Ⅰ:Spectral Models[J]。ASCE,1997，123(4).[6] Davies A M,Jones J E, and Xing

J. Review of Reccent Developments in Tidal Hydrodynamic Modeling Ⅱ:Turbulence Energy Models[J]。ASCE,1997，123(4).

[7] Ebrahim Jabbari，Jean Berlamont. Finite Element Modeling of Cohesive Sediment Transport in the Western Scheldt Estuary[C]。Advances in Hydro-Science and-Engineering, Volume Ⅱ,Beijing, China,1995.

[8] 窦国仁，董凤舞，窦希萍。河口海岸泥沙数学模型研究[J]。中国科学(A辑)，1995，25(9).

[9] 匡翠萍。长江口拦门沙冲淤及悬沙沉降规律研究和水流盐度泥沙数学模型[D]。南京水利科学研究院博士学位论文，1993.

[10] 曹祖德，王运洪。水动力泥沙数值模拟[M]。天津：天津大学出版社，1994.

[11] 黄胜，卢启苗。河口动力学[M]。北京：水利电力出版社，1995.

[12] 王尚毅，等。河口工程泥沙数学模型[M]。北京：海洋出版社，1990.

[13] 林秉南，等。潮汐水流泥沙输移与河床变形的二维数学模型[J]。泥沙研究，1988，(2).

[14] 窦希萍，罗肇森。波浪潮流共同作用下的二维泥沙数学模型[J]。水利水运科学研究，1992，(4).

[15] 李蓓，张征。珠海电厂煤港港区海域二维波浪、潮汐泥沙淤积数学模型[C]。第七届全国海岸工程学术讨论会论文集，北京：海洋出版社，1993.

[16] 徐建益，袁建中。长江河口局部区域内泥沙冲淤计算方法及应用[J]。水利学报，1992，(12).

[17] 史英标。潮汐河口平面二维悬沙输移及河床变形的数值模拟[J]。河口与海岸工程，1997，(1).

[18] 张世奇。二维动边界潮汐输沙及河床变形计算[J]。泥沙研究，1988，(4).

[19] 张世奇。黄河口输沙及冲淤变形研究[J]。水利学报，1990，(1).

[20] 张世奇。黄河口及三角洲冲淤演变计算原理及方法[J]。泥沙研究，1997，(2).

[21] 何少苓，等。隐式破开算子法在二维潮流计算中的应用[J]。海洋学报，1985，7(2).

[22] 何少苓，王连祥。窄缝法在二维边界变动水域计算中的应用[J]。水利学报，1986，(12).

[23] 虞志英，等。连云港吹泥区岸滩自然冲淤及吹泥条件下海滩演变的观测分析[J]。海洋与湖沼，1986，17(4).

[24] 黄建维。粘性泥沙在静水中沉降特性的试验研究[J]。泥沙研究，1981，(2).

[25] 曹文洪。黄河河口海岸泥沙输移规律和演变机理及湿地变迁研究[D]。中国水利水电科学研究院博士学位论文，1999.

[26] Bijker E. W. Littoral Drift as Function of Waves and Current[C]。Proc. 11th Conf. Coastal Engineering,1968.

[27] 韩其为。悬移质不平衡输沙的初步研究[C]。第一届河流泥沙国际学术讨论会论文集，1980.