二维水沙数学模型的建立和应用

1、整理ppt1整理ppt2摘摘 要要 二十世纪70年代以来，随着电子计算机的普及与计算机性能的不断提高，人们在研究河流泥沙问题时，越来越多的采用了数学模型。它以数值方法和计算机技术为手段，通过对河流的数值模拟计算，解决河流工程所关心的问题。对二维水流泥沙数学模型的探索和研究已在水利、环境、港口航道、防洪等工程领域中取得了一定成果，但水流泥沙运动机理仍是工程问题中的重要研究课题。 本文结合河北省水利水电勘察设计研究院河北省庙宫水库汛期调度运用方式项目，研究了平面二维水流泥沙数学模型的理论、方法和实际应用问题，主要工作如下： 1.系统地推导了二维非恒定流、二维泥沙连续方程，应用有限体积法思想，采用分

2、类简化离散方程的方法，降低了计算难度、节省了计算时间,为建立平面二维水沙数学模型奠定了基础。 2.在已有的泥沙研究成果的基础上，提出并推导了新的泥沙连续方程形式，结合适合庙宫水库情况的爱因斯坦水流挟沙公式，建立了平面二维水沙数学模型。 整理ppt3 3.根据离散方程，编制Fortran计算程序，对庙宫水库库区内的水流和泥沙冲淤进行计算并与已有的实测资料对比，吻合效果较好，说明模型可以用于该水库的河道变形模拟。对上游来水、来沙情况和水库排沙控制水位序列13年水沙组合过程进行计算模拟，给出了全河段的流场分布、含沙量分布、河床冲淤演变分布的计算模拟成果。 4.通过综合对比，定量分析了不同计算组合下床

3、面冲淤演变情况，计算出排沙比,提出了水库现状泄流条件下最优清淤方案。 研究成果不仅为理论研究提供了新思路，而且也为实际工程提供了可行的数值模拟方法,具有一定的理论创新和应用价值。 关键词关键词：平面二维水流泥沙数学模型；庙宫水库；泥沙连续方程；平面二维水流泥沙数学模型；庙宫水库；泥沙连续方程；冲淤演变冲淤演变整理ppt4目 录第一章 绪论 1.1泥沙数学模型发展简况 1.2水沙数学模型基本方程及计算方法 1.3数值计算方法 1.4水沙基本理论存在的问题 1.5淤积上延现象第二章 河道复杂边界的处理技术 2.1 斜对角笛卡尔方法的基本原理 2.2斜对角笛卡尔方法的优点 2.3边界条件及边界网格单

4、元类型处理 2.4 本章小节第三章 平面二维水流泥沙数学模型的建立 3.1 二维水沙数学模型的基本方程 3.2 初值条件及边界条件的处理 3.3 网格划分 3.4 水流挟沙能力第四章 平面二维水沙数学模型的数值模拟 4.1 水流控制方程的离散 4.2 泥沙连续方程的离散 4.3 河床变形方程的离散 4.4 模型求解过程整理ppt5第五章 水流泥沙数学模型验证及应 用 5.1 模型的基本特征 5.2 网格划分 5.3 模型验证 5.4 模型预报 5.5 排沙量及冲淤分析第六章 工作与展望 6.1 工作 6.2 展望整理ppt6目前研究河流泥沙及河床演变问题，主要有三种手段，即原型资料分析法、物理

5、模型方法、数学模型模拟方法。这三种方法各有千秋，视具体工程问题，可独立运用或联合运用。物理模型方法是将原型的尺寸以及各种水力要素在遵循相似原则的基础上按一定比尺缩小为模型。物理模型方法具有直观、准确、易于操作等特点，是早期研究河床演变的主要方法。近年来随着测量技术的进步，物理模型的精度得到了很大的提高，但是与数学模型相比，物理模型耗资大，对外界条件要求高，且完成实验方案的周期长。河流泥沙数学模型能够模拟水流泥沙运动过程，河床冲淤量、河床形态及河床泥沙组成的变化。它的研究和应用可以追溯到上世纪50年代，由于受到当时计算工具的限制，这种计算仅限于一些简单的状况，而且精度较差、工作量大，因此，未能得

6、到普遍推广。二十世纪70年代以来，随着电子计算机的普及与计算机性能的不断提高，这些困难已逐步克服，数学模型的优越性越来越突出。经过近二、三十年的发展，水流泥沙数学模型经历了由一维到二维、到二 第一章第一章 绪论绪论1.1 泥沙数学模型发展概括泥沙数学模型发展概括整理ppt7 维嵌套、到三维、到三维嵌套，由原先在概化水文、泥沙及河床条件下数学模型的建立、率定和验证，发展到与物理模型相互配合共同回答工程水流和泥沙问题(复合模型)的阶段。1.2 水沙数学模型 1.水流控制方程 2. 泥沙控制方程0yxqqZtxy222227 322()xxyyxxxxxxtqqqqqqqqqqZghgnthxhyx

7、hxy222227322()yxyyyyyyyxtq qqqqq qqqqZghgnth xh yyhxy整理ppt8*2222()()()()()()kxkykskskskkkksxsyhSq Sq SSStxyhShSxy *1()Nbs ks kkkkZSSt整理ppt91.3 数值计算方法概述1.3.1 数值离散方法描述河流水沙运动的控制方程一般是非线性的，对于这些控制方程，一般很难求得解析解，只能通过离散控制方程、求解代数方程组来得到近似数值离散解。数值计算方法的精度及速度依赖于控制方程的离散方法、代数方程组的求解方法、网格的划分及边界条件的处理等。目前常用的数值计算方法有：特征线法

