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

文章编号:1005-9865(2006)01-0139-16海岸河口泥沙数学模型研究进展李孟国(天津水运工程科学研究所,天津　300456)摘　要:对海岸河口泥沙运动的基本方程、数值方法、边界条件、参数选取等进行了归纳总结和评述,以期对本学科的发展起到一定的引导和促进作用。

关键词:海岸河口;泥沙;数学模型;挟沙力;悬移质(悬沙);推移质(底沙)中图分类号:TV148 文献标识码:AA review on mathematical models of sediment in coastal and estuarine watersLI Meng -guo(Tianjin Research Institute of Water Transport Engineering ,Tianjin 300456,China )A bstract :This paper makes s ummaries and com ments on various mathematical models (basic equations ,numerical methods ,boundary condi -tions and parameter determination ,etc .)of sediment movement in coastal and estuarine waters ,in the hope of guiding and promoting the stud y of the subject .Key words :coast and estuary ;sediment ;mathematical model ;sediment -carrying capacity ;suspended load ;bed load收稿日期:2004-10-25基金项目:中国博士后科学基金资助项目(2003034327)作者简介:李孟国(1964-),男,天津武清人,研究员,博士,主要从事海岸河口水动力泥沙数值模拟研究。

黄河小浪底调水调沙问题数学建模
黄河小浪底调水调沙问题是指通过调整黄河水流的水量和输沙量来解决黄河小浪底河道的淤积问题。

数学建模可以帮助我们分析和预测黄河小浪底的水流和沙传输规律，从而提出合理的调水和调沙方案。

以下是数学建模中可能涉及的一些步骤和方法：
1. 数据收集和处理：收集黄河小浪底相关的水文数据、地质资料和历史数据，对数据进行整理和处理，建立合适的数据模型。

2. 建立水流模型：通过流体力学理论和水流实验数据，建立黄河小浪底水流的数学模型，包括水流速度、水动力和水力调控方面的参数。

3. 建立沙传输模型：根据黄河小浪底河道的地质特征和沙传输规律，建立沙传输的数学模型，包括输沙通道的沙动力和沙质输运规律方面的参数。

4. 模型验证和参数拟合：利用已有的观测数据和实验数据验证建立的水流和沙传输模型，并通过参数拟合来优化模型的准确性和适用性。

5. 模拟预测和优化调控：利用建立的数学模型，进行水流和沙传输的模拟预测，通过调整输水和输沙量来优化黄河小浪底的调水和调沙方案，以达到降低淤积和维护航道的目的。

数学建模可以辅助相关专业的研究人员和决策者做出科学的决策，使调水和调沙方案更加合理和有效，减少淤积和保护黄河流域的生态环境。

河口海岸泥沙数学模型研究河口海岸是地球上一种独特而重要的地理环境，具有复杂的动力和物质输运过程。

其中，泥沙输运是河口海岸过程的重要部分，它影响着河口海岸的形态、地貌和生态系统的功能。

为了更好地理解和预测河口海岸的行为，我们构建并研究了一个新型的泥沙数学模型。

我们的模型基于以下假设：河口海岸的泥沙输运主要受到水文条件、地形和海洋环境的影响。

我们用一系列偏微分方程来表达这个系统，包括水流速度、泥沙浓度、地形变化等。

我们还考虑了泥沙的沉积和侵蚀，以及与周围环境的相互作用。

我们选取了一个具体的河口海岸作为案例，将我们的模型应用于此，以检验其有效性和准确性。

通过与实地观测数据进行比较，我们的模型在预测泥沙输运、沉积和侵蚀方面表现出良好的性能。

这表明我们的模型可以有效地应用于实际问题的解决。

我们的模型具有几个主要的优点。

它考虑了多种影响因素，如水流、泥沙浓度、地形等。

我们的模型具有良好的灵活性，可以适用于不同的河口海岸环境。

然而，我们的模型还有一些局限性，例如在处理一些极端环境条件时，可能需要更复杂的物理机制和更精确的参数设定。

我们的河口海岸泥沙数学模型提供了一种有效的工具，可以帮助我们理解和预测河口海岸的行为。

尽管还有改进的空间，但这个模型已经展示出其在研究和应用中的重要价值。

希望我们的工作能为未来河口海岸研究提供有价值的参考和启示。

我们将继续研究和改进我们的数学模型，以更好地理解和预测河口海岸的行为。

我们将以下几个方面：一是提高模型的精度和适应性，以应对更复杂的环境条件和需求；二是将模型与其他相关模型进行集成，形成更完整的河口海岸系统模型；三是加强模型的验证和测试，以确保其准确性和可靠性。

我们也将利用先进的计算技术和算法，提高模型的计算效率和性能。

这将使我们能够更有效地解决实际问题，并为河口海岸的研究和管理提供更强大的支持。

河口海岸泥沙数学模型研究是一项富有挑战性和实用性的工作。

通过建立和应用数学模型，我们可以更好地理解和预测河口海岸的行为，为相关研究和应用提供有力的支持。

黄河三角洲海岸冲淤及泥沙输运模型数值模拟分析王奎峰;张太平;宋新强;许国辉;尚桂勇;郑建国【摘要】黄河三角洲是陆海交界地带,陆地海洋相互作用显著,泥沙所引起的岸滩演变、港口航道的淤积、水动力环境的改变等问题比较复杂.该文采用三维HEM-3D 数值模型对黄河三角洲海域流场变化、盐度、悬浮泥沙浓度及海岸冲淤分布进行了潮流周期内的数值模拟分析.结果表明,黄河三角洲海域的悬浮泥沙浓度分布与潮流场变化和河口泥沙输入有密切的关系,在三角洲北部受五号桩外强潮流区的影响,近岸海底的泥沙发生明显的再悬浮,并在涨潮流向南输送,含沙量达1.5 g/L左右.在三角洲南部(现行河口区域和莱州湾区域),受现行河口入海泥沙扩散的影响显著.海域年冲淤分布,在北部废弃三角洲区域,由于海洋动力作用强烈,浅水冲刷,海底侵蚀显著,形成了明显呈沿岸展布的侵蚀中心,在侵蚀中心以外,侵蚀快速减弱.在现行河口区域,以淤积为主,在羽状流扩散的控制下,泥沙入海后向南输运,至莱州湾区域逐渐减弱.数值模拟的结果与卫星遥感解译的岸线变化基本一致,效果良好.【期刊名称】《山东国土资源》【年(卷),期】2018(034)011【总页数】10页(P22-31)【关键词】三维数值模型;海域流场;悬浮泥沙;海岸冲淤;黄河三角洲【作者】王奎峰;张太平;宋新强;许国辉;尚桂勇;郑建国【作者单位】山东省地质科学研究院,国土资源部金矿成矿过程与资源利用重点实验室,山东省金属矿产成矿地质过程与资源利用重点实验室,山东济南 250013;山东大学土建与水利学院,山东济南 250061;山东省沉积成矿作用与沉积矿产重点实验室,山东青岛 266590;山东省地质科学研究院,国土资源部金矿成矿过程与资源利用重点实验室,山东省金属矿产成矿地质过程与资源利用重点实验室,山东济南250013;东营市国土资源局河口分局,山东东营257100;中国海洋大学环境科学与工程学院,山东青岛 266100;东营市公路管理局河口分局,山东东营 257200;中国海洋大学环境科学与工程学院,山东青岛 266100【正文语种】中文【中图分类】P737河口是河流沉积物向海域传输的通道，河口泥沙传输过程相当复杂，河口泥沙输运过程对于研究河口陆海相互作用及近岸地质环境演化有重要的参考意义[1-2]。

波浪作用下我国河口海岸地区的一种输沙计算马晶河海大学海洋学院，南京（210098）Email: sxyzmj@摘要：本文讨论了我国河口海岸泥沙的一般特点，对现有波浪作用下的非粘性沙输移的计算方法作了介绍，同时结合前人的结论和数学模型，提出了一个波浪底部边界层流速计算的模型，并在此基础上讨论了泥沙输移的计算。

关键词：粉沙；波浪边界层；紊动边界层；k-ε模型1.前言河口海岸泥沙问题是河口海岸动力研究过程中的一个很重要的组成部分。

波浪、潮流的作用是河口海岸地区泥沙运动的主要动力因素。

波流共同作用下的泥沙运动是海岸演变的直接原因[1]。

泥沙运动的结果会产生岸滩冲刷、淤积、岸线变形、港口航道淤积、防波堤等建筑物根底部淘刷等问题，而这些问题是海岸河口地区海岸工程、水利工程、港口航道工程等必须要考虑和解决的问题[2]。

