（应用数学专业论文）河流河道二维数学模型.pdf

摘要 论文题目：河流河道二维数学模型 学科专业：应用数学 研究生：王阳 指导教师：王秋萍副教授 秦毅教授 摘要 随着计算机技术的飞速发展，人们在研究河流泥沙问题时，越来越多的采用了数 学模型，研究对象也由一维问题逐步发展到二维以及三维问题。水沙数学模型以其方 便高效、费用低廉的特性，己经成为研究河道水沙运动规律和河床演变过程的一种重 要手段之一。二维水流泥沙模型的研究与运用己有一段的时间，许多成果己经应用到 水利、环境、港口航道、防洪等工程领域中，为工程的规划设计、施工及水域的可持 续发展提供了重要的指导作用。 本文结合陕西高校省级重点实验室重点科研项目渭河下游河槽萎缩机理及其对 防洪形势的影响研究，主要针对河流河道二维数学模型进行了研究，做了如下一些 工作： 1 在平面二维水流泥沙数学模型中，挟沙力的计算采用水流挟沙力双值公式，同 时考虑贴边淤积现象，依据质量守恒定律和动量守恒定律，对已有的平面二维水流泥 沙数学模型进行改进，建立新的平面二维水流泥沙数学模型。 2 本论文采用斜对角笛卡尔方法对不规则的复杂边界进行模拟，模型为非耦合输 沙模型。数值计算方法采用a d i 方法，离散时采用迎风格式，计算网格采用正交均匀 网格。 3 利用m a t l a b 编制程序对平面二维水流泥沙数学模型进行数值计算，计算水槽 实验数据，修正模型。 关键词：斜对角笛卡尔方法；挟沙力双值公式；平面二维水流泥沙数学模型； a a t l a b ； 数值模拟 本研究得到陕西高校省级重点实验室重点科研项目( 0 4 j s l 6 ) 的资助 a b s t r a c t t i t l e ：t w o d l m e n s i o n a lm a t h e m a t i c a lm o d e lo f f l u v i a lw a t e r c o u r s e m a j o r ：a p p l i e dm a t h e m a t i c s n a m e ：y a n gw a n g s u p e r v i s o r ：a s s o c i a t ep r o f q i u p i n gw a n g p r o f y iq i n a b s t r a c t s i g n a t u r e ：酶鳖 s i g n a t u 陀：峰雩坠丝篮 s i g n a t u r e ：盥荔 w i t h r a p i dd e v e l o p m e n to f c o m p u t e rt e c h n o l o g y , m a t h e m a t i c a lm o d e li sa p p l i e dm o r ea n d m o r ei nt h es t u d yo f p r o b l e ma s s o c i a t e dw i t hs e d i m e n tt r a n s p o r to f ar i v e r a n dt h et y p eo f t h e m o d e lh a sb e e ne x p a n d e df r o mo n e - d i m e n s i o nm o d e lt ot w o - d i m e n s i o n , e v e nt h r e e - d i m e n s i o n o n e t h ef l o w - s e d i m e n tm a t h e m a t i c a lm o d e lh a sp l a ya ni m p o r t a n tr u l ei ns t u d y i n gt h e m o v e m e n tl a w so fw a t e rf l o wa n ds e d i m e n tm o t i o na n dr i v e r b e de v o l u t i o no far i v e rd u et oi t s c h a r a c t e r i s t i c so fc o n v e n i e n c e ，h i g he f f e c t i v e n e s sa n d l o wc o s t i th a sb e e na l o n gh i s t o r yf o r t h i sk i n do ft h em o d e lt ob ea d o p t e da n dr e s e a r c h e d m a n yr e s u l t sf i o mt h em o d e lh a v eb e e n a p p l i e dt op r a c t i c a le n g i n e e r i n g , s u c ha sw a t e rc o n s e r v a n c y , e n v i r o n m e n t , p o r ta n dc h a n n e l ， n o o dd e f e n s e a n de t c ；a n dt a k e na sa ni m p o r t a n tg u i d a n c ef o rt h ep l a n n i n g , c o n t i n u o u s d e v e l o p i n ga n du t i l i z i n go f w a t e r r e s o u r c e s a s s o c i a t e dw i t ht h ep r o j e c t ( 0 nt h em e c h a n i s mo f t h es h r i n k a g eo f d o w n - s t r e a mc h a n n e lo f w c ir i v e ra n dt h er e s p o n s eo ff l o o dc o n 拄o lt ot h es h r i n k a g e s u p p o r t e db yt h ee d u c a t i o n d e p a r l m e n to f s h a u n x i t k st h e s i sa i m st or e s e a r c ht h ec h a n n e lt w o - d i m e n s i o n a lf l o w - s e d i m e n t m a t h e m a t i c a dm o d e lf o f i o w sa r ct h ew o r k sd o n e ： 1 f r o mt h el a w so fc o n s e r v a t i o no fm a s sm o m e n t u m , an e wp l a n et w o - d i m e n s i o n a l f l o w - s e d i m e n tm a t h e m a t i c a lm o d e lw a sb e e ne s t a b l i s h e do nt h eb a s i so f t h em o d i f i c a t i o no f t h e c x i s t i n gm o d e l 1 1 1 en e wm o d e la d o p t st h ed o u b l ev a l u es e d i m e n tc a r r y i n gc a p a c i t yf o r m u l aa n d t a k e st h eb o u n d a r yg l o b u l es i l t a t i o ni n t oa l la c c o u n t , w h i c hw a sn e v e rb e e nd o n eb e f o r e 2 t h et h e s i su s e dt h ed i a g o n a lc a r t e s i a nm e t h o dt oc a n - yo u tt h es i m u l a t i o nt ot h ei r r e g u l a r c o m p l e xb o u n d a r y 1 1 p r o p o s e dm o d e li san o n - c o u p l i n gs e d i m e n tt r a n s p o r tm o d e l 1 1 硷a d i m e t h o dw a sf o rt h en u m e r i c a lc o m p u t a t i o nm e t h o d n l eu p s t r e a ms c h e m ew a su s e di nt h ec a s e o f d i s c r e t e ，a n dt h eo r t h o g o n a le v e n 鲥d i sf o r t h ec o m p u t a t i o ng r i d 3 1 km a t l a bc o m p u t e rp r o g r a m m i n gf o rt h em a t h e m a t i c a lm o d e lw i t ht h ed a t af i o m e x p e r i m e n th a dc a r r i e do u tt h ev a l u ec o m p u t a t i o n a n df i n a l l yt h em o d e l w a sb e e nm o d i f i e d k e yw o r d s ：d i a g o n a lc a r t e s i a nm e t h o d ；d o u b l ev a l u es e d i m e n tc a r r y i n gc a p a c i t yf o r m u l a ； p l a n e2 - d f l o w s e d i m e n tm a t h e m a t i c a lm o d e l ；m a t l a b ；n u m e r i c a li m i t a t i o n 1 1 1 et h e s i si s s u p p o r t e db yt h ep r o j e c to fk e yl a bc o n s t r u c t i o ni s s u e db yt h ee d u c a t i o n d e p a r t m e n to f s h a a n x ip r o v i n c e ( 0 4 j s l 6 ) n 主要符号表 主要符号表 水深 水位 地形高程 节点的冲淤厚度 分别为工、y 方向的单宽流量 垂线平均流速沿工，j ，方向的分量 糙率 流水的紊动粘性系数 分别为第k 粒径组泥沙的含沙量 分别为第k 粒径组泥沙的有效沉速、恢复饱和系数 非均匀沙的加权平均沉速 重力加速度 分别为x 、j ，方向的泥沙扩散系数 床沙干密度 泥沙和浑水容重 泥沙起动速度 水流速度 水力学半径 摩阻流速 玑在x 、y 方向上的分量 非均匀悬移质泥沙分组数 以 b k 以 矗z 互嵋” 以 咋卧蚧g扑y盂以 独创性声明 秉承祖国优良道德传统和学校的严谨学风郑重申明：本人所呈交的学位论文是我个 人在导师指导下进行的研究工作及取得的成果。尽我所知，除特别加以标注和致谢的地 方外，论文中不包含其他人的研究成果。与我一同工作的同志对本文所论述的工作和成 果的任何贡献均已在论文中作了明确的说明并已致谢。 本论文及其相关资料若有不实之处，由本人承担一切相关责任 论文作者签名：至鱼坦7 年；月堀日 学位论文使用授权声明 本人垒垫在导师的指导下创作完成毕业沦文。本人已通过沦文的答辩，并 已经在西安理工大学申请博士硕士学位。本人作为学位论文著作权拥有者，同意授权 西安理工大学拥有学位论文的部分使用权，即：1 ) 已获学位的研究生按学校规定提交 印刷版和电子版学位论文，学校可以采用影印、缩印或其他复制手段保存研究生上交的 学位论文，可以将学位论文的全部或部分内容编人有关数据库进行检索；2 ) 为教学和 科研耳的，学校可以将公开的学位论文或解密后的学位论文作为资料在图书馆、资料室 等场所或在校园网上供校内师生阅读、浏览。 本人学位论文全部或部分内容的公布( 包括刊登) 授权西安理工大学研究生部办 理。 ( 保密的学位论文在解密后，适用本授权说明) 论文作者签名：垄l 叠导师签名：! 塞聋。勰妒每；月z 2 日 第一章绪论 1 绪论 1 1 研究背景及意义 随着社会的发展，人们对河流的利用日趋多样，对河流的认识也日益深入。模拟河流 的方法有物理模拟和数值模拟两种类型。物理模拟即比尺模拟，以其直观性和物理概念明 确一直受到水利工程界的重视，但其模拟程度往往受模拟试验手段和比尺关系的限制，而 且费用比较高，试验时间、周期也相对较长。河流数值模拟是以挟沙水流运动和河床演变 控制方程为基础，以数值方法和计算机技术为手段，通过对河流的数值模拟计算，解决河 流工程所关心的问题。随着计算机科学的迅猛发展，使数值模拟有了强大的生命力，它以 其省时、高效、无比尺效应等优点，在众多河流工程问题上被广泛应用，取得了良好的效 果。如在水力学计算中，用一维数学模型或二维数学模型进行河流、河网、河口的数值模 拟已很普遍，模拟技术已较成熟，并且已在流域规划、水环境保护、河网开发、河口整治 等方面取得了巨大的经济效益，因此，计算机数值模拟己成为一种水科学研究的重要手段。 渭河是黄河的最大支流，自宝鸡峡出峡后，横贯作为陕西粮仓的关中平原，沿岸的西 安、咸阳、宝鸡、渭南等城市均依赖渭河滋润灌溉，渭河的安危对陕西省一线两带经济发 展战略有着重要影响。渣关断面位于黄河与渭河汇流区下游，由于受三门峡水库回水淤积 和近期黄河来水减少等因素的影响，致使近年来潼关高程持续上升且居高不下，渭河下游 河道也淤积抬高，常有洪泛发生，给渭河下游防洪带来了较大危害。因此，长期以来，寻 求降低渭河下游高程的有效途径成为一项十分重要的任务，有必要对渭河的冲淤与影响作 数值模拟。 以计算水动力学、泥沙运动力学与河床演变学等学科为基础的水流泥沙数学模型近 年来得到了很大发展，它是预测冲积河道内水流泥沙运动及河床变形的重要工具之一，而 且与物理模型相辅相成，是解决实际工程问题的重要手段。水动力数值模拟是水流泥沙数 值模拟的重要组成部分，它为泥沙的运动迁移提供了最基本的环境背景。泥沙的任何行为 都离不开水体本身这个介质，通过水动力的数值模拟，可了解水质点在某一时刻的流速、 流向、水量的输送，湍流的强度以及涡旋、环流状况等。在水动力数值模拟基本方程的基 础上，再耦合泥沙运动方程和河床变形方程以及一些辅助方程就构成了水流泥沙数学模型 的基本方程组“1 本论文结合陕西高校省级重点实验室重点科研项目渭河下游河槽萎缩机理及其对 防洪形势的影响研究，主要针对平面二维水流泥沙数学模型的数值模拟，在已有的研究 成果的基础上，采用适合渭河情况的水流挟沙力双值公式，同时考虑高含沙水流河床切应 力的影响、贴边淤积现象，对已有的平面二维水沙数学模型进行改进，建立新的平面二维 水沙数学模型，利用斜对角笛卡尔方法对不规则的复杂边界进行模拟，用a d i 法求解数值 西安理工大学硕士擘位论文 解。研究成果不仅为理论工作者进一步深入研究提供可行的数值模型，而且为实际工作者 具体应用提供简单易行的数值算法，具有一定的理论价值和应用价值。 1 2 国内外水沙数学模型研究 早在二十世纪5 0 年代，苏联、美国、西欧和我国就开展了水流泥沙数学模型的研究 和应用，主要是用来计算大型水库的淤积和水坝下游的沿程冲刷。由于受到当时计算工具 的限制，这种计算仅限于一些简单的状况，而且精度较差，工作量大，因此，未能得到普 遍推广。二十世纪7 0 年代以来，随着电子计算机的普及与计算机性能的不断提高，这些困 难已逐步克服，数学模型的优越性越来越突出。数学模型投资相对较少、运行周期短和效 率高，具有通用性，重复模拟时具有理想的抗干扰性，可严格控制并可方便修改边界条件 及其他条件，同时具有模拟真实条件和模拟理想条件的能力，因而日益成为研究水流、泥 沙等运动规律的一个强有力的手段，普遍受到人们的重视“1 。经过近二、三十年的发展， 水流泥沙数学模型经历了由一维到二维、到二维嵌套、到三维，到三维嵌套，由原先在概 化水文、泥沙及河床条件下数学模型的建立、率定和验证，发展到与物理模型相互配合共 同回答工程水流和泥沙闯题( 复合模型) 的阶段“1 。 一维水流泥沙数学模型模拟的是变量在断面上的平均值，基本上能够满足实际工程 的需要，它是至今使用最广泛的一种模型。一维模型在理论及实践上都比较成熟，模型计 算省时，可快速方便地进行长河段、长时间的洪水和河床演变预报，因而国内外使用较为 普遍。已出现很多比较成熟的一维水流泥沙数学模型，如美国陆军工程兵团开发的h e c - 6 模型0 f e l d m a n , 1 9 8 1 ) ，能用来计算河道和水库的冲淤过程“1 韩其为等0 9 8 7 ) 根据泥沙运 动统计理论，建立了一维非均匀悬移质泥沙的不平衡输沙模型，该模型运用范围较广，并 通过多个长河段实测资料的验证仆，。还有杨国录和吴卫民开发的s u s b e d - i 模型( 1 9 9 4 ) ， 该模型为一维恒定平衡与不平衡输沙模式嵌套计算的非均匀沙模型，可用于计算水库及河 道内的水沙变化与河床变形“1 李毓湘等利用一维非恒定流基本方程组、河网节点连结 方程及边点方程建立了珠江三角洲河网区水动力模型m 1 一维水流泥沙数学模型很多， 可根据特点分为三大类：一是以适宜水温相关分析为基础的模型；二是以水动力学和泥沙 运动力学为基础的模型；三是介于上述两类之间的模型。 由于一维水流泥沙数学模型研究较早，故而也较早应用于实际工程中但是，一维 模型无法给出各物理量在断面上的分布，因而在模拟河床细部变形、河口和港湾等水域的 流动和冲淤问题时，必须用二维模型甚至三维模型才能解决。 为了克服一维水流泥沙数学模型无法模拟水流沿横向( 河宽方向) 、垂向( 水深方向) 变 化的不足，二维水流泥沙数学模型便应运而生，并且发展迅速，目前在工程中应用广泛。 二维水流泥沙数学模型按照横向和垂向分为平面二维水流泥沙数学模型和立面二维水流 数学模型( 或垂向二维水流数学模型) 。在水平尺度远大于垂向尺度的河流及河口地区，各 2 第一章绪论 物理量在垂向的变化可以使用沿水深平均的量来表示，这就构成了平面二维水流数学模 型。这种模型基本上可以给出工程上所需要的一些水力参数，因此在工程上应用也较为广 泛。另一种是立面二维水流泥沙数学模型。在宽度沿程变化不是很明显的航道地区，以及 流态比较顺直的河流水域，为了研究各物理量沿水深方向的变化情况，就要用到立面二维 水流泥沙数学模型。它以在水流中截取的纵剖面上的水流因素作为研究对象，主要研究水 流垂向运动问题。与一维水流泥沙数学模型相比，二维水流泥沙数学模型不管是在建模的 基本方程上，还是使方程组封闭的辅助方程上，以及针对某一具体工程的实际应用中，都 要复杂和困难的多。二维水流泥沙数学模型早期应用于河口、海湾、湖泊等宽阔水域的数 值模拟，经过十多年的发展已经广泛应用到河流中。 平面二维水流泥沙数学模型近年来发展较快，主要用于研究短时间内重点或局部河 段的河床冲淤变形，经实测资料检验，已能近似反映实际情况。国外从2 0 世纪7 0 年代后期 开始推广运用，而国内在近几十年来发展较快。国内现已有多个平面二维水流泥沙数学模 型，在实际工程中用得较为成功( 李义天，1 9 8 7 ”洞建军等，1 9 9 3n 1 ；窦国仁等，1 9 8 7 “，1 9 9 5 1 ) 。此外，国外也已开发了许多平面二维水流泥沙数学模型，如德国k a r l s r u h e 大学的 f a s t 2 d 模型、荷“- d e l f t 的v a n r i j n 模型、t r i m 2 d 模型、i o w a 水利所的m o b e d 2 模型、美 1 虱m i s s i p s i p p i 大学水科学计算中一d c c h e 2 d 模型、丹麦水利所的m i k e 2 1 模型等。 局部问题的三维水流泥沙数学模型，由于计算工作量大，目前还很难应用于实际工 程( s h i m i z u , 1 9 9 0 ) t t o f a n g 和w a n g ( 2 0 0 0 ) 采用非正交曲线坐标下的三维水流泥沙数学模型， 模拟了三峡库区近坝段的冲淤过程“”。三维计算还不成熟，且工作量繁琐巨大，而平面 二维计算既能满足实际需要，又相对方便经济，在实际工程中广泛采用。平面二维模型己 经用于研究河口海岸泥沙运动及河道泥沙运动问题上。目前对平面二维泥沙模型的研究工 作主要在两方面进行：一是对平面二维水流、泥沙运动基本规律的研究；二是对计算方法 的研究。前者关系到数学模型是否真实地反映描述物理现象的基础，后者关系到是否可以 得到数值解以及所得到的数值解是否可靠，各家数学模型的主要区别就在于这两个方面问 题的处理上“ 1 3 水沙数学模型基本方程及计算方法 1 3 1 一维水沙数学模型 尽管目前有多家一维水流泥沙数学模型，但各种水流泥沙数学模型其理论基础都基 于四大控制方程，即水流连续方程，水流运动方程，泥沙连续方程和河床变形方程，其主 要区别在于求解方程时，各家采用的数值计算方法、水流挟沙力计算方法、动床阻力计算 方法等有所不同。 一维水沙数学模型基本方程( 谢鉴衡，1 9 9 0 ) n ”： 1 水流连续方程式： 西安理工大学硕士学位论文 望+ 曰丝玑 ( 1 1 ) a x o t 。 2 水流运动方程式： 詈+ ( 鲥一口第) 罢+ 2 罢罢一l q 彳：2o 缸a 儿、i + 哟叫吼 c t 3 泥沙连续方程： 掣+ 掣饥聃一) & 吼 ( 1 3 ) 4 由悬移质引起的河床变形方程： p 訾一冀口t 丑 一) ( 1 4 ) 式中：妫流程；g 为流量；吼为旁侧流入量；坼为旁侧流入主流流速分量；z 为水位； b 为河宽；沩时间：为过水断面面积；g 为重力加速度；，为能坡；& 、s 。、& 分别为 第女粒径组泥沙的含沙量、水流挟沙力、旁侧入流含沙量：n ，t 、口。分别为第t 粒径组泥沙 的有效沉速、恢复饱和系数；4 为总的河床变形面积；p 为床沙干密度；为非均匀沙 分组数。 一维水流泥沙数学模型的计算方法可以分为两大类：一是将水流和泥沙方程直接联立 求解，称为耦合解，适用于河床变形比较剧烈的情况；二是先解水流方程式求出有关水利 学要素后，再解泥沙方程式，推求河床冲淤变化，如此交替进行，称为非耦合解，适用于 河床变形比较缓和的情况。对长河段、长时间的河床变形计算，多采用非耦合的一维恒定 水流泥沙数学模型，这种方法简单实用( 谢鉴衡，1 9 9 0 ) “对于动床条件下的洪水演进 过程，常用一位非恒定水流泥沙数学模型，采用p r e i s s m a n n 四点偏心隐格式离散水流方 程，结合线形化技术，并用追赶法求解个变量( 杨国录，1 9 9 3 ) ”1 。 一维水沙数学模型还存在着一些关键性问题需要解决，如动床阻力变化的问题、水流 挟沙力与床沙级配调整的计算问题、断面冲淤面积的横向分配问题以及恢复饱和系数的取 值问题等“。 1 3 2 平面二维水沙数学模型 4 平面二维水沙数学模型基本方程 1 水流连续方程式： 丝+ 堕4 - 盟0( 1 5 )+ = = 1 5 ) m瓠 却 2 ，水流运动方程： 鲁+ 警警+ 等詈+ 劝罢+ 秽2 丁q 斯+ q 2 y v r ( 誓+ 争) n s ， 第一幸绪论 鲁+ 鲁詈+ 等詈+ 劝芳+ 明2 学- v , f 誓+ 争】， 3 悬移质泥沙扩散方程： 砉( 蝇) + 去( 以) + ( q y s t ) + 口止( 墨一) - 号笋+ 了a 2 ( 舾d ( 1 8 ) 4 悬移质河床变形方程： p 导罗口t 蛾一) ( 1 9 ) 脚 镯 以上各式中：吼、g ，分别为x 、y 方向的单宽流量；z 为水位；h 为水深；n 为糙率。 坼为流水的紊动粘性系数；、s k 分别为第k 粒径组泥沙的含沙量、水流挟沙力；g 为 重力加速度；n ka 。分别为第k 粒径组泥沙的有效沉速、恢复饱和系数；f 。、s 。分别 为x 、y 方向的泥沙扩散系数；p 为床沙干密度；n 为非均匀悬移质泥沙分组数；a 乙为 节点的冲淤厚度。 平面二维水流泥沙数学模型在计算中，常用的方法仍是有限差分法，如采用交替方向 的隐格式差分格式，即a d i 法，即能保持计算结果的稳定，又能节省计算时间( 郑帮民 等, 2 0 0 1 ) “。部分研究者采用破开算子法计算，该方法将基本方程按不同的物理意义破 开成对流、扩散、传播三个子方程，然后分别求解( 周建军等，1 9 9 1 ) n ”。破开算予法因 占用内存小、计算简便等优点从而得到广泛应用，但是有时由于计算区域采用方形网格剖 分，这种算法往往对边界拟合不够理想，经常在计算区域的角点处出现奇点，这样势必影 响计算结果的准确性。用有限差分法的计算网格一般为矩形，计算量虽小，但一般不能模 拟复杂边界的流动问题，采用三角形网格可克服这一缺点( 窦国仁等，1 9 9 5 ) n ”。也可采 用正交曲线坐标系下的控制方程来求解，虽然计算量便大，但能适用于不规则的计算边界 ( m o l l s ，1 9 9 5 ) 1 “。平面二维水流泥沙模型的计算，也可采用有限元法( 韦直林，1 9 9 0 ) 和有限体积法( a l e x a n d e r , 2 0 0 0 ) 汹1 ，这两种方法能较好地模拟不规则岸边界。但有限元 方法存在离散格式较复杂、大型矩阵较难求等困难，计算量和存储量远比有限差分法大得 多，因而在现有的河床变形计算模型中应用很少。有限体积法逐渐在二维浅水流动的模拟 中开始应用，该方法的误差主要来自界面通量的估算。它严格遵循守恒律，不存在水量、 动量不平衡的守恒误差。该方法需要解决的难题是如何自动生成网格、如何估计界面通量 等( 谭维炎等，1 9 9 9 ) 2 “。 在一维模型中存在的某些问题，在二维模型中同样存在，而且还更复杂。这些问题主 要为在一维条件下获得的计算公式，如阻力公式、挟沙力公式及有关参数能否用于二维问 题还需作深入研究。在二维模型中，值得注意的是二维控制方程的数值解法问题，在计算 中应该考虑计算的稳定性与收敛性，而且应加快运算速度，减少计算时间( 曹叔尤等，1 9 9 2 ) 西安理工大学硕士擘住论文 1 4 数值计算方法概述 描述河流水沙运动的控制方程一般是非线性的，对于这些控制方程，一般很难求得解 析解，只能通过离散控制方程、求解代数方程组来得到近似数值离散解。从数学模型研究 开始至今，寻找一种集稳定性、相容性、守恒性、精确性、健全性、通用性和经济性于一 身的“优化”的数值计算格式，一直是理论界和工程界的一个活跃课题。通过数十年来的 努力，到目前为止，问世的各种数值计算方法不下百种，但还没有一种能够同时具有这么 多的优良品质 1 4 1 数值离散方法 数值模拟中传统的离散方法旧1 有：有限差分法，有限元法，有限体积法，有限分析 法。 有限差分法是数值解法中的经典方法，它将求解域划分为差分网格，用有限个网格 节点( 即离散点) 代替连续点的求解域，以泰勒级数展开为工具，将偏微分方程的导数用差 商代替，推导出在离散点上含有有限个未知数的差分方程组。求解差分方程组，其解就作 为微分方程定解问题的数值近似解。有限差分法数学概念直观，表达简单，其解的存在性、 稳定性及收敛性的研究成果较完善。 目前在平面二维水沙数学模型中应用较多的有限差分法主要是交替方向隐格式法 ( a d om 1 和破开算子法“”。a d i 是p e a e e m a n ，r c h f o r d 和d o u g l a s 于1 9 5 5 年提出的，是一 种比较流行的古典差分法。a d i 法是一种显隐格式交替使用的有限差分格式，具有显式和 隐式两种差格式的一些优点。它的表达式简单，不像完全隐式格式那样，每一步都需要解 一个大型代数方程组并且要进行迭代，因而所需内存少，计算量也相应减小；在计算过程 中没有显格式常常出现的波动现象，具有较好的计算稳定性和精度。破开算子法也称算子 分裂法或分步法。它是将多维问题分解成一系列的低维问题，通常是化成一系列的一维问 题求解，或者将原始的包含若干物理性质不同的复合过程分解成若干性质单一的物理过 程。虽然可用破开算子法灵活多变地将一个复杂的数学物理问题化成较为简单的闯题来处 理，在数值处理和计算上带来了许多便当之处，但仍存在不少问题，如物理现象被解体等， 目前许多学者对破开算子法提出疑义，因此该方法应用不是很广泛。 人们曾经认为有限差分法只能应用于均匀的网格，当边界形状与网格线不一致时便 产生一些困难。但自1 9 7 4 年汤普森( t h o m p s o n ) 等人发展了与边界适配的曲线坐标技术后， 可将不规则的物理区域映射到均匀的计算区域，在计算区域上求解椭圆方程得出物理平面 上正交网格的位置，再在均匀的计算区域上求解变换后的差分方程。该技术使有限差分法 重新获得优势，不论是处理简单边界还是复杂边界，都具有很强的灵活性，而且离散方程 系数的物理意义明确，计算程序易于编制，计算方法成熟且具有良好的稳定性。现在有限 差分法可以在曲线坐标系中求解，一种方法是，直接用曲线坐标系来计算，对河岸线可用 6 第一章绪论 一定的曲线来描述的计算域来说，计算比较简便；另一种方法是进行坐标变换，如用拉普 拉斯方程将曲边域变成矩形域，对控制水沙运动的方程也做相应变换，在变换后的计算平 面进行求解，然后将计算结果再变换回原来的曲边域去，这种方法计算十分复杂，计算工 作量大。 有限元法是将一个连续的求解区域任意分成适当形状( 三角形或四边形) 的许多微小 单元，并在各小单元分片构造插值函数，然后根据极值原理( 变分或加权余量法) ，将问题 的控制微分方程化为控制所有单元的有限元方程，把总体的极值作为各单元极值之和，即 将局部单元总体合成，形成嵌入了指定边界条件的代数方程组，求解该方程组就得到各节 点上待求的函数值。 有限元的求解方法常见的有直接法、变分法、加权余量法和能量平衡法等。其优点是 网格划分灵活，拟合复杂河岸边界容易，网格节点可局部加密，稳定性好，精度高，适合 于几何、物理条件复杂的问题。对于隐式f e m ，其精度较高，但数学推导繁杂，计算量 和储存量较大，而且在误差估计、收敛性和稳定性等方面的理论研究与有限差分法相比还 显得不够成熟和完善。尤其在多维计算中，由于有限元法贮存量大，直接影响着计算速度。 有限元在用于非恒定流计算时，每个时间步都要求解一个大型线性方程组，耗时多。另外 g a l e r k i n 有限元在数学上适于求解椭圆方程的边值问题，其性能类似中心差分格式，缺乏 足够的耗散，不适于计算间断问题，要加入人工粘性。因此，在河道计算方面的应用不如 有限差分法广泛。 有限体积法将所计算的区域分成一系列控制体积，每个控制体积都用一个节点作代 表。通过将守恒型的控制方程对控制体积做积分来导出离散方程。在导出离散方程的过程 中，需要对界面上的被求函数本身及一阶导数的构成做出假定。 有限体积法在计算出通过每个控制体边界沿法向输入( 出) 的流量和动量通量后，对每 个控制体分别进行水量和动量平衡计算，便得到计算时段末各控制体平均水深和流速。因 此，有限体积法正是对于推导原始微分方程所用控制体途径的回归。因为跨控制体间界面 输运的通量，对相邻控制体来说大小相等、方向相反，故对整个计算域而言，沿所有内部 边界的通量相互抵消。对由一个或多个控制体组成的任意区域，以至整个计算区域，都严 格满足物理守恒律，不存在守恒误差，并且能正确计算间断。由于采用守恒型的微分方程 并对每一计算单元进行质量和动量守恒形式的离散，使得微分方程包含的守恒性质在每一 个控制容积上都得到满足，若保持各单元界面两侧相邻控制体的计算输运通量相等，那么 整个计算区域上都能保持守恒。 有限分析法是由美国华裔科学家陈景仁教授在1 9 8 1 年提出的。它克服了有限差分法在 求解不可压粘性流体时大雷诺数下的困难，避免了差分近似中的数值效应，是求解大雷诺 数下的各种流体力学问题行之有效的方法。有限分析法的基本思想是；有限分析法把求解 域划分为多个有限小的单元，在每个单元内求解基本方程的解析解。单元内的解析解把该 单元的一个内点上的参变量与该域边界上八个点的参变量联系起来，以代数方程表示。由 7 西安理工大学硕士学位论文 于单元边界上的点即为其邻域内的点，而每个内点有一个代数方程，因此，联立这些代数 方程组可解出原方程的数值解。求解域边界的参变量为求解联立代数方程组时嵌入的边界 条件。 陈景仁通过局部单元的线性化方程，利用原方程的分块局部分析解，将其相互关联起 来而建立各个离散节点上函数值的关系式，以一定的插值近似边界条件下的分析解来近似 原定解问题的精确解，从而在整体上保持了控制方程的非线性性质，同时使这种格式可以 免除将微分方程离散化而带来的各种误差和虚假效应。有限分析法由计算网格内的分析解 和整个计算区域内的数值解相结合，不存在对流扩散项离散格式选择问题，数值解的精 度高。事实证明，用这种方式所建立起来的有限分析格式确实能比较好的保持原有问题的 数学物理性质，它能够准确的模拟对流扩散效应，也不存在迎风差分格式和有限元方法中 的伪数值扩散现象以及有限体积法中存在的数值耗散和数值频散问题，计算稳定性好，收 敛较快，特别比较明显的是格式系数自动具有迎风的数值特性。同时，有限分析法可以直 接对有限分析解求导，这样可以避免数值求导的误差。 有限元法适合于几何、物理条件比较复杂的问题；有限体积法导出的离散方程对整个 计算区域可以保证具有守恒特性，而且离散方程系数的物理意义明确，但在导出离散方程 的过程中，需要对界面上的被求函数本身及一阶导数的构成做出假定；有限分析法计算精 度高，并且具有自动迎风特性，能准确地模拟对流项，计算稳定性好，收敛较快，但由于 有限分析系数中含有无穷级数，给实际计算带来了一些困难；有限差分法是建立在经典的 数学逼近理论的基础上，其解在线性前提下的存在性、收敛性、稳定性早己有研究成果， 是目前较成熟且应用最广的方法。 1 4 2 不规则边界的处理 在水利工程中，水流的边界往往是很复杂的。而在数学模型计算中如何处理不规则的 复杂边界将直接影响到整个模型计算结果的可靠性、精度和计算时间，是数值模拟的一大 难点，选择复杂边界的处理方法与算法的选择有着同等重要性。因此在研究江、河、湖、 泊的水流数值计算过程中，学者们提出了许多处理边界的方法，例如阶梯近似法，任意网 格差分法，通度系数法，不同坐标系组合法，等参元函数法及贴体坐标法等。计算网格选 择的好坏不仅影响解的稳定性，事实上，网格系统的繁简程度已成为计算流体力学技术普 及应用的主要障碍。 在这些方法中，常用的不规则复杂边界处理方法主要有：美国学者t h o m p s o n 嘲1 提出 的贴体坐标法，因其能使计算域边界上的坐标线与边界线密切贴合，而获得了广泛的研究 与应用。贴体坐标法的主要优点是用曲线坐标逼近复杂边界，在边界模拟上具有较高的精 度。但该方法在较复杂的边界上不易形成正交贴体网格“，如果是非正交曲线坐标系， 则会产生大量的交叉导数项，给求解带来很大的困难。 8 第一章绪论 非结构性网格使用三角形和四边形计算单元，使得在复杂区域形成计算网格更加灵活 和方便。但非结构性网格往往要建立一套标明网格节点之间关系的复杂系统，因此在编程 计算处理时比较麻烦，特别是针对自由面流动和干湿边界变化时，这种困难更加突出。多 重网格法则由于网格系统之间的相对独立性，难以做到对系统守恒性的保证。 传统的笛卡尔坐标方法被应用的历史最长，由于其简单、易形成计算网格、正交等 特点，在工程界应用得最为普遍。但传统的笛卡尔坐标方法在边界处形成锯齿状边界，导 致边界的拟合精度不高，因此该方法在大尺度计算中应用较多。 斜对角笛卡尔坐标方法哺1 不仅保留了锯齿法正交、易形成计算网格、计算效率高等 特点，同时又能有效地改善锯齿法拟合边界精度差的缺点，更能贴近实际，使边界流场更 加合理。 在国内余利仁，许卫新，张书农啊1 使用在一对相邻边界上可任意选定边界网格点的 数值方法，通过逐次求解拉普拉斯方程组，成功的生成了非对称渠槽中的正交贴体坐标系。 河海大学的何国建，汪德憎等以求解拉普拉斯方程组为基础，不合并或聚集任何控制函 数，生成了与边界正交性良好的网格线，并可随意控制网格的疏密程度。清华大学的揭冠 周，介玉新，李广信”1 提出用基于适体坐标变换的有限差分法求解复杂边界的自由面渗 流。西安理工大学的魏文礼咖1 等采用p o i s s o n 方程进行曲线网格生成，并提出了一种新的 构造p ，q 函数的方法，该方法直接利用边界网格节点的信息来控制区域内部网格节点的信 息，利用动边界技术力求实现边界处网格的正交。董耀华n 1 1 基于河势概念和h e n n i t e = 次 插值函数，提出了河势贴体河道平面二维正交四边形网格的生成方法。清华大学的黄炳彬， 方红卫，刘斌”1 在直角坐标系下发展了一套处理复杂边界的斜对角笛卡尔方法。西安理 工大学的刘玉玲。1 将曲线拟合坐标变换技术应用于三角形网格生成中，形成一种方便、 有效的三角形网格生成方法。 1 4 3 本文采用的计算方法介绍 在本文的平面二维水沙数学模型的求解中，考虑到模型的修正方便，数值计算采用 成熟的a d i 法离散控制方程，离散时采用迎风格式，计算网格采用正交均匀网格。a d i 法， 即交替方向隐式格式法，是由p e a c e m a n - r a c h f o r e l 1 和d o u g l a s “1 首先提出的一种专解二 维问题的特殊分步法。其实质是把时间步长分成两个半步，前半步仅对一个空间方向用隐 格式( 如笛卡儿坐标系中的x 方向) ，而对另一个空间方向用显格式方向) ；后半步则相反， 彷向用隐格式，x t y 向用显格式。如此就可以得到两个三对角系数矩阵，可以用追赶法求 解，也可以用迭代法求解。 在不规则边界的处理上，采用斜对角笛卡尔方法对不规则的复杂边界进行模拟清 华大学的黄炳彬等( 2 0 0 3 ) e m 在直角坐标系下发展了一套处理复杂边界的斜对角笛卡尔方 法。斜对角笛卡尔方法基本原则就是按就近逼近原则，把距实际水边线较近的网格节点连 西安理工大学硕士学位论文 线以拟合实际水边线。斜对角笛卡尔方法的优点：更接近实际水边线，处理简单，附加计 算量小，用斜对角笛卡尔法进行边界处理，不仅保留了离散简便，计算效率高的优点，同 时又使边界拟合更加精细，避免了锯齿状边界对流场计算的影响。 总而言之，水沙模型的数值模拟是一个复杂的过程，它的发展是一个不断完善、逐 渐精确的过程。尽管现在所使用的一些计算公式带有经验或者半经验性质，而且一些公式 的通用性还有待商榷，但伴随着计算机的发展，人们对水流泥沙的运动机理不断的认识与 提高，数值模拟技术也将更加完善。 1 5 贴边淤积现象 贴边淤积现象对过洪断面影响巨大，近十几年来，随着渭河下游河槽淤积的加重， 洪水过程线由尖瘦型变成矮胖型，峰线时刻明显推后，洪水历时显著增加，峰量关系也发 生了变化，于防洪形势不利。其原因不仅有河道边界条件变化的影响，更主要的是贴边淤 积影响。这些影响因素的存在，致使洪水位居高不下，因而酿成小水大灾的后果 t i l l 贴边淤积现象如图1 - 1 ，图中虚线为初始水面线，实线为贴边淤积后水面线，右侧发 生贴边淤积。可明显看出贴边淤积直接造成水面线抬升，在冲淤总量相对不变时，水面线 抬高比一般淤积高。 图1 - 1 贴边淤积现象 f i g l i t h ep h e n o m e n o no f b o u n d a r yg l o b u l es i l t a t i o n 1 6 本课题的主要目的及内容 综上所述，目前很多学者致力于河道平面二维问题的研究工作，而且得出了许多实 用的数学模型。但同时也发现由于天然河流情况复杂多变，针对具体河流需要不同的处理 方法。本课题主要目的是在已有的模型基础上，采用适合渭河情况的水流挟沙力双值公式， 考虑贴边淤积现象，建立河流河道平面二维水流泥沙数学模型，并利用有限差分法求数值 解。具体研究内容包括： ( 1 ) 简述水流泥沙模型问题的研究现状。 1 0 第一章绪论 ( 2 ) 采用斜对角笛卡尔方法，来克服由于复杂边界而引起的计算困难，同时在边界处 理时考虑贴边淤积现象的影响。 ( 3 ) 在已有模型的基础上，将适合渭河情况的水流挟沙力双值公式用在具体的模型 中，建立河流河道平面二维水流泥沙数学模型。 ( 4 ) 采用有限差分法对平面二维水流泥沙数学模型控制方程进行离散，利用a d i 法， 实现在边界处正交的曲线网格系统下，对方程求解。 ( 5 ) 应用得出的数学模型和计算方法，编制计算程序，对水槽实验进行数值模拟并与 实验数据进行比较，验正模型。 1 7 小结 本章对本课题的研究背景、研究现状及意义，并综述了国内外水流泥沙数学模型研 究现状，水流沙数学模型的计算方法做一些介绍。描述了贴边淤积现象及其危害。还对本 文将采用的计算方法作了具体的介绍，阐述了本课题的研究目的和研究内容。 西安理工大学硕士学位论文 2 斜对角笛卡尔方法对边界的拟合 传统的笛卡尔坐标方法对边界的拟合被应用的历史很长，由于其简单易形成计算网 格、正交等特点，在工程界应用得最为普遍。但传统的笛卡尔坐标方法在边界处形成锯齿 状边界，导致边界的拟合精度不高，因此该方法在大尺度计算中应用较多m 。但利用斜 对角笛卡尔方法后，则可显著改善计算的精度m ”。 2 1 斜对