8、、有限差分法、有限元法 、有限体积法、 有限分析法等。 1.特征线法特征线法是河流模拟数值计算在计算机尚未普遍之时所采用的一种方法，主要采用手工计算，其最初思路是在平面上绘制特征线。其优点是能反映问题中信息沿特征传播的性质，算法符合水流运动的物理机制，稳定性好、计算精度高，但对于求解周期短、变化急剧的问题计算效果较差，因此目前很少直接用于数值计算。但是其原理仍很重要，经常用于作为了解其他数值方法的基础。2. 有限差分方法有限差分方法(FDM)是计算机数值模拟最早采用的方法，至今仍被广泛运用，FDM以Taylor级数展开等方法，把控制方程中的微商用差商代替进行离散，从而建立代数方程组来求解该方法

9、数学概念直观、表达简单，其解的存在性、收敛性和稳定性早已有较完善的研究成果，是比较成熟的数值方法，目前应用广泛。 整理ppt10 3.有限元方法有限元法是根据极值原理(变分或加权余量法)，将问题的控制微分方程化为控制所有单元的有限元方程，把总体的极值作为各单元极值之和，即将局部单元总体合成，形成嵌入了指定边界条件的代数方程组，求解该方程组就得到各节点上待求的函数值。有限元的求解方法常见的有直接法、变分法、加权余量法和能量平衡法等。其优点是网格划分灵活、拟合复杂河岸边界容易、网格节点可局部加密、稳定性好、精度高,适合于几何、物理条件复杂的问题。对于隐式FEM，其精度较高，但数学推导繁杂、计算量和

10、储存量较大，而且在误差估计、收敛性和稳定性等方面的理论研究与有限差分法相比还显得不够成熟和完善。尤其在多维计算中，由于有限元法贮存量大、会直接影响计算速度。整理ppt11 4.有限体积法有限体积法是计算出通过每个控制体边界沿法向输入(出)的流量和动量通量后，对每个控制体分别进行水量和动量平衡计算，便得到计算时段末各控制体平均水深和流速。因此，有限体积法正是对于推导原始微分方程所用控制体的回归。因为控制体间界面输运的通量，对相邻控制体来说大小相等、方向相反，故对整个计算域而言，沿所有内部边界的通量相互抵消。对由一个或多个控制体组成的任意区域，以至整个计算区域，都严格满足物理守恒律，不存在守恒误差

11、，并且能正确计算间断。由于采用守恒型的微分方程并对每一计算单元进行质量和动量守恒形式的离散，使得微分方程包含的守恒性质在每一个控制容积上都得到满足，若保持各单元界面两侧相邻控制体的计算输运通量相等，那么整个计算区域上都能保持守恒。整理ppt12在数值模拟的过程中，虽然采用的离散求解方法不同，但都有相同的特点，即首先把计算区域划分成许多控制体或网格，然后在这些小控制体上把微分方程离散成代数方程，再把小控制体上的代数方程汇合成总体代数方程组，最后在一定的初始条件下求解此方程组，从而求得计算区域内各节点的物理量。所以数值模拟的正确性和精度取决于网格的划分、方程的离散、初始值条件、代数方程组的求解以及

12、所建模型的物理理论依据是否正确合理等几个因素。各种方法均有其自身的优点和适应性，在实际计算时选择什么数值方法应根据所研究问题的特点和计算精度要求，以及研究者的习惯而定。整理ppt131.3.2边界处理技术1.不规则边界的处理在水利工程中，水流的边界往往是很复杂的。而在数学模型计算中如何处理不规则的复杂边界将直接影响到整个模型计算结果的可靠性、精度和计算时间，是数值模拟的一大难点，选择复杂边界的处理方法与算法的选择有着同等重要性。常用的不规则复杂边界处理方法主要有：美国学者Thompson16提出的贴体坐标法、非结构性网格以及多重网格法。贴体坐标法因其能使计算域边界上的坐标线与边界线密切贴合，而

13、获得了广泛的研究与应用。其主要优点是用曲线坐标逼近复杂边界，在边界模拟上具有较高的精度。但该方法在较复杂的边界上不易形成正交贴体网格，如果是非正交曲线坐标系，则会产生大量的交叉导数项，给求解带来很大的困难。非结构性网格使用三角形和四边形计算单元，使得在复杂区域形成计算网格更加灵活和方便。但非结构性网格往往要建立一套标明网格节点之间关系的复杂系统，因此在编程计算处理时比较麻烦，特别是针对自由面流动和干湿边界变化时，这种困难更加突出。多重网格法则由于网格系统之间的相对独立性，难以做到对系统守恒性的保证。另外对于复杂边界的处理问题，俞国良于1987年提出了镜像法，周建军于1988年提出了边界对称点法

14、。这两种方法都是根据对称原理，假设流场外有一定的虚拟流动，从而导出计算域外节点各物理量的计算式。该方法简单易行，提高了不规则边界上运用差分方法的计算精度，但由于域内外对称点一般不在节点上，仍需一定插值过程，会带来一定误差。整理ppt14 2.动边界的处理由于河床冲淤变化和水位的波动，在河岸及洲滩等区域都会遇到动边界问题。动边界数值模拟的困难主要在于沿动边界法向流动不同于明渠均匀流，其非恒定和非均匀性强，常用的曼宁摩阻公式等在形式上难以套用，再者是因为水深很小，对离散格式的要求很高，要求数值解不产生假振，保证水深总是大于零。常见的处理动边界方法有两种：一种是对实际边界的数值模拟，即时刻追踪动边界

