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

黄河河口海岸二维非恒定水流泥沙数学模型曹文洪，何少苓，方春明(中国水利水电科学研究院泥沙研究所)摘要：针对黄河河口海岸岸线变化剧烈和含沙量变幅大的特点，开发和建立了适合黄河河口海岸应用的平面二维动边界非恒定水流泥沙数学模型。

验证表明，本模型可以较好地模拟黄河河口海岸泥沙输移和冲淤变化，为研究和解决多沙河口海岸的泥沙问题提供技术手段。

关键词：黄河口；挟沙能力；窄缝法；非恒定流；数学模型收稿日期：2000-01-06基金项目：国家重点基础研究发展规划项目(G1*******)资助作者简介：曹文洪(1963-)，男，(满族)，黑龙江省人，中国水利水电科学研究院教授级高工，博士。

自本世纪七十年代以来，由于计算机技术的迅猛发展，国内外相继出现了众多的河口海岸泥沙数学模型[1-7]，有力地促进了河口海岸的泥沙研究的发展。

然而，已有的河口海岸数学模型大多是模拟含沙量较低的河口海岸的泥沙运动，而能够模拟多沙和岸线延伸剧烈的河口海岸泥沙数学模型还极为少见。

近年来，已有个别学者尝试用泥沙数学模型模拟黄河河口海岸的泥沙运动，如张世奇开发了一套黄河口平面二维泥沙冲淤数学模型，得到了较好的效果[18～20]。

为了全面系统地反映黄河三角洲海陆动态交互影响机理和泥沙运动与湿地演替关系，本文开发和建立了径流、潮流和波浪作用下的黄河河口海岸平面二维动边界非恒定流非均匀沙不平衡输沙数学模型。

1 模型结构1.1 水流运动基本方程(1)(2)(3)为谢才式中：U、V分别为潮流速在x及y方向的垂线平均值分量；Z为潮位；Cf系数；F为柯氏系数，F=2ωsinφ，式中ω为自转角速度，φ为地理纬度；h为水深；Z为海底起始高程。

b1.2 潮流和波浪共同作用下泥沙运动方程在河口海岸地区，潮流和波浪是泥沙运动的最主要动力。

“波浪掀沙和潮流输沙”使河口海岸地区的泥沙运动极为活跃，但也更为复杂。

因此，在计算和预报河口海岸地区的泥沙运动和冲淤变化时，仅仅考虑潮流的作用是不全面的，还应考虑波浪的作用。

二维非恒定流不平衡情况下悬移质含沙量随空间和时间的变化规律的偏微分方程式(考虑泥沙粒径不均匀性时，以下各方程式中含沙量和沉速是相应于某一粒径组的)(4)式中：l表示非均匀沙中的第l组泥沙，h为水深；U、V分别x、y方向垂线平均流速；Sl 为第l组泥沙垂线平均含沙量；S*l为第l组泥沙水流挟沙力；ωl为第l组泥沙沉速；εs为泥沙扩散系数；α为恢复饱和系数。

为方便起见，在以下分析中省略下角标l.设Uc 、Uw和Vc、Vw分别表示潮流速和波浪质点速度在x和y方向上的分量，则在流场任意点x和y方向上瞬时流速为：U=Uc +Uw(5)V=Vc +Vw(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.28Hb(16)破波带内波高衰减过程可由下式计算[23]：H5/2=0.455hb 5/2-1.138/β+1hb1.5-β(hbβ+1-hβ+1) (17)式中：β为表征破波能耗散特征的经验常数，对淤泥质海滩β≈1.6.1.5 絮凝在河口地区，细颗粒泥沙在盐水中易于发生絮凝，因此，必须考虑絮凝的作用。

如以ωs和ω分别代表絮凝团粒及单颗泥沙在水中的平均沉速，两者之比值F称为絮凝因子，采用黄建维根据我国河口淤泥沉降试验建立的关系式[24]：F=ωs /ω=7×10-4d50-1.9(18)式中：d50为泥沙中值粒径，以mm计。

该式的适用范围为d50<0.02mm，大于0.02mm的泥沙基本不发生絮凝现象。

对于细颗粒泥沙在清水中的沉速一般采用Stockes公式计算：ω=1/18γs-γ/γgd2/ν(19)1.6 挟沙能力黄河口地区的泥沙主要来源于径流，并在潮流和波浪的共同作用下向深海输移，因此，单向流(包括径流和潮流)与波浪共同作用下的挟沙能力公式是否符合实际，是泥沙数学模型能否准确模拟河口海岸泥沙运动规律和冲淤变化的关键所在。

本数学模型采用笔者[25]建立的挟沙能力公式：S *=ρsδ/νT+BdSVmu4*/ωsu2*k[(1+/Cfk)J1-/CfkJ2]/(1-δb/δb)Z(20)式中：δ为单位床面面积内的平均猝发面积；T+B为近壁区低流速带无因次时间间隔，大量的水槽试验表明：T+B ≈100，δ≈0.016；d为泥沙粒径；Svm为单位床面层的极限体积含沙量；v为水的运动粘滞系数；K为卡尔曼常数；b为近底某处距床面距离，一般取δb=b/h=0.01～0.05；;Z=ωs /kωs. u*为波浪和单向流共同作用下的摩阻流速，采用Bijker[26]的研究成果:u *=u*c[1+1/2(δum/U)2]1/2 (21)其中，u*c为单向流的摩阻流速；U为垂线平均流速。

δ=akC f/(22) Bijker给出a=0.45，本文根据1984～1994年黄河口实测资料得到a=0.18.um为床面波浪质点运动最大水平分速：um=πH/Tsinh2πh/L (23)u*k为临界摩阻流速，根据谢尔兹(Shields)曲线确定。

目前，人们对水流分组挟沙能力或者挟沙力级配的机理尚不十分明了，各家采用的方法也各不相同。

本次计算采用韩其为的处理方法[27]，即:S *l =plS*(24)(25)式中:pl为悬沙级配。

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.60.007～0.01 1.74 0.40 1.34 77.00.01～0.025 6.08 5.46 0.62 10.20.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.83.3 新生湿地演变预测黄河三角洲是世界上新生湿地增长最快的地区，湿地演替与泥沙运动密切相关，湿地内动植物分布及生态演替过程与湿地形成的大小、高程、部位的关系亦息息相关。