对于河道内的单向水流作用下的泥沙运动研究已比较成熟，但是在河口海岸地区，由于受波浪，潮流等海洋动力因素的影响，泥沙运动变得更为复杂。

故搞清楚河口海岸泥沙的计算是很有价值的。

物理模型实验和数学模型研究是泥沙问题研究的两种重要手段。

2.泥沙情况2.1我国海岸河口地区的泥沙的一般特点我国海岸河口地区的泥沙以细沙、粉沙为主，具有如下一般特点[2]：泥沙运动具有非恒定性；泥沙运动具有非饱和性(不平衡输沙)，即水体含沙量不一定等于水体挟沙力；泥沙颗粒很细，悬移质(悬沙) 一般呈絮凝沉降；由于泥沙粒径较细，泥沙运动通常属于悬移质运动范畴，在潮汐河口及其附近海域，底床的冲淤变化主要由悬沙造成；海岸河口水体大多属于低含沙水体，悬移质的存在可视为不影响水流动力；泥沙粒径变化不大，可视为均匀沙。

因此在海岸河口地区使用均匀沙模型是合适的。

由以上几点可知在我国河口海岸地区主要关注的泥沙类型是细沙、粉沙，主要运动形式是悬移质运动，考虑泥沙沉降时要注意絮凝问题。

由于悬移质的存在可视为不影响水流动力，所以水流计算（流场）可以独立进行。