15、的准确位置。另一种方法是：把整个计算区域都按照常进行计算，在计算过程中取一定的方法对无水区域进行处理。但是由于第一种方法在编制程序中处理较复杂，因此目前经常采用的方法是第二种方法，在计算过程中不追踪动边界的准确位置，把整个计算域都参与运算，对于干河床部分采用一些处理技巧进行处理，比如何少苓提出的“窄缝法”，其思想是在床面上人为设置窄缝，由于缝足够深，总是有水流动，而缝很窄，当河床露出水面时由窄缝引起的水量守恒误差很小。此外，还有周建军提出的渗透边界法，程文辉、王船海等运用的“冻结”技巧。整理ppt151.4水沙基本理论中存在的问题1.阻力问题 阻力计算是阻力系数(糙率及涡流粘性系数)的确定问题

16、。处理阻力的方法有两种，一种是采用水力半径分割法或能坡分割法划分阻力单元，分别计算各单元的阻力，然后计算总阻力；另一种方法是根据实测资料直接计算总阻力。以上两种方法都只能用于一维水流糙率的计算。多数的二维数学模型均是直接采用一维阻力系数进行计算，这种处理方式显然比较粗略。李义天在考察一维情况下断面综合糙率沿河宽变化分布的情况下，提出了二维糙率计算公式。实际应用表明，该方法比直接采用一维方法确定的糙率更合理，能比较好地反映实际情况，但由于因素比较复杂，由此得到的推论在定量上尚未能被实测资料所证实。2.水流挟沙力问题水流挟沙力的计算一般有两种方法，一种是将含沙量沿垂线分布公式积分求垂线平均含沙量，

17、但河底泥沙含沙量很难确定，不同的研究学者给出了不同的计算方法，所得的计算结果与实际情况差别较大。另一种是半经验半理论方法，这样的公式很多，公式中包含两类参数，一类是反映水力条件的参数，如流速水深等；另一类是反映泥沙特性的参数，如沉速、泥沙粒径等。这种方法是建立在实测资料的基础上得，因而结果比较接近实际情况，得到了广泛的应用。现在挟沙力的研究成果多基于一维问题，对二维问题研究的较少，而大多数的二维问题也是直接采用一维挟沙力公式进行计算的，但是李义天通过经验比较，发现二维挟沙力与一维挟沙力有很大的差异。整理ppt16 3.床沙悬移质沿程变化的问题目前在数学模型计算中存在的尚未很好地解决的问题是床沙

18、质及床沙级配的沿程变化问题。对于冲淤甚小的处于自然状态下的河流，在许多情况下床沙及其级配在冲淤过程中可以假定基本不变，而冲泄质级配则视来水来沙条件而定；对于冲淤幅度甚大的河流，对床沙质及床沙级配的沿程变化不能不进行考虑。现在大多数模型将推移质和悬移质分开，将悬移质中床沙质与冲泄质分开的做法基本相同，至于将床沙质进一步分级来求不同粒径的冲淤变化，一般都比较粗糙。韩其为23提出了一套考虑泥沙级配变化的计算模式，其基本思想是不区分床沙质与冲泄质，假定水流挟沙力的级配和实际输移的级配一致，而悬移质和床沙质级配在每一时段内的变化都看成本时段内河段冲淤的直接后果，计算断面每个粒径组的含沙量须反复试算才能确

19、定，该方法对泥沙输移的物理考虑比较细致但计算过于复杂。整理ppt171.5淤积上延现象当水流进入到水库回水区后，由于水流条件的改变，流速减小、泥沙落淤，而落淤的结果反过来又促成水流条件的进一步改变，使得回水曲线在淤积区及其上下游一定范围内普遍抬高，为泥沙继续在原淤积区落淤并同时向上下游发展提供条件。这种淤积自回水末端(通常指正常蓄水位回水末端)向上游发展的现象，称为水库淤积的上延现象，也就是通常所说的“翘尾巴”现象。按淤积中剖面形态的类型，可以分为三角洲淤积、带状淤积、和锥体淤积见。整理ppt1805101520253035010203012345相对高程(米)距 坝 里 程 (公 里 )1-

20、尾部段；2-顶坡段；3-前坡段；4-异重流淤积段；5-坝前淤积段0204060801001200102030405060相对高程(米)距坝里程(公里)1-原平均河底高程；2-淤积后平均河底高程12图 1-1 官厅水库淤积纵剖面图 1-2 丰满水库淤积纵剖面整理ppt192130740750760高程(米)距坝里程(公里)211-原河床；2-淤积后滩面图1-3 陕西黑松林水库淤积纵剖面整理ppt20 前的3.74变缓为1.8，淤积末端距大坝13.7km，上延系数1.55。随着淤积末端不断上延，上游河床逐年抬高，距大坝9.6km的四合永公路桥处河床淤高3m多，使河床由原来的地面以下到现在与两岸持平

21、，甚至高出两岸耕地及村基，产生严重浸没。浸没区内数千亩耕地盐渍化、沼泽化，不能正常耕种，近千户居民住房地基沉陷、墙垣断裂甚至倒塌，给群众生产带来严重损失，以至出现影响社会稳定的问题，并严重威胁四合永镇的安全，受到省、部领导和有关部门的高度重视。 为延缓庙宫水库淤积浸没的发展，自1987年至1992年连续6年进行了空库排沙原型实验。试验中对不同水沙条件，空库排沙时间，不同排沙方式等进行了全面的试验，总结分析了水库现有泄流排沙设施条件下的减缓水库淤积浸没的可能性和运用效果。1993年由水利部天津院完成的庙宫水库空库排沙试验物理模型研究和同年由天津大学泥沙研究室完成的庙宫水库排沙洞前期工程泥沙数学模

