（精密仪器及机械专业论文）系统仿真技术在黄河调水调沙控制中的研究与应用(精密仪器及机械专业优秀论文).pdf

哈尔滨工程大学硕士学位论文 摘要 系统仿真技术作为一项非常重要的技术，越来越多的引起了人们的重 视，它不仅可以节约资金，而且又可以节省项目从开始到实验完成的周期。 黄河作为目前世界上独一无二的高含沙河流，9 0 年代以来，慢慢进入枯 水期，水量在不断减少。因此，下游河道泥沙的淤积非常严重，因为流速缓 慢水位上涨，危险情况频频发生。为了改善下游河道，进行了三次调水调沙 实验，运输泥沙2 6 亿吨入大海，使黄河下游的河道流量能力增加6 0 0 万 m 3 s 。本论文围绕如何能够正确预报黄河下游河道的流量、水位、泥沙量、 淤积量，详细的进行了研究。 论文主要包含三方面的工作：首先在设计系统仿真平台的基础上，分析 了各个系统仿真平台的组成部分及其各个部分的功能，然后设计了仿真数据 库并研究数据库建立要考虑的内容，最后设计了仿真控制系统的结构框图： 确定了一维模型作为系统的仿真模型，并对模型进行了合理的简化，然后对 各个模型进行了离散和线性化，在算法上进行了分析与研究：最后，研究了 各个模型的水力参数，根据仿真模型，编制了程序，进行了仿真实验，因为 结果误差很大，所以本文分析了误差产生的原因。 关键词：系统仿真；黄河；一维水沙控制方程；p r e i s s m a n n 哈尔滨工程大学硕士学位论文 要约 系统傲真技衍c 土巴- e0 大切吞o 寸、多i 于扎浮多0 1 旺p 、人重祝 芒扎| l 亡、资金窑溅少卡5 茫i 于- c 、项目始曲女、弓终扣 t 周期苗 短缩l 寸。 黄河c i 今 t 、世界t 一。烂i ，高泥砂o t 寸、9 0 年代以降、徐l 二渴水 期【二入9 、水量c i 溅少l 【o 、，亡。圣。亡曲、下流。河道cc 土泥砂口) 沈殿 深刻化l 、流扎滞，亡：t ， ，岛、水位抹上_ 、危险各状沉：频繁c 二 生b 为上j 忙畚，亡。下流。河道力 改善寸为止曲c 二、3 回水调整砂调整j 麦颤啻行l 、针2 6 德p 二，o 泥砂七海 逢搬芒寸、黄河下流。河道流量能 力它6 0 0 万秒【二 t 向上芒世el 、右。：o 谕文娃下流河道忙流量、水位、 泥砂量、滞稹量啻正l 兄达扎5 允曲l a 、群l 研究l ，- eo 、 寸。 ：o 谕文；i 三? 力方面。仕事含曲与扎| l 盎。主矿、系统傲真弑骏台 童彀射b 亡基e 、褚试酸台。中t 部分褚弑黢台。劾能七分析b - ek 、b 亡。干lt 、s i m u l a t i o n d a t a b a s e 老彀卦l ，( 、圣。他。内容窑群l 研 究l p ；即使在含沙量较高的情况下， 也采用简化处理，只是在系数方面作些修正9 ”。用( 3 1 0 ) 式计算含沙量的模 型称为非饱和输沙模型，国内发展的模型基本属于此类；不用( 3 1 0 ) 式求解 含沙量，而认为含沙量近似等于挟沙力，这种模型称为饱和输沙模型，国外 早期发展的模型多属此类。实际应用中常常需要考虑泥沙颗粒级配不均匀的 影响，可以将泥沙按粒径分组，每一组内以一个平均粒径代表，当作是均匀 的，于是又可以建立分组泥沙连续方程和河床变形方程【3 4 ”5 1 ，其形式与前述 相同，但含沙量、挟沙力及沉速等意义均指某一粒径组。河床变形计算往往 历时较长，为了简化计算，通常将非恒定流过程概化为梯级式的恒定流过程 求解，即取罢= 0 ，罢：o 及娑：0 ；此外含沙量、水深随时间变化所引起 u |o io l 的水体中的沙量变化通常也忽略不计，即取皇竺：0 ；在岸线比较顺直，没 u t 有回流产生的条件下，水流运动方程中的紊动剪力项和泥沙连续方程中的泥 沙扩散项也可忽略不计。近几年，非恒定流模型研究发展较快【3 6 】【3 7 1 。 3 1 3 2 仿真算法 平面二维泥沙数学模型的计算比较复杂，计算量很大，所采用的数值方法 是否合理，直接影响到它的适用性。常用的方法有有限差分法、有限元法 ( f e m ) 、有限分析法( f a m ) 、边界元法( b e m ) 等。 3 1 4 剖面二维及三维数学模型 3 1 4 1 基本控制方程 下面给出在低含沙量条件下三维模型的控制方程【4 0 】 哈尔滨工程大学硕士学位论文 水流方程 丝+ 坐+ 业：o( 3 1 6 ) l 一+ + 一= u l j 一 积咖瑟 詈+ “警+ v 考+ w 警= 一吉罢q 窘+ 窘+ 害 c s 一， 百栅面+ v 瓦+ w 西5 一i 言+ v r 矿+ 矿+ 矿 - 1 宴+ ”宴+ v 宴+ w 堡：一土望+ v 。窑+ 窑+ 垂 ( 3 一1 8 ) 瓦棚瓦+ v 瓦+ w a z2 一i 亩w t 萨+ 矿+ 萨 - 1 警+ 群塑o x + v 考+ w 警= 一吉謇+ q 窑a x + 害+ 窑i g z 睁，。， 8 t 却 a z p 瓠 i 巩1 2 泥沙连续方程 署+ “篆+ v 雾+ w 暑一哮= 墨- 扩j s + 窘+ 害 c s z 。， 瓦瓦+ v 万+ w 西叫瓦“s+ 矿+ 萨 3 2 0 河床变形方程 儿掣+ 掣+ 掣= f l c o ( s o - ( 3 _ z - ) 式中z 为铅垂方向坐标，疋、s ，分别为悬移质近底含沙量和挟沙力。 对于剖面二维模型，种是将三维模型沿河宽方向平均，得到所谓侧向 平均二维模型另一种是直接将控制方程中与y 有关的项去掉，前者适应于河 宽较小的情况，后者适应于河道宽浅的情况。 3 1 4 2 仿真算法 由于剖面二维和三维问题和一维及平面二维问题有很大的不同，它涉及 到的关键问题都极为复杂，目前研究成果也不多，但由于它有广阔的应用前 景，因此随着泥沙运动理论的进一步研究和泥沙工作者高度的重视，三维泥 沙数学模型也一定会逐渐成熟起来。仿真算法除了上述二维所用到的算法外， 还有下述算法【4 。制体积法( f ) 、l a t t i c eb o l z t m a n n 法( l b m ) 。 1 9 哈尔滨工程大学硕士学位论文 3 2 模型方案的确定 3 2 1 一维模型的适用条件 一维数学模型通常将计算区域划分成若干小区域，计算各断面的平均水 力、泥沙要素以及上下两断面之间的平均冲淤厚度的沿程变化情况，研究对 象比较简单。但由于挟沙水流与动河床的相互作用十分复杂，即使一维问题 至今也未能彻底解决。一维模型一般用来研究来水来沙条件发生巨大变化所 引起的河床变形，这种巨大变化主要是由修建水利枢纽等人类活动所引起的， 例如修建水库等。 3 2 2 准二维模型的适用条件 一维数学模型只能给出各小河段的平均冲淤情况，对于需要了解河床细 部变形的工程建设问题，它就显得无能为力了。例如坝区引航道的淤积问题、 弯道冲淤问题等等，这些问题只有通过平面二维计算( 甚至要进行三维计算) 才能解决。高含沙水库的非恒定流泥沙数学模型基本方程由水流控制方程和 泥沙输移方程组成，由于未知物理量的个数大干基本方程数，因此为封闭基 本方程，必须引入一些补充关系式。 3 2 3 二维模型的适用条件 平面二维泥沙数学模型的研究工作早在5 0 年代初期就已开始。但由于当 时计算手段的限制，因而不得不作一些简化，应用范围受到很大限制。随着 计算机的发展二维泥沙数学模型得n t 迅猛发展，特别是在河! z i 海湾的潮流 输沙计算方面取得了较多成果，在河道的浅滩挖槽的冲淤计算等方面也有较 大进展。平面二维泥沙模型作为重要的研究与预测手段之一，它的优势在于 能够计算出全计算域内任意一点的水流流态及河床冲淤变化，与一维或准二 维泥沙数学模型相比功能上有了很大程度的增强。相对于河工模型试验而言， 平面二维泥沙数学模型不受缩尺及变态的影响，对于长河段、长时段和不同 水沙组合及不同边界条件下的水沙演进及河床变形预测，具有周期短、投资 少的优势。 哈尔滨工程大学硕士学位论文 3 2 4 剖面二维及三维模型的适用条件 剖面二维及三维数学模型主要是用于研究沿水深方向水沙运动变化对河 床冲淤影响较大的问题，前者的研究工作在4 0 年代就已开始，由于研究手段 的限制以及紊动结构难以测定，多数工作仅限于均匀流问题的研究：后者是近 2 0 年来才发展起来的，还存在很多问题有待解决 3 2 5 模型方案的确定 上述中，虽然模型的维数越高对物理系统描述得越精确，但会给仿真所 需要的时间和解算的方法带来很大的工作量。而且由于维数的增加就会遇到 很多目前理论无法解决的问题，所以就会给模型带来很大的误差。比较各个 不同维数数学模型适用条件，从既能完成仿真的需要，又可以节省不必要的 仿真过程中的工作量，我们选择非恒定流的一维模型作为本研究的方案。算 法采用先求解水流方程，得到有关水力要素后，再解泥沙方程求出河床冲淤 变形，即所谓的非耦合解法。数值计算方法采用p r e i s s m s n n 隐格式方法，这 样使计算的精度问题得到了保证。详见下章。 3 3 本章小结 本章在分析仿真模型建立必须注意的三方面问题基础上，介绍了目前应 用于黄河水沙仿真方面的四类数学模型，并分析了它们在水沙仿真中的算法 及其各个模型的适用条件。通过综合这几类数学模型，确立了本文的仿真模 型及其算法。 2 l 哈尔滨工程大学硕士学位论文 第4 章非恒定一维水沙数学模型 4 1 数学模型及其定解问题 对于狭长的库区，或者相对于库区长度来说，宽、深特征量都是微小量， 可以在整个过水断面上取所有因变量的平均值，这样在空间上就可以只考虑 这些平均变量沿流程的变化，从而问题就变成一维了。由于内陆性多沙河流 的水库泥沙一般以悬移质为主，因此，本章主要针对低含沙浓度建立悬移质 水沙数学模型，基本方程采用( 3 - 1 ) 一( 3 - 4 ) 式及( 3 - 6 ) 式，并根据模型需要作 进一步简化，对上述方程组的形式也作些改变，以流量q 、水位z 的形式写 出，河床变形公式也改写成另外一种形式。 一维非恒定不平衡输沙模型的基本方程【3 4 1 为 水流连续方程 罢+ b 祟：g ( 4 - 1 ) 水流运动方程 詈+ 芸c 手，+ 鲥塞+ g a 气势= 鼋q c a z ， 泥沙连续方程 了o ( a s ) + 昙舀+ f l c o b ( s 一鼠) ：0 ( 4 - 3 ) 水流挟沙力方程 s 。= s ( 甜，h ，) ( 4 - 4 ) 河床变形方程 以孕：p c o ( s 一) ( 4 5 ) 边界条件为 q ( x ，f ) = q o ( x o ，f ) s ( x ，t ) = s o ( x o ，f ) ( 4 6 ) ( 4 - 7 ) 哈尔滨工程大学硕士学位论文 z ( x ，t ) = z j v ( x ，f ) ( 4 8 ) 初始条件为 渊芒；茹 汁。， ( x ，f ) i 瑚= z o ( 功 以上各式中q 为流量，a 为断面面积，b 为断面宽度，z 为水位，q 为旁 侧入流，为旁侧入流流速在水流方向上的分量，这里q 、k 近似为零。s 为断面平均含沙量，为断面平均挟沙力，为河床高程，g 为重力加速度， 口为动量修正系数，反映河道断面流速分布均匀性系数。为泥沙非平衡 恢复饱和系数，甜为泥沙颗粒沉速，以为泥沙干容重，k 为流量模数， 一 1 ， k = a c r ，其中c = 三r y 6 ，r 为水力半径，n 为曼宁糙率系数。 4 2 参数的选择 4 2 1 泥沙沉速的计算 泥沙沉速的大小反映着水流紊动强度与泥沙重力二者之间的相对作用， 泥沙沉速越大，则泥沙落淤的倾向越强。这样泥沙的淤积也越多。 黄河下游泥沙的运动过程比较复杂? 泥沙的形状、絮凝水流、含沙量等 这些都是影响泥沙沉速的重要因素。泥沙沉速是影响水流挟沙力的一个重要 因素，泥沙沉速公式的选取直接影响水流挟带泥沙能力的准确性。舒安平的 分析结果表明，张红武水流挟沙力公式中引用费祥俊的公式计算沉速，得出 的挟沙力与天然河流实际比较接近。另外，该泥沙沉速计算公式也反映了黄 河下游含沙量越大，泥沙沉速越小的输沙特性。 本模型采用既能反映泥沙的形状、絮凝水流、含沙量等因素的影响，又 可满足水流挟沙力计算要求的费祥俊沉速计算方法。该方法以临界浓度s 。 来划分含沙水流，当来沙量小于临界浓度时，泥沙颗粒沉速采用水电部1 9 7 5 年水文测验试行规范中的沉速公式。仅当s 。 s ，o 时，才可能出现宾汉体。 含沙量变化的影响由下式反映1 4 2 j ： 鼠o = 1 2 6 s m 3 2 ( 4 - 1 0 ) 哈尔滨工程大学硕士学位论文 q n n ，九口丫1 2 、(41192 ) s v m _ 0 - o 2 ( 1 9 子) 一 妒生笋 ( 4 - 1 2 ) 只= p h h i + p h h 2 + 1 p h h 3 + p h h 一4 + p h h 5 ( 4 - 1 3 ) 式中：s 加为牛顿体转变为非牛顿体的临界浓度；。为悬浮液粘滞系数接近 无穷大时的极限浓度；只某一粒径的级配；研某粒径的平均直径；d 。、d 。i 。 为某一粒径组的最大最小粒径范围；将表中数值代入得到 s ，。= o 9 2 0 2 ( 1 9 7 6 2 2 6 3 4 ) = 0 5 4 3 6 ，所以s ，o = o 1 7 9 2 ，当且仅当 s 、， s 。o = 0 1 7 9 2 时候才可能由牛顿体转化为宾汉体。用混合含沙量表示为 s = 几s ，= 2 6 5 0 0 1 7 9 2 = 4 7 4 8 8 k g m 。所以只有大于此临界浓度数值时候 才可能由牛顿体转化为宾汉体。而本模型中最大含沙量都没有达到这这个临 界浓度，所以运动体都是牛顿体或者牛顿体和宾汉体的过渡体，这样就给计 算带来了方便。 单个泥沙在清水中沉速其大小都可用张瑞瑾公式【4 j j ： 0 3 0 2 厣事耳 1 3 9 5v o( 4 1 4 ) d 国s = 国o s v ” ( 4 1 5 ) 式中：为单个泥沙在清水中的沉速；。为群体沉速；s ，为体积比表示的 含沙量；i l l 为经验指数，取值为8 ；v 0 为清水水流运动黏性系数此处 v o = 1 0 0 3 1 0 。m 2 s ，d 为泥沙平均粒径；p 。沙的密度此值与形成沙子的 岩石种类有很大关系，一般都在2 5 5 0 k g m 3 2 7 5 0 k g m 3 之间变化，经查 资料黄河通常取风= 2 6 5 0 k g 删。来计算；p o 为清水密度p o = 1 0 0 0 堙m 3 。 当来沙含沙量大于临界浓度时，认为此时浆液已由牛顿体变为宾汉体。 这时先以包含全部粗细颗粒的悬浮液计算流变参数，用下式计算雷诺数，最 后用阻力系数一雷诺数关系曲线试算各级颗粒的沉速【4 引。 耻t 筹矗 阵m ) g l 叩五+ 面j 斗m 甜( 学 4 玎 件m 7 b = 9 8 1 0 2e x p ( b 占+ 1 5 ) ( 4 1 8 ) 占：墨二坠(4-19) 式中：r 为雷诺数；y 。为悬浮液的容重；叮为悬浮液的相对粘滞系数： 为悬浮液宾汉体极限剪切力；b 为系数，取8 4 5 ；非均匀沙的代表沉速采 一n ( 4 2 0 ) 。= 最 、 式中：最为第k 组泥沙( 相应沉速吼) 占全沙的重量比；n 非均匀沙的分组 4 2 2 水流阻力 其糙率一般都比较大，随着水位的上升和流量、流速的加大，曼宁系数逐渐 减小，在一定的水流条件下，曼宁系数达到一个最低值；滩地糙率也是十分 复杂的，当滩地被淹没不多、农作物未遭水淹时，随着水深的增加和相对糙 率的减小，曼宁系数逐渐减小，最后接近一个常数。如漫滩较广、农田被淹 着水流条件、流态的改变，河流的阻力也在不断的变化，因此，进行阻力计 算时合理的做法应该是根据计算过程中每一时刻的水流条件来确定河道滩、 哈尔滨工程大学硕士学位论文 槽对水流的阻力。 由于冲淤河流的河床为动床，床沙阻力主要是由沙粒组成，它随着河道 冲淤在不断地发生变化。当河道发生淤积时，河道的床沙被细化，床面阻力 减小；反之，当河道发生冲涮时，河道的床沙被粗化，床面阻力增大。因此， 模型对河床阻力随河道冲淤不断变化这一特点进行了初步考虑。 开展黄河下游洪水模拟计算，沿程水位是最重要的内容之一。若仅仅为 保证河床冲淤量与原型相似，而对水力摩阻特性方面要求不高，可使水位出 现较大的偏差。因此，专门引入一套新的糙率计算公式，不仅可以描述水力 泥沙因子的变化对摩阻特性的影响，也能反映天然河道中各种附加糙率的影 响，( 黄河河道宽浅，r 。h ) 即式l “j ： 珂= 生 0 4 9 审0 7 7 + 等c ，一和n c 铲 5 ) 1 汁z ， 4 9 h 6 该式中河滩上的摩阻厚度最，即为当量粗糙度，可根据潍地植被情况等， 由水力学计算手册查得，而在主槽内，黄河沙坡尺度及沙波波速对摩阻特性 有较大的影响。对这一复杂的影响过程，根据动床模型试验资料，建立了摩 阻厚度反与弗汝德数( f r ：兰) 等因子之间的经验关系 2 2 】，即 、劝 茂= d ( 1 + 1 0 8 1 - 1 3 f r 。 5 ( 1 - f r 3 ) 1 ( 4 2 2 ) 式中：职。为床沙中值粒径 冲积河流局部阻力一部分是由河流的平面形态特征引起；另一部分是由 孤立的特大粗糙所引起。由于黄河下游河段河宽变化引起的局部水头损失。 至于对河流中不多见的形态特征如孤石、小岛、浅滩等，因资料有限，模型 阻力计算中未加考虑。 4 2 3 挟沙力的计算 水流挟沙力是表征一定来水来沙条件下河床处于冲淤平衡状态时的水流 挟带泥沙的综合指标，也是研究水流模型相似率的不可缺少的概念。特别黄 河泥沙数学模型中水流挟沙力的确定更是关于仿真结果的关键所在。越来越 哈尔滨工程大学硕士学位论文 多的研究结果表明，水流挟带的泥沙影响水流的物理性质和絮动结构，从而 又影响其能量损失、流速和含沙量分布，因此，多沙河流的挟沙力所反映的 冲淤平衡是相对的，它不仅与水流条件等常见的因子有关，而且还同本河段 含沙浓度等因素有密切关系。对于黄河下游，由于受上游水库“畜清排浑” 运用的影响，黄河下游既会出现高含沙水流，也会出现低含沙水流，因此含 沙水、瀛特性要求水流挟沙力公式既适甩于高含沙水流，又能适用于低含沙水 流；另外，由于泥沙含量高和细颗粒的存在，改变了挟沙水流的流变、流动 和输沙特性，河道冲淤计算也要求选取的水流挟沙力公式能反映其“多来多 排”的输沙特性。 关于水流挟沙力的研究，长期以来，国内外工程和学术界的许多专家学 者或理论出发，或根据不同的河渠测验资料和实验室资料，提出了不少半理 论的、半经验的或经验性的水流挟沙力公式。远在两千年前，我国西汉时期 的张戎就开始对泥沙左右的角度分析河床冲淤规律，由此提出以水排沙的治 河主张，并对水流挟沙能力随流速大小而变化的关系有所认识；1 9 1 4 年西方 学者g i l b e r t 通过室内水槽试验，对水流输沙能力进行了系统的研究；2 0 世 纪4 0 年代末前苏联出现了著名的3 a m a p uh 公式，主要用于渠系挟沙力的 计算，在5 0 年代和6 0 年代前期，e i n s t i n 通过求单宽悬移质输沙率的方法 分析水流挟沙力问题，b enhnah0b 从重力理论出发，认为浑水在单位 时间的能量损失除用于克服阻力做功外，还用于悬浮泥沙做功。我国早期以 张瑞瑾1 4 1 和沙玉清【5 1 为代表，后来曹如轩“、杨志达”i 、张红武8 1 等人都相 继投入研究，并取得一定的成果。由于考虑到低含沙量和各公式的适用条件 及精度问题，本模型的水流挟沙力选用张红武公式【2 2 】： ，、r = 2 5 1 喜( 0 0 0 2 2 + s v ) v 3i n f 土1l ( 4 2 3 ) i i 姓咖国。t , 6 d 5 0j i l。 。 国j = 脚0 s v “ y 。= 1 0 0 0 + 0 6 2 2 s 。 ( 4 2 4 ) ( 4 - 2 5 ) 哈尔滨工程大学硕士学位论文 ， 一 、 k = k o l l 4 2 c ，( o 3 6 5 一s ，) ) ( 4 2 6 ) 式中：d ；。为床沙中值粒径；。为泥沙群体沉速；品为体积含沙量；0 9 。清 水沉速；d ；。为悬移质泥沙的中值粒径；浑水容重；k 为卡门常数；k o 为 清水时候卡门常数此处为0 4 。上式采用k g 、m 、s 制。 可以看出，张红武公式反映了含沙量及泥沙群体沉速对水流挟沙力的影 响，并通过卡门常数考虑了浑水流速分布的影响，还通过引入床沙中径d ，。， 使悬沙的挟沙能力与河床粗度联系起来。在模型的调试过程中发现，在不同 的河床条件下，该公式的适用性较强。 4 3 相关水力要素的预测 4 3 1 拉格朗日( l a g r a n g e ) 二次插值 圆为断面水力要素a 、b 、r 是水位z 的函数，所以计算过程中需要建立 各断面水力要素与水位的关系，来预测下一时刻不同断面的水力要素值，通 常采用拉格朗日( l a g r a n g e ) 二次插值公式。即内插法和外推法确定水沙模型 中相关的参数l ”j 。 有时，仅知道函数f ( x ) 在一系列点x l ，x 2 ，x 3 x n 上的数值( 假如 x l x n ) ，但并没有f ( x ) 的解析式来求得任意一点的数值。例如，f ( x 1 ) 可 能来源于一些物理量值或来源于不能用简单函数形式表示的冗长的数值运 算。通常x i 是等间距的，但不是必须如此。现在的任务是，画一条通过( 可 能越过) x i 的平滑曲线，来对任意x 估算f ( x ) 的值。如果所求的x 处于x i 的最大值和最小值之间，则称为“内插法”，如果x 超出范围则称“外推法”， 后者冒险性更大些。内插法和外推法是在已知点之间或越过已知点，用某些 可能的函数形式模拟此函数。函数的形式应该是很常见的，而且能够近似于 应用中可能遇到的几大类函数。函数中最常用的形式为多项式。有理函数也 被证明是非常有用的。拉格朗日( l a g r a n g e ) 二次插值理论公式如下式： y = y t ( x x 2 ) ( x x 3 ) ( x l x 2 ) ( x l x 3 ) 】+ y 2 ( x x 1 ) ( z x 3 ) ，、 ( x 2 一x o ( x 2 一x 3 ) 】+ y 3 ( x x o ( x x 2 ) 【( x 3 x 1 ) ( x 3 一x 2 ) 】 哈尔滨工程大学硕士学位论文 其中x 。、x ：、x 。为插值基点，y ，、y 。、y 。为插值基点对应的函数值，它们都是 已知的。利用上式可求插值点为x 的插值函数值y 。本模型中采用以往来水 来沙条件下，同一断面的不同数据进行内插。这样保证了计算的精度要求。 4 3 2 断面水力动要素的插值函数 将水位z 当作x ，断面水力要素如a 、b 、r 等作为y ，就可求得相应于断 面不同水位的水力要素值。 下面是得到的黄河小浪底河段h h i 断面的几组水位z 与断面面积a r e a 的数据，a ( 2 3 0 4 7 ，9 1 3 3 3 ) ，b ( 2 3 1 3 7 ，8 2 3 0 8 ) ，c ( 2 3 3 0 3 ，6 8 1 8 ) ， d ( 2 3 3 7 6 ，6 8 6 9 6 ) ，选择a 、c 、e 三点按照上式( 4 - 2 7 ) 用拉格朗曰插值得 到的曲线，如图4 i 。 图4 1 断面面积与河流平均水位高程曲线 由内插法求得的曲线得到b ( 2 3 1 3 7 ，7 8 7 6 5 ) ，与实际测量的值绝对误 差也仅在4 3 可以看出断面面积与水位高程具有很好的相关性。 下面也是给出断面的宽度w i d t h 与水位z 之间的关系，采用的点有四组 a ( 2 3 0 4 7 ，1 9 1 1 ) 、b ( 2 3 1 3 7 ，1 7 1 8 3 ) 、c ( 2 3 3 0 3 ，1 4 1 8 ) 、d ( 2 3 3 7 6 ，1 4 2 2 ) 随意的选取a 、c 、d 三组数据作为拉格朗日插值的基点，按照上式( 4 2 7 ) 求得z w i d t h 的函数曲线。如图4 2 ： 哈尔滨工程大学硕士学位论文 图4 2 断面宽度与河流平均水位高程曲线 可以通过曲线求得b ( 2 3 1 3 7 ，1 6 4 7 7 4 1 ) ，与真实值点误差仅为2 3 3 ， 几乎完全和插值曲线吻合，证明具有相关性。 下面也是给出断面的平均水力半径r 与水位z 之间的关系如图4 3 ，采 用的点：a ( 2 3 0 4 7 ，7 4 3 ) 、b ( 2 3 1 _ 3 7 ，1 5 8 5 ) 、c ( 2 3 3 0 3 ，2 0 8 5 ) 、d ( 2 3 3 7 6 ， 1 6 1 4 ) 如图4 3 断面的截断面水力半径与平均水位高程曲线 通过曲线求得水位在b ( 2 3 1 3 7 ，1 7 4 5 ) 与实际测得到的误差在1 0 0 9 。 根据河流预测规范误差在2 0 d a 内，所以误差完全符合仿真的要求。通 过求得的相对误差数据可以用它来预测不同平均水位高程下的截断面积、河 槽宽度、截断面水力半径。 4 4 数值计算方法 天然河流的边界条件十分复杂，为了便于求解，我们做如下假定： 哈尔滨工程大学硕士学位论文 1 、假定计算中，河流为顺直河道，在不考虑流量的渗透及其蒸发，以及 沿程的引用和农工业用水。 2 、定河床发生冲淤过程中，在每一个短时段内河床变形对水流条件影响 不大，即可采用非耦合解法进行计算，只要限制每个计算时段内冲淤量不太 大就可满足该假定条件。 3 、有限差分法是将偏微分方程离散成代数方程，进而求得离散节点上的 近似解的方法，在理论上可以获得较为可靠的推导和证明，因而在数学模型 计算中应用广泛。根据差分格式的不同，人们建立了各式各样的计算方法， 如t v d 格式、l a x w e n d r o f f 格式、a b b o t t 隐格式、p r e i s s m a n n 隐格式等等。 本文就采用p r e i s s m a n n 格式能够通过选择加权因子以获得较高的精度，可以 避免数值扰动现象，并且计算过程中只涉及四个网格点，使得内外边界的处 理十分方便，得到广泛应用。一维非恒定泥沙数学模型的计算采用非耦合方 法，首先求解水流连续方程和运动方程，然后求解水流挟沙力泥沙连续方程 和河床变形方程具体求解过程如下。 4 4 1 水流方程的离散 水流方程属于一维圣维南方程组，求解一般采用有限差分法，由于一般 显格式差分格式关于时问的导数项取前差或后差，只有一阶计算精度。这样， 对于纯对流方程时间项取中差会导致计算格式的不稳定。所以在研究非恒定 问题p r e i s s m a n n 格式( g e n e r a lp r e i s s m a n ns c h e m e ) 既吸取了中差具有二阶 精度的优点，同时又可避免计算的不稳定的优点。其基本思路为：在空间导 数取平均加权数的基础上对时间导数项也取加权平均，适当地选取两种权因 子，可以达到较高的计算精度，而又能保证计算稳定性【4 5 1 。p i c h t n y e r 首先 于1 9 5 7 年使用隐格式研究热传播问题，后来隐式方法在水流波动问题上得以 广泛应用，本文用p r e i s s m a n n 格式1 4 6 j 因为该格式在处一维非恒定水流转播 和泥沙问题时精度很高。 p r e i s s m a n n 使用四点偏心差分格式，此格式的网格形式如图4 4 所示， 在每一个结点上同时求出流量和水位。因变量及其导函数的差分形式1 4 7 】为 哈尔滨工程大学硕士学位论文 ；i i i ；i ；i ；i i i ；i ；i i i ；i i 图4 4p r e i s s m a n n 叫点倔心格式 ，( m ) = 【9 1 ”1 + ( 1 - o ) f j ”】+ ( 1 一) 暾”1 + ( 1 一e ) f j ”】 ( 4 2 8 ) a f ( m ) ：8 墨! ! ：二互：= ! + ( 1 一们五! ! ：二五：( 4 - 2 9 ) 咖xx 型玉曼坌：! ! 五! ! ：二互! ! ：! ! ! 二兰坚五：二二五：! ( 4 3 0 ) 一 u v 1 、当目0 5 和0 5 时候，差分格式精度具有一阶的： 2 、0 = = 0 5 格式，y = 坛，l 具有二阶精度； 3 、又当0 = y = 0 5 和y = ，= 1 ，差分格式具有更高的精度，至少具有四 阶精度。 此处= 0 5 ，将差分格式代入圣维南方程可以得出p r e i s s m a n n 四点隐 式格式下的差分方程，同时利用线性替换，在q 、z 之间建立线性关系， 结合初始条件和边界条件，利用“追”“赶”法求出下一时刻的水位和流量。 由于隐式解法最后往往需要通过解一个矩阵卡求出计算结果，四点偏心所提 出的算法则避免了对矩阵的直接求解，而通过线性交换关系把矩阵的迭带计 算过程与水力要素的中间处理结合起来，形成一个“追”和“赶”的过程， 使得计算过程中可以对水力要素的变化进行控制，比如在处理河道有侧流加 入时p r e i s s m a n n 算法的优越性就得到了体现。并且这样一种处理没有影响到 计算的精度，因而受到广泛的重视，在上面的推导中，空间权重因子的引入 实际上是为了增加人工阻尼效应来抑制数值弥散，理论和实践证明这样的方 法达到了预期的目的。 哈尔滨工程大学硕士学位论文 把方程( 4 1 ) 变为： 盟+ 鱼：o b o xo t 按照以上三式对方程上式离散得： 2 ，+ 1 “1 一。h ”+ z ，”1 一z ，”2 【o ( o ”1 + l 一9 ，”1 ) + ( 1 一目) ( q ”，+ 1 一o ”朋n 二。=二一i-一i， 2 a t a x o ( b j + l ”干1 + b ”) + ( i o ) ( b j + l “+ b ”) ( 4 3 1 ) 同理对( 4 2 ) 离散可得到如下： 芝1 尘塑n ：兰! 二堕+ 旦i 堕鱼掣一塑掣l + 坐l 塑曼一旦! l + 2 肿 缸la j + l “一j 缸l 铂”4 ”j 十丛乎+ ( 1 一。笙号互l 在蔷丛+ ( 1 一句生字卜 啦黟+ 镨p 心笋+ 等p ( 4 3 2 ) 以上两式中j = 1 ，2 ，3 ，n _ 1 当0 5 扑剃 洚 南2 南2 可1 玎井刳 而2 雨可2 而等r 帝r 刮 ( q i ”1 ) 2 = ( q j ”+ a q j ) 2 ( q j ”) 2 + 2 q j ”q j g ”1j 9 ”1p g “i g “i + 2 i g ”i 9 3 5 ( 4 4 9 ) ( 4 5 0 ) ( 4 - 5 1 ) 哈尔滨工程大学硕士学位论文 鲋，= d d ，如，”( & ) = b j ”业j ( 4 5 2 ) a k j = d k ”，i d z j ( a z ，) 把( 4 4 8 ) ( 4 - 5 3 ) 及( 4 3 1 ) 式代入( 4 3 2 ) 式整理得 a l j 扯i n + b ! i a q j n + c 1 j i 七d 2 j q i + e 2 j = 0 式中系数为 毛：面a s t 弋2 a ( f j j 丽+ 1 ) 2 b 广 j + j + 鲋_ + ，) + 靠”川 p 川等剡 ( 4 5 3 ) ( 4 - 5 4 ) = ”肌卜锚掣。， 驴- + 警( 鬻) + 2 9 触裂 = 芸p 铲一g 口_ - 十毋，+ 如n j + _ z n j 册 + 驰拶 b j 2 码a j 矧 小警( 酬恤融r 锣掣 ( 4 5 6 ) ( 4 5 7 ) ( 4 - 5 8 ) = 尝眸一警+ 卿圳阶一力) + 刎警治跏 其中系数式( 4 4 3 ) ( 4 4 5 ) ( 4 5 5 )( 4 5 7 ) 中还包括有关z 的导数项。d b 、坚 舷d z 一般使用差分方法离散上述关系式可得到 d b ? 一b ? 一b r - ”一 d z ”2 ”一。j ”一1 ( 4 - 6 0 ) 哈尔滨工程大学硕士学位论文 k ”k ，”一1 ( 4 6 1 ) 这样，通过( 4 3 6 ) 将原来拟线性方程组( 4 3 3 ) 、( 4 3 4 ) 作了线性化处理，新 方程组( 4 - 4 2 ) 、( 4 - 5 4 ) 中的系数a 、b , j 、c q 、d 矿e o ( f _ 1 、2 ) 不再包含有变量 z 和a q ，而仅仅依赖于n a t 时层上的结点函数值，因此是一个线性方程组。 但是方程( 4 - 4 2 ) 、( 4 - 5 4 ) 只能看成非线性方程( 4 3 3 ) 、( 4 3 4 ) 的近似，要使 线性方程( 4 4 2 ) 、( 4 5 4 ) 能较好的逼近方程( 4 3 3 ) 、( 4 3 4 ) ，必须满足 鲈f o ，所以 1 g 1 一只止l ( 4 7 7 ) ( 2 ) q = 0 ( t ) 给定： 因为q 1 = q l ”1 一q 1 ”，并设e = 0 ，据( 4 7 5 ) 式可得 g l = q l ( 4 7 8 ) ( 3 ) q = q ( z ) 给定： 因为姨= 9 ”1 9 l ”，q 1 ”为己算出的值，从给定的边界条件有 q 1 “：她”) + 掣止。，( 4 - 7 9 ) 2 z i 所以有 q l ：q ( 毛) 一q l ”+ 皇罂z ( 4 8 0 ) ( d z l 上式对照( 4 7 5 ) 式可得 曩：掣( 4 _ 8 1 ) o - z 1 g i = q ( z l ”) - q , “ ( 4 8 2 ) 在下游边界同样也有三种不同的形式用于确定z 。： ( 1 ) z = z ( t ) 给出： 直接离散& 。可得 a z = 知”1 2 ”( 4 8 3 ) 实践应用表明，计算时取计算值知”较佳。 哈尔滨工程大学硕士学位论文 ( 2 ) q = q ( t ) 给定： 由于玑= q 。”1 一q ”，并由( 4 6 2 ) 式得 血。：旦 q ? 一g n f n ( 3 ) q = q ( z ) 给定： 由于q 一：q o 。一) + 掣止。，可得 “n a q n ：_ d o k n 代入( 4 6 2 ) 式可以得到 ( 4 8 4 ) ( 4 8 5 ) 位= 巧 ( 4 8 6 ) 尝一目 本文上下游边界条件都是利用q = o ( t ) 给定来求得的结果。 4 4 4 泥沙连续方程求解 天然河流泥沙运动种类依性质主要可分为推移质、悬移质及冲淤质，其中 推移质及悬移质合称为河床质。总输沙量为与河床变动有关的泥沙输移量， 冲淤质亦是以悬浮型态输送，但存在底床之量极小，可视为不参与底床改变 量，一般忽略不考虑。 在一般的动床模式中，多采用河床质公式来计算输沙量，在将所计算之 输沙量带入泥沙连续方程中，来求得底床之变动量，这种做法是建立在推移 质与悬移质淤积运动过程中，对流动均有同步反应大假设上，但这和泥沙运 动之基本特性并非完全一致。因此为了更确实反应泥沙运动的基本特性，常 将推移质及悬移质分开来计算，推移质利用推移质计算公式计算，悬移质则 直接解析悬移质方程计算，并将分开计算所得之推移质输沙量及悬移质浓度 带入泥沙连续方程式以求解河床变形情形。 由于推移质运动的时间尺度远大于悬移质运动的时间尺度，因此对于短 延时的模拟，般将推移质忽略不计，总输沙量只考虑悬移质。本文就是如 此忽略推移质只考虑悬移质。 4 1 哈尔滨工程大学硕士学位论文 悬移质泥沙连续方程一般采取有限差分的迎风格式离散求解 1 4 】【15 1 ，但由 于该方程对数值解法的要求较高【2 5 1 ，采用一般的差分格式，精度不能得到保 证，含沙量的计算结果有时可能会出现不合理的数值，甚至为负值；采用高 精度的差分格式，不仅计算复杂，而且计算量太大，在长时间、长河段的河 床冲淤计算中应用比较困难。为了避免含沙量计算结果的不合理现象，同时 不致使计算量过大，本文参考文献 4 9 ，采取在相邻时层之间用差分法求解， 而在同一时层上用分析法求解。这种方法的实质与恒定流状态下求解悬移质 泥沙运动方程式的思路是一致的，将式( 4 3 ) 中的时间导数项用隐式差分格 式离散。 首先将( 4 3 ) 式两边同除以q 可改写成 三箜+ 箜+ 4 c o ( s 一1 ：、0 ( 4 8 7 ) “拼出 q 式中u 为断面平均流速，q 为单宽流量。将上式中的时间导数项用隐式差分 格式离散【5 0 得 击譬+ c 耖1 + c 争c s k + l _ s k + l m s s ， 由于式( 4 8 8 ) 仅随自变量x 而变化，因此可改写为常微分方程 c 扩+ c 争+ 南s k + 1 = c 争川击( 4 - 8 9 ) 式中u 、q 、s 斛1 、s 、卢及国均为随x 变化的已知函数，可以看出上式 为一阶非线性齐次微分方程，类似于方程：6 + p ( x ) y = q ( z ) 其通解为y = c e j p 3 出+ 8 一jp 。出f q ( x ) e j 。p 。冲d x ，在此处 p ( 州争“+ 赢，q ( x ) _ ( 争驴1 + 瓦1 s 。 4 2 哈尔滨3 - 程大学硕士学位论文 p ，i 扭钞+ 南出 i 争k + l s , k + l + 志州“钞+ 南出出十q ( 4 9 0 ) 由( 4 9 0 ) 式可知，要进一步积分必须确定1 1 、q 、s 。“1 、s 的变化规律， 按照恒定流不饱和输沙模型的近似计算方法，可以利用这些量在短河段内的 平均值近似代替这些量，也可以假定这些量在短河段内呈线性变化 1 。为节 省计算工作量，本文采用前一种近似方法，即认为上述变量在短河段内为常 数，并用其平均值代替。考虑到在进口断面s “1 = 置“1 ，据此可确定积分常 数c ，并由( 4 - 9 0 ) 式得 艮川：牡吨秒+ 击+ 譬善笺群里卜。 【( 钞+ 击m ， ( 4 - 9 1 ) 式中i = l ，2 ，3 ，一n - 1 ，地小鼠及s 分别为缸。库段内的平均流速、单宽流 量、挟沙力及含沙量。 4 4 5 河床冲淤量计算 在一定的水流和泥沙综合条件下，水流能够携带的临界含沙量，当水流 悬移质的含沙量大于临界值时，水流处于饱和状态，河床将发生淤积：反之 当不足这一临界值时，水流处于未饱和状态，水流将向床面寻求补给，河床 将发生冲涮，通过淤积或冲涮，使悬移质中的泥沙含量恢复临界值，达到不 冲不淤积的新平衡状态。因此透过水流的含沙量与水流挟沙力之间的关系， 配合上泥沙连续方程式则可计算河床冲淤量并建立输沙动床模型。 按照上述理论将( 4 5 ) 式写成差分方程各断面淤积厚度为 ：螳( s ，川一s , j n + i ) ( 4 - 9 2 ) 然后利用式 2 0 ，“+ 1 = 。o l ”+ a z o j ( j = 1 ，2 ，n 一1 ) 即可计算出淤积后河床高程而第j 库段的淤积量为 a w j 。( a z o j b j + a z o j 十1 b 川) 子a x ( 4 9 3 ) ( 4 9 4 ) 哈尔滨工程大学硕士学位论文 恢复饱和系数口理论上应为河底含沙量与垂线平均含沙量的比值，是一 个综合系数，一般计算中都根据现有的计算成果，经过本人大量的翻阅资料， 发现应用于黄河中的恢复饱和系数要比一般的计算情况下要小得多，所以黄 河中的恢复饱和系数当淤积的时候取为0 0 1 ，冲涮时候取0 0 3 。 4 4 6 河床淤积泥沙级配的调整 床沙级配一般是计算水流挟沙力的基本条件，因此在验算中要确定起始 床沙级配情况。床沙级配一般测次较少，精度也不高，由于资料有限，本论 文资料是在2 0 0 2 年汛前的床沙级配为依据，并参考2 0 0 2 年汛后的资料进行 编排求其平均值( 因2 0 0 2 年小浪底干流调水调沙后的级配突然变细不符合实 际情况) ，小浪底水库的级配是根据以下几个水文部的级配插补得到的水文 站中间断面级配则是由小浪底水库干流和五个基本水文站的级配值编排见表 4 2 。验证算法正确性时候，采用淤积或冲涮泥沙级配为p b 8 ，如表4 2 。 表4 2 泥沙级配表 站 小于某一粒径的沙重百分数 名粒径级，、m m ) 0 0 0 40 0 0 8o 0 1 60 0 3 20 0 6 2o 1 2 5o 5l h 50 0 7 7 o 1 1 40 2 0 3 0 5 0 50 8 3 9 0 9 7 7 0 9 9 81 o h 90 0 6 50 1 2 30 2 1 2 0 5 0 8 0 8 4 00 9 7 01 io 1 0 h 1 7 0 0 7 8 0