2000年9月水 利 学 报SHU IL I XU EBAO 第9期收稿日期:1999209227作者简介:江春波(1960-),男,吉林镇赉人,副教授,博士.文章编号:055929350(2000)0920020205河道立面二维非恒定水温及污染物分布预报模型江春波1,张庆海1,高忠信2(11清华大学水利系,北京　100084;21中国水利水电科学研究院机电所,北京　100038)摘　要:通过建立一种立面二维模型,考虑自由水面的变化,采用全显式的有限体积法模拟了河道流动的水温及悬浮污染物质分布.该模型适用于大尺度水域中较长期的流动、水温及悬浮物质迁移问题.通过优化取水口位置高程将水体中污染物质残存总量控制到最少,对水库实际流动、水温、悬浮物质分布的模拟结果与实际观测符合良好.关键词:温度分层流;立面二维模型;水污染控制中图分类号:X143 文献标识码:A由于人类活动及工农业生产,使河流两岸的植被遭到破坏,大量泥沙及其悬浮物质流入江河、湖泊等水域,造成河流、水库水质的严重污染.随着人们对环境保护的日益重视,已经开展了大量的对水污染治理的研究,并且提出了相应的工程措施.其中,建立一个预报河道水域中流动、水温和污染物质分布的数学模型,是水污染控制方法之一.模拟河道水流流动特性的数学模型有二维与三维[1,2]之分.与三维数学模型相比,二维模型输入数据少、计算效率高.对于流速、温度以及浓度等在水深方向上分布基本均匀的情况,用沿水深方向积分后的数学模型[3]就可以得到满意的结果.而对于河道比较窄,流速、水温或污染物在深度方向上分布具有明显差异的情况,如垂向有回流旋涡、密度分层等,用立面二维模型[4,5]可以充分反映流动特征,与完全三维模型相比还可以节省计算工作量.对于简单均匀流动(垂向流速近似为零),文献[4]在假定悬浮泥沙颗粒为单一粒径且不考虑颗粒本身沉降速度的情况下,采用剖面二维模型模拟了泥沙颗粒在河道中的迁移过程;对于考虑复杂流态的立面二维模型,文献[5]采用k 2ε双方程湍流模型计算了一个立面二维流动,实际算例为一个定常的悬浮泥沙分布问题.以上两种模型均未考虑自由水面的变化.本研究提出了一种预测河道型水库中流速、温度和悬浮污染物质分布的立面二维数学模型.模型的特点为:(1)数值格式全部采用显式,具有节省机器内存、计算速度快的优点;(2)可以模拟较长时期(如几个月,几年)的流态、水温和悬浮污染物质分布;(3)考虑河道宽度的变化以及汛期自由水面的变化,适合于非定常问题的模拟;(4)考虑悬浮物质粒径不均匀性及自身在水中的沉降速度,并通过取水口布置对水库中悬浮物质分布及残存总量进行控制.1　数学模型111　基本方程　从Navier 2Stokes 方程和对流扩散方程出发,将方程沿河宽方向积分,通过时均化处理并引入湍流模型,就得到本文将要求解的立面二维方程组.—02—55x (B u )+55z(B w )=0(1)55t (B u )+55x (B uu )+55z (B w u )=-B ρ05p 5x +55x B v x 5u 5x +55z B v z 5u 5z -τw x ρ0(2)55z(p )=-ρg 或p =-ρg (h +ζ)(3)55t (B T )+55x (B u T )+55z (B w T )=55x B k x 5T 5x +55z B k z 5T 5z +B H ρ0c p (4)55t (Bc )+55x (B uc )+55z [Bc (w -w s )]=55x B λx 5c 5x +55z B λz 5c 5z(5)式中:u ,w 为水平和铅直方向的流速分量;ζ为水面相对高度;w s 为悬浮物质的沉降速度;T 为水温;C 为悬浮物质的浓度;v x ,v z 为涡粘性系数分量;k x ,k z 为水温扩散系数;λx ,λz 为污染物质的湍动扩散系数;p 为压力;ρ为水的密度;ρ0为水的参考密度;g 为重力加速度;h 为水深;c p 为水的比热;τw x 为边壁阻力;H 为太阳辐射热源项.H =55z(1-β)<0exp [-ξ(ζ-z )](6)式中:<0为太阳辐射量;β为热辐射的反射系数;ξ为常系数.对于直径在10～100μm 的悬浮质,其沉降速度w s 可由如下的斯托克斯公式确定:w s =g 18v ρs ρ-1d 2(7)式中:ρs 为悬浮物质的密度;d为悬浮物质的直径;v 为水的分子粘性系数.水的密度ρ是水温T 和污染物浓度c 的函数,由如下公式确定:ρ=ρ(T )+c 1-ρ(T )ρs (8) 从式(1)～式(8)中可见,模型考虑了河道宽度B 的变化、自由水面与大气的热交换、悬浮物质的沉降速度.对于悬浮颗粒浓度的计算,模型根据粒径的变化将其分为不同的组次分别利用式(5)、式(8)进行计算,然后将得到的不同粒径的悬浮质浓度按百分比进行加权平均,就可得到悬浮颗粒的浓度.模型的着重点在于掌握悬浮污染物质分布与水流速分布的关系,而不考虑悬浮污染物的生物及其化学反应过程.图1　计算区域和边界条件112　边界条件　计算区域如图1所示.(1)上游入口条件:水平方向流速—根据上游来流流量Q in 确定(垂向流速为零);入口处水温T in —根据气温的变动情况而确定;悬浮质浓度C in —根据与流量Q in 的相关分析而确定.(2)取水口(出口)边界条件:取水口处的流速根据发电和下游用水量确定;温度和浓度的法向梯度5T/5n ,5c/5n 为零.(3)水面边界条件:在水面由于水温和气温不同,有热量交换产生,其值<L 可由下式计算<L =<e +<c +<ra(9)式中:<e ,<c ,<ra 分别表示在水面上由于蒸发、传导和辐射而引起的热量交换.另外,水面上还给定悬浮物质的沉降通量w s c =0及扩散通量5c/5n =0.(4)底部河床边界条件:给定无滑移流速条件;温度和浓度的扩散通量为零.由于悬浮物质具有沉降速度,沉到河床底部的悬浮质的数量为:{w su c u B u -w sd c d B d -(B u -B d )w s c}Δx i (10)式中:下脚标u ,d 分别表示控制体的上、下表面;w su ,w sd 分别为在控制体的上下表面的沉降速度;Δx i 为控制体在水平方向的尺寸.在式(10)中,第三项是沉降到河床底部的悬浮质的数量.2　数值格式图2　控制体和变量的位置示意采用控制体积方法对基本方程(1)～(6)进行空间离散,根据交错网格方法,变量u ,w 以及标量(c ,T )分别布置在不同的空间位置上.如图2所示,求解水平流速分量u 的控制体用虚线表示,求解垂直流速分量w 和标量c ,T 的控制体用实线表示,两种控制体在水平位置上相差半个网格.若第k 时间步的各变量数值已知,则第k +1时间步的浓度值c k +1(i ,j )可表示为ΔV (i ,j )Δtc k +1(i ,j )-c k (i ,j )+B r u k (i +1,j )c k r -B L u k (i ,j )c k L Δz j +B u (w k (i ,j +1)-w s (i ,j +1))c k u -B d (w k (i ,j )-w s (i ,j ))c kd Δx i =λx B r c k (i +1,j )-c k (i ,j )Δx i +1/2-λx B L c k (i ,j )-c k (i -1,j )Δx i -1/2Δz j +λz B u c k (i ,j +1)-c k(i ,j )Δz j +1/2-λz B d c k (i ,j)-c k (i ,j -1)Δz j -1/2Δx i (11)式中:ΔV (i ,j )为关于浓度的控制体体积;下脚标(i ,j )为变量对应于网格点的位置;u r ,u L ,u u 和u d 为流速在控制体的右、左、上和下各边中点的流速值;B r ,B L ,B u 和B d 为各边中点处的河宽;Δx i ,Δz j 为网格尺寸;Δt 为时间步长.温度方程的离散形式与式(11)完全相同,在此略去.计算第k +1时间步的水平流速分量u k +1(i ,j )的格式为ΔV ′(i ,j )Δtu k +1(i ,j )-u k (i ,j )+B ′r u k r u k r -B ′L u k L u k L Δz j +B ′u (w k u u k u -B ′d w k d u k d Δx i -1/2=-1ρ0[p k (i ,j )-p k (i -1,j )]B ′c Δz j -1ρ0τk w x ΔV ′(i ,j )+v x B ′r [u k (i +1,j )-u k (i ,j )]Δx i -v x B ′L [u k (i ,j )-u k (i -1,j )]Δx i -1ΔZ j +v z B ′u [u k (i ,j +1)-u k (i ,j )]Δz j +1/2-v zB ′d [u k (i ,j )-u k (i ,j -1)]Δz j -1/2Δx i -1/2(12)式中:ΔV ′(i ,j )为水平速度分量的控制体体积;相应的河宽记为B ′r ,B ′L ,B ′u,B ′d .在水平流速分量u k +1(i ,j )求出之后.利用连续方程可以直接求出垂直方向的流速分量w k +1(i ,j +1).普通的交错网格法在求解垂直流速时,采用与水平流速和浓度均不同的控制体.本研究考虑到河道立面二维模型的特点,采用与求解标量(温度、浓度)相重合的控制体求解垂直流速分量.计算先从河道底部的格子开始,利用河道底部垂直流速为零的条件和已知的水平流速,直接可以得到格子上边中点的垂直流速,这样从下到上直到自由水面.对于自由水面上的格子,利用连续方程可以直接得到水面的升降变化值.在求出新的水面位置后,就可以实现网格的自动调整,满足变动水位情况下的非定常模拟.详细离散公式如下:w k +1(i ,j +1)=w k +1(i ,j )B d B u-u k +1(i +1,j )B r -u k +1(i ,j )B L Δz j /(B u Δx i )(13)Δζk +1i /Δt =w k +1(i ,j )B d B u -u k +1(i +1,j )B r -u k +1(i ,j )B L Δz j /(B uΔx i )(14) 从上述空间离散格式中可见,由于在变量的布置上考虑了立面二维河道的特点,采用显式的方法求解各变量,并能够模拟自由水面的变化过程,因此适用于河道流速、水温以及污染物分布的中长期的非定常计算.由于水库中的流速比较小,在处理对流项时采用了中心格式.对于网格雷诺数大于2—22—的强对流情况,根据水平来流方向采用迎风格式.3　数值模拟结果实际河道由主流和支流组成,在应用立面二维模型求解河道中温度分层流动及悬浮污染物的分布时,应在几条河流交汇处考虑质量、动量守恒条件.在主流河段布置34×80个立面有限体网格,在支流布置15×80个立面网格.网格在流动方向尺寸约为120～200m ,在水深方向上为1m ,在河床高程以下的网格在计算中不加考虑.根据相关分析来确定来流的悬浮质浓度,其公式为:c in =αQβ(16)式中:系数α=01135,β=1146,是利用观测数据通过最小二乘法确定的.其他计算参数的取值为:v x =10m 2/s ;v z =0101m 2/s ;k z =110×10-5m 2/s ;ρs =2650kg/m 3;c p =411784kJ /(kg 1℃);λx =015m 2/s ;ρ0=2650kg/m 3;k x =10m 2/s ;λz =110×10-4m 2/s ;g =918m/s 2.其中粘性系数和扩散系数的选取,是经过原形观测数据检验确定的.计算时段为某年9月4日～11月30日,在这期间有原型观测水温和悬浮物质浓度分布可以用来验证模拟结果.图3(a )为计算开始后96h (9月8日)的计算结果和实际观测结果的比较,图3(b )是在计算开始后600h (9月29日)的结果.图中:浊度为残留悬浮物质浓度.(a )(b )图3　计算与观测结果的比较 从这些结果中可见,由于太阳辐射的影响,在水面附近水温较高,数值模拟能够再现温度分层现象.计算得到的水温值与观测结果符合良好.图4　残留污染控制总量随时间的变化在洪水期间,水库中悬浮物质浓度较大,如何尽快清除悬浮污染物质及其掌握水域中的水温分布是水资源保护的一个重要课题.为了控制水库中的悬浮质分布和残留总量,根据实际工程需要对以下几种取水口位置高程方案进行了模拟和比较:(1)现状:取水口保持在430m 高程;(2)表面取水:取水口随着水面的位置而变化;(3)高浓度取水A (高A ):在440～480m 高程内,取水口随着坝前最高浓度的位置而变化;(4)高浓度取水B (高B ):在430～480m 高程内,取水口随着坝前最高浓度的位置而变化.不同工况下水库中残留污染控制总量随时间的变化如图4所示,从中可见在现状运行条件下有大量悬浮物质停留在自由水面附近;采用表面取水,可以保持水面清洁,但是水库中悬浮物质总量增多;采用高A和高B运行方案,可以使得悬浮物质总量减小,并且保持较厚的清水层,相比之下高B方案在取水口上方还有少量悬浮质存在.同时,采用高A和高B运行方案可以减小水库中污染物的残留量,并且能够保持较厚的清水层存在.4　结论(1)应用立面二维流动模型、显式的数值格式模拟较长时期、较大流域范围内的水温、悬浮污染物的迁移过程,计算效率高.(2)模型考虑了自由水面的变化和河道宽度的变化,适用于模拟河道非定常流动问题.计算结果与实测资料符合较好.(3)水域中的悬浮物质残留总量可以通过优化取水口位置高程的运行方案进行控制.比较4种运行方案(现状取水、表面取水、高A和高B取水)数值模拟结果,表明高浓度取水方案最优.参　考　文　献:[1]　ASCE Task Committee on Turbulence Models in H ydraulic Computations.Turbulent Modeling of Surface WaterFlow and Transport:Part1[J].Journal of Hydraulic Engineering.1988,114(9):971-991.[2]　Van Riji L C.Field verification of22D and32D suspended sediment models[J].Journal of H ydraulic Engineer2ing.ASCE,1990,116(10):1270-1288.[3]　Falconer R A,Owens P H.Numerical modeling of suspended sediment fluxes in estuarine waters[J].estuarine,Coastal and Shelf Science,1990,(31):745-762.[4]　张耀新,吴卫民.剖面二维非恒定悬移质泥沙扩散方程的数值求解[J].泥沙研究,1999,(2):40-45.[5]　雒文生,周志军.水库垂直二维湍流与水温水质耦合模型[J].水电能源科学,1997,15(3):1-7.A22D unsteady flow model for predicting temperature and pollutantdistribution in vertical cross section of a riverJ IAN G Chun2bo1,ZHAN G Qing2hai1,G AO Zhong2xin2(11Tsi nghua U niversity,Beiji ng　100084,Chi na;21Chi na Instit ute of W ater Resources andHydropower Research,Beiji ng　100038,Chi na)Abstract:A model for predicting the thermal stratified flow and turbidity transportation in the vertical section of a river is established.Since the scheme is explicit for velocity components,pres2 sure and scalar quantities,the model is efficiency and able to simulate long period transportation process of suspending materials in reservoirs and rivers.The model can be used to forecast and control of pollution.The validity of this model is verified by in site observation data.K ey w ords:thermal stratified flow;vertical22D flow model;explicit scheme;pollution control ——42。

河道二维非恒定流场计算方法的研
究
“河道二维非恒定流场计算方法的研究”是水文学中的一个重要课题，其在研究河流的水文特征、水位变化和水流运动方面具有重要意义。

它涉及多种数学工具，如矩阵分解、图论、函数逼近、微分方程等，并且有助于理解河流的传输过程、水文特性和水力学过程。

河道二维非恒定流场计算方法是指在河道内部运用数学模型，用来模拟河道水文特征及水流运动的方法。

其基础是对河流形状、河床材料以及河床高程的详细分析，然后根据河道的实际情况建立模型，用来计算河道的水力学参数，如水力坡降、流量、流速等。

河道二维非恒定流场计算方法有两种基本形式：一种是基于河流形状和水力学参数的数学模型，另一种是基于水文参数的大气模型。

前者是指根据河流形状、河床材料以及河床高程等物理参数，来建立模型，用来研究河流水文特征及水流运动，而后者是指根据河流水文参数，如水位变化、流量、流速等，来建立大气模型，用来模拟河道的水力学特性。

河道二维非恒定流场计算方法的研究包括：首先，根据河流的实际情况，建立模型，使模型能够准确地反映河流的水文特性，并且模型的精度随时间变化；其次，利用矩阵分解、图论、函数逼近、微分方程等数学工具，求解河流二维非恒定流场的数学模型；最后，根据计算结果，探讨河流水文特性及水力学过程，以此来提高河流水文预报的准确性。

河道二维非恒定流场计算方法的研究是水文学中的一个重要研究课题，其能够帮助我们更深入地理解河流的传输过程、水文特性及水力学过程，从而有助于提高河流水文预报的准确性。

非恒定流泥沙数学模型原理及其应用泥沙研究1999（4）梁国亭，高懿堂，梁跃平，龚坚（黄委会黄河水利科学研究院; 三门峡水利枢纽管理局）摘要：本文介绍了一维非恒定泥沙数学模型的原理和计算方法，并且对于泥沙数学模型中的一些关键技术等问题进行了较详细的讨论。

利用已建立的一维非恒定泥沙数学模型，对黄河小北干流1981年洪水资料进行了模拟计算，计算结果表明：黄河干流、渭河、北洛河计算水位、流量过程线与实测值符合良好，可用于黄河的洪水演进计算。

关键词：非恒定流；泥沙数学模型；河床变形基金项目：国家自然科学基金和水利部联合资助重大项目(59890200)作者简介：梁国亭(1963-)，男，河南省淅川县人，黄河水利科学研究院高级工程师。

1 泥沙数学模型基本方程明渠或天然河流常被考虑作为一维流动，根据洪水波运动的圣维南方程、泥沙连续方程和泥沙扩散方程，可以简化推导出一维非恒定流泥沙数学模型的基本方程，其形式为水流连续方程(1)水流动量方程(2) 泥沙连续方程(3) 不平衡输沙方程(4)挟沙力方程S=f(Q,A,B,ω,S) (5)*其中 Q为流量，A为断面面积，B为断面宽度，Z为水位，K为流量模数，S为断面平均含沙量，S*为断面平均挟沙力，g为重力加速度，α为泥沙非平衡恢复饱和系数，ω为泥沙颗粒沉速，Ad 为断面冲淤面积，γ′s为泥沙干容重。

2 计算方法简介一维非恒定泥沙数学模型的计算采用非耦合方法，首先求解水流连续方程和动量方程，然后求解水流挟沙力、泥沙不平衡输沙方程和泥沙连续方程，具体求解过程如下。

2.1 水流方程的求解首先利用Preissmann隐式差分格式将水流连续方程和动量方程转化为差分方程，然后对差分方程进行线性化，在线性化过程中，略去增量的乘积项，最后得到以下线性方程组A 1j ΔQj+B1jΔZj+C1jΔQj+1+D1jΔZj+1=E1j(6)A 2j ΔQj+B2jΔZj+C2jΔQj+1+D2jΔZj+1=E2j(7)其中 Aij 、Bij、Cij、Dij、Eij（i=1,2）为第j单元河段差分方程的系数（j=1,2,......,N-1,其中N为断面个数）。

黄河水沙数学建模黄河作为中国第二长河流，其水沙问题一直备受关注。

水沙是指河水中所携带的泥沙颗粒的含量，它对河流的水质、水量、泥沙输移等方面都有重要影响。

因此，对黄河水沙进行数学建模是研究黄河流域水资源管理和生态环境保护的重要手段。

黄河水沙数学建模是通过对黄河流域的水沙数据进行统计分析，建立数学模型来描述黄河水沙变化的规律。

首先，我们需要收集大量的实测数据，包括黄河各测站的水位、流量、泥沙含量等信息。

然后，根据这些数据，我们可以进行统计分析，研究黄河水沙的时空分布特征。

在建立数学模型时，我们可以采用回归分析、时间序列分析等方法，通过对水沙数据的拟合和预测，来揭示水沙变化的规律。

例如，我们可以建立水位与泥沙含量之间的数学关系，以及流量与泥沙输移速率之间的数学关系。

这些数学模型可以帮助我们预测未来的水沙变化趋势，为黄河流域的水资源规划和生态环境保护提供科学依据。

除了建立数学模型，黄河水沙数学建模还可以通过地理信息系统（GIS）技术来进行空间分析。

利用GIS技术，我们可以将黄河流域的水沙数据进行空间插值，得到水沙分布的空间图像。

这些图像可以直观地展示黄河水沙的空间变化规律，为水资源管理和环境保护提供决策支持。

黄河水沙数学建模的研究还可以结合其他相关领域的知识，如水文学、地质学、生态学等。

通过综合分析不同学科的知识，我们可以更深入地研究黄河水沙问题，探索黄河水沙变化的影响因素和机制。

黄河水沙数学建模是一项复杂而重要的研究工作。

通过建立数学模型和利用GIS技术，我们可以揭示黄河水沙变化的规律，为黄河流域的水资源管理和生态环境保护提供科学依据。

这项研究的成果将对黄河流域的可持续发展和生态安全具有重要意义。

2001年10月水 利 学 报SH UI LI X UE BAO 第10期收稿日期:2000208230基金项目:国家自然科学基金及水利部联合资助重大项目(59890200).作者简介:张修忠(1972-),男,山东临沂人,博士生.文章编号:055929350(2001)1020082206河道及河口一维及二维嵌套泥沙数学模型张修忠1,王光谦1(11清华大学水沙科学教育部重点实验室,北京　100084)摘要:建立了一种河道及河口一、二维嵌套的泥沙数学模型,对基本的控制方程、方程的离散和求解方法、嵌套连接条件以及非均匀沙的处理等问题进行了研究.以非恒定非均匀不平衡输沙理论作为本文建模的基础,为方便处理二维计算域的不规则边界,采用有限元数值离散格式.验证算例对河道做一维简化,对口外海域做二维处理,通过交界面的水位、流量和含沙量等的传递,在每一迭代步内进行耦合计算.数值模拟结果与实测资料吻合较好,且计算省时,表明本文建立的嵌套模型是一种解决某些实际工程问题的可靠的和高效的工具.关键词:河口;泥沙输运;嵌套连接;有限元离散中图号:T V149 文献标识码:A泥沙数学模型作为研究和解决河流、水库和近海等水域的水流运动和泥沙冲淤问题的有效工具,已得到了较为普遍的应用.一维模型计算省时,可快速方便地进行长河段、长时期的洪水和河床演变预报,但无法给出各物理量在平面范围的分布,因而在模拟河床细部变形、河口和港湾等水域的流动和冲淤问题时,显得无能为力.水深积分的二维模型克服了一维模型的缺陷,但因计算量剧增,模拟长河段、长系列、平面大范围的水流运动和河床演变问题时很不经济,即使是短时期问题也不易做到实时预报.因此,将一维和平面二维模型嵌套连接,发挥其各自的优势,对于解决许多生产问题是必要的和有意义的.文献[1]在这方面做了比较细致的研究工作,文献[2]应用一、二维嵌套技术成功的模拟了黄河口的演变.1　水流泥沙数学模型及其求解方法111　河道一维非恒定流水沙方程　河道水流运动的圣维南方程:5A 5t +5Q 5x=0(1)5Q 5t +55x (Q 2A )=-gA 5ζ5x -gA Q 2K 2(2)悬移质不平衡输运方程及河床变形方程:5(AS k )5t +5(QS k )5x=-αωk B (S k -S 3k )(3)γ′5A sk 5t=αωk B (S k -S 3k )(4)式中:A 、B 、Q 、ζ分别为河道的过水面积、河宽、流量和水位;K 为流量模数,由谢才公式计算;S k 、S 3k 、ωk 、A sk 分别为第k 粒径组泥沙的含沙量、挟沙力、沉速及冲淤面积;α为恢复饱和系数;γ′为淤积物干容量;x 、t 为空间和时间变量.112　口外平面二维水沙运动基本方程　对于平面大范围的自由表面流动,由于水深尺度一般远小于水面尺度,可以引入浅水假定以简化基本守恒方程.假定压力沿水深服从静压分布,对基本方程(N 2S 方程)沿水深积分,可得到如下守恒型的浅水方程:55t (h )+55x j (q j)=0(5)55t (q i )+55x j (u j q i )=f δij q j -gh 5ζ5x i -λq i +55x j (νt 5q i 5x j)(6)悬移质不平衡输沙方程和海床变形方程:55t (Ψ)+55x j (u j Ψ)=αωk s 3k -βΨ+55x j (εs 5Ψ5x j)(7)γ′5Z b 5t =∑N s k =1αωk (s k -s 3k )(8)式中:u j 、q j 为x j 方向的平均流速和单宽流量;f 为柯氏力系数;δ为系数矩阵,除δ12=1和δ21=-1外其余元素均为0;λ=g u j u j Π(C 2h );C 为谢才阻力系数,可由曼宁公式计算;β=αωk Πh ;Ψ=hs k ;涡粘性系数νt 由νt =κu 3h Π6.计算,κ为卡门常数,u 3为摩阻流速;泥沙紊动扩散系数εs 假定与水流涡粘性系数相等;h 表示水深;水位函数ζ由水深和床底高程确定,即ζ(x j ,t )=h (x j ,t )+Z b (x j ,t );i ,j =1,2.113　水流挟沙力　潮汐河口挟沙力可由下式表示[3]:s 3=K V 2gh (9)在风、浪和潮流联合作用下,流速应该是风、浪和潮的合成流速,即:V =| V T + V b |+|V w |(10)式中:V T 为潮流速度;V b 为风吹流的平均速度;V w 为波流的平均速度;V b =0102W ,W 为平均风速;V w =012ch c Πh ,c 为波速;h c 为波高;K 为率定系数.114　基本方程的有限元离散11411　河道单元的离散　河道单元的流动守恒方程和泥沙输运方程可写成如下统一形式的对流方程:5φ5t +5(U φ)5x=F (11)式中:φ=[A ,Q ,AS k ]T ,F =[0,-gA (5ξ5x +Q 2K2),-αωk B (S k -S 3k )]T ,U =Q ΠA .对流方程的有限元离散可写成:M 5φi 5t =C ij φj +F i (12)式中:M 表示集中质量矩阵,M =∫ΩN i N j d Ω;C ij 表示对流矩阵,C ij =-∫Ωw i5UN j 5xd Ω;F i 代表源项,F i =∫Ωw i F d Ω.11412　口外平面二维单元的离散　有限元在本质上属于非结构化网格离散方法,便于处理复杂边界问题.因此,本文对控制方程采用有限元法离散,方程(5)～(7)的弱解形式经分部积分后可得如下的空间半离散方程:M 5h i 5t=C ij h j (13)M 5<i 5t=C ij <j +D ij <j +F 1+M ・F 2<i (14)<=[q x ,q y ,hs ]T ,F 1=[F x ,F y ,F z ]T ,F 2=[-λ,-λ,-αωΠh ]T式中:Cij为对流矩阵;D ij为扩散矩阵;由下列各式表示:C ij=-∫Ωw i5uN j5x+5νN j5y dΩD ij=-∫Ωνt5N i5x5N j5x+5N i5y5N j5y dΩF x=-∫Ωgh5ζ5x w i dΩ　F y=-∫Ωgh5ζ5y w i dΩ　F z=∫Ω(αws3)w i dΩ式中:N、wi分别为插值函数与权函数.若上述离散中的权函数与插值函数相等,则构成经典的G alerkin有限元法.对于对流占优问题, G alerkin法等价于中心差分格式,因缺乏足够的耗散,往往导致数值振荡.为此,本文采用高分辨率格式对流项进行重构,即通过引入几乎相等的扩散与反扩散以保证格式的高精度,同时利用限制因子保证影响系数的非负性及解的保单调性[4].115　离散方程的求解　为使计算收敛或加快收敛,离散中对源项M・F2<进行负坡线性比.离散后的方程(13)、(14)为常微分方程,可采用多种显式或隐式方法求解.本文对时间导数项采用C2N格式离散,对离散后的代数方程采用QMR[5]方法迭代求解,该方法具有节省内存,收敛快的优点.本文顺序求解离散后的方程.对于二维海域,先由二维对流扩散输运方程(6)计算流速,由对流输运方程(5)计算水深,由悬沙对流扩散方程(7)计算含沙量,最后由河床变形方程(8)计算节点的冲淤深度.对于河道一维计算,运动方程的单元离散转化为求解流量Q的方程,连续方程的单元离散方程则转化为求解水位ζ的方程.2　嵌套连接条件一、二维嵌套模型是通过交界面连接的,由于一维模型只给出物理量的断面平均值,二维模型给出节点的水深平均值,因此在交界面上存在各物理量断面平均值和节点垂线平均值的相互转化和衔接问题.水沙运动在交界面上的连续性是模型嵌套连接的基本原则,因此,一、二维嵌套的连接条件是:水位相等,即:ZB=∫B0z d y.式中:Z为断面平均水位;z为节点的水位;B表示交界断面的水面宽度.流量相等,即:Q=∫B0uh d y.式中:Q为通过交界面的流量:u、h分别是节点的垂线平均流速和水深.悬移质输沙量相等,即:SQ=∫B0uhs d y.式中:S、s表示交界面的断面平均含沙量和节点垂线平均含沙量.此外,还有阻力、挟沙力、河床变形等连续条件.本文口门水位由二维控制,流量由河道一维计算给出,含沙量进行相互传递.3　非均匀沙水流挟沙力级配及床沙级配计算311　挟沙力级配计算　鉴于水流中的泥沙源于上游水流挟带和床沙紊动扩散进入,因此由水流条件和床沙组成推求非均匀沙的分组挟沙力的做法是较为合理的.本文采用李义天通过建立输沙平衡状态下的床沙质级配和床沙级配间的关系以及垂线平均悬沙浓度和河底悬沙浓度之间的关系得到的挟沙力级配公式[6].312　床沙级配计算　床沙级配随河床冲淤而变化,对阻力、输沙率及河床冲淤影响显著.若已知各粒径组泥沙的冲淤厚度ΔH sk及总的冲淤厚度ΔH s,则床沙级配调整计算可分为以下两种情况.(1)ΔH s>0,即发生淤积的情况,床沙活动层级配由下式计算:ΔP t+Δt bk =[ΔHs k +ΔP t bk (H t m -ΔHs )]ΠH t +Δtm 式中:ΔP t bk 、ΔP t+Δt bk 分别为t 时刻和t +Δt 时刻的床沙活动层级配;H tm 、H t+Δt m 为相应时刻的床沙活动层的厚度.(2)ΔH s <0,即发生冲刷的情况,床沙活动层级配由下式计算:ΔP t+Δt bk =[ΔHs k +ΔP t bk H t m +|ΔHs |ΔP remk ]ΠH t +Δt m式中:ΔP t bk 、ΔP t+Δt bk 、H tm 、H t+Δt m 同前,ΔP remk 为若干个记忆层内的床沙平均级配.床沙活动层是指河床发生冲淤变化过程中,河床表层参与河床冲淤变形的那一层床沙,水流挟带的泥沙与床沙的交换在这里发生,河床冲淤变形也在这一层里发生[7].床沙活动层厚度是指一个冲淤计算时段内感受到水流作用并且泥沙组成发生变化的床沙厚度[7,8].受河床变形和来水来沙条件的影响,活动层的厚度和组成不断变化.由于问题的复杂,要从数学上严格定义和表达活动层厚度目前还比较困难.尽管有许多学者对这一问题进行了研究,给出了一些计算方法.但这一问题不能考虑过细,一方面缺乏床沙级配沿垂向变化的实测资料;另一方面考虑过细不一定能提高精度.故本模型采用较常用的处理方法,即取一固定值2m.313　床沙级配的分层记忆模式　为了模拟床沙组成的变化过程,将床沙划分为床沙活动层及其下面的记忆层[8]两部分.记忆层可根据实际情况分n 层.计算中,当河床发生淤积时,记忆层层数增加,增加层的级配为t 时刻的床沙活动层级配ΔP t bk .当河床发生冲刷时,根据冲刷量的大小,记忆层数相应减少,且级配作相应的调整.4　模型的验证411　计算条件　漳卫新河是海河流域南系的一条尾闾河道,担负着漳河、卫河的泄洪排涝任务.自1973年扩大治理以来,由于入海径流少,辛集闸闸下河道被潮汐动力所控制,源源不断的海相来沙使河道严重淤积.据94年实测地形资料分析[9],淤积河道长达26km ,淤积总量达到1262万m 3,河道行洪能力下降47%.本文验证计算的河道一维计算域取自漳卫新河的辛集闸至河口,长3716km ;口外海域的下边界至-20m 等深线,纵向长30km ,横向宽20km.河道地形资料采用94年大断面资料,口外地形采用1∶50000海图.412　边界条件　河道进口给定流量、含沙量过程,口外开边界条件采用潮位控制,岸边界采用水流无滑移条件.口外各角点水位由实测潮位根据潮波传播相位差推延得到,再根据域内测点流速过程验证情况稍作调整,以90年实测大潮潮型概化计算潮型.一、二维连接断面采用流速边界,并按曼宁公式进行分配.413　有关参数的选取[10]　在现有的认识条件下,河口水沙预测的关键是选取可靠的基本参数,如糙率、挟沙力系数和泥沙恢复饱和系数等.为此,需对河口现状水流泥沙条件进行验证,它一方面是对数学模型本身的检验,另一方面也是率定水流泥沙基本参数,为各方案科学预报提供依据.河道糙率采用01025,河口二维海域糙率采用0102;水流挟沙系数采用海河口数据K=100;淤积物干容重取0165t Πm 3;根据验证计算确定泥沙恢复饱和系数α冲=011,α淤=0125;波高取大口河测波站平均波高.图1　计算与实测潮位过程对比图2　计算与实测流速过程对比414　验证计算结果　图1～3给出了94年8月26～27日河口处的潮位、流速和含沙量计算与实测的对比,图中零时刻对应于26日14时.可以看出计算与实测潮位、流速吻合良好,表明水流计算参数的选取是合理的,计算方法也是可靠的.含沙量过程计算与实测有一定差别,主要是由于在潮流和波浪共同作用下泥沙参数的选取还有待进一步改进.图4给出了83年～93年河道累计淤积量计算与实测的对比,全河段累计淤积量计算值约590万m3.图5～6给出了河口局部涨急和落急流场,可以看出,涨潮流速明显大于落潮流速,与实测资料一致.这也是涨、落潮输沙不平衡,河道淤积的一个重要原因.受资料限制,口外海床变形未作验证.图3　计算与实测含沙量过程对比图4　计算与实测河道累计淤积量对比图5　涨急局部流速矢量场图6　落急局部流速矢量场5　结语对口外海域进行平面二维计算,对河道采用一维模拟;或者对河道流动复杂段应用二维模型,流动简单或顺直河段应用一维模型的一、二维耦合算法,既具有一维模型的快速方便,又能获得局部河段或平面大范围的细部信息.这样可以用较少的机时复演和预测长河段的河床变形及其重点段的细部变形,是一种解决某些实际工程问题的有效方法.致谢:论文得到大连理工大学土木系金生教授的指导和帮助,在此表示衷心的感谢.参　考　文　献:[1]　Wu W M,Li Y T.One2and T w o2Dimensional nesting mathematical m odel for river flow and sedimentation[C],5thInternational sym posium on river sedimentation,1992,K arlsruhe547-554.[2]　Zhang S Q.One2D and T w o2D combined m odel for estuary sedimentation[J],Int.J.Sediment Research,1999,14(1):37-45.[3]　刘家驹,张镜潮.淤泥质海岸航道、港地、淤积计算方法及其应用推广[J].水利水运科学研究,1993(4).[4]　张修忠,王光谦,金生.浅水流动有限元分析及其高分辨率格式[J].长江科学院院报,2001(1).[5]　Freund R W,Nachtigal N M.An im plementation of the QMR method based on coupled tw o2term recurrence[J].SI2AM.J.Sci C om put.,1994,15(2):313-337.[6]　李义天.冲淤平衡状态下的床沙质级配初探[J].泥沙研究,1987,(3).[7]　李义天,胡海明.床沙混合活动层计算方法探讨[J].泥沙研究,1994,(1):64-71.[8]　吴卫民,等.河床床沙组成数值模拟方法[J].武汉水利电力大学学报,1994,27(3):320-327.[9]　王文治,梁永立.漳卫新河冲淤变化及发展趋势的分析[R].水利部天津勘察设计研究院,1996.[10]　金生.漳卫新河河口泥沙冲淤计算[R].大连理工大学,2000.12D and22D nesting sediment transport model for rivers and estuariesZH ANG X iu2zhong1,W ANG G uang2qian1(11T singhua Univer sity,Beijing　100084,China)Abstract:A12D and22D combined sediment transport m odel for rivers and estuaries is presented.The basic equation,the numerical method,the coupling conditions and the treatments for non2uniform sediment are stud2 ied.The m odel is based on unsteady non2uniform and non2equilibrium sediment transport theory.The finite el2 ement method is abopted to solve the g overning equations for its capability of accepting com plex geometry.The m odel is verified by the simulations of the flow and sediment transport in the estuaries of the Zhangweixin River,in which the river area is treated as12D and the sea area is treated as22D.By trans ferring the water level,dis2 charge and sediment concentration at the interface,the coupling calculations are conducted in each iterative step.The results are in g ood agreement w ith the measured data and a lot of CPU time is saved,which shows that the proposed m odel is reliable and high efficiency in solving practical engineering problems.K ey w ords:estuary;sediment transport;nesting linking;finite element discretization。

河道平面二维水沙数学模型的有限元方法内容摘要：摘要：采用有限元方法建立起一套河道平面二维水流泥沙数学模型。

在前人研究的基础上，采用了质量集中的处理方法，提出了压缩存储的方法，从而大大减少了计算存储量。

针对有限元法时间步长需取得较短问题，采用了“预报-校正-迭代”的算法，提出了“非恒定-恒定-非恒定流”的算法，既能解决工程实际问题，又大大减少了计算量。

以下荆江监利河段为例进行泥沙冲淤计算，计算结果与实测值符合较好，从而证明了模型的可靠性。

关键词：水流泥沙有限元模型验证三峡工程建成后，水库将拦蓄大量泥沙，下泄水流含沙量减小，对三峡工程坝下游河道将产生以冲刷为主的影响，包括对荆江河段的河势及荆江大堤带来影响。

为研究坝下游重点河段的河床冲淤分布、河势变化、近岸流速变化等问题，一维模型显得无能为力，但可采用平面二维模型来解决。

有限元方法可采用无结构化网格，能很好地模拟不规则的几何形状，因此很适合于对天然河道的模拟。

然而，正如其它方法一样，有限元法也有它的缺点，主要是计算存储量和运算量较大。

为扬长避短，使有限元方法能运用到对天然河道的模拟上来，本模型运用质量集中［4］的方法将系数矩阵转化为三对角矩阵，并提出了紧凑的分块压缩存储方法，从而大大减少了计算存储量，使得计算能在一般微机上进行。

采用质量集中方法的不足之处是时间步长需取得较短，且在河道模拟中尤为突出(因河道比较窄长，网格需划分很细，而该法的稳定性要求时间步长与网格尺度成正比)。

针对该问题，笔者采用了“预报-校正-迭代［5］”的算法，该法可加大时间步长，同时有效避免了数值震荡。

针对长系列水沙条件下计算量较大问题，又提出了“非恒定-恒定-非恒定流”的算法，该算法既能解决工程实际问题，又大大减少了计算量，使有限元方法能够很好地运用于河道水流泥沙问题的实际计算。

1基本方程平面二维水流方程(1)(2)(3)悬移质泥沙扩散方程(4)推移质不平衡输移方程［6］(5)河床变形方程由悬移质引起的河床变形方程为<![endif]>(6)由推移质引起的河床变形方程为(7)以上各式中U，V分别为垂线平均流速在x,y方向上的分量；Zs、Zb和H分别为水位、河底高程和水深；g为重力加速度；vt为水流紊动粘性系数；ρ为水的密度；τx、τy、舄瓂分别为底部切应力在x和y方(τx、τy)=，向上的分量：C 为谢才系数，常用曼宁公式计算：C=H1/6/n；S和S*分别为垂线平均含沙量和挟沙力；N和N*分别为推移质输沙量和推移质输沙能力折算成全水深的泥沙浓度；εs为泥沙紊动扩散系数；ω为泥沙沉速；γ′为床沙干容重；α为悬移质泥沙恢复饱和系数，淤积时取0.25，冲刷时取1.0；β为推移质泥沙恢复饱和系数，取0.25。

黄河口泥沙异重流基本控制参数的数值试验摘要：为了定量地评价泥沙异重流基本控制参数对泥沙异重流流动的影响，本文设计了一系列数值试验来展示各种参数与泥沙异重流流动特征之间的定量关系。

采用ADI法求解三维方程经垂向积分后的平面二维泥沙异重流方程。

所选取的基本控制参数包括水下斜坡坡角，总阻力系数，水力挟带系数，泥沙饱和浓度和泥沙沉速。

通过数值试验得出水下斜坡坡角θ、泥沙饱和浓度Es越大，总阻力系数cd、泥沙沉速ω越小，相应地泥沙异重流流速、流动厚度以及悬沙浓度越大。

其中斜坡坡角对泥沙异重流影响最大。

当泥沙异重流流动处于超临界流状态，水力挟带系数Ew大于零，这时由于有一定水量的卷入，所以流动厚度沿程增加。

关键词：黄河口泥沙异重流控制参数数值试验 ADI法泥沙异重流及形成的地质体在沉积学和石油地质中有重要意义。

在黄河口，泥沙异重流是泥沙搬运的重要方式，是河口区沉积动力的重要过程，是认识河口三角洲发育演化的重要途径。

为使研究深入一步，我们进行了泥沙异重流的数值模拟。

近年来，国内对异重流数模研究较少，仅见方春明等(1997)关于泥沙异重流潜入时的立面二维数值模拟以及邱晨霞(1995)对盐水异重流的二维两层数值计算。

1 平面二维泥沙异重流控制方程及数值模拟目前，包含三个基本控制方程(流体质量、动量守恒方程和泥沙质量守恒方程)的泥沙异重流数值模拟被广泛应用，本文把一维方程推广到平面二维的情形，成为四方程模型，如下所示[1](1)(2)(3)(4)其中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为柯氏力参量。

黄河下游河道准二维泥沙数学模型研究摘要：本文首先简要回顾了以往黄河河道泥沙数学模型研究成果，再以水流运动方程及经过作者修正的泥沙运动方程为基础，同时引入与实测资料相符合的水流挟沙力、动床阻力、泥沙级配等计算公式作为补充方程，构造出黄河下游河道准二维泥沙数学模型。

然后，采用1986年11月～1996年10月这10年长系列实测资料，开展了验证计算。

其结果表明，该模型不仅能计算黄河下游河道一般洪水引起的河床冲淤变形，还能成功地模拟出大沙年下游处于强烈淤积时的规律。

关键词：黄河下游河道准二维泥沙数学模型1 黄河河道数学模型研究的简要回顾黄河数学模型的研究，是与流域规划、工程建设和管理运用等生产紧密结合的。

早在1955年编制黄河综合利用规划技术经济报告中，就曾在黄河三门峡水库规划阶段用初级的一维恒定平衡输沙模型对水库淤积和下游河道的河床冲刷变形进行计算[1]。

当时的计算结果认为，在桃花峪下游冲刷9年后河床刷深27m，冲刷量及冲刷速度显然比实际夸大甚多。

三门峡水库投入运用后，库区严重淤积，并迅速向上游延伸，与原来的计算结果截然不同。

为研究改建方案，通过实测资料建立了许多泥沙冲淤量与水文要素之间的经验关系式。

麦乔威、赵业安、潘贤娣等学者在分析下游河道冲淤变化及挟沙能力变化规律的基础上，提出了三门峡水库下游河床冲淤计算的方法[1]。

在1975年～1977年进行治黄规划过程中，为研究各种规划方案对黄河下游的减淤作用，黄河水利委员会及其有关科研、设计部门的科技人员，使用两种河道冲淤计算的经验模型开展了大量的计算工作，麦乔威、李保如两位学者进行了总结。

其中，在规划方案的初步比较阶段，使用的是李保如提出的框算模型，只考虑影响输沙能力的最主要因素（来水、来沙）。

对于问题比较复杂和初选的方案，则使用较详细的经验模型[2]。

20世纪80年代，刘月兰、韩少发、吴知等学者在上述经验模型的基础上，根据不同河段汛期、非汛期输沙经验公式，引入水流连续方程、动量方程及泥沙平衡方程式，建立了一套颇具特色的黄河下游河道冲淤计算模型[3]。

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

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

关键词：黄河口；挟沙能力；窄缝法；非恒定流；数学模型收稿日期：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 潮流和波浪共同作用下泥沙运动方程在河口海岸地区，潮流和波浪是泥沙运动的最主要动力。

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

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

1998年4月水 利 学 报SH UILI XUEBAO第5期黄河下游河道水沙运动仿真模型的开发研究程晓陶　黄金池(中国水利水电科学研究院)薛云鹏(黄河水利委员会河务局)摘　要　本文介绍的黄河下游宽浅河段汛期水沙运动仿真模型建立在二维不恒定流与不平衡输沙理论的基础上,采用无结构不规则网络,结合了有限差分法和有限体积法的优点,能较充分地反映滩区、生产堤与滩槽冲淤变化对洪水演进的显著影响,并具有实时动态图象显示和人机对话功能.关键词　二维不恒定流,水沙运动仿真.黄河下游滩宽槽浅,汛期来水来沙变幅很大.大堤间滩地面积占河道面积83%,有居民近140余万.滩区生产堤纵横,在纵比降为1/6000的河道内,滩面横比降可达1/2000;平滩流量下降到约4000m3/s.洪水一旦破堤漫滩;顶冲大堤,顺堤行洪,甚至倒灌,对滩区人民生命财产及黄河大堤的威胁更大.漫滩行洪和滩槽水沙激烈交换可引起洪峰严重变形,且河槽冲淤幅度与水位涨落达到同等量级.大幅快速的河床调整又反过来影响水流的挟沙能力与河道的行洪能力,引起水位-流量关系的复杂变化,使洪水预报异常困难.对这种大范围的、在计算域内部受堤防溃决等突发性因素影响十分显著、水沙交换激烈的动边界、动床洪水现象,一般一维或一维扩展模型很难模拟.已有的洪水预报模型或河势演变模型对漫滩洪水大都作了简化,难以充分体现生产堤溃决、滩地行洪与滩槽水沙交换对洪水传播变形的影响.新研建的黄河下游河道水沙运动仿真模型,既要求计算精度高,以便为滩区洪水风险分析提供科学依据;又要求运算速度快,在大洪水发生时能及时准确为滩区居民发布洪水预警报信息.如此矛盾的要求增加了模型设计的难度.目前,黄河下游夹河滩至高村(包括封丘倒灌区)防洪险要河段的水沙运动仿真模型已通过鉴定,其主要特点包括:①以一、二维不恒定流与不平衡输沙理论为基础,采用无结构不规则网格,可以较好地概化河道主槽、滩地的地形特征,充分反映生产堤、渠堤等对水沙运动的影响.②泥沙计算在黄金池方法的基础上作了改进,引入了张红武的一套挟沙能力公式[1],使黄河下游主槽冲淤变化对洪水演进的影响能够较好体现出来,并能适应高含沙洪水的计算.③模型具有在运算过程中实时动态显示流量、水位、洪水淹没范围或泥沙冲淤计算成果的功能.④模型中设计了人机对话的窗口,对生产堤在洪水期间是否发生决口这种带有随机性的问题,模型使用者有机会根据实际情况或经验进行选择.⑤模型运算速度显著提高,对1427.8km2的计算域模拟6天洪水过程,仅需十余分钟机时,可满足为洪水预报服务的要求.⑥开发了模型的菜单对话式图形信息管理系统,使模型所需大量基本信息的前处理与计算成果的后处理工作大为简便.⑦计算结果与损失评估模型相结合,可以快速算出相应洪水条件下的受灾人口及可能造成的经济损失.1　水沙运动的计算方法1.1　基本控制方程　从黄河水沙运动特点出发,模型中选垂向平均单宽流量Q与水深H作为基本状态变量,其深度积分形式的控制方程组可记作如下的质量与动量守恒形式:本文于1996年7月1日收到.水流连续方程 H t + M x + N y=q ,(1)水流动量方程 M t + (uM ) x + (v M ) y +gH Z x +g n 2u u 2+v 2H 1/3=0,(2)N t + (uN ) x + (vN ) y +gH Z y +g n 2v u 2+v 2H 1/3=0,(3)式中:H —水深;Z —水位;q —源汇项,本模型中代表有效降雨强度,不计降雨情况下q =0;M ,N —分别为x 、y 方向的垂向平均单宽流量;u ,v —分别为垂向平均流速在x 、y 方向的分量;n —曼宁糙率系数;g —重力加速度.基本状态变量与其他交量之间有以下转换关系:Z =H +B ,(4)u =M /H ,(5)v =N /h ,(6)式(4)中B 为地面高程.在考虑泥沙冲淤计算时,基本控制方程组中要增加泥沙运动方程与河床变形方程.根据不平衡输沙理论,这两个方程通常可表达为如下形式:泥沙运动方程(HS ) t + (HuS ) x + (HvS ) y -e s 2 x 2(HS )+ 2 y 2(HS )高含=a ω(S ＊-S ),(7)河床变形方程 r ′ B t=αω(S ＊-S ),(8)式中:S ,S ＊—分别为深度平均含沙量与水流挟沙力;ω—悬移质泥沙深度平均沉速;α—悬移质泥沙恢复饱和系数;e s —悬移质泥沙扩散系数;γ′—淤积物干容重;其余符号与水流模型意义相同.1.2　水流控制方程的简化与离散方式　黄河滩区计算域大,河势复杂多变,控导工程、生产堤、渠堤密布其间.在规则网格情况下,若为省计算时间而将网格空间步长放大,则会引起严重的地形失真,难以满足洪水风险评价对计算精度的要求.若加密网格以提高对地形的分辩率,又必然会增大计算量,难以满足缩短运算机时、为洪水预报服务的要求.为此,本模型在开发中保留了体积积分形式的有限差分法的优点[2],如在网格周边通道上计算流量、泥沙通量,在网格的形心计算水深、含沙量、挟沙力与底床的冲淤变化;同时又吸收了有限体积法中非结构性不规则网格的优点[3].通道上定义的单宽流量,不限于正交的X 或Y 的方向,而可以是任意的通道法线方向的流量,同时,水深与流量在时间轴上分层布置,交替求解,物理意义清楚,并有利于提高计算稳定性,如图1所示.a .H 与Q 的空间布置方式b .H 与Q 的时间交错计算方式图1　基本状态变量的离散化布置方式 将式(5),(6)代入式(1),则水流连续方程式改写后对计算域进行面积分,再利用高斯定得,可得∫A H t d A +∮(H u ·n )d l =∫A q d A ,(9)式中:u 为计算域边缘上任一点的流速矢量,n 为该点的外法线方向单位向量.令Q =H u ·n ,则Q 为任意n 方向垂向平均单宽流量,当n 取为x 、y 方向的单位向量时,Q 即为M 、N .由于黄河下游河道比降较小,滩区地势平坦开阔,洪峰流量靠河宽变化调整,水位变幅不大,可以假定水深H 随时间的变化在一个有限的网格内差异很小,降雨在一个网格内也是均匀分布的,则式(9)简化成:A H t +∮LQ d l =Qq .(10) 对任一K 边形网格,上式左端第2项的线积分可记为:∮L Q d l =∑K k =1Q k L k ,(11)式中,A 为网格面积,L 为通道长度,下标k 为对K 边形网格逆时针方向的通道编号.设对任意网格,Q 流入为正,流出为负,则连续方程对i 网格的显式离散化形式可以写成H T +2D Ti =H T i +2DT A i ∑K k =1Q T +D T ik L ik +2D Tq T +D T ,(12)式中下标i 为网格的编号.式(12)表明,在已知T 时刻水位的情况下,为求出网格中T +2D T 时刻的水位,只要能合理确定各通道上T +D T 时刻水流的单宽通量Q 既可.为此,针对黄河特点,将通道分为河道型、滩地型、有堤型等,采用不同假设进行简化计算.河道型通道设在主槽中.对于黄河下游平原型河道,动量方程中仅忽略了对流项,即假定第j 通道上的垂向平均单宽流量Q j 不仅与相邻两网格的水位Z j 1,Z j 2有关,而且与前一时刻流量有关,则动量方程的离散化形式简化为:Q T +DT j =Q T -D T j -2D T ×gh j Z T j 2-Z T j 1DL j -2D T ×g n 2Q Q h 7/3j,(13)式中h j 为j 通道上的平均水深,一般由下式计算:h j =DC j 1H T j 2+DC j 2H T j 1有限/DL j ,(14)式中DL j =DC j 1+DC j 2,为空间步长,如图2(a )所示.考虑到黄河河道宽浅,有时需将一个横断面分成若干网格,而网格间高差较大,成阶梯状,式(14)计算的h j 可能过大.因此,判断河底比降大于0.0006时,h j 改由下式插值给出h j =DC j 1Z T j 2-B m g +DC j 2Z T j 1-B m c k p t g a /DL j ,(15)式中B m 为相邻网格中较高的河底高程.滩地型通道设在陆地上.漫滩洪水传播主要受重力与阻力的作用,可按扩散波运动进行处理.从动量方程中完全略去惯性项,其离散化形式简化为:Q T +D T j =sign Z T j 1-Z T j t o 2h 5/3j Z T j 2-Z T j 1DL j ,1/2·1n ,(16)式中DL 的定义为相临网格形心距在通道法线方向的投影,如图2(b )所示.当滩地水深达到一定深度后,顺滩行洪状况型等同于河道主槽.模型计算中将根据给定的判别条件由式(16)自动转换为式(13)的主槽算法.对有阻水建筑物的通道,可分别概化成连续堤防、有缺口堤防及有桥涵的堤防等几种情况,采用适当的堰流或孔口出流等经验公式进行计算. 根据模拟对象的特点,在网格类型上也区分了主槽、嫩滩、老滩.网格的非结构性设计,使其既便于针对主槽的游荡性进行计算域的调整,又便于今后计算域向花园口至陶城铺河段扩张;网格的不规则性,使地形概化更接近实际,能更好反映复杂的河势变化与工程条件对洪水演进的影响;网格分类上的灵活性,使其能针对计算对象的基本特点采取不同的简化方式,既抓住主要矛盾,又节省计算时间.a.河道型通道b.滩地型通道图2　不规则网格空间步长的分类处理 为适合黄河的特点,模型开发中对水流与泥沙计算的结合、网格形心流速的算法、滩唇的处理方法、出口断面无水位-流量关系及滩区横比降引起的数值效应等问题均作了有益的探讨,受篇幅限制,详细计算与处理方法将另文介绍,或请参见文献[4].2　模型的验证在模型开发阶段,先后以1982年与1973年两场洪水资料对模型进行了验证.前者是三门峡水库建成后黄河下游发生的最大一场洪水,代表了大流量、小含沙的情况;后者则代表了小流量、高含沙的情况.模型经多种验证均得到令人满意的结果[4,5].高含沙洪水的计算及河床冲淤成果的详细分析,笔者将另文介绍,本文仅以1982年洪水验证为例.模型验证范围取黄河下游夹河滩至高村河段,河道长83km,大堤间距最宽达20km,平均约11km,面积908km2,包括长垣、东明两个大滩,面积均在200km2以上.滩区生产堤完整,两岸生产堤堤距约3km左右.验证结果表明:①洪水淹没范围与实际相当吻合,而计算图更清楚地给出了整个计算域上的最大水深分布.计算成果可对任意时刻绘制出洪水平面二维淹没水深分布图与相应的流速矢量分布图[4];②合理模拟了洪水破堤漫滩后顶冲大堤,顺堤行洪、甚至倒灌的现象.比较图3、图4,前者较好模拟出了洪峰传播变形的特点;而后者由于洪水不能进滩,四个断面的洪水波形几乎不变,后期出口断面水位得不到滩地归槽洪水的补充,明显偏低;③计算最高水位在整个域中分图3(a)　1982年高村站计算与实测流量过程比较图3(b)　1982年高村站计算与实测水位过程比较布合理(图5);④滩槽冲淤分布合理,充分体现了大洪水期间“淤滩刷槽”的现象,主槽中出现了一条清晰的连续的冲刷槽,反映了黄河汛期“大水出好河”的特点.根据实测资料,“82.8”洪水期间,泥沙淤积主要分布在生产堤间嫩滩上,约占滩地淤积量的95%,而生产堤外滩地的淤积量很少,计算结果表明泥沙淤积也是绝大多数分布在生产堤间的滩地上.主槽与嫩滩冲淤幅度的量级是合理的,主槽冲刷深度约为1.0—2.5m ,嫩滩淤积厚度约为0.2—1.0m .图4(a )　假定生产堤不溃的流量过程模拟结果图4(b )　假定生产堤不溃高村站水位过程的比较图5(a )　1982年洪水黄河北岸大堤计算最高水面线与实测资料比较图5(b )　1982年洪水黄河南岸大堤计算最高水面线与实测资料比较3　模型的实时动态显示功能以往的数值计算模型总是先“盲算”,只有待计算结束后,再进行数据分析,才能知道计算结果.而黄河水沙运行模型计算域广、数据量大,为了改变模型调试、运行效率低的状态,并适应计算中人-机对话的需要,专门开发了一整套绘图子程序,与计算程序一起运行,使各种计算结果能够实时动态显示出来,包括:①以不同颜色表示的水深或泥沙冲淤深度的平面分布.两者之间在计算中可以随时用指定的功能键进行往复切换,想看水深分布还是冲淤分布,随意选择.②以不同颜色线段表示的若干指定河道断面的流量变化过程与若干指定地点的水位变化过程.在有实测资料的情况下,还可以显示计算结果与实测过程的比较。