22、型初步研究，为当时增设新排沙洞位置、洞径和底高程的确定提供了理论依据。 近年来随着水环境关注意识的提高和水资源高效利用思路的确立，如何在水库保水和水库排沙间选择优化的水库运行模式是需要进行研究和解决的重要问题。根据庙宫水库设计正常蓄水位778m(大沽高程，以下不加说明处均采用大沽高程)和库区死水位768m的要求，优选即能泄水排沙又能保水兴利的调水调沙方案是本次研究的主要目的。研究过程拟采用库区限蓄水位768m、770m、772m、774m、776m、778m六级水位标准，来水来沙过程采用19581959、19711975、19851989和1990的四个系列三种组合来模拟不同年份水沙过程，在1

23、8个计算系列中进行方案优选。整理ppt21第二章 河道复杂边界的处理技术自然界中，水流的边界往往是十分复杂的，而在水流数学模型计算中如何处理不规则的复杂边界将直接影响到整个模型计算结果的可靠性、精度和计算时间，是数值模拟的一大难点。复杂边界的处理方法与算法的选择有着同等的重要性，因此在研究江河湖泊的水流数值计算过程中，必须要对边界作出一个适当的处理方法。2.1 斜对角笛卡尔方法的基本原理斜对角笛卡尔方法基本原则就是按就近逼近原则，把距实际水边线较近的网格节点连线以拟合实际水边线。生成网格时，根据边界是否固定，可采取不同的判定方法。当边界为固定边界时，可把固定边界离散为若干离散节点，确定各离散节

24、点坐标，按距离最小原则选定固定边界上各离散节点的逼近网格节点，然后依序连接各网格节点，即形成拟合边界；对于动边界，由于边界线变化不定，则可借助网格节点间水深的干湿变化和水深值，参照相应网格节点间的地形高差，即可判定距离实际水边线较近的网格节点，合并纵横向的拟合网格节点，即可形成拟合瞬时水边线。整理ppt22图2-1 斜对角笛卡尔法边界拟合示例整理ppt232.2斜对角笛卡尔方法的优点斜对角笛卡尔方法的优点 1.边界拟合效果良好，更接近实际边界。由图2-2可以看出，计算水边线与实际水边线拟合良好，其边界长度误差一般小于6%，仅为锯齿法长度误差的，且随网格间距的缩小而更为精细。图2-3为按锯齿法(

25、也以距离最近为原则)确定的计算水边线，对比图2-2、图2-3可知斜对角笛卡尔法所拟合的水边线比锯齿法拟合的水边线更为贴近实际水边线。事实上，斜对角法就是在锯齿法基础上通过斜对角连线添加一个三角形网格单元，使边界形态更接近真实形态，且避免了锯齿状边界的产生，使流场计算更精确。2.处理简单，附加计算量小在矩形网格水域边界上运用斜对角笛卡尔法，计算时，只须对三角形边界网格单元进行处理，而三角形边界网格单元若把斜边中点视为一角点，则仍可视为四边形网格，可用统一方程求解，与内部计算网格无异，并不像贴体坐标拟合一样须对基本方程进行坐标变换而增添许多交叉导数项，因此，用斜对角笛卡尔法进行边界处理，不仅保留了

26、离散简便、计算效率高的优点，同时又使边界拟合更加精细，避免了锯齿状边界对流场计算的影响。整理ppt24实 际 边 界 线斜 对 角 法 拟 合 线实 际 边 界 线锯 齿 法 拟 合 线图2-2斜对角笛卡尔法边界拟合 图2-3锯齿法边界拟合整理ppt252.3边界条件及边界网格单元类型处理根据斜对角笛卡尔法的逼近原则，在平面二维计算中拟合水边线上形成的计算边界网格单元存在八种基本类型，见下图。 图2-4八种网格单元基本类型整理ppt26图2-5 网格类型及差分点分布整理ppt27斜对角笛卡尔方法不仅可以应用于有限差分法中，也可应用于有限体积法等其它离散格式中，它对复杂边界条件和复杂水流条件的适

27、应性和计算精确性将使复杂边界的处理变得简单高效而且准确，为复杂边界问题和近岸工程问题的计算提供了一种新思路。整理ppt28第三章第三章 平面二维水流泥沙数学模型的建立平面二维水流泥沙数学模型的建立天然河道水沙数学运动一般都属于三维运动，运动要素既沿程变化，又沿水深和河宽发生变化28。由于三维水流运动比较复杂，难于计算，常将运动要素沿水深方向平均，把三维问题转化为平面二维问题。平面二维非恒定流方程，也就是浅水方程，它描述的对象是符合以下条件的均匀流体：(1) 具有自由表面。(2) 以重力为主要的作用力，以水流与固定边界之间及水流内部的摩阻力为主要耗散力，有时还存在水面气压场风应力以及地球柯氏力等

28、的作用。(3) 水平流速沿垂线近似均匀分布，不考虑实际存在的对数或指数等形式的垂线流速分布(4) 水平运动尺度远大于垂直运动尺度，垂向流速及垂向加速度可以忽略，从而水压强接近于静水压强分布。整理ppt29在实际问题中，通常可以作为非恒定流处理的具有自由表面的实际水流，有如下情况：(1)水深相对较浅。可以根据水深h和波长L之比区别，当0.4时称为浅水长波。(2) 水底坡度较缓。设底坡倾角为，判断缓坡的条件是。此时，可以忽略底坡引起的垂直速度和垂直环流，也不必考虑垂直加速度及由此产生的动水压力。(3) 水面渐变且坡度较缓。既使满足以上两个条件，在特定的情况下，由于受到天然或人工控制，有时水面比降仍

