（水利水电工程专业论文）梯级水库非恒定流研究.pdf

华中科技大学硕士学位论文 摘要 明渠非恒定流是水力学中最常见的工程实际问题之一，在很多工程中都需要对 非恒定流进行预报。在水电站系统中，非恒定流预报的结果对水电站工程的规划、设 计和运行等都会产生很大的影响，因此研究明渠非恒定流问题是十分必要的。刁7 本文研究了国内外非恒定流的现状和发展方向结合“清江梯级水电站调度仿真 系统”项目的实际，以圣维南方程和有限差分法为理论基础，建立了描述梯级水库非 恒定流的数值模型，并探讨了它的数值解法。从有效利用水资源和提高水电站经济效 益角度出发，本文根据非恒定流的特性与上下边界条件以及初始水位和流量之间的关 系，在满足航运和电力系统调峰任务这两个相互制约的条件下，建立了一个“弃水量” 最小的优化模型，为此提出了一种动态调节法，该方法主要从两方面入手满足优化模 型的要求，一是调节作为“反调节”水库的下游电站的日放水过程，另一方面是选择 一个最佳的初始条件。 最后通过清江高坝洲下游河道的计算实例验证了本文建立的数学模型和计算方 法，通过三峡一一葛洲坝一一枝江的计算实例验证了本文针对梯级水库提出的“弃水 量”最小的优化模型和相应的动态调节法。通过分析两个实例的计算结果可知，本文 的解决方案是合理的也是可行的，计算和预测的结果对水电站的运行有很重要的参考 价值。 关键词：非恒定流圣维南方程组有限差分 梯级水库水位约束动态调节 i 华中科技大学硕士学位论文 a b s t r a c t u n s t e a d yf l o wi no p e nc h a n n e l si su s u a l l ym e ti n a c t u a l p r o j e c to fh y d r a u l i c s a n d n e e d e dt op r e d i c ti t sc h a r a c t e r i s t i c ss u c ha sw a t e ra n dr u n o mi nh y d r o p o w e rs y s t e m t h e p r e d i c t i o nr e s u l to fu n s t e a d yf l o wm a yg r e a t l yi n f l u e n c ei t sl a y o u t ，d e s i g na n dr u n n i n g ，s oi t i sn e c e s s a r ya n ds i g n i f i c a n c et or e s e a r c ht h eq u e s t i o n t h i sp a p e rr e v i e w e dt h er e s e a r c hs t a t u sa n dd e v e l o p m e n td i r e c t i o no f u n s t e a d yf l o wa t h o m ea n da b r o a d ，b u i l tam a t h e m a t i c sm o d e lt od e s c r i b et h eu n s t e a d yf l o wi nc a s c a d e r e s e r v o i rb a s e do nt h et h e o r yo fs a i n t - v e n a n te q u a t i o n sa n df i n i t ed i f f e r e n c ea s s o c i a t e d w i t h p r o j e c to f c o m b i n e d o p e r a t i o no fq i n g j i a n gc a s c a d eh y d r o p o w e r s i m u l a t i o ns y s t e m a n dd i s c u s s e dd e e p l yt h en u m e r i c a ls o l u t i o n f r o mt h ep o i n to f e f f i c i e n c yu s i n g t h ew a t e r r e s o u r c ea n di n c r e a s e dt h ee c o n o m i cb e n e f i t so f h y d r o p o w e r s ，a n o p t i m i z em o d e lw h i c h i s t ow a s t ea sl i t t l ew a t e ra sp o s s i b l eu n d e rt h ec o n d i t i o n so f m e e t i n g t h es i t u a t i o no f s h i pa n d r e q u e s to f p e e ko p e r a t i o ni np o w e rs y s t e mw a se s t a b l i s h e d i nt h i sp a p e r t h i s p a p e r o f f e r e d ad y n a m i co p e r a t i o nm e t h o dt os e t t l et h ec o n f l i c tb e t w e e nt h er e s t r i c t i o n sa n d t a r g e tf r o m t w oa s p e c t s ，o n ei st oa d j u s tt h ed o w nb a n dc o n d i t i o nw h i c hi st h ed a yr u n o f fc o u r s eo f r e v e r s e r e g u l a t i v eh y d r o p o w e ri nd o w n s t r e a m ，t h eo t h e ri s t oc h o s ei n i t i a lc o n d i t i o n s a c c o r d i n gs o m e r o l e s f i n a l l y , t h ep a p e rp r o v e di t si d e ab yt w oc a l c u l a t i o ne x a m p l e s t h er e s u l to f g a o b a z h o u c h a n 翻i a n g k o uo f q i n g j i a n gv a l i d a t e dt h em o d e ld e s c r i b i n gt h eu n s t e a d yf l o ww a sr i g h t t h er e s u l t a n a l y s i so fs a n x i a g e z h o u b a - z h i j i a n gr i v e rs h o w e dt h ed y n a m i cr e g u l a t i v e m e t h o dw a sr e a s o n a b l ea n du s e a b l ea n dh y d r o p o w e rr u n n i n gc a nc o n s u l tt h ep r e d i c t i o n r e s u l t k e yw o r d s ：u n s t e a d yf l o w , s a i n t v e n a n te q u a t i o n s ，f i n i t e - - d i f f e r e n c e c a s c a d er e s e r v o i r , w a t e rr e s t r i c t i o n ，d y n a m i ca d j u s t i n g 华中科技大学硕士学位论文 1绪论 1 1 课题简介 清江是长江出三峡后的第一条大支流，该流域水能资源丰富，主要集中在干流 中下游，约占全流域可能开发量的8 5 8 ，占干流可开发量的9 7 4 。目前中下游隔 河岩、高坝洲水电站已建成，隔河岩装机容量1 2 0 0 m w ，年电量3 0 4 亿k w h ，下游 高坝洲为反调节水库，装机容量2 5 2 m w ，年电量8 9 亿k w h ，上游水布垭水电站已 完成设计，在施工准备阶段，水布垭装机容量2 0 0 0 m w ，年电量3 9 2 亿k w h 。清江 梯级三大水电站建成后，总装机容量3 4 5 2 m w ，联合运行年发电量8 1 6 亿k w h ，为 湖北省电能供应基地之。 梯级水电综合调度不仅是一个极富挑战性的前沿研究课题，而且亦是国际学术及 工程界研究的热点问题，具有深远的社会效益和巨大的经济效益，同时其研究成果对 推动我国水电行业的科学技术进步也有重大的科学意义和工程实用价值。自然条件 下，清江只能局部河段通航，梯级建成后恩施以下河道基本上渠化，年通过能力大为 提高，为库区沿岸运输和物资交流带来了便利条件，必将对促进当地资源的开发和经 济的发展带来积极的影响，而电站日调节造成的非恒定流会对航运会产生很多不利影 响，因此必须研究满足航运要求的水电站日运行方式。 “可视化清江梯级水电联合调度仿真系统”研究项目主要是以清江流域防洪、发 电、航运、环保以及经济可持续发展等综合效益为总目标，充分发挥清江水电工程的 综合效益，并提供最大的调峰发电效益，为国内外其它大型水电厂梯级水电综合调度 决策过程的自动化、智能化工程起到示范作用，研究成果达到国际领先水平。该项目 由以下几大模块组成：梯级系统长期优化调度、梯级系统日优化调度、径流预报、负 荷预报、非恒定流计算和梯级水电预演仿真。 本课题结合“可视化清江梯级水电联合调度仿真系统”项目展开，与该项目中非 恒定流计算模块联系紧密，非恒定流计算主要是研究各种日运行方式下任一时刻各断 面的水位流量关系以及它们其对航运的影响，在满足航运要求的条件下，制定梯级 华中科技大学硕士学位论文 系统的最佳调峰运行方式，提高电力系统的运行稳定性和梯级系统的效益。 1 2 明渠非恒定流的描述 一维明渠非恒定流指的是流场包含一个重力自由水面且流场沿着流动方向的尺 度远远大于其它两个方向的尺度的流动，这种流动的断面平均流速及水深都随着时间 过程而发生变化”1 。 在工程实践中经常要遇到明渠非恒定流问题，有时人们要求正确估计过水断面上 的水位过程与流量过程，这时解明渠非恒定流就具有直接的工程意义：而对于有些情 况，像确定污染物质在河流中的输移、扩散规律等问题时，虽然明渠非恒定流场本身 并不是直接关心的对象，但由于所要解决的问题同流动特性有紧密的联系，这时解明 渠非恒定流也不可避免。 一般水电站都会担任一定的调峰任务，在负荷日调节过程中，电站日出力变化 较大，这样就会造成下游河道水位陡涨陡落，水面坡降较大，对有航运要求的河道来 说，将会严重影响通航，因此必须研究这种情况下非恒定流的流态，以便更好地解决 发电和通航之间的矛盾，提高电站的运行效益。水电站负荷的变化、船闸的运转或拦 河闸的运用也会造成河道中的非恒定流，这些情况的水流状态对于水电站工程、船闸 工程或拦河闸工程等本身的规划、设计和运行方案的选择等都有很重要的影响，因此 也需要对相应的明渠非恒定流动进行研究。 当河道上游出现洪水时，为了给下游防洪提供依据，要迅速演算出下游可能的洪 水情况，这就是短期洪水预报，它给下游提供水位过程和洪量过程。在设计一个重要 的水利工程时，通常都要考虑拦河坝溃决的极端情况对下游的影响，这就要求进行溃 坝洪水的估计，预报溃坝可能造成的洪水情况，这也都属于非恒定流的研究范围。 现在在工程建设中人们越来越重视根据预测的工程效益来优化方案，再加上明渠 非恒定流动本身以及明渠非恒定流动中质量、热量及能量的输移、扩散规律的工程背 景的广泛性与重要性，明渠非恒定流动的预报问题已成为很多工程中一个带有共同性 的问题。有关明渠非恒定流的理论研究很多，但针对各种工程实际情况，实用的算法 并得到实践验证的方法还在探索中，由此可见，深入研究明渠非恒定流的计算方法是 十分必要也是很有意义的。 华中科技大学硕士学位论文 1 3 明渠非恒定流的研究现状 描述一维非恒定流的圣维南方程在1 8 7 1 年提出“1 ，明渠流曼宁公式“1 诞生也有 一百多年，在相当一段时期内，人们对非恒定流的研究都处在理论分析、现场观察及 室内水力学模型并举的局面，在解决工程实践问题时又以后两者为主，这些途径都存 在一定的局限性。圣维南方程为一阶拟线性双曲型偏微分方程组，只能在十分有限和 简化的条件下才能求其解析解，所得的解一般不能应用到实际问题当中。自2 0 世纪 5 0 年代电子计算机被广泛应用以来，数值计算方法的在解圣维南方程中得到快速的发 展6 1 。随着高速度、大容量计算机的问世和水利建设的不断发展，现在非恒定流计算 已广泛应用于防洪、除涝、灌溉、海涂围垦、航运、发电、防沙和环境保护等各项工 程。研究对象有天然河流、人工渠道、水库、湖泊、潮汐河口、海湾、城市下水道系 统及径流形成过程等。研究的问题也由一维发展到二维甚至三维。目前应用数值解法 大致可以分成三大类：有限差分法_ n 、有限元法m 卜n 6 1 、有限分析法“7 卜。“。 有限差分法是数值计算方法中最古老的方法，但到现在这种方法的应用仍最广 泛。在求解流体运动问题的数值方法中，有限差分仍占主导地位。有限差分法的基本 精神是把求解域划分成许多矩形网格，将因变量予以离散化，然后用离散化后应变量 的差商作为应变量导数的近似，即对水流运动微分方程中的导数项用差分式来逼近， 将原来的微分方程简化为代数方程求解。 有限元法从7 0 年代起开始应用于计算水力学中。有限元法的基本原理是将连续 介质划分成许多具有合适形状的微小单元，如三角形、四边形等等，在单元内选定一 些适当的点称为节点，把单元内因变量表达成该单元节点变量或其导数的差值函数， 借助于变分原来求加权剩余法，将控制微分方程转换成控制所有孤立单元的“有限元 方程”，最后将所有局部单元的方程汇集成总体的微分方程组或代数方程组求解。 有限分析法是1 9 7 7 年陈景仁提出。“。它不像有限差分法那样用差商近似控制方 程中的微商，也不像有限元法那样用插值近似来使控制方程通过积分形式在不同意义 下得到满足“”。该法的基本思想是把求解域划分成有限多个小的子域，在每个子域内 求基本方程的分析解，对非线性的微分方程，须将方程在子域范围内局部线性化，子 域内的分析解把该子域的一个内点上的变量与该子域边界上各点的参变量联系起来， 华中科技大学硕士学位论文 这种联系崩- 代数方程表示，由于每个子域边界上的点即为领域的内点，而每个内点 有个代数方程，因此全场所有内点的代数方程式组成了。个联立的代数方程组，求 解域边界的参变量为求解联立代数方程组时嵌入的边界条件。 有限元法用于边界复杂，水面宽阔的二维水流，如河口海湾的潮汐计算，有它一 定的优点，对于一维水流，使用有限元法和有限分析法就没有什么优点了。近2 0 年 来有限体积法”“在计算水力学中的优越性也受到国内外研究者的重视，对于非恒 定流的计算方法，国内外研究者仍在不断的进行改进。 1 4 本文的内容安排 本文在研究国内外有关非恒定流的研究现状和发展方向的基础上，结合梯级水电 站系统中水库非恒定流的特点，建立了梯级水电站系统中非恒定流的计算模型，采用 p r e i s s m a n n 隐格式的有限差分法和追赶法求得其数值解。从节约水资源和提高水电站 经济效益的角度出发，本文充分利用非恒定流预报的结果，建立了一个减少弃水的优 化模型，根据非恒定流的预测结果动态地从两方面进行调节，一方面调节作为反调节 水库下游水电站的日放水过程，另一方面调节各个断面的初始状态，通过计算机仿真 出一种满足优化模型的水流状态，使上下游电站日放水过程和非恒定流这两个相互制 约的条件，在满足航运和上游担任调峰负荷的要求下尽量达到充分利用水资源的目 的。本文最后通过清江的高坝洲一长江口河段以及三峡一葛洲坝一枝江河道的非恒定 流进行了模拟计算，通过分析计算结果证明本文提出的方案是可行的。全文各章内容 大致安排如下： 1 第一章介绍了国内外非恒定流的研究现状和发展趋势，课题的背景和意义： 2 第二章介绍了明渠非恒的数学模型圣维南方程建立的前提条件和推导过程： 3 第三章介绍了有限差分法的理论基础，比较了不同格式差分形式的优劣，重点介 绍了p r e i s s m a n n 隐格式离散微分方程的过程和相应的追赶法求解数值解的步骤，对三 种常见的边界条件处理也分别进行了详细的阐述； 4 在前面章节介绍的理论基础上，第四章讨论了梯级水库非恒定流的计算与调节， 从充分利用水资源和提高电站经济效益的角度出发，建立了减少弃水量的优化模型， 详细介绍了调节反调节库下游目放水过程和初始条件来满足目标函数的过程； 华中科技大学项士学位论文 5 第五章通过清江的高坝洲长江口的河段和长江三峡葛洲坝枝江两个 实际例子对本文针对梯级水库非恒定流提出的方案进行了验证和分析； 6 第六章对全文进行了总结，并提出了需要进一步研究的工作。 华中科技大学硕士学位论文 2 明渠非恒定流的数学模型 2 1 圣维南方程的推导 明渠非恒定流的数学处理是个重要而困难的问题，前人根据数学或试验观察 得到很多非恒定流的近似公式，其中被大家普遍接受的是圣维南方程，近百年来的水 力学研究也证明了它的正确性。尽管很多科学家企图修正、改进圣维南方程，但实质 上都没有很大的变化，有些方程可能要完善一些，当应用到实际中时都要对其进行简 化，简化结果得到的仍是圣维南方程。 天然河道的河形极为复杂，严格讲，天然河道是不存在一维水流运动的，因此圣 维南方程在描述非恒定流的过程中作了一些假设，虽然这些假设在物理意义上很难做 到，但从数值计算角度讲是可以近似做到的。圣维南方程的基本假设。”有以下几点： 1 波面是渐变的，即沿垂线压力呈静水压力分布和垂直加速度为小量； 2 不恒定流的摩阻损失与恒定流的摩阻损失无显著差别： 3 断面流速分布对波的传播无实质性的影响； 4 波的运动可认为是二维的，自由横断面中最终可能出现的水位差造成豹影响 可忽略不计； 5 渠道底部的平均坡度很小，s i na 可以用t a na 代替，c o sa 用1 代替，其中a 是渠 底与水平线的夹角。 浅水理论的基本假定是压力在垂线上按静水压力变化，即 | p = p g ( h z 、 z h( 2 1 ) 其中：z 一垂向坐标：h 一水深；p 一流体密度；g 一重力加速度；如图2 1 。下面 以平均速度v ( x ，) 和水深h 为因变量，纵向距离( 流程) x 和时间，为自变量，推导一维 方程，由于采用垂向平均值而除去了垂向尺度，坐标系统并非完全的直角坐标，其中 x 轴沿河槽地线，z 轴是铅直的。当假定河槽底坡底的余弦接近于1 时，可近似地认 为是直角坐标，下面就分别推导描述一维水流的圣维南方程。 华中科技大学硕士学位论文 1 连续方程 图2 1 一维水流x 、z 和v 符号说明 如图2 2 ，设基元水体在y 方向( 垂直于水流的水平坐标) 的宽度为l ，x 坐标 的位置在基元体的中点，基元体的高度h 是有限的，令z 方向的长度缸趋于零。有限 基元水体的质量守恒表述为：进入基元水体的净流量等于其体积的变化率。这样，质 量守恒可表示为： ( v - 是等坳一夏o h 等) - ( 矿+ 罢争( + 警等) = 百a h 缸 ( z z ) ( 流入)( 流出)( 体积变化) 将上式展开整理，除以缸，并取缸专o ，即可得到连续方程： 旦竺+ l f y h 1 ：0( 2 3 ) 图2 - 2一维水流的有限基元水体 2 动量方程 动量守恒方程表述为：进入基元体的净动量率加上作用于单元体上的诸力之和等 于基元体动量的变化率。流体的动量为质量流体与速度之积，从而有： 一 华中科技大学硕士学位论文 进入基元体的动量通量。p ( 砌) 一芸 矿( 砌) 1 了a x 流出基元体的动量通量2 p 矿( 砌) + 芸 矿( 砌) 】a f x ) = 肛= 偌如一z ) d z = 础2 ( 2 5 ) 尸 ，v 2 h - - 基c v 砌2 h ) 尹a x p v 2 哆内，争+ 艄缸 ， + 圭偌鼢2 一警丁a x 卜( 矗2 + 豢等) 】_ 彬，缸= 昙( p 砌) 血 。 v h ) + 昙( y 2 + 詈昙( 2 ) = 劝( s x - 曲) ( 2 8 ) 百o v + 矿詈+ g 罢- g ( 沪_ ) ( 2 9 ) 华中科技大学硕士学位论文 则需要更为普遍的方程。与前面的推导类似可得： 连续方程： 丝+ 旦( 蹦) ：0 o to x 、 7 动量方程： _a(va)+呈(矿2一)+厶=ga(sdto x ，一s ，)rj ( 2 1 0 ) ( 2 1 1 ) 其中厶对应于基元水体所受的静水压力项，作用于基元体端面上的压力为： ：k ( 矗一z ) b ( z ) d z ( 2 1 2 ) = 一吉鲁= 一昙c 一z ，b c z ，彪- 9 7 0 b c z ，龙瓦o h = 一鲥芸c z 。， 所以动量方程为： 丁a ( v a ) + 鲁( y 2 彳) + 鲥石o h = 鲥( 墨一墨) ( 2 1 4 ) 由于流量q ( x ，f ) = v ( x ，f ) a ( x ，t ，) ，水面线z = z ( x ，f ) = z 。o ) + h ( x ，f ) ，其中z a ( x ) 为河 槽底线则丝：一o h ， 丝：掣+ o h ：o h 一只。oto to x缸o xo x 。 2 流量和水位为因变量的方程 如果以o 、z 为因变量，则圣维南方程又可写为： 连续方程： 挈+ b _ o z ：0( 2 1 5 ) 戚优 动量方程为： 署+ 鲁( 争+ ( o 苏z + s s ) = 。 ( 2 1 6 ) 9 华中科技大学硕士学位论文 3 沿程有旁流q 流入或流出 实际应用中，有许多问题应考虑流体沿水流方向增加或减少，例如，地下水分散 地流入河道，支流集中地流入河道，降水或地表径流流入河道，灌溉系统引入的集中 出流及由渗流损失产生的分散出流。在这里，设旁流汇入河道时q 取下号，当旁流流 出河道时，q 取负号。考虑沿程有均匀旁流( 单位距离) 时，则连续方程变为： 望+ b 丝：口( 2 1 7 ) 融a 考虑旁流对动量方程式的影响时，还应附加一项由于旁流进出该水体带进的动量，一 般旁流附加项对动量方程式的影响很小，所以常常忽略不计，因此只考虑旁流对连续 方程的影响。 2 3 摩阻力公式 要使圣维南方程组中的两个方程定量地确定流量和水位的变化过程，方程中的 s ，必须确定。假定在非恒定流情况下摩阻损失与恒定流情况下差别不大，从河道水 力学可知，较常用的经验公式。”有如下几种： 曼宁公式 谢才公式 史觉克公式： s j ：而( n v ) 2 ：黑：掣 ( 2 18)ra2 rk 4 34 ，32一 aq 2 6 f - - c2 a 2 r ( 2 t 9 ) 驴斋斋 亿z 。， 三 其中记k ：二坠称为流量模数，r ：! 一水力半径，p 为湿周，a 为断面面积，n 一 ”月 糙率，c 一谢才系数，k ，一史觉克系数。 上述三个阻力公式中，由于引入了流量的平方，故不论水流方向如何鼠均为正号。 但对于潮水河道流向是变化的，阻力s ，的方向亦应随流向改变而改变，考虑流向的 1 0 华中科技大学硕士学位论文 变化，上面三个阻力公式应改成下面三个更为普遍适用的形式： 彭= 帮= 辔 眨z ， s ，= c q ：彳iq _ 盟2 r i ( 2 2 2 ) 曲= 署备 ( 2 z s ， 比较可知，这三个公式是非常相似的，其差别仅仅表现在各式中的系数上。在模型中 重要的不在于公式的选择而在于合理确定公式中的阻力参数，这些公式可以通过系数 拷拯同时信用存周一橙型中太寸存建奇檬利时引用的县寻宁绎赊公式。 华中科技大学硕士学位论文 3 数值解法 苹维南方程为双曲线性偏微分方程，存一股情况f ，目前无法求出年维南方程组 的解析解，只能采用数值解法或简化解法求其近似解。其中简化解法因其精度1 i 高而 较少使用，数值解法般分为婀步，第。步将其离散化，目j 把微分力程连续的定解离 散到定解域i t 的一些网格点卜，把偏微分方程转化为一组代数方程；第二步求解这组 代数方程。 3 1 有限差分法 何限差分法的基本思路是把描述连续变量( 如流量、过水面积、水f 证等) 的微 分力程，在讨论域内化为有限差分力程( 通常为代数打程) 求近似解的办法，也就是在 有限个网格节点卜求微分方程近似解的一种方法“。 如图3 一l 所示，取水平轴表示距离x ，垂直轴表示时间，用一系列平行r x 轴和 ，轴的直线将x t 平面划分成许多小矩形，称为嘲格。平行于，轴的垂直线表示1 i 同 的断丽位置，平 j ：于x 轴的水平线表示不同的时刻。平 j 二线之间的距离x 、，分别 称为计算的距离步长和时臼j 步长，定义为a x = 工，一工，_ l 、a t = 一f 一j 和n 为网格 的编号，称为序号，并以次表示网格函数在x 一，甲面卜- 的位置。各氲线之间的交点 称为甘点，节点处函数值今后记为9 7 、z ? ，其中下标表示断而位冒，上标表示时间。 3a ， ， 图3 1 网格和节点示意图 数值求解的任务就是根据“叫刻各节点卜已知的起始函数值( 即初始条件) 及给 华中科技大学硕士学位论文 定的上、下边界条件，计算f i 时刻各节点上的函数值。根据计算得到的 时刻值又可 计算下一时亥l j t 2 各节点上的函数值，依次类推，就可求得网格图上各节点的函数值。 由此可知，在计算过程中除了边界条件外，总是把前一时刻各节点上的函数值作为己 知，而把下一时刻各节点上的函数值作为未知值。 由于所有计算均在节点上进行，故网格的划分十分重要，若时问步长相同，距离 步长也相同，这种网格称为均匀网格。它的优点是计算方便，程序容易编制，要求机 器容量小，但缺点是： 1 断面位置受均匀分格限制，往往不是实际工程中要求计算的断面位置； 2 不能很好地根据河道断面变化情况来划分河段； 3 不能根据水力要素变率的大小来划分河段； 例如当在接近上边界断面处水力要素变化比较剧烈时，断面位置应设得密些，即) c 小些，而靠近下边界条件水力要素变化较缓慢，因此衄可取得较大。同样当所给的 上边界条件水力要素变化剧烈时，出可取得小些，即时间间隔短些，而水力要素变化 比较缓慢时，甜可取得大些，所以在解决天然河道非恒定流的计算中，往往较多采用 不均匀网格，特别是缸的不均匀风格，网格的划分除了考虑上述因素外，还应考虑 计算稳定性和精度要求。 明渠非恒定流的有限差分格式按照其推导的出发方程不同可以分为特征差分格 式。”、矩阵方程差分格式。”和直接差分格式。 所谓特征方程差分格式，就是先将偏微分形式的基本方程化为沿两条特征线上具 有全微分形式的特征方程，然后再对特征线方程和特征方程进行差分化得到一组差分 方程，它包括自由格式和固定格式两种情形。自由格式是以离散的初始条件和边界条 件作为出发点，通过求特征线的交点确定x t 平面上的计算点，同时确定计算点上 的流量g 和水位z 值，并依次作为新的出发点逐步往后时刻推进计算。这种格式的缺 点是不能得到某一断面的流量和水位的时间过程，也不能直接得到某一时刻流量和水 位沿河道的分布情况，因此对计算成果的后处理工作量比较大。固定格式则是通过规 l ， 华中科技大学硕士学位论文 定计筲时刻和计算断面预先确定x l 平面上的计算点，通过计算点向前时刻作特征 线，利j j 日0时刻的l 知 条件求出汁算点 ：的流量和水位，依次逐时段进行，该格式 在明渠非，晦定流的计算一i ，比较成助。 矩阵方程差分格，足先将基奉方程化为甜。，阵方程，然后将适合丁这一方程的各 种格式，包折特征中心差分格式、特自e 偏心差分格式等，推广到矩融：方程中然后得到 组等分力程，称为矩阵蓐分格式利用矩阵乃程莽分格式可以由已知的条件直接求 卜一刚刻沿计算河段的流晕：和水位。在矩阵方程差分格式t h 一阶精度的特征偏心 差分格式也称为逆风格式具有很强的扩散作片j ，它的稳定性比较好。 直接差分格式可以分成最示和隐式两人类，每一类中又包括许多不同形式的差 分格式。其中显格式包括菱形格式、跳点格式、l a x - - w e n d r o f f 格式1 等，它上要是 指任意节点卜待求冈，叟氐在新的时间层( 如在n + l 层) 的值可以通过早先时问层( 如 n 层n 一1 层等) 上变量值显式解出，由丁二早先时间层的变量值都是已知的，所以， 当时h j 向前推进时，新时问层上的待求变量值可逐点求出。丽隐式格式则小同，新时 j 日j 层上待求量的值不可能由甲先时闸层的已知值求出，还必须用同时问层( 新时问层) 相邻节点值( 未知的) 作为信息，这样，常常在一个代数方程中包含着几个未知数， 凶此，新时f 日j 层的待求变量的值是不町以用显式解出来的，往往需要借助于迭代求近 似解的方法爿能获得解。“。，隐式差分格式丰要包括p r e i s s m a n n 格式、v a s i l i e u 格式、 c r a n k - - n i c o | s o n 格式等。”。隐式格式的最大优点是对问步长可以取较大，且稳定性 好，对急缓，变化问题都适用，缺点是计算上作量大。显式格式则不同，时间步长不能 选得过大，稳定性限制较强，一般显示格式都是条件稳定的，但具有计算简单，工作 量小等优点。在实际问题中，选用哪一种格式要根据具体情况综合考虑决定。 3 2p r e is s m a n n 隐式格式 这是个四点偏心隐式格式，是酱菜上曼( p r e i s s m a n n ) 1 9 6 1 年提出，它土要是 对变量及其一阶偏微商在柏邻点和相邻时m 层采用加权平均进行离敞，即对时间， 的偏微商分别墩点，和，+ l 上的差商半均值，刈+ 审问工的偏微上商则是分别墩n ，和 ( n + 1 ) ，时问层的差商加权平均值，变量厂采用同一网格周围四个相邻点的加权半 华中科技大学硕士学位论文 均值进行逼近，在这里详细介绍用p r e i s s m a n n 格式进行求解的过程。 1 把微分方程离散为矩阵方程组 方程组： 丝+ 一io _ q ：0 ( 3 1 ) 詈+ 昙 譬 + 鲥 豢+ k i z g l ，、= 。 。， 首先在x f 平面上，将计算域离散成如下图的矩形网格，计算万= 厂( x ，。) ，- ，= l ，2 ，门2 1 ，2 n 。用( ，= 1 ，2 ，j ，疗= 1 ，2 n 。) 作为微分方程的解厂( x ，) 的近 似， 对其中任一网格如图3 - - 3 ，因变量及其导数的差分形式为： 厂( m ) ：厂( x ，f ) = 导( 塌+ 万+ t ) + 半( 最。+ 疗) ( 3 3 ) 要z 曰笔+ ( 1 帼盟a x a ， 缸缸 、7 望。垡二墨立笪：二笠 t ( 3 5 ) 幽3 - - 2 计算域禹散网格图3 - - 3p r e is s r l t a r l f l 网格布置 式中，0 为加权系数0 s 0 1 根据( 3 3 ) 、( 3 4 ) 和( 3 5 ) 的差分格式，连续方程 和动量方程可化为： c ( z ? ，z 二，饼，鲸，z ? “，z 剁，g “，彰n + + 1 ) = 0 ( 3 6 ) d ( 刁，乃。，鲜，鲸。，刁“，z 爿，饼“，纠n + + 1 ) = 0 ( 3 7 ) 当a t 充分小时，“与刀较接近，如果记，”1 = 厂“+ 4 厂，并略去上标，f = 厂”则 华中科技大学硕士学位论文 f ( x , t ) 2 i 0 ( 删+ ) + 丁1 - 0 ( 片l + l o ) 望：口型二： 8 xx 可 = 西 州卅华= ，篆 ；( 够一m ) + 去一，) ( 3 8 ) ：臼( 竺! ：! 二堑) + 丛王( 3 9 ) a x ，a x ， 丛。二墨! 互：! 二笪：笠：! 二竺 ( 3 1 0 ) 2 a t2 a t 利用( 3 8 ) 、( 3 9 ) 和( 3 1 0 ) 的差分式，可得连续方程的差分形式： 竺! ：! 竺! 。 2 ， o ( z x b ，+ 舳+ i ) + ( b j + b ，+ 1 ) 0 二二_ 一r 动量方程的差分形式为 g + 1 a q ， 缸， + 坠二鸟：o( 3 1 1 ) 血 。 堕2 a 塑t + 寿c 鬻一警，+ 詈a xr 警一等， 缸爿x i + 鲋川一：+ 鲋， 。 一鼻i； 1 + 譬( 鲋“+ 鲋，) + 詈( 以+ 一专( a z 一一z ，) + 寿( ，一z 铷 ( 3 2 ) 拦2 【堕些篆卷a k 紫堑型+ 墨螋2 【 ( k i l +川) 2 爿0 q 二j q 二 ( k - ) 2 利用下面的关系式将( 3 1 1 ) 和( 3 1 2 ) 两式线性化 + 笔k 掣m( ：) 2 。 再1 2 南2 扣争 ( 3 1 3 ) j一2磊1(k；+akj)2 z 寿”z 争 ( 即2 ( 1 + 管) 2 2 。巧 “7 ( q ? + g ) 2z ( 研) 2 + 2 q ，n g ( g + a q j ) i 鲜+ a q j 睁饼i 饼i + 2 i 彰l q j ( 3 1 5 ) ( 3 1 6 ) da?dk“db“ 幽t 。茁2 i 2 b j 2 i 丝i 2 i 2 | 柚j 。云a 2 j 在线性化的过程中，略去增量的乘积项，例如，爿q ，等，线性化后的式( 3 11 ) 可以写成： 华中科技大学硕士学位论文 4 1 ，a q ，+ 占1 ，z - i - c 1 a q ，+ 1 + d 1 j a z ，+ l = 止、j 3 l ，j 其中各系数表达式如下： 铲一硒4 0 丽a t 一，小筹等一= 菇 n ，小鬻筹确，一锄c 盱研，。7 x ( b ? + 口1 1 ) 2c 亿备l”r ( b ；+ b 盖1 ) 、圹 一 对于梯形断面罢：p ，即为边坡系数。 线性化后的式( 3 1 2 ) 可以写成： 一2 ，q ，+ b 2 j z j + c 2 j q + l + d 2 j z + l = e 2 ( 3 1 8 ) 其中的系黼铲- 一i 4 0 a tl 虿q t ) + 2 船目筹 = 芸【警叫叫) + g ( 碥n 叫舛】+ g 触鬻孵一筹】 c = j = l + 4 峨0 a t q ；，+ l + 2 船口哿 驴争萼争吲钳钞孵蚓引删，裂一器】 = 争警+ 等叫锅州抄晖_ 舻酬警+ 等等， 如果乙、q j 、z 川、易。已知，因函数a ，、k ，、日是水位z ，的单值函数，而且其 形式( 多项式表示或用表格表示的函数) 为已知，所以在任一时刻f 。就可计算出上式 中的系数爿l ，、a 2 、b 。、b 2 j 、c i 、c 2 、d l 、d 2 、e l 、e 2 j 。 有上面两式对解出z ，、g 、z 川、g + 。而言，它们还不够，因为相应于四个 未知数，仅有两个方程可利用，但上式可针对任一计算点( ，_ ，+ 1 ) ，如在模型里有j 个 壬卜笪点就能对2 j 个蛮量写m2 ( ，一1 ) 个汶样的方稗。网一凶然有两个沩界条件可 华中科技大学硕士学位论文 供利用，实际上就获得包含2 j 个未知数( a z ，q ，) 的2 ( j 一1 ) + 2 = 2 j 个代数方 程组，因此这个系统对任意时i 瑚步长出都可求解。 应当指出，仅当i a f f f 1 且充分小时，原差分方程的线性化才f 确，当i a f f 1 较 大时，就不能认为线性化后的方程组是非线性方程组的良好近似，不过，如对一个选 定的出， 具有与，同样的数量级，那么，就总是或几乎是能采用减小时间步长f 来减小比值f 4 厂卅，从而使j 4 ， 卅充分小。 2 求解矩阵方程组 对于计算周期中每一个时间步长的所有计算点，必须解出前面的( 3 1 7 ) 和( 3 1 8 ) 两个方程式，它们构成一个线性的代数方程系统，如果边界条件对z 和a q 也是线性 化，那么任何标准的求解法都是可以采用的。这一求解过程是程序中最耗费时间的部 分，在这里详细介绍用于解线性系统的追赶法( 也称双扫描法) 。对于( 3 1 7 ) 和3 1 8 ) 两式组成的线性系统，假定对于点_ ，存在线性关系： 厶q i = f i 监j + g i ( 3 1 9 ) 就能证明在_ ，+ 1 点也存在线性关系：9 + = + a z + 。+ g 川 事实上，将( 3 1 9 ) 代入( 3 1 7 ) 可得 a 1 j ( fa z j + g j ) + b i z j + c 1 j g + l 十d l ，z 卅= e l j ( 3 - 2 0 ) a z j = 一赤q ，r 击心s + l l 一巧e i i - - 而a i j g z t ， ( 3 2 1 ) 式可写成： a z j = 日j g + l + ，z ，“+ i ，j ( 3 2 2 ) 其中：一= 磊，一们o 峨l ，以= 巧e l i - - 而a i j g j 将( 3 1 9 ) 和( 3 2 2 ) 两式代入( 3 1 8 ) 得。 1 8 华中科技大学硕士学位论文 爿2 ，【一( 日，q ，+ i 十，a z ，+ i + j ，) + g ，i + 君2 ，( 抒，q + l + ，j z + 1 + ，) + c 2 ，q ，+ i + d 2 ，a z = e 2 ， 整理后得： 啄一一筹黼+ 搿群z 。， 若记：a2 4 ，c + 岛，则：c “= 一三 苦，q “= 墨群式 ( 3 2 3 ) 即为：纺+ = + i a z j + l + g s + l ，这就证明了在后继点也存在类似于( 3 1 9 ) 的关系，并且此关系中的系数巧_ rg j 一可由系数a f 、岛、q 、b 、e i = l ，2 ) 和 吩g ，通过四则运算直接计算出。这些递推关系启示了求a z ，、a q ，及 刁“= 刁+ z ! l z 、g “= 纺+ g 的计算方法，当上边界条件线性化后，可得关系 a q , = 曩z l + g l ，由递推关系，可依次求出关系式：z ，= h q ，+ l + j j z 川+ j ， 9 + t = c + i 乙+ ，+ q “( ，。1 ，2 j - - 1 ) 这一过程形象地称为追a 在下边界处，如果计算得z ，代入上面所得关系式中，可 以得出q ，求z ，、q ，( _ ，= d l 2 ，1 ) 这个过程称为赶。由以上过程可以 看出，用追赶法求解系统所需得的初等运算的数目与结点数j 成正比，而矩阵求逆的 标准方法所需的运算次数与，3 成正比，由此可见追赶法的优越性。 3 p r e i s s m a n n 格式的稳定条件和精度 按照不同格式将微分方程离散成代数方程，用差分格式代替微分格式，将求微分 方程的解转化为求代数方程组的离散解，这只是求解问题的第一步，在得到代数方程 后，还必须考虑以下一些问题： ( 1 ) 用什么方法解代数方程组，这是保证其一定可解的必要工具； ( 2 ) 收敛性：若求得的差分方程精确解为，在所考虑区域内的任意点( j ，n ) 上， 当离散的网格步长( z k x ，t ) 趋于零时，差分方程精确解是否趋近于微分方程的解， 华中科技大学硕士学位论文 这就是收敛性问题； ( 3 ) 稳定性：求代数方程解是各层逐点计算的，在某一时问层上引入的误差( 包含 舍入误差问题和初始误差) 必然会影响到以后各层，如果这种误差随时问的推进而不 断累积，最终误差分方程的解被引入的误差所“淹没”，则这种数值解是不稳定的。 第一个问题讨论的是代数方程解的方法，对相同的微分方程采用不同的数值格式，解 代数方程的方法是有所不同的，收敛性问题是研究代数方程组解精确解与微分方程解 的关系问题，而稳定性问题则是研究代数方程组近似解与代数方程组精确解的关系问 题。 通过p r e i s s m a n n 格式得到的差分方程，对于系数0 的任意值，具有对a x 的一阶 精度，而0 = o 5 时，则具有二阶精度，当0 口0 5 ，它总是不稳定的；当0 5 饼+ g = f ( z ? + z 抄( 刃) + 砉心l ( 3 2 7 ) 比较后可取互= 誓t ：f ，g 。= 厂( z ? ) 一研，显然应有g = 。，但是由于具有截断和 2 1 华中科技大学硕士学位论文 其它误差，它不总是i # 女- d b 零，因此g 是一个修正因子。当作下游边界条件时，以z 为未知量，则： 又有： 则 q = q ? + a , 2 ，= q ? + f a z ，+ g ， ( 3 2 8 ) q a 。= 厂( z ：+ 哆) z 厂( z ；) + 羞j 。彪，( 3 2 9 ) 嵋= 警 一荔b ( 3 3 0 ) 在用追赶法解此方程的过程中( 假设j 为河段的断面数) ，对于任意时段， 当上边 界条件线性化后，在追的过程中求得h ，、j ，、j ，( j 2 l - j1 ) ，、g ( j 一2 j ) ， 再线性化下边界条件，然后利用式( 3 1