（水力学及河流动力学专业论文）河道洪水数值模拟理论的研究与应用.pdf

华北水利水屯学院硕十学位论文 河道洪水数值模拟理论的研究与应用 摘要 论文结合“黄河下游二级悬河问题数值模拟研究”与“甬江清水浦 码头防洪影响评价”两个工程实例，对河道洪水数值模拟的理论与计算方法 进行了较系统的研究。 首先，在洪水数值模拟的基本理论方面进行了三项研究：i ) 通过水力学 方法分析得出洪水数值模拟中四种水流边界条件的合理性与洪水流场流态相 关，并通过一维水流模型验证了在上下游边界条件都给流量关系时，如果不 附加边界条件，便会影响模拟流场的唯一性。2 ) 结合有限差分法和控制体积 法，提出了一维水流控制方程在交错网格系下的一种新的离散格式：时间项 采用偏心差分，空间项采用控制体积法离散，对流项采用迎风差分格式；3 ) 建立了一种无重叠网格、能同时保证质量、能量和动量守恒的新型一、二维 耦合水流模型，同时给出模型中相关问题的解决方法。 其次，在河流数值模拟的前后处理和程序设计上做了两方面的研究工作： 1 ) 采用a u t o c a d 和m a t l a b 相结合的方法进行河流数值模拟的前后处理，解决 了混合编程中数据传输问题，可以大大提高河流数值模拟前处理工作效率， 实现了数值模拟成果快捷方便的查询和编辑，以及数值模拟成果的直观可视 化。2 ) 应用面向对象程序设计方法进行河流数值模拟程序设计，探讨了团队 开发日益复杂的河流模拟程序的解决方案。 再次，通过两个实际工程问题的河流数值模拟实践，取得一些基本认识： 应用平面二维水沙数学模型研究黄河下游“二级悬河”问题，分析模拟 成果得到了二级悬差的临界控制条件，不同治理方案模拟成果分析表明，采 取“增强主槽过流能力”和“滩地治理”相结合的方式对“二级悬河”治理， 效果较好。 应用一、二维耦合水流模型研究了甬江清水浦码头的防洪影响，较好模 拟了码头修建前后的流场变化，解决了近海河口长河段工程数值模拟研究中， 精细局部流场的嵌套模拟问题。 关键词：数值模拟，离散格式，水流耦合模型，前后处理，二级悬河， 潮汐河段 a b s t r a c t t h es t u d ya n da p p l i c a t i o no fr i v e rf l o o dm a t h e m a t i c a lm o d e l i nt h ep a p e r ，s o m es t u d yw o r k sa r es h o w na n da p p l i e do nt w op r o j e c t s ，n u m e r i c a l s i m u l a t i o nr e s e a r c ho f s u s p e n d e dr i v e r ”i nl o w e ry e l l o wr i v e ra n dn u m e r i c a ls i m u l a t i o n r e s e a r c ho ff l o o de v a l u a t i o nf o rq i n g s h u i p ud o c ki ny o n g j i a n gr i v e r ，b ys t u d y i n gn u m e r i c a l s i m u l a t i o nt h e o r ya n dc a l c u l a t i o nm e t h o do nr i v e tf l o o ds i m u l a t i o n f i r s t l y , t h r e es t u d yw o r k sh a v eb e e nd o n eo nr i v e rf l o o dn u m e r i c a ls i m u l a t i o n ，w h i c ha r e e x p r e s s e da sf o l l o w ： 1 b va n a l y z i n gi nh y d r a u l i cm e t h o d ，t h ec o r r e c t n e s so ff o u rk i n d so fh y d r a u l i cb o u n d a r yi s m l a t e dt of l o wp a f l e m i ft h e r ei sn oa d d i t i o n a lc o n d i c t i o n ，t h es t e a d yf l o ws i m u l a t i o n r e s u l t sw o u l db ee f f e c t e db yt h ei n i t i a lc o n d i c t i o nw h e nt h ed i s c h a r g ec o n d i c t i o n sa r es e t i nt h eb o t hb o u n d a r i e s 2 c o m b i n e dw i t hf i n i t ed i f f e r e n e em e t h o da n df i n i t ev o l u m em e t h o d ，an e wd i s c m tf o r m a t h a sb e e nf o u n df o ro n ed i m e n s i o nf l o ws i m u l a t i o n i nt h en e wd i s c r e tf o r m a t ，t h et i m e r p a r t i a ld e r i v a t i v ei ss o l u t e db yf i n i t ed i f f e r e n c ew i t hw e i g h t i n gc o e f f i c i e n t ，b u tt h es p a c e p a r t i a ld e r i v a t i v ei ss o l u t e db yf i n i t ev o l u m em e t h o da n dc o n v e n t i o ni t e mi ss o l u t e db y u p w i n df o r m a t 3 an e w1 d 2 dc o u p l ef l o wm o d e l ，w h i c hn e e d n ts u p e r p o s eg r i d sa n dc a nm a k et h em a s s ， e n e r g ya n dm o m e n t u mb ec o n s e r v e d ，h a sb e e nf o u n d s o l u t i o n sf o rs o m ep r o b l e m si n c o u p l em o d e l h a v eb e e ng i v e no u t s e c o n d l y , t w os t u d yw o r k sh a v eb e e nd o n eo nn u m e r i c ls i m u l a t i o np r o g r a md e s i g na n d p r e - p o s td e a l tw a y ，w h i c ha r ee x p r e s s e d a sf o l l o w ： 1 c o m b i n e da u t o c a da n dm a u a b as o f t w a r eh a sb e e nd o n ef o rd e a l i n gt h ep r e - w o r ka n d p o s t - w o r ki nr i v e rf l o o dn u m e r i c a ls i m u l a t i o n n l es o l u t i o nf o rt r a n s m i t t i n gd a t ab e t w e e n m a t l a ba n dv i s u a lc + + h a sb e e ng i v e no u t t h es o f t w a r eh a sb e e nf o u n dt h a ti tc o u l d p e r f e c t l yi m p r o v ee f f i c i e n c yo fp r e - w o r k ，e a s i l yl o o ku pa n de d i ts i m u l a t i o nr e s u l t s ，a n d v i s u a l i z es i m u l a t i o nr e s u l t s 2 n l ep r o g r a md e s i g nm e t h o do f o b j e c t - o r i e n t e dp r o g r a m m i n gh a sb e e na p p l i e di nt h er i v e r f l o o dn m n e r i c a ls i m u l a t i o 玑a n di sab e t t e rs o l u t i o nf o rd e v e l o p i n gm o r ea n dm o r e c o m p l e xn u m e f i c a ls i m u l a t i o np r o g r a mb yt e a mw o r k t h i r d l y , s o m ei d e a sh a v eb e e ng a i n e di nt h er i v e rf l o o dn u m e r i c a ls i m u l a t i o np r a c t i c eo f t w op r o j e c t s b ya p p l y i n g2 df l o wa n ds a n dm a t h e m a t i cm o d e lt or e s e a r c h “s u s p e n d e df i v e r i nl o w e r y e l l o wr i v e r ，ac r i t i c a lc o n d i c t i o no fs u s p e n d i n gh e i g h tv e r t i c a l l yf r o mt h eb o r o m l a n dn e a rb y r i v e r b a n kt oc h a n n e la l i g n m e n tt o ph a sb e e nf o u n d b a s e do nt h es i m u l a t i o nr e s u l t s t h e c o m p o u n dm e t h o do fe n l a r g i n gb a n k f u ld i s c h a r g ea n dd o i n gt r a i nt ob o t t o m l a n dw o u l dg e ta b e t t e rr e s u l tf o rt r a i n i n gs u s p e n d e dr i v e r b ya p p l y i n g1 d - 2 dc o u p l em o d e lt od of l o o de v a l u a t i o nf o rq i n g s h u i p ud o c ki ny o n g j i a n g r i v e r ，p r e c i s i o nf l o wf i e l d sh a v eb e e ng a i n e dw h e nt h ed o c kb u i l tb e f o r ea n da f t e r as o l u t i o n f o re l a b o r a t ef i e l ds i m u l a t i o ni na l o n gr e a c hi nt i d er i v e rh a sb e e ng i v e no u t k e yw o r d s ：n u m e r i c a ls i m u l a t i o n 。d i s c r e t ef o r m a lc o u p l ef l o wm o d e p r e p o s tw o r k s u s p e n d e dr i v e r , t i d a lr i v e r 华北水利水电学院碗上学位论文 卜垂线水深，m ； f l _ 一断面平均的水深，m ； 一水位，m ； “x 方向的垂线平均流速分量，m s ； v y 方向的垂线平均流速分量，m s ； 卜垂线平均含沙量，k g m 3 ； 一垂线平均挟沙力，k g m 3 ； 岛一水流饱和挟沙力，k g m 3 ； s 一紊动运动粘性系数： g b ，。呵方向单宽推移质输沙率，k g s m g b 厂7 方向单宽推移质输沙率，k g s m 疗糙率系数； p 一泥沙干密度，k g m 3 ： 一垂线泥沙沉速，m s ； o 断面平均泥沙沉速，m s ； 岛；一悬移质引起的河床冲淤厚度，m ； z 。一河床总冲淤厚度，m ； 卜底坡： “r 一摩阻流速，m s ： c _ 一谢j 才系数，m 2 s ： 符号说明 a 一悬移质泥沙恢复饱和系数。 矗一泥沙粒径，m i l l ； 以一泥沙容重，k g m 3 ； ，一清水容重，k g m 3 ； ，1 一床沙干容重，k g m 3 ： u 断面平均流速，m s ： u 一泥沙起动流速，m s ； 9 流量，m 3 s ； 4 过流面积，m 2 ； 肛一水力半径，m ； 卜垂线水面坡降； 山一断面平均水面坡降； g _ 一重力加速度，9 8 1 n k g ； 妒一单宽流量分布函数 西一流速分布函数 函一水位分布函数 一水深分布函数 独立完成与诚信声明 本人的学位论文是在导师指导下独立撰写并完成 的，学位论文没有剽窃、抄袭等违反学术道德、学术 规范的侵权行为，否则本人愿意承担由此产生的一切 法律责任和法律后果，特此郑重声明。 学位论文作者( 签字) ：庳j 、蔻 芦硐年月廖日 华北水利水电学院硕上学位论文 1 1 洪水演进数值模拟概述 第一章绪论 洪水是1 种不呵抗拒的自然现象，给人们的生命财产带来了巨大的灾难。因此，人类 为了存，与洪水的斗争一刻也没停止过。在古代，由于人口稀少，经济不发达，洪水威 胁并不十分明显，但随着社会的高度发展和人口的急剧增长及人类对洪泛区的频繁开发， 洪水灾害呈现逐年增长的趋势，从而使人们对洪水的认识不断加深。但洪水的成因大大超 出人类的控制能力，人们对洪水发生的规律，也只能作出部分预测和预报i ”。 目前，人们对洪水的分析和预测，主要有两种方法：资料分析和洪水模拟【2 】。水文学 主要采用统计学方法对已发生洪水进行资料分析和对可能发生的洪水进行预测【3 j ，而水动 力学则主要采用水动力学的基本理论对洪水进行模拟来研究洪水传播时间、影响范围和影 响程序。 水动力学模拟一般有物理模拟与数值模拟两种类型。物理模拟即比尺模拟。该方法以 其直观性和物理概念明确一直受到水利工程界的重视，然而其模拟程度往往受模拟试验手 段和比尺关系的限制，费用比较高昂，试验周期相对较长【4 】。本课题论及的是后一种模拟 的方法，即河流工程的水动力数值模拟。它具有省时、高效、无比尺影响等优点，目前正 以强劲的势头发展，并在众多河流工程问题上被广泛应用，取得了良好的效果。 1 1 1 水动力数值模拟的发展 自2 0 世纪5 0 年代以来，随着计算机的出现和现代高速电子计算机的发展，数值模拟 和系统模拟在科学技术的许多领域中已广泛应用。在水利工程方面，因计算机的发展和应 用，诞生了一个新的分支学科计算水力学。这一学科的发展促进了求解流体力学方程 数学方法的发展，从而使流动现象的数值模拟达到了更高水平。2 0 世纪7 0 年代，在计算 水力学的基础上，河流泥沙工程领域的数值模拟开始发展起来，并以其省时、省力和灵活 方便等优势在研究水沙运动规律和河床变形等方面的研究中起着愈来愈重要的作用u 扣1 。 在计算机与水沙运动计算相结合的现代水动力数值模拟之前，2 0 世纪5 0 年代在中国、 美国、前苏联，以及西欧，一些国家曾经使用计算方法研究水库淤积、水库下游河段的冲刷 和潮汐河口等方面的问题忉。现代水动力泥沙数值模拟是在2 0 世纪7 0 年代后期才发展起 来的，比许多学科的数值模拟晚1 0 2 0 年。其主要原因有三：首先是人类对泥沙问题严 重性的认识不足，由此决定其在国民经济发展中的地位不高；再者是泥沙数学模型的发展 必须建立在水流数学模型的基础上，所以其发展必然是滞后于水流模型的发展且受到水流 模型发展的制约；另外一个原因是由于泥沙问题本身的复杂性，使得泥沙数学模型的发展 受到很大的约束i 。 2 0 世纪7 0 年代后期以来，许多国家在发展冲积河流数值模拟方面作了大量的工作。 竺苎查型查皇兰堕塑主兰些兰塞一 一维数学模型是最简单，也是发展最早、最完善的模型，它主要用于研究长时期内长河段 的水流及河床变形情况。一维数值模拟是以断面平均的河床、水流及泥沙因素作为研究对 象，只能给出各河段的平均流速及其冲淤情况f 7 l 。然而，在一些水利工程中，需要对河床 细部变形进行了解，这就需要二维甚至三维模型进行研究。 11 2 永动力数值模拟的特点 水动力数值模拟实质上是采用数值近似计算方法来模拟研究对象的自然现象和动力 地貌演变过程。 水动力数值模拟的理论基础是水力学、泥沙动力学和数值计算方法。在水力学和泥沙 动力学中，实际自然现象又都是非线性的，描述这些现象的微分方程也是非线性的，因此， 对这些微分方程很难求得解析解，一般只能求得近似解或进行数值模拟。在求解这些非线 性微分方程时，许多问题( 如稳定性、收敛性、误差估计等) 至今尚无合理的理论分析， 在数值模拟中，一般只能以线性理论为基础，参考一定的经验，结台一。定的分析，边试验 边修改；在试验和修改过程中，还要有丰富的经验和一定的技巧，并且以现场资料为依据， 最后才能得出可靠的解答心j 。 通常人们认为数值模拟仅是一种模拟手段，只能复演一些已知现象，不熊作为一种研 究新手段。事实卜，数值模拟不仅可以复演自然界中各种已知现象和过程，还能发现一些 以前没有被发现的新现象并对其进行预测和预报。如c a m p b e l l 和m u e l l e r ( 1 9 6 8 ) 通过数 值模拟发现了亚音速斜坡引起的分离现象1 9 ；后来，他们在风洞试验中也证实了这种现象。 又如，近年来在数值模拟中发现了水流以小角度穿越航道时的水流偏转和速度加大的现 象，这种情况与后来水池和物理模型试验所获得的结果相一致。因而，应该把数值模拟的 方法作为一种试验或研究及预测方法来考虑。 数值模拟与物理模型相比较，各有其优缺点【射。 数值模拟的优点： ( 1 ) 试验费用少。数值模拟的费用通常比物理模型的费用要少，而且物理模型的费 用还有逐渐提高的趋势，而数值模拟的费用却因计算机技术的发展而有减少的趋势。 ( 2 ) 试验速度快、试验周期短。研究人员可以在一天内进行几个方案的模拟计算， 还可以在几台计算机上同时进行几组方案的试验或者进行同一个方案的并行运算：但同样 的工作量在物理模型中却要花费几倍、几十倍的时间。 ( 3 ) 可以模拟多种因素相互作用的复杂物理过程。如可以模拟波浪、水( 潮) 流、 风、柯氏力等多种因素共同作用下多种泥沙运动形态及地形演变的复杂过程。 ( 4 ) 不占用大片模型试验场地。 ( 5 ) 可以完全控制流体的物理性质( 如密度、容重、粘度、含沙量等) ，也不存在试 验探头对水流扰动的影响。 ( 6 ) 不存在模型比尺的影响。 ( 7 ) 模型建成以后，可以长期使用。修改方案容易。 2 华北水利水电学院硕士学位论文 物理模型则恰好和数值模拟相反，数值模拟所具有的优点正是物理模型的缺点；而数 值模拟存在的无法模拟微分方程不能描绘的物理现象、结果不直观、物理概念不明确等缺 点，恰好是物理模型的优点。因此，物理模型和数值模拟作为两种不同的模拟手段，目前 在水力、海洋、环境工程中得到了快速发展和应用，任何一方想完全取代另一方都是不可 能的，两者只有互相补充，互相验证，取长补短，共同提高，才是正确的方法。最近发展 的复合模拟就是吸收了这两种模拟技术的优点而使其有机结合的一种新兴模拟技术。 1 1 3 水动力数学模型的分类 目前划分数学模型的方法有很多种，根据不同的划分标准数学模型也可以分为许多种 类【5 1 。根据所研究问题的维数来划分，可分为一维、二维、三维数学模型；根据对来水来 沙过程处理方法的不同来划分，可分为恒定水沙数学模型和非恒定水沙数学模型；根据模 型水流及泥沙计算方法划分，可分为耦合解及非耦合解模型【7 】。 另外，按照所模拟的泥沙运动状态进行分类，可分为仅模拟悬移质运动的悬移质模型， 仅模拟推移质的推移质模型，以及同时模拟悬移质和推移质运动的全沙模型。按照计算含 沙量的方法进行分类，又分为悬移质饱和输沙模型和非饱和输沙模型( 或不平衡输沙模 型) 。非饱和输沙模型通常仅限于模拟悬移质运动，这是因为推移质运动达到饱和输沙状 态速度比较快，可以采用饱和输沙法计算。所以在悬移质饱和或非饱和模型基础上加入推 移质运动模型，又可以分为全沙饱和输沙模型和全沙非饱和输沙模型 7 l 。 此外，按照各类模型不同的组合也可以分为不同的模型，如耦合解恒定饱和输沙模型、 耦合解恒定非饱和输沙模型、耦合解非恒定非饱和输沙模型、非耦台解恒定饱和输沙模型、 非耦合解非恒定饱和输沙模型及非耦合解非恒定非饱和输沙模型等等。 一般的河流水沙数学模型，为了简化计算，大多采用非耦合解数学模型。在工程实际 应用中具体采用什么模型，应对所要研究的实际问题进行分析，再确定采用哪一种数学模 犁。 1 1 4 一、二维水动力数值模拟的研究对象 ( 1 ) 一维水动力数值模拟的研究对象 一维水动力数值模拟是以断面平均水流和输运物质( 如泥沙、b o d 等) 作为研究对象， 研究断面水力要素和输运物质的时空变化情况。 如前所述，一维数学模型是发展最早、最完善的水动力模型，所以常用于研究洪水 和污染物的传播、沿流程的水位和流量变化以及长时期的河床变形。近年来，为了研究河 流工程对地区防洪与河流的影响，一维数学模型常与二维数学模型结合，为二维模型提供 一个比较合理的进出口边界条件。本文的第二个应用例子就属于该情况。 ( 2 ) 平面二维水动力数值模拟的研究对象h , 1 1 平面二维水动力数值模拟是以垂线平均的水流及泥沙因素作为研究对象，研究它们在 平面上的变化情况，主要用来进行平面流速分布和河床细部变化问题的研究。 华北水利水电学院硕j j 学位论文 在大部分河道中，水流流动具有以下特点：具有自由表面，为明渠水流；重力为水流 流动的主要驱动力，水流内部及水流与固体边界的摩阻力为水流能量的主要耗散力；水流 流速沿垂线近似均匀分布，不必考虑实际存在的对数或指数等形式的垂线流速分布；水平 运动尺度远大于垂向运动尺度，垂向流速及垂向加速度可忽略；水压力接近静压分布等。 在上述水流流动情况下，对于需要研究水流或河床细部变化的问题，水沙因素可以沿垂线 积分，采用平面j 二维水动力数值模拟对水沙运动情况进行模拟研究。 目前，平面二维水动力数值模拟研究的问题主要有：水工及河工建筑物附近的流速分 布与河床变形，如坝区引航道淤积、桥渡的冲刷、浅滩挖槽的回淤；自然河流泥沙成型堆 积体消长、运动所引起的河床变形；分汉河段主要支汉的交替发展；交错边滩向下游运行； 弯道冲淤；浅滩演变；河道流量及洪水预报；蓄滞洪区洪水演进等l ” 1 4 1 。 1 2 论文主要理论基础和研究任务 1 2 1 主要理论基础 水动力泥沙数值模拟是一门综合性很强的交叉学科，涉及到水力学、流体力学、河流 动力学、数理方程、数值分析、矩阵论、计算机科学等学科，数学模型的选择、建立和数 据计算及其整理、分析，都要以此为基础。 ( 1 ) 数理方程、数值分析和矩阵论 这三门学科是基础学科，为水沙数值模拟提供了必需的数学工具，包括水沙数学模型 控制方程的推理和离散格式，各种计算格式的稳定性、相容性与收敛性，以及代数方程组 的求解方法等。 ( 2 ) 水力学、流体力学与河流动力学 这三门课程为专业课程，为水沙数值模拟提供了必须的专业知识，如各种过流建筑物 过流能力的计算方法、水沙基本控制方程以及泥沙输运计算理论。 ( 3 ) 计算机科学 计算机科学足水沙数值模拟搭建的平台。地形数据的读取与插值、平面网格的生成、 各种差分方法的程序实现、各种水力要素分布的场图及其动态演示都是以计算机图形学、 数据可视化、程序设计为基础的，与计算机科学密不可分。 1 2 2 主要研究任务 本论文圭要研究任务分两个方面：一是建立平面二维水沙数学模型，以黄河下游“二 级悬河”问题研究为应用实例，研究“二级悬河”对黄河下游防洪和大堤安全的影响；m - - 是建立一、二二维嵌套水流模型，以“甬江清水浦码头的防洪评价研究”为应用实例，研究 拟建码头对甬江沿岸防洪的影响。其中的主要研究任务和研究重点有以下几个方面； ( 1 ) 数学模型的选择 在实际河流工程数值模拟中，由于研究问题侧重点的差异、计算区域大小的差别、以 4 华北水利水电学院硕士学位论文 及河流形态与计算精度要求的影响，需要根据具体情况选择合适的数学模型。 黄河下游河道滩面宽阔，水流含沙量高，洪水期水域范围变化较大。为了能合理模拟 “二级悬河”对防洪的影响，宜采用平面二维水沙数学模型。 甬江是宁波市的一条海口河流，河道宽度为5 0 0 m 左右，航道占断面的绝大部分过水 面积，滩地面积比较小，而且大部分河段已经渠化，因此一维模型可以模拟洪水在该河段 的传播情况。但一维模型不能模拟码头修建后，其周围水流结构的变化以及这些变化对周 围河床的影响，因此需要采用二维或三维数学模型。鉴于二维数学模型比三维数学模型简 单，而且基本上能解答上述问题，因此本论文采用一、二维嵌套数学模型来研究甬江清水 浦码头的防洪影响。 ( 2 ) 一、二维模型的连接方法和计算格式 一、二维耦合模型已有一些研究成果：导师孙东坡教授和彭文启博士于2 0 0 0 年建立 了一、二维耦合模型研究近海尾闾河段1 1 5 1 ；熊涛等利用特征线法对一、二维模型进行了 耦合【l6 】；赖易军等利用重叠一投影法建立了耦合模型i l7 】：程杭平等针对热电厂之间区域 建立了热污染的一、二维耦合模型1 1 8 ；e d b a t e s 等通过将h e c r a s 与t e l e m a c 一2 d 等 模型耦合，得到了一个组合模型【l 9 j 。 本文通过研究，采用合适的、二维控制方程离散格式和求鳃程式，通过连接一、二 维代数方程的循环系数，建立满足质量、能量和动量守恒的一、二维显隐交替的耦合水流 模型。 ( 3 ) 河流模拟前后处理和面向对象程序设计 河流模拟需要的资料包括地形地物、过流建筑物和水文资料，数据量一般比较大，因 此数值模拟的前处理工作一般要占到整个模拟工作的一半以上。为了提高河流模拟的效 率，进行河流数值模拟前处理研究是非常必要的。此外，为了减少计算程序对数据正确性 判断，提高程序运行效率，也是有必要在前处理过程对数据进行较核和加工处理。 随着大河网和大区域河流模拟的进展，不管是一维模拟还是二维模拟，模拟计算输出 的数据都是巨大的。而我们所关心的数据成果往往并不是很多，而且有很多数据并不是数 值模拟输出的直接数据，而是综合的数据成果，因此需要研究从海量数据提取和总结有用 数据的方法。此外，从河流模拟成果的应用看，这些模拟成果不仅用于理论分析，还有用 于生产实践，因此模拟成果不仅需要能直观表现，还要便于设计和施工部门使用。 为了实现河流模拟的前后处理，需要编写相关的辅助程序或软件。在过去，f o r t r a n 语占凭借强大数值计算能力占据数值模拟的大部分领域，但随着数值模拟越来越复杂，需 要处理的事务剧增，新的数模程序一般需要团队开发才能完成，因此结构化程序设计的 f o r t r a n 语占面临新的挑战。而今，由于c + + 标准化协会的努力，c + + 的计算速度可以与 f o r t r a n 相当，而且c + + 是一个面向对象的程序设计+ 适合团队开发和复杂问题处理。 本文结合自己三年多的数值模拟实践，研究和提出了基于a u t o c a d 和m a t l a b 的 河流模拟前后处理方法，采用面向对象程序设计方法，编写了将基础数据转为数学模型计 华北水利水电学院硕上学位论文 算数据的前处理程序，开发了数值模拟成果可视化和c a d 化的相关工具，较好地满足当 前河流模型前后处理的一般要求。 ( 4 ) 工程实例的计算与分析 本论文进行了两个实例的应用研究。 第一个实例是应用平面二：维水沙数学模型和前后处理辅助软件，对黄河下游“二：级悬 河”问题进行数值模拟研究。本论文选择黄河下游“二级悬河”典型河段( 东坝头至高村) 作为研究河段，研究该河段在“二级悬河”的三种可能发展趋势( 维持现状、加剧发展和 治理) 条件下发生典型洪水时，洪水的演进、漫滩及河势变化状况，分析研究三种发展状 况的滩槽分流比、滩槽流速分布，顺堤流速，堤沟河流量及滩槽冲淤状况。最后对“二级 悬河”对黄河下游防洪和大堤安全的影响作出评价，并对其治理提出合理的建议。 第二个实例应用一、二维嵌套水流模型和前后处理辅助软件，对宁波市甬江清水浦码 头的防洪评价进行数值模拟研究。本课题研究在三种频率洪水与大小潮流遭遇时，码头修 建前后沿程的水位变化，局部河势的变化，以及码头附近的冲淤等问题。最后，根据研究 成果，对拟建码头的防洪影响作出评价，并为减小码头的不利影响提出合理的建议。 6 望! ! 查型查皇笺堕堡兰兰垡兰塞 第二章一、二维河流数学模型的基本理论 2 1 一维河流水沙数学模型基本控制方程和辅助方程 数学模型包括水流模型和非耦合的泥沙输移模型两部分。水流模型主要有两种形式： 一是恒定流模型，二是非恒定流模型。由于泥沙输移计算的大部分理论还是以恒定流输沙 为主，在非恒定流输移上还没有重大突破，因此采用以上两种计算方法对泥沙模型来说没 有什么本质差异。 关于明渠非恒定流问题，1 8 7 1 年法国圣维南( s a i n t v e n a n t ) 根据b o u s s i n e q 提出的缓变 流定义而建立的非恒定流方程，即圣维南水力方程组，至今仍然被工程界所采用，成为一 维非恒定流数值模拟的基础1 4 j 。 一维河流数学模型常用于研究河道的长时期的河床演变问题【7 】，因此泥沙的输运和沉 降计算是一维水沙数学模型的一个重点内容。本论文没有对该问题进行深入研究，只是对 已有研究成果作一些分析和资料整理，并根据河流的特点选用比较合适的泥沙输运和沉降 计算公式，采用比较合理的计算处理方法。 2 1 1 一维非恒定非平衡水沙数学模型的基本控制方程 本论文水流模型采用的水流基本控制方程是圣维南水力方程组的一种变形，因此本模 型的适用条件必须满足圣维南水力方程组的基本条件。为此，下面首先简单介绍圣维南水 力方程组的基本条件。 b o u s s i n e q 在引入缓变流的定义时，曾对缓流作了以下几点解释 4 1 ： ( 1 ) 明渠中的水流具有小坡度和微弯曲的轴线； ( 2 ) 过水断面可以为任意形状，但渠底与自由水砥和水平轴线间的倾角较小； ( 3 ) 水流参变数( 如流速、水面坡等) 变化缓慢，随时间而变化的速率很小： ( 4 ) 水流速度几乎和渠轴平行，流层的曲率也很小。 根据上述定义，1 8 7 1 年法国圣维南进一步作了以下几点假设： ( 1 ) 水流是一维的，即过水断面的流速呈均匀分布，渠道曲率所产生的离心力可以 忽略不计： ( 2 ) 水力压力分布与静水压力分布一致，即垂直方向的加速度可以忽略不计： ( 3 ) 水流表面在横向是水平的； ( 4 ) 流体密度是均匀不变的： ( 5 ) 边界阻力及湍流效应均作为阻力处理。 根据以上假设，圣维南建立了缓变流的两个方程： 7 华北水利水电学院硕士学位论文 连续方程 丝+ 塑：0( 2 1 ) 研苏 运动方程 百a u + 昙陋引叫筹 c ：吲 式中，a 为过水面积，q 为过水断面平均流量，【，为断面平均流速，亭为水位，c 为谢 才阻力系数，r 为过水断面的水力半径，g 为重力加速度。 上面所给方程，a 和q 是守恒物理量，但u 不是守恒物理量，因此上述方程不是守 恒形式，虽然可以保证质量守恒，但不能保证动量守恒。由于天然河道过流面积沿流程的 变化比较大，而水位的变化比较小，可以将连续方程的触研改为b a 孝，o t ，虽然这样的 形式不是守恒形式，但仍然能保证质量守恒 。基于上面的分析，模型采用的基本控制 方稗为： 连续方程 运动方程 丝+ 土塑：o 0 tb 良 署+ 鲥一占( 等) 2 豢+ 去( 等) = 一g 嬲 ( 2 3 ) ( 2 4 ) 式中，a 为过水面积，q 为过水断面平均流量，亭为水位，b 为水而宽度，c 为谢才 阻力系数，r 为过水断面的水力半径，g 为重力加速度。 本论文的维水沙数学模型只考虑悬移质的输运，因此泥沙输运和沉降计算的基本控 制方程为f 2 】： 泥沙连续性方程： 河床变形方程： 塑+挈：喇b(s一)at盘 、 厂訾= a a , b ( s 删 ( 2 5 ) ( 2 6 ) 式中，s 为泥沙混合浓度，为悬沙沉速，n 为泥沙恢复饱和系数，& 为水流挟沙力， 也s 为由悬移质引起的断面冲淤面积，7 为冲淤部分床沙干容重。上述各物理量的单位都 采用s i 制。 2 1 2 一维非恒定非平衡水沙数学模型的辅助方程 上面水沙数学模型的基本控制方程只有四个，因此只能选择4 个基本物理量。基本物 理量可以选择原始物理量、守恒量和特征量。在本论文建立的水沙数学模型中，选择水位、 流量、含沙浓度和断面冲淤面积作基本物理量，其他物理量必须借助一些辅助方程求解。 ( 1 ) 谢才系数俨卅 谢才系数c 是反映明渠水流阻力的系数，它直接影响明渠的过流能力。它的计算常 华北水利水电学院硕士学位论文 用满宁公式和巴浦洛夫斯基公式，在本文的水流模型中采用满宁公式，即：c = 二r 。 以 ( 2 ) 水流挟沙力公式 水流挟沙力的定义为：在一定的水流泥沙综合条件下，水流能够携带的悬移质中的床 沙质部分的临界含沙量【2 ”。当悬移质中的床沙质含量超过这一临界数值时，水流处于超 饱和状态，河床将发生淤积；反之，当悬移质中的床沙质含沙量不足这一临界数值时，水 流处于次饱和状态，河床将发生冲刷。目前，水流挟沙力的计算公式有很多，如张瑞瑾公 式、沙玉清公式、窦国仁公式、曹如轩公式和张红武公式等【2 m 4 1 。这些公式大多为根据 经验公式，具有1 定的适用条件和范围。为了使水沙数学模型具有比较广的应用范围，水 流挟沙力的计算不应该采用单一计算公式，而应该将几种常用的计算公式都纳入水沙数学 模型，让使用者根据具体模拟的河道选择比较合适的计算方法。考虑到张瑞瑾对泥沙输运 和沉降有其一套比较完善的理论，所以张瑞瑾公式为本论文水沙模型首选公式。根据文献 阱】的介绍，张红武公式具有适用范围广、精度高的优点，因此张红武公式也被本论文水 沙模型选用。由于各种因素，本论文水沙数学模型还没有将其他水流挟沙力公式纳入备选 公式。 张瑞瑾于2 0 世纪5 0 年在整理了大量的长江、黄河及若干水库、渠道和室内水槽试验 资料后，通过对前人的计算公式修正，得到如下著名的水流挟沙力公式i 划： “ ( 2 7 ) 式中：s + 、u 、h 、分别为断面平均的挟沙力、流速、水深和泥沙沉速；聊、j | 为指 数和系数，一般由实测资料确定。对于一些缺乏实测资料的河道，可以参考张瑞瑾给出的 r ， 、七、m 与的关系图。 g n m 为了进行黄河下游洪水的输沙计算，张红武等对高含沙水流进行了深入研究，得到了 如下的包括全部悬移质泥沙的水流挟沙力计算公式【2 5 】： 鼠= 2 5琴嘧r 沪鼬 j e 式单位k g m s 制，其中跏为体积含沙浓度，矿为流速，h 为水深， 为浑水容重 ( 以= y + o 6 2 3 s ) ，以泥沙容重，为卡门常数( r = o 郇一4 2 瓦( o 3 6 5 一& ) 】) ，d 5 0 为 泥沙中值粒径，q 为泥沙群体沉降速度。 9 华北水利水电学院硕l ：学位论文 ( 3 ) 泥沙沉降速度 在泥沙的基本控制方程和水流挟沙力计算公式都涉及到泥沙的沉降速度。泥沙的沉降 速度受泥沙颗粒形状、温度、絮凝和泥沙浓度等因素影响。此外，在泥沙输运计算中，关 心的不是单个泥沙的沉速，而是泥沙群体或某一级配泥沙的沉速。目前泥沙的输运计算研 究向非均匀沙非平衡输运，而由于本人能力和时间有限，目前的水沙数学模型只考虑非粘 性均匀沙的非平衡输运。 泥沙沉速的计算公式已有比较多。为了与水流挟沙力计算公式相对应，选用张瑞瑾泥 沙沉降计算公式和文献f 2 5 】推荐的泥沙群体沉速公式。 张瑞瑾泥沙沉速计算公式是根据阻力叠加原理和分析实测资料公式得到的过渡区泥 沙的沉速公式1 2 ，具体公式如下： 瓜扎c ：等鲥 ( 2 9 ) 式中q = 1 3 9 5 ，c 2 - 1 0 9 ；d 为泥沙粒径；，。为泥沙容重；y 为清水容重。 由丁张瑞瑾公式计算的沉速是单个泥沙的沉速，因此需要换算为泥沙群体的沉速，本 论文采用钱宁总结的非粘性均匀沙的沉速公式【2 4 】： 堕：( 1 一鼠) ” ( 2 1 0 ) 式中：峨代表含沙量为s ，时的均匀沙沉速，国为清水中单颗泥沙沉速，s 。体积百分 比含沙量，n 为与沙粒雷诺数有关的指数。 文献( 2 5 l 在介绍上述张红武公式时，给出了如下的群体沉速计算公式： 一 ( 1 - 毒卜乃鼠) 1 沪 式中：哦代表含沙量为s v 时的均匀沙沉速，o 为清水中单颗泥沙沉速，为体积百 分比含沙量，以。为泥沙中值粒径。 ( 4 ) 泥沙的起动流速和止动流速 泥沙恢复饱和系数q 是一个与河床的冲淤状态相关的系数。河床的冲淤状态一般可 以根据水流速度与泥沙的起动流速和止动流速的关系确定口6 1 。当流速大于起动流速时， 可以认为河床处于冲刷；当流速小于止动流速时，可以认为河床处于淤积；而当流速处于 起动流速和止动流速之间时，认为河床处于平衡状态。 本模型的起动流速公式采用张瑞瑾提出来的包括散粒体和粘细沙在内的统一的起动 流速公式【2 4 】： 虬= 肛磊霹c ( 2 - - 1 2 ) 0 兰! ! 查型查皇堂堕堡主兰焦笙苎 式中： 为水深；a 为浮容重系数；d 为泥沙粒径；单位采用1 1 一s 制。 泥沙由运动转为静止状态的临界水流流速称为止动流速( 也) 。止动流速一般根据它 与起动流速的关系确定，目前常用公式( 2 1 3 ) 计算阱】。 u ：= ( o 7 1 0 8 3 ) u ( 2 1 3 ) 2 1 3 断面水力参数处理方法 为了能真实反映河床的阻力影响，本论文一维模型没有对断面形态进行概化，而是 采用建立水面宽度、面积、湿周、断面综合糙率与水位的关系来计算这些断面水力参数。 断面综合糙率的常用 = f 算方法有水力半径分割法和能坡分割法。本模型采用水力半径 分割法，具体的计算公式为哪l ： 胛：宇岛芦 智z ( 2 1 4 ) 式中，胛为综合糙率，? c 为断面总湿周，j c f 为第i 个子断面湿周，r l i 第f 个子断面的糙率， ， l 为断面划分的子断面个数。 虽然本论文没有对断面形态做限制，但水流在所给断面形态下必须是缓流，否则可能 会计算失败。 2 2 平面二维河流水沙数学模型的基本控制方程和辅助方程 对于平原河流、河口、海岸、湖泊、大型水库等广阔水域地区，水平尺度远大于垂向 尺度，水力参数在垂直方向的变化较之沿水平方向的变化要小得多，假定沿水深方向的动 水压强分布符合静水压强分布，可以略去水力参数沿垂线的变化，用沿水深积分的平均流 动量来表示，因此可采用平面二维水动力数值模拟技术进行模拟。 2 2 1 平面二维河流水沙数学模型的基本控制方程 目前河, 流s t z 面二维数学模型一般采用非饱和、非耦合、浑水模型，控制方程包括：水 流连续方程、水流运动方程、泥沙连续方程及河床变形方程。水流连续方程和水流运动方 程根据三维时均雷诺方程( r e n o l d z ) 沿水深积分平均得到，雷诺应力采用布辛涅斯克 ( j b o u s s i n e s q ) 的零方程模型和勃兰特动量传递及掺长假设处理，使雷诺方程组得以封 闭。根据质量守恒定律及分子扩散定律( 费克定律) 得到三维泥沙连续性微分方程，继而 沿水深积分平均可得平面二维泥沙连续性方程 2 s , 2 9 1 。忽略风应力和柯氏力的影响，采用中 值粒径作为悬移质和推移质的代表粒径，基本控制方程采用如下的守恒形式： 水流连续方程： 善+ 掣+ 警：0 ( 2 1 5 ) 华北水利水电学院硕十学位论史 x 方向水流运动方程； 鲁+ 警+ 警+ 咖箬+ g 孕硝f 塑o x 2 芬) 协 百+ i + 百+ 驴蔷增歹一邵圳一+ 矿j u 。 y 方向水流运动方程： 百a h v + i a h u v + 可a h v v + 曲雾+ g 丁v 拓2 + v 2 = 砌f 窘+ 窘 c z 椰， 百+ i + 可憎苗增丁可_ 丽+ 矿j 坦0 ” 泥沙连续方程： 警+ 警+ 等+ 鲫p s ) = 瓦a ( 4 - 罢) + 言( 巧 等) c ：郴) 西 良却 、 7 知、1 出7 却、却7 悬移质河床变形方程： ，鲁= 鲫( s 一) ( 2 - 1 9 ) 式中：h 为水深；l 为水位；“、v 分别为x 、y 方向的垂线平均流速分量；s 、分别为按 均匀沙计算的垂线