29、很大，水深沿程变化很大的急变流，水压力不能用静水压强分布来近似。(4) 无明显垂直环流。由于外力改变或地质环境的改变在垂直平面内产生二次流，会给浅水模型带来误差，需要加以改正。此时方程不能沿垂向积分，要按不可压流来处理。密度沿垂向变化较大，改变了静压分布，需要作为密度流处理。因此，通常的浅水流动指在重力作用下密度均匀，具有自由表面，流动近似水平的长波传播现象。 相比一维水沙数学模型，平面二维水流泥沙数学模型对沿水深平均的水沙各运动要素(垂线平均流速、悬移质含沙浓度、推移质输沙率)以及河床变形等随流程沿河宽两个方向的变化情况均可以进行比较细致的模拟。整理ppt301.水流控制方程3.1二维水沙数

30、学模型的基本方程0yxqqZtxy222227 322()xxyyxxxxxxtqqqqqqqqqqZghgnthxhyxhxy222227 322()yxyyyyyyyxtqqqqqqqqqqZghgnthxhyyhxy整理ppt312.泥沙控制方程*2222()()()()()()kxkykskskskkkksxsyhSq Sq SSStxyhShSxy *1()NbsskskkkkZSSt 整理ppt323.关于悬沙连续方程的一些探讨泥沙在水流中分布就是悬沙的含沙量，而分布的变化即为悬沙含沙量的变化，它是可以用泥沙连续方程来描述的。悬沙含沙量的变化不外乎是由对流项(上游输沙)，和本地的紊

31、动扩散项引起的，由此可以把对流项，紊动项引起的悬沙含沙量变化用两个连续方程来表示。主要思想如下：由三维泥沙扩散方程：其中：是单元体内的含沙量，可以这样认为单元体内含沙量是由上游来沙和本地床面紊动扬沙引起的。我们把上游来沙影响用来表示，称为输沙含沙量。本地床面紊动扬沙用来表示，称为扬沙含沙量。这样，可以把单元体总的含沙量用如下式子表示： ()()()xyzSSSSSSSuvtxyzxxyyzz SSS 整理ppt33分别根据时均，脉动的概念，可以得到两个方程： 可见，输沙含沙量和扬沙含沙量都是满足三维泥沙扩散方程的。三维方程在物理意义描述上是十分准确的，但是在实际应用中，因为计算条件在河道发生冲

32、淤变化时，悬移质底层及床面扬沙含沙量的变化是引起床面冲淤的主要因素，而输沙含沙量只是上游来沙在本地的一种扩散效应，在当时刻冲淤变化过程中并不会起到举足轻重的作用。我们把含沙量区分为输沙含沙量与扬沙含沙量也是为了在计算冲淤变化过程中体现这种关系，以使结果更符合实际情况。()()()xyzSSSSSSSuvtxyzxxyyzz()()()xyzSSSSSSSuvtxyzxxyyzz整理ppt34对输沙含沙量方程和扬沙含沙量二维化()()xySSShShShuhvhtxyxxyy *()()xySS SSSShShShhuhvtxyxxyy 整理ppt35S水 面Z0Z0+haSaSbSAB输 沙

33、含 沙 量扬 沙 含 沙 量整理ppt36整理ppt37如上所述，悬移质中造床质所占比例是要受条件限制的，同样的，床面泥沙能否起动，以及起动后是以何种形态运动，这在模型计算中，也是要予以考虑的。 Shield(1936)对各种泥沙颗粒进行了临界起动实验，实测得到了无量纲临界起动剪切应力与颗粒雷诺数的关系曲线， 整理ppt38泥沙颗粒的悬浮指标是通过悬浮表达式来计算， *Zu (3-111) 式中：k是卡门系数，是泥沙颗粒的沉速，*u是摩阻流速。 根据Einstein的观点，对于给定的水流和泥沙条件，可取5z 作为泥沙是否进入悬浮状态的临界判断值。 整理ppt39整理ppt40222222222

34、222220SSSSSSuvuvuvtxyx ty tx y (3-117) 设式子(3-108)的解为1sin()kiiiiSAatb xc y，代入式子(3-108)，每一项均可以写为(为简单起见去掉下标i)： 22222sin()(222)0Aatbxcy ab uc vabuacvbcuv (3-118) 如果式(3-108)解存在，则上式(3-109)要满足如下条件： 222222220ab uc vabuacvbcuv (3-119) 即 2()0abucv (3-120) 则 abucv (3-121) 其中解sin()SAatbxcy中，, ,A b c为待定系数，由初始条件和

35、边界条件确定，为定性给出解S的变化，给定特定参数vucbaA,，可以画出床面定性变化趋势。 整理ppt41式(3-108)是椭圆方程，通解肯定是正弦余弦函数的线性组合，为了定性显示解得形式，并取0t 时刻，取式(3-113)、(3-114)为S的解。 2sin(23 )2cos(2 )Sxyxy (3-122) 2sin(75 )4cos(46 )Sxyxy (3-123) 图(3-2)、图(3-3)分别显示式(3-122)、式(3-123)值的变化以及等值线的分布。为清楚起见，图(3-2)中Z放大比例为20，, x y空间步长为0.2，, x y取值范围-10到+10，图(3-3)中Z比例为

36、10，, x y都以步长0.05变化，, x y取值范围-5到+5。 整理ppt42-10 -8-6-4-20246810-10-8-6-4-20246810(a)(b)空间分布等值线分布图 3-2 式(3-122)值的变化及等值线的分布 -5-4-3-2-1012345-5-4-3-2-1012345(a)(b)空间分布等值线分布图 3-3 式(3-123)值的变化及等值线的分布 整理ppt43结果表示具有波动性，表明床面有可能出现沙纹或沙波的变化形式，同时床面的变化形式与扬沙含沙量的变化相关。图中表示的床面比较规则，但对复杂的水动力条件和床面初始形态，床面变化具有拟周期性，同时变化程度非常

37、复杂不规则。由于泥沙问题的复杂性，目前对泥沙运动机理还没有完全认识清楚。现有的各家对基本问题的认识及处理方法也不尽相同，甚至相差甚远。在20世纪70年代开始出现了大量根据实测资料，利用计算机进行回归方法计算得出的阻力和输沙力经验公式，由于有大量的实验及天然资料做基础，各种变量的覆盖范围大，这样得出的一些公式精度很高，甚至超出了许多理论公式，其原因在于现有的理论体系往往不能够把所有的泥沙的输移过程完整地、正确地反映出来。本节内容从泥沙扩散的物理角度出发，基于现有的泥沙理论推导出两个可分别计算输沙和扬沙的公式，在理论扩展、公式推导上做了一新的改进，从公式结构上看与现有悬移质泥沙连续方程结构相同，如

38、将输沙连续方程与扬沙连续方程相叠加得到的全沙连续方程与现有悬移质泥沙连续方程具有完全相同的形式，说明方程推导正确，只是在划分输移泥沙和扬动泥沙标准上可能还存在进一步研究的空间。整理ppt443.2 初值条件及边界条件的处理 初值条件可以由初始时刻的实测资料给出， 但所需条件在实际计算中一般难以全部获得，只能通过估算加以补足，当然估算得越接近实际越好。 对计算中所用的初始水流条件，一般以下游水位为起点，根据平均河床比降给出各点水位的初值，并令初始流速u、v均为0。含沙量初值有实测值或取为042。 通常对水流而言，进口边界给出的流量过程线，出口边界给定水位流量过程线。对悬移质而言，进口一般给定含沙

39、量过程及相应泥沙级配，出口认为/0kSx 。 对于岸边界条件，一般采用水流无滑移条件，即取u、v均为0。对悬移质而言，岸边界处取/0kSn。(n为边界法线方向) 上游进口边界： 000( )( )( )uu tvv tSS t (3-124) 下游出口边界： 00( )( )/0kuu tvv tSx (3-125) 整理ppt453.3网格划分网格的确定对模拟实际流动的可靠性同样重要。建网基本要求是：符合流动特点；易于建立；比较光滑和规则；满足精度和计算的稳定要求；便于组成节约、高效的数据结构，必要时可随时依解的梯度作适应性调整。本文为了更好的模拟计算域内的实际地形、地貌，采用的是无结构有限

40、体积网格划分计算区域，使模型地形、地物与实际情况基本一致。一般而言，无结构网格划分采用三角形网格，因其有较为成熟的自动生成技术，故自动化程度较高，但是针对本文设计的二维河道情况，采用的是不规则四边形组成的网格，主要考虑到：一方面若完全采用三角形网格作均匀划分，反而会难以客观反映小水量时水流沿主槽流动的特点，有可能造成计算结果失真；另一方面，网格呈四边形格式，局部辅以三角形网格，这样不但可以很好地拟合边界形状，而且格子数及其边数，要比完全采用三角形网格时少，同时实际空间步长也可以取大些。整理ppt46按照有限体积法布置二维网格的方式如图 3-4 所示。它取单元网格为控制体，在网格中心处计算水位H

41、，在网格周边通道的中点处计算流量Q。即在平衡计算时，沿控制体每一边的法向通量用该边中点处的通量作代表，乘以边长即为通量沿该边的积分。中点的通量可用中心格式(如取相邻两格子形心处通量的平均)或逆风格式确定。 计算所采用的方程改写成矢量形式， 按照有限体积法，将其在控制体内进行积分，对水位H和流量Q按时间交错方式计算(如图 3-5)。 节 点通 道网 格通 道网 格T-2dtT+2dtT+dtT-dtT 图 3-4 H和Q的空间布置方式 图 3-5 H和Q的时间交错计算方式 整理ppt473.4 水流挟沙能力1.悬移质输沙率公式的选取本文中采用的积分求输沙率的方法，选用了爱因斯坦挟沙力理论。*12

42、30.211.62.303lg()sSvashgu aSIIk (3-131) 其中： 11110.216(1)ZZZAAIdA (3-132) 11210.216 ln(1)ZZZAAIdA (3-133) 整理ppt481I和2I均是A及Z的积分函数， 可通过数值积分求解预先做成图标， 如图3-7所示。 如果没有床面形态时， 上述公式中的沙粒阻力对应的剪切流速*u就等于*u。 爱因斯坦假定推移质的运动速度为*11.6u，推移质厚度为二倍粒径。若求解悬移质输沙率时把沿垂线积分的下限取为2aD，记该处泥沙的体积比参考浓度vaS，则 */(11.62)vabsSguD (3-134) 整理ppt

43、492.泥沙计算辅助公式 (1)Manning公式48 21321uR Jn (3-135) 式中：n为Manning系数，也称为糙率，R是水力半径，J是水力坡降。 (2)Manning-Strickler公式 1 66524()uRDggR J (3-136) 式中：u是断面平均流速，R是沙粒阻力对应的水力半径，J同上，65D是泥沙级配中65%对应的泥沙粒径。 (3)泥沙颗粒统一沉速公式 29(9)sgDDD (3-137) 式中：v是水流的运动粘性，D是泥沙颗粒粒径。 整理ppt50整理ppt51(7)起动临界流速 参考沙莫夫公式49 1/61.47()chUgDD (3-141) (8)

44、止动临界流速 参考沙玉清公式 3/40.011 0.39sUD (3-142) 1 21 21()()bssggD 式中：为推移质强度。 整理ppt52第四章第四章 平面二维水流泥沙数学模型的数值模拟平面二维水流泥沙数学模型的数值模拟本文采用有限体积法离散方程， 它将计算区域划分为网格，并使每个网格点周围有一个互不重复的控制体；将待解微分方程对每一个控制体积分，从而得到一组离散方程，求解差分方程组，其求解就作为微分方程定解问题的近似解。 4.1 水流控制方程的离散2122kTdtTTdtiiikikinidthhQ LdtqA (4-9) 式中，iA为第i个网格单元面积，ikL为i号网格的第K

45、号通道长度；ikQ为i号网格的第k号通道的单宽流量。 整理ppt53动量方程的求解，采用分类简化处理的方法，把它概括为地面型通道、河网型通道和特殊通道来分类处理。假定同一时段内同一网格水位变化不大，同时由于平面浅水中对流项作用较小，在计算中略去不计，此外，有根据堰流流量公式计算的连续堤或缺口通道。整理ppt541.河道型通道，即通道两侧网格均为河道型网格，动量方程中保留局地加速度项、重力项和阻力项，利用差分方法离散得到河道型通道的动量离散方程整理ppt55 2.地面型通道，即通道两侧单元为陆地地面，且通道上没有堤防等阻水建筑物。考虑到库区内的地形变化，地面洪水演进主要受到重力和阻力的作用，利用

46、差分方法离散得到地面型通道的动量离散方程 3.对于高于地面的阻水建筑物，如堤防、铁路、公路等，可以概化成连续堤或缺口堤通道，其流量采用宽顶堰溢流公式来计算，离散后得到 泥沙连续方程离散方法与水流连续方程离散方法相似，此处不再赘述。整理ppt564.2 泥沙控制方程的离散 泥沙连续方程离散方法与水流连续方程离散方法相似，此处不再赘述。整理ppt574.3 河床纵向变形方程的离散 对式子(3-104)的离散如下 22*TdtTbibidtSZZ (4-32) 式中：bZ是河底高程。 在计算河床变形要注意以下情况 SSSSSSSSSSStZb|, 0 |, 0 , 0 0, 0 且如果且如果但不满足

47、泥沙起动条件如果且满足泥沙起动条件如果 (4-33) 整理ppt584.4 模型求解过程 模型求解的过程具体步骤如下： (1)给定计算所需要的基本参数，比如, ,xyn 等，代入初始条件和边界条件； (2)对整个计算区域赋初值； (3)在时间步长之内，求解速度, u v，水深h； (4)根据 Einstein 悬移质理论计算vaS；应用输沙含沙量方程、扬沙含沙量分别计算,S S、此后应用悬浮指标z及王尚毅的悬浮功理论判断悬沙、底沙的交换比例，最后计算河床纵向变形bZ。 (5)至此一个时间步长求解完毕，所求得的结果作为新的初始值，进入下一个步长的计算过程，重复第(3)(5)步骤，直至时间步长循环

48、到最终结果。 (6)输出所需的计算结果。 整理ppt59第五章第五章 水流泥沙数学模型验证及应用水流泥沙数学模型验证及应用本模型的建立主要用于庙宫水利枢纽库区泥沙变化的研究，选择了库区河道为模拟对象，在分析库区河道水流泥沙运动特点的基础上，建立了二维水流泥沙数学模型，并用实测水文资料对模型进行验证分析，在此基础之上，进行了水沙组合的计算。整理ppt60图1-1庙宫水库位置示意图滦 平张 百 弯波 罗 诺风 山滦河白 虎 沟成 子燕 上道 坝 子大 唤 起棋 盘 山龙 头 山围 场蚁蚂吐河隆 化承 德庙 宫 水 库钓 鱼 台 水 库武烈河三 家荒 地大 庙伊 逊 河 流 域 平 面 位 置 图整

49、理ppt61庙宫水库是河北省大型水库中淤积浸没最严重、问题最突出、社会影响最大的水库， 在海河流域乃至全国大中型水库中也是淤积浸没最严重的水库之一。该水库位于滦河的多沙支流伊逊河上，由于建库设计时资料短缺，上游来沙量计算偏小，对泥沙淤积的严重性估计不足，因此在水库泄水建筑物的布置和运用方式的拟定上，没有充分考虑排沙的要求，致使大量泥沙淤积在库区。水库淤积形态为锥体淤积，坝前淤积高程 768.5m，超过死水位 0.5m。库区河道比降由建库前的 3.74变缓为 1.8，淤积末端距大坝 13.7km，上延系数 1.55。随着淤积末端不断上延，上游河床逐年抬高，距大坝 9.6km 的四合永公路桥处河床

50、淤高 3m 多，使河床由原来的地面以下到现在与两岸持平，甚至高出两岸耕地及村基，产生严重浸没。 整理ppt625.1模型基本特征年际间输沙量变化悬殊，庙宫水库入库最大沙量1360万吨(1959年)，最小沙量61万吨(1981年)，两者相差22.3倍。蚁蚂吐河下河南站1959年输沙量为1080万吨，1966年输沙量仅32万吨，变幅达33.8倍。水沙特点表现为水沙过程的不完全对应一致，有时沙峰先于洪峰，有时反之；有时大水大沙、小水小沙、也有时大水小沙、小水大沙，主要是由于暴雨强度、过程和中心位置的差异以及各支流水土流失程度不同等原因造成的。整理ppt635.2网格划分庙宫水库数学模型自水库大坝向上

51、游延伸，在四合永镇北河道岔口分为两支，一支沿伊逊河向西北至掌字水文站，另一支沿不澄河向东北至边墙山水文站。模型主体南北走向，全长约12.6km，模型最宽河段约1km。模型内包含27个实测断面，公路桥、铁路桥各一座。模型网格划分以实测断面为基础，采用无结构有限体积方法划分网格，考虑主河槽走向，划分后得到断面194个，单元3451个，结点3626个，通道7076个。断面平均间距为65m，断面内空间步长为50m，单元平均面积为2124m2，模型覆盖水域面积约为8km2，模型内实测断面分布见图5-3，模型划分网格见图5-4，网格单元编号、通道编号、节点号的分布形式见图5-5。整理ppt64整理ppt6

52、5整理ppt66 123456789101112131415161718192021222324252627282930313233234567891011121314151617181920212223242526272829303132333435363738123456789101112131415单 元 号通 道 号节 点 号：红 色黑 色紫 色整理ppt67 5.3模型验证 1.水流模型验证 验证水流模型上游边界采用伊逊河掌字水文站和不澄河边墙山水文站1990年7月14日至15日历时6小时流量过程（见图5-5），下游边界采用庙宫水库水位流量关系过程（见表5-1），验证计算中没有考虑排

53、沙洞泄洪影响。假定大沽基面与国家高程基准面之间的换算关系为：假定大沽基面-16.8681985国家高程基准。整理ppt68整理ppt69整理ppt70水流数学模型验证伊逊河掌水文字站和不澄河边墙山水文站1990年7月14日至15日历时6小时水位过程，验证结果见图5-8，验证计算流场见图5-9。整理ppt71整理ppt72从水位验证过程来看，计算值与实测值变化趋势一致，数值偏差不大，基本能保证模型预报计算精度。从不同时刻流场上看，当上游边界产生大流量时，公路桥附近收缩断面出现较大流速，下游大坝附近还没有达到溢流水位；当上游边界流量回落时，公路桥附近断面流速减弱，上游流量已传递到大坝附近，达到溢流

54、水位，在坝前和库区弯道附近产生较大流速，基本符合水流运动规律，说明水流计算结果合理，水流模型可用于模型预报计算。从水位验证过程来看，计算值与实测值变化趋势一致，数值偏差不大，基本能保证模型预报计算精度。从不同时刻流场上看，当上游边界产生大流量时，公路桥附近收缩断面出现较大流速，下游大坝附近还没有达到溢流水位；当上游边界流量回落时，公路桥附近断面流速减弱，上游流量已传递到大坝附近，达到溢流水位，在坝前和库区弯道附近产生较大流速，基本符合水流运动规律，说明水流计算结果合理，水流模型可用于模型预报计算整理ppt732.泥沙模型验证泥沙数学模型验证采用计算模式验证的方法进行，验证1990年汛期前后平均

55、滩面的变化，汛期上游来流过程见图5-6，上游边界含沙量过程见图5-10验证结果见图5-11。整理ppt7420:0021:0022:0023:000:001:002:0090年7月14日7月15日时间（h）含沙量（kg/m ）305010015020020004000600080001000077077578078590年汛前平均床面90年汛后实测平均床面90年汛后计算平均床面距坝里程（m）高程（m）34 567891011 12 13 14 1516 17 18 19 20 212223实测断面序号图 5-10 汛期上游边界含沙量过程图5-11 1990年汛期前后水库平均滩面变化整理ppt7

56、55.4 模型预报近年来随着水环境关注意识的提高和水资源高效利用思路的确立，如何在水库保水和水库排沙间选择优化的水库运行模式是需要进行研究和解决的重要问题。根据庙宫水库设计正常蓄水位778m(大沽高程，以下不加说明处均采用大沽高程)和库区死水位768m的要求，优选即能泄水排沙又能保水兴利的调水调沙方案是本次应用的主要目的。研究过程拟采用库区限蓄水位768m、770m、772m、774m、776m、778m六级水位标准，来水来沙过程采用19581959、19711975、19851989和1990的四个系列三种组合来模拟不同年份水沙过程，在18个计算系列中进行方案优选。下游排沙控制水位分别选取7

57、68m、770m、772m、774m、776m和778m，与三种水文类型年组合系列再组合形成18种计算方案，方案命名以水文类型年组合序号为先导字符，用横线连接排沙控制水位值，计算方案见表6-2。整理ppt76整理ppt77整理ppt78整理ppt79整理ppt80整理ppt81整理ppt825.5排沙量和淤积量分析根据方案的计算结果，得到沿庙宫水库上游至下游各过流计算断面的床面平均变形量曲线，便于比较将同一水文系列的情况在同一图中给出(图5-33至图5-35)。整理ppt83从图5-33至图5-35中可以看出，三种不同类型的来水来沙条件所塑造的河床变化形态相差较小，变化趋势完全一致，只是下游不

58、同的排沙控制水位对河床冲淤趋势和冲淤量影响较大，故只给出水文条件下的最低排沙控制水位768以及最高排沙控制水位778情况下的库区含沙量、冲淤变化分布、流量沿程变化(见图5-36至图5-45)。其他情况与此类似，应当注意的是，输沙含沙量与挟沙能力都是正值，而扬沙含沙量、冲淤量都存在正负，对于扬沙含沙量正值代表冲刷，负值代表淤积；冲淤量正值代表淤积，负值代表冲刷。整理ppt84整理ppt85整理ppt86整理ppt87整理ppt88整理ppt89整理ppt90整理ppt91整理ppt92整理ppt93整理ppt94整理ppt95整理ppt96从图5-33至5-35可以看出，不同的排沙控制水位，对于离坝距离大于5200m(铁路桥附近)的地方，影响很小，而影响范围，主要集中在大坝到铁路桥这一区段，变化趋势是，随着控制水位的不断增大，排沙洞到大坝这一段淤积增加；排沙洞附近冲刷能力加强；排沙洞上游淤积增加。总之就是，排沙控制水位高，离排沙洞较远的地方淤积就会加重，而排沙洞附近的冲刷反而会提高。进一步比较各方案之间在整体河道上的排沙效果，在表5-6中给出了庙宫水库在三种不同类型的来水来沙条件和不同排沙控制水位下各计算方案比较，13年