（水文学及水资源专业论文）河道二维水流模拟高精度算法研究.pdf

摘要 随着计算机软硬件技术的迅猛发展，人们对水流数值模拟的广度和深度提 出了更高的要求。本文在前人研究的基础上，从河道二维水流数值模拟的实际应 用需求出发，对j 下交曲线坐标系下的河道二维水流数值模拟方法进行了进一步研 究。在传统有限差分离散的基础上，本文对数值模拟的性能和精度提高做了一些 探索，研究了的河道二维数值模拟的时间项外推、算子迭代和并行计算的求解方 法。 时间项外推法是在同一时刻对两倍时间步长的数值解与一倍时白j 步长的数 值解进行外推组合，消除离散的低阶时问项误差，得到更高精度的结果。 算子迭代是在同一时间步长内，用算子迭代求解离散的方程，逐步消除离 散格式引入的误差，最终达到同全隐非分步算法同等的精度。 河道二维水流并行计算方法是对河道沿程进行区域分解，重叠相邻区域部 分网格，通过不同计算单元对各分段河道单独的迭代计算，以重叠区域计算结果 是否相同为收敛准则，实现河道二维水流的并行计算。 关键词：河道二维水流、数值模拟、外推法、算子迭代、并行计算 a b s t r a c t w i t ht h er a p i dd e v e l o p m e n to fs o f t w a r ea n dh a r d w a r et e c h n o l o g yi nc o m p u t e r , p e o p l er a i s ean e w n e e di nb i a d t ha n dd e p t hf o rn u m e r i c a ls i m u l a t i o no f w a t e rf l o w b a s e do nt h er e s u l t so fp r e v i o u ss t u d i e s ，t h i sp a p e rs a t i s f i e st h ea c t u a la p p l i c a t i o n n e e df o rt h e2 dw a t e rf l o wn u m e r i c a ls i m u l a t i o ni nc h a n n e la n dd o e st h ef u r t h e rs t u d y f o rt h a ti nt h eo r t h o g o n a lc u r v i l i n e a rc o o r d i n a t e s t h i sp a p e rd o e ss o m ee x p l o r i n ga t t h ea s p e c to fi m p r o v i n gp e r f o r m a n c ea n da c c u r a c yo fn u m e r i c a ls i m u l a t i o no l lt h e b a s i so ft r a d i t i o n a lf i n i t ed i f f e r e n c 七m e t h o d i th a s3p a r t s ：t i m ei t e me x t r a p o l a t i o n m e t h o d ，a r i t h m e t i co p e r a t o ri t e r a t i v em e t h o da n dp a r a l l e lc o m p u t i n g t h et i m ei t e me x t r a p o l a t i o nm e t h o di st or e c a l c u l a t e st h er e s u l to fl i n e a r c o m b i n a t i o nf o rd o u b l ea n do n et i m ep a c ev a l u ea tt h es a m et i m e t h ec o m b i n a t i o n r e s u l th a se l i m i n a t e dt h ee r r o ro f l o wr a n dt i m ei t e ma n dh a sm o r ea c c u r a t ev a l u e t h ea r i t h m e t i co p e r a t o ri t e r a t i v em e t h o ds o l v e st h ed i s p e r s ee q u a t i o nb yi t e r a t i v e m e t h o di nt h es a l n et i m es t e pc i r c l ea n dg r a d u a l l ye l i m i n a t et h ee r r o ro fd i s p e r s i n g m e t h o d a st h er e s u l t , i th a st h es a m ea c c u r a c ya sa l li m p l i c i tn o s p l i td i s p e r s e m e t h o d 。 t h ep a r a l l e lc o m p u t i n gm e t h o df o rc h a n n e l2 dw a t e rf l o wi st os p l i tt h ew h o l e r i v e rc h a n n e li n t os o m ep a r t sa l o n gt h er i v e r a tt h en e i g h b o rr e g i o n ，p a r tm e s h e sh a s b e e no v e r l a p p e da sac o u p l i n gj o i n t i nt h ec a l c u l a t i o np r o c e s s ，t h ed i f f e r e n tp a r t s h a v eb e e ns o l v e di nt h ep a r a l l e lw a yb yd i f f e r e n tc o m p u t i n gu n i t s t h es i n g l es t e p s o l u t i o ni st e r m i n a t e dw h e nt h ed i f f e r e n c ei sl e s st h a nt h es p e c i f i e dv a l u ea to v e r l a p r e g i o n b yt h i sm e a n s ，t h ep a r a l l e lc o m p u t i n gm e t h o dh a sb e e nr e a l i z e d k e y w o r d s ：c h a n n e l2 dw a t e rf l o w 、n u m e r i c a ls i m u l a t i o n 、e x t r a p o l a t i o n m e t h o d 、a r i t h m e t i co p e r a t o ri t e r a t i v em e t h o d 、p a r a l l e lc o m p u t i n g 学位论文独创性声明： 本人所呈交的学位论文是我个人在导师指导下进行的研究工作 及取得的研究成果。尽我所知，除了文中特别加以标注和致谢的地方 外，论文中不包含其他人已经发表或撰写过的研究成果。与我一同工 作的同事对本研究所做的任何贡献均已在论文中作了明确的说明并 表示了谢意。如不实，本人负全部责任。 论文作者( 签名) ：壁惶玺 学位论文使用授权说明 擗凶月山e l ， 河海大学、中国科学技术信息研究所、国家图书馆、中国学术期 刊( 光盘版) 电子杂志社有权保留本人所送交学位论文的复印件或电 子文档，可以采用影印、缩印或其他复制手段保存论文。本人电子文 档的内容和纸质论文的内容相一致。除在保密期内的保密论文外，允 许论文被查阅和借阅。论文全部或部分内容的公布( 包括刊登) 授权河 海大学研究生院办理。 论文作者( 签名) ：塑 1 1 研究的目的和意义 第一章绪论 随着社会的发展，人们对河流的利用日趋多样，修建了很多水利工程，这 些水利工程与河道水流相互影响，形成一些水利、环境、生态等方面的问题。 因此许多与河流相关工程项目的设计及规划都需要建立相应河道的水流模型来 研究兴建工程对河道水流的影响，以及评估河道水流变化对其它工程如护岸、 桥梁等的影响。河道工程对于河道水流的影响分析主要采用模型模拟，常采用 物理模型或数学模型来模拟实际的河道水流。 自1 9 世纪弗劳德( f r o u d e ) 奠定物理模型的理论基础以来，物理模型在 很长一段时间几乎是验证流体力学理论和研究工程水力学问题的唯一手段i l 】。 物理模型是根据研究问题的目的，以及研究问题关注的重点不同，采取不同的 相似准则设计比尺模型，构建实体模型对实际河道的水流进行模拟。物理模型 由于其直观性和物理概念明确，一直受到水利工程界的重视，但其模拟程度往 往受模拟试验手段以及比尺关系的限制，不能同时满足所有的力学相似准则。 另外，物理模型还具有研究费用高，试验周期长的缺点。 近几十年来，由于计算机科学的迅猛发展，计算机对水流运动的模拟能力 和模拟精确程度也不断地得到提高，河道水流数值模拟逐渐成为人们认识水流 运动、研究水流运动的种重要的方法，通过计算机程序模拟真实、直观地反 映了河道水流的运动情况。与此同时，水动力学、数值计算方法和计算机技术 相结合，形成了一门新的交叉学科计算水力学。从理论到实际应用，河道水流 的数值模拟技术已经曰趋完善。相对物理实体模型而言，数值模拟试验具有容 易改变方案、初边值条件以及各种模型参数，可以重复进行，不受试验环境和 场地的影响等优点，因而在水利、航运等工程的规划设计中应用也越来越多。 在实际科研项目中，由于数学模型容易与其它研究手段有机结合，相互验证， 越来越多的科研人员把数值模拟的结果作为基本资料与其它各种模型进行耦 合，用来研究更复杂的相关问题。在空气动力学界，航空航天器设计所需的二 维及三维气动力学数值模拟正在逐步取代费时费钱的风洞试验，可以预见水流 数值模拟在实际水利工程项目和科学研究中将会得到越来越多的应用1 2 j 。 词道= 堆水流模拟离精度算法研究 常见的河道二维水流数值模拟应用有：分滞蓄泻区的防洪规划、洪水漫滩、 河道整治、河床演变分析、渍坝决堤模拟、闸堰及枢纽布置、河流输沙计算、 污染物输运扩散模拟、河口潮汐计算、盐水入侵研究等。对这些相关问题实现 精确数值模拟的前提是首先实现水流的精确模拟，因此水流数值模拟的精度非 常重要，是其它功能数值模型实现精确模拟的基础。此外，河道二维水流模型 的计算量都较大，高精度的算法一般都需要更大的计算量和更复杂的计算过程， 因而对河道水流模型求解效率和精度的提高会促进水质、泥沙、温排水等其它 专业模型求解效率和精度的提高。 随着社会经济发展、人们对问题关注程度以及数值模拟需求不断提高，河 道水流数值模拟在广度和深度上不断发展，从原来的一维模拟发展到二维、三 维模拟，从原来的单一河道模拟发展到河网、流域模拟，这些需求都会带来计 算机模拟的计算量成倍增加。与此同时，计算机架构也向超线程、多c p u 、多 核方向发展【3 】1 4 】【5 l ，这也需要研究适应新型计算机硬件系统结构的数值求解算 法，充分利用计算机系统资源，更高效地实现对水流的数值模拟。 1 2 河道二维数值模拟方法介绍 1 2 1 浅水方程 实际的天然水体流动在空间上表现是三维的，描述三维水流运动的数学方 程是n s ( n a v i o - s t o k e s ) 方程。一般来讲，直接数值求解n s 方程比较困难， 在有些情况下甚至是不可能的。为了满足实际工程应用需要，人们一般对三维 的n s 方程进行简化得到浅水方程，进行数值求解。 天然河道水流满足浅水流动的四个条件【l 】：( 1 ) 有自由表面；( 2 ) 以重力 为主要驱动力，以水流与固体边界之间及水流内部的摩阻力为主要耗散力，有 时还存在水面气压场、风应力及地转柯氏力等的作用；( 3 ) 水平流速沿垂线近 似均匀分布，不必考虑实际存在的对数或指数形式的垂线流速分布：( 4 ) 水平 运动尺度远大于垂直运动尺度，垂向流速及垂向加速度可忽略，从而水压力接 近静压分布。 由于一般的河道水体流动平面尺度远大于垂直尺度，因此对三维的n s 方 第一i 爿r 论 程在垂直方向积分取平均，得到沿水深平均的二维水流方程，再利用上述4 个 条件，进行简化得到的浅水方程就可以模拟实际河道的平面二维水流运动。 1 2 2 常用数值求解方法 对浅水方程的数值求解根据差分方程逼近微分方程方式不同，有不同的求 解方法。以下简要介绍水流数值模拟中最常用的的三种求解方法：有限差分法、 有限单元法和有限体积法。 1 2 2 1 有限差分法( f o m ) 有限差分法是直接用差分方程逼近微分方程。首先建立规则的差分网格， 然后求解网格节点上的水位、流速等未知变量的数值，利用求出的网格节点上 的水位、流速等来代替整个求解域内的水位、流速等。根据离散求解方法的不 同，有限差分法又可以分为以下几种方法： ( 1 ) 直接差分法 直接差分法就是直接用差商代替微分方程中的导数，主要形式有前差分、 后差分和中心差分三种。根据求解又可以分为显式差分格式和隐式差分格式。 显式差分格式可以通过初始值直接求出未知节点上的值，不需要联立方程组求 解，然而显式差分格式收敛性和稳定性较差，时间步长不能过大。隐式差分格 式稳定性和收敛性好，可以取较长的时阔步长，求解精度高，在数值计算中有 广泛应用。1 9 9 1 年，王船海、程文辉在文献7 1 推导出全隐非分步算法，消除了 分步误差，当c o u r a n t 数达到4 0 0 时仍能很好模拟水流流态变化较剧烈区域。 2 0 0 2 年，王船海等在文献【1 4 】采用迭代求解，将有限差分求解推广n - - - 维树状河 网水流模拟，取得了良好的应用效果。2 0 0 6 年，向小华在文献【15 】中将一维环状 河网追赶法推广n - - 维，提出了通用河网二维水流模拟方法。实际的数值计算 应用表明，这种全隐非分步算法可以很好的解决河道类水流的模拟。 ( 2 ) 特征线法 特征线法是离散求解基于沿特征线成立的特征方程，而不是求解利用普通 空间坐标中的描述水流运动的原始方程。特征方程反映了双曲型问题中信息沿 河遭= 维水漶模拟离精度算法研究 特征线传播的性质，因而算法符合水流的物理机制【i j 。不足之处是，特征方程 常为非守恒形式，用差分法离散特征方程时会带来守恒误差，另外处理多维问 题时特征方程非常复杂。韦春霞等人采用算子分裂法将二维水流问题分裂为沿 x 和y 方向的一维流动，并推导出相应的特征线方程和相容方程l l6 】，对瞬间部 分溃坝时二维洪水波的演进过程进行了数值模拟，得到的数值模拟结果表明， 这种方法有较强的溃坝涌波捕捉能力，展示了溃坝波的运动特性，为溃坝引起 的灾害预测和分析及防灾减灾等工程问题提供一种有效的数值模拟方法。 ( 3 ) 分步法 分步法是将包含具有几个不同物理机制的复合过程分解成机制单一的物 理过程组合，不同的物理过程对应不同的微分算子，根据不同微分算子的不同 性质，分别选取适宜的计算方法，最后得出综合的结果。分步算法优点是物理 概念清楚，对于兼含对流项、扩散项和源项的方程，可以分别采用适宜的方法 求解。缺点是在地形复杂、流态变化剧烈的地方会出现较大的分步误差【“。魏 文礼、孙广才采用物理分步法【7 引，将原方程拆分成三部分：即对流过程、扩 散过程和传播过程，对横门出海航道局部区域流场的数值模拟结果表明，计算 结果与实测资料符合较好，这种算法适宜计算河口、湖泊、海洋等宽浅一类水 域的大体积水流数值模拟。 ( 4 ) a d i 法 a d i 法用于求解二维水流问题是一种特殊的分步法，它将一个时间步长分 成两个半步长，前半步长在x ，y 方向的偏导项分别用隐、显格式离散，后半 步长在x ，y 方向的偏导项分别用显、隐格式离散。采用a d i 离散的优点是离 散求解的方程比较简单，计算工作量比较小，用于二维计算是无条件稳定的。 a d i 法收敛快，在一条网格线上，无论有多少节点，来自网格末端的边界条件 信息可以通过t d m a 法一次传入内部区域。a d i 法其时间步长不宜取得太大， 否则其模拟效果会违背实际的水流规律。研究表明r 玎，当时间步长超过5 倍的 c o u r a n t 数时会出现严重的轴化现象( p o l i a r i z a t i o n ) 或称为一维化效应。 有限差分法建立在经典的数学逼近理论的基础上，理论成熟且简单易为人 们接受，处理效率较高，是发展比较成熟的方法，在水流数值模拟方面得到广 泛的应用。有限差分法的缺点是，差分方程组必须建立在规则的矩形网格上， 4 第一幸铺论 而规则的矩形网格不能很好地适应复杂的河道边界。在实际应用中，工程界广 泛采用正交曲线坐标变换1 9 h 矧，将实际复杂的区域边界变换到规则的矩形计算 区域边界，有效解决了边界拟合问题。本文第二章将对正交曲线坐标系下的贴 体网格生成技术做具体的介绍，这罩从略。 1 2 2 2 限单元法( f e m ) 有限单元法将计算区域划分为许多互不重叠的单元，一般为三角形或四边 形，在各单元内部选定适当的节点，将微分方程中的变量和其导数表示成由节 点值与其所选用的插值函数组成的线性表达式，应用变分法或加权余量法，生 成有限元离散方程，将各单元的离散方程合成为总体方程组，结合初值、边值 条件得到近似解。有限元方法在数学上适于求解椭圆型方程的边值问题，不适 于求解以对流为主的输运问题。在二维浅水方程中扩散过程属于椭圆型方程， 而对流过程属于双曲型方程，使用传统的有限元方法容易产生较大的数值弥散 和数值振荡而使计算失稳或者产生非物理解。1 9 8 2 年d o u g l a s r u s s e l l 提出特征 有限元方法【】，可以解决有限元方法求解对流占优的过程中的数值弥散问题， 成为解决有限元方法求解对流占优的过程中的数值弥散问题的有效方法。这种 方法在引入动坐标系后，将传统的欧拉方法和拉格朗日方法有机地结合起来， 通过控制动坐标系的移动速度，使对流扩散方程中的对流项不出现或者使之变 得非常小，此时在动坐标系中看到的将始终是扩散过程，从而使得对流作用和 扩散作用可以分别进行计算。采用隐式有限元求解最终会形成一个大型非线性 方程组，求解比较困难而且计算量大【1 2 】。为了解决运行效率的问题，国内外许 多学者提出了不同的显式有限元方法，通过一定的技术处理使隐式有限元的代 数方程组的系数矩阵转变为对角矩阵，使每个网格节点都存在一个独立方程， 无需联立求解方程组。从而大大减少了计算量。但显式有限元由于受稳定性限 制，时间步长较隐式有限元要小很多，在公式推导、程序编制等方面也比较复 杂【2 5 l 。周杰在文献【2 6 1 中采用分步有限元方法离散二维浅水方程，并初步探索了 行指标矩阵压缩存储稀疏矩阵的方法和效率以及预条件双共扼梯度法在求解有 限元方法中形成的大型线性方程组中的效率，取得了满意的结果，在定程度 上缓解了有限元在非恒定流计算中存在的耗内存大、耗时多的问题。 河建二维才0 流模拟商精度算法研究 1 2 2 3 有限体积法( f v l v l ) 有限体积法，又称为有限容积法、控制体积法，其以守恒型的浅水方程为 出发点，将计算区域划分为一系列不重复的单元( 称为控制体积) ，以网格点上 的因变量为未知数，并假设其在网格点之间的分布规律，对每个控制体的控制 方程进行积分离散来构造离散方程。有限体积法有两种导出方式，一种是控制 容积积分法，另一种是控制容积平衡法。不管采用哪种方式导出的离散化方程， 都描述了各个有限控制容积内物理量的守恒，所以有限体积法是守恒定律的一 种最自然的表现形式。有限体积法适用于任意类型的单元网格，无论网格尺度 大小以及网格的形状，只要在单元边上与相邻单元估计的通量是一致的，就能 保证方法的守恒性。有限体积法各项近似都含有明确的物理意义，是吸收了有 限元分片近似的思想以及有限差分方法的思想而形成的高精度算法，是目前最 为流行的水流数值模拟方法，但对于非结构网格有限体积法的离散求解以显式 居多，稳定性与时间步长等都受到较大限制【l3 1 。 1 2 3 网格技术 在计算流体力学中，从有限差分法、有限元法到有限体积法，所有的方法 都离不开网格。网格质量的好坏直接影响数值结果的精度，甚至影响数值计算 的成败。数值计算所采用的网格主要有两种：结构网格和非结构网格。 1 2 3 1 结构网格 结构网格，指的是一类具有一定的分布特征，可以用相应的行列关系来顺 序描述的网格。常用的结构网格有：矩形网格、曲线网格以及块结构网格等。 矩形网格是f d m 或f v m 最常用的网格，它不对原始的控制方程( 直角坐标) 作任何坐标变换，对河道的岸边线采用简单的阶梯形网格来近似复杂的河岸边 界，边界拟合性比较差，而且还会带来很大的内存浪费，在实际河道水流数值 模拟中应用较少。正交曲线网格以其良好的贴体性、较高的计算精度以及比较 节约内存而在水流数值模拟中得到了广泛的应用。目前常用的生成网格的坐标 系大致可分为2 类：一类是代数坐标系( a c s ) ，即直接应用已知的代数方程， 第一章爿障分 如折线、二阶曲线等，网格生成较方便，但贴体性差、人为因素较大，仅适用 于规则边界问题，对于实际流场的计算，有其不适之处；另一类是微分坐标系 ( d c s ) ，即通过求解特定的微分方程来生成坐标系，其坐标生成较前类复杂， 但可以密切贴合边界，比较合理地分配网格，在某些情况下能做到正交，这不 仅能简化工程问题的控制方程，还可以减少误差来源【2 9 】。t h o m p s o n 提出的利 用p o i s s o n 方程生成贴体曲线网格的t h o m p s o n 方法【2 7 】应用最广泛。文献【2 8 】【捌 采用该方法通过调整源项控制网格的分布，在梯度大的区域人工或自动加密网 格，很方便地实现了满足精度的正交曲线网格。 1 2 3 2 非结构网格 非结构网格可以采用任意形状的网格单元，单元边的数目也没有限制，弥 补了结构化网格对复杂区域及连通区域的网格剖分困难的欠缺。非结构网格能 很好地模拟自然边界及水下地形，有利于对复杂边界调节的实现；便于控制指 定区域网格的疏密，易作修改及适应性调整。但是非结构网格应用也带来一系 列需要解决的问题：单元网格排列不规则，须建立相应的数据结构存储单元格 信息：控制方程离散形成的代数方程系数矩阵是高度稀疏的非对角型矩阵，需 要寻求合适的存储方式及高效的求解方法；对于非结构网格实现隐式求解比较 困难。非结构网格原则上可应用于包括f d m 在内的任何类型的数值方法，但 非结构网格在f e m 和f v m 中算法更成熟，应用更广泛。非结构的网格生成方 法主要有d e l a u n a y 三角化方法【3 1 3 2 1 和阵面推进法【3 3 l ，可以采用开源软件 e 勰y m e s h 或商业软件g a m b i t 、s m s 等生成。 1 2 。4 动边界处理 在实际数值模拟中，宽浅河流中的滩地和江心洲的水边线会随着水位的变 化不断改变，形成滩地和江心洲上水边线的移动边界。目前对动边界处理的方 法主要利6 】【2 l 】：冻结法、窄缝法、切削法、干湿法等。冻结法，根据网格单元 中心处水深判断网格单元是否露出水面，对露滩单元取糙率系数为一极大的数， 使单元四周的流速都为趋于零的微小量，使露滩单元水位在计算时被冻结不变。 河道= 维水流模拟高精度算法研究 窄缝法，假想在岸滩的每个网格上存在一条很窄的缝隙，它的深度和岸滩前的 水深一致，其宽度是水位相对于单元中心高程距离的函数，根据水量平衡原理 将窄缝内的水量平铺到岸滩上，成为化引水深，相当于把岸滩前的水域延伸到 岸滩内，从而可以把计算网格点布置在岸滩上。切削法，与冻结法一样判断单 元是否露出水面，但并不冻结水位，而是引入一个富裕水深来保证计算过程的 完整和稳定，相当于将原始地形切削降低，而一旦判断实际水深大于富裕水深 时，即恢复原始地形。干湿法，根据计算点的水深及其相邻网格处的水深和水 位值判别该计算点的干湿，而使露滩节点流速沿此方向为零。这些动边界处理 方法大同小异，基本都可以很好地处理计算过程中河道水流边界变化的问题。 1 3 本文主要研究内容 目前对河道水流数值模拟常用的方法有：a d i 法、分步法、全隐非分步法。 尽管采用全隐非分步法求解精度高，没有分步误差，但求解每一时问步需要做 大量系数矩阵的求逆运算，计算量较大，所需的时间较长，程序实现也相当复 杂。因此本文研究了时间项外推、算子迭代以及区域分解的数值计算方法，在 数值方法的求解精度、加速收敛方面做了一些有益的尝试，通过数值试验，分 析了求解加速及提高离散精度的效果。 l 、采用时间网格的外推算法，对变量q ”1 分别采用时间步长为f 和出2 离散求解，将计算的结果对应到同一时刻，按q ”1 = 2 q 阻2 y “一q ( f y “进行 线性组合，实现计算的q ”1 达到较高的d ( f 2 ) 精度。 2 、对动量方程和连续方程采用迭代求解。在同一时间步长内将对流项、 摩阻项等采用不同的离散格式进行算子迭代，通过迭代消除离散格式引入的误 差，在给定收敛准则下，实现算子迭代算法和全隐非分步算法达到一致的精度。 3 、采用区域分解的并行方法，将河道按上下游分段，在连接处重叠部分 节点，上下游河段段由不同的线程计算。在同一时间步长内，上游段的下边界 由下游段的计算结果给定，下游段的上边界由上游段的计算结果给定，反复迭 代直至重叠部分上下游两段的计算值相等，从而实现对整个河段的并行求解算 法。 第二章河道二维水流控制方程 2 1 河道二维水流控制方程 描述河道平面二维水流运动是用浅水方程，浅水方程有守恒型和非守恒型 两种形式，方程的变量的表示也有采用原始变量表示和采用守恒变量表示两种。 本文采用的方程都是用原始变量表示的浅水方程。 2 1 1 非守恒型方程 直角坐标系下非守恒型二维浅水方程的为： 鱼+ l i l f 塑+ 鱼1 + ”o h + v o h ：o o t l 反砂j 魂 砂 扣罢+y面au+百2+v2+g瓦oh咖三hf(型ax2+可02(uh)gu ( 2 - )石棚夏面+ 历厂+ g 瓦一少2 一一十可j 旺j 罢+ v 塞+ 百2 + v 2 + g 面o h 咖哥掣+ 掣西+ “瓦+ ”石+ 矿+ g 面+ 扣2 i 1 1 y + 1 芦j 式中，“、v 分别为石、y 方向的垂线平均流速：z 为水位；h 为水深， h = z z d ，z d 为河底高程；厂为柯氏力系数；c 为谢才系数；占为运动粘滞系 数。 2 1 2 守恒型方程 直角坐标系下守恒型二维浅水方程【2 l 】的为： 鱼+ 型+ 型：o o t o x 却 掣+ 掣+ 掣+ 孕拓u 2 + v 2 + 曲塞一钟= 拶+ 拶 ( 2 2 ， 掣+ 掣+ 掣+ 孕3 厅u 2 + v 2 + 磅嘣= 拶+ 等) 式中，“、1 分别为瓢j ，方向的垂线平均流速；z 为水位；h 为水深， 9 河海大掌硕士掌位论文词道= 维水漶相撤离精度算法研究 h = z z d ，z d 为河底高程；，为柯氏力系数；c 为谢才系数；f 为运动粘滞系 数。 从流体微元的角度看，守恒型和非守恒型控制方程是等价的，都是物理守 恒定律的数学表示。但对于一个具有一定大小的计算水体单元，这两种形式的 控制方程在数值离散上是有区别的，非守恒型控制方程便于对采用不同格式离 散形成的方程进行稳定性等理论分析，而守恒型控制方程更能保持控制体内物 理量守恒的性质，便于处理对流项非线性引起的问题，而且也方便采用非矩形 网格离散，如在非结构网格下的方程离散。在数值求解中一般采用守恒型的浅 水方程。 2 2 正交贴体坐标变换 2 2 1 曲线坐标生成 对直角坐标系下的浅水方程采用有限差分法求解，由于采用矩形网格计算 区域，计算区域的边界往往不能贴合实际河道边界，严重地影响了模拟精度。 同时，由于实际河道的方向不一定与直角坐标系中的坐标轴平行，使得计算网 格区域远大于河道地形区域，这也增加了很多无用的存储量及计算量。 t h o m p s o n 等人在1 9 7 4 年提出了通过求解椭圆型偏微分方程组来生成曲线 坐标系的网格生成方法l l9 1 ，把平面直角坐标系下的坐标转换到曲线坐标系下的 贴体网格坐标。它以在计算区域边界上的坐标线必须与边界线密切贴合为坐标 系的建立原则，同时网格的实际尺寸可按需要加以变化，从而获得了广泛的应 用。由于t h o m p s o n 等人所提出的网格生成方法，不能保证所生成的网格正交。 而正交的曲线族出现于多种物理现象中，如果以这些物理现象为依托建立坐标 网格，则可以保证所建立的网格正交。w i l l e m s e 等人根据流体力学势流理论中 流线与势线正交的原则【2 0 】，导出下列一组方程： + + o 0 l o ( 2 3 ) = = 、，、j 西一砂铆一砂 一k一置 ，。，。、 a y a y 骘面却| 遗一k一k 甄吖虱 ：l l = - l t 河道= 照水漶控制方程 式中：髟= 止乏了霹刀丽，_ 、儿、y , 1 分别是工、y 对乒叩的偏 导数。如图2 1 所示， 工、y 是物理平面坐标，参r 是计算平面坐标。 f 一徽影坐标平面 图2 i 坐标变换示意图1 5 1 由( 2 3 ) 可以进一步写成 对应t h o m p s o n 等人提出的椭圆型偏微分方程： a 2 掌c a 2 亭 紊+ 矿。尸 垒+ 塑：d 玉2 矿。 比较( 2 4 ) 与( 2 5 ) 得： ( 2 4 ) ( 2 5 ) 砷础恤仳 。一砂a一砂管一砂塑砂 一 一 曲砷啦舢 a一苏a一缸鸳瓦塑缸 = = 孵矿塑驴 + + 丝酽塑酽 河递= 维水流橇拟高期- 度算法研究 ：一芒晏( 1 n k ) 一譬昙o n k ) 缸缸、7 却却、 7 = 一粤兰( 1 n k ) 一譬晏( 1 n k ) 缸缸、7 却却、 ( 2 6 ) 式中p 、q 为调节网格疏密的调节因子。( 2 4 ) 式的逆变换方程为 + ：哆+ 瓯2 ( 2 7 ) 【a y 菇+ + ，2 慨+ 纵) = o 。 式中，口= + y ：，= + ，= 鼍) ，一y f ，p q 分别由式( 2 6 ) 给出。在计算过程中，如果计算域四周边界线上网格点选择不当，直接求解( 2 7 ) 式生成的曲线网格的正交性精度不高。特别对于河道，网格的长宽比相差很大， 单靠收缩因子p ，q 是很难使曲线网格点达到足够精度的正交性。文献采用不 把四个边界上的网格点都固定死，而是允许其中一条岸边线上的网格点根据正 交性沿边界滑动，这样不仅提高了正交性的精度，迭代收敛速度也得到加快。 2 2 2 曲线坐标系下的水流方程 在正交贴体坐标系下，用虬，u 分别表示沿f 和刁方向的流速m s ) ，则原 直角坐标系下的流速“、v 与、v 之间的关系可以表示为： ：! ：兰量二兰墨! 兰 y x r , 一y q x v = 竺：墨兰2 二兰兰! 兰 x y q x n y 式中，g f 、岛分别为曲线网格的长度b ) 和宽度如) 。 ( 2 8 ) 考虑f 和玎方向粘滞系数不同，并用曼宁公式c ：h l 1 6 代替谢才系数，将 甩 ( 2 8 ) 代入( 2 1 ) ，化简后就得到正交曲线坐标系下平面二维水流控制方程【1 4 1 ： 第= 河道= 擅水漶控制方程 害+ 爿喜k “。曲+ 南k u 明= 。 堡+尝詈+昙茜+竖等一了2虿09,atj + 学一a + 舌g 妻=g ：粥g h8 n a q j8 0 | 必 去毒幢爿) 一言南幢动 心。 鲁+尝等+毒等+坠鼍一芷j玺orl+gn2_v,届厂+vfj+ 旭+ 言高=a g ；粥g | 8 n 越 移一g ，a q 去南幢一) 一去南眩功 其中， 一= 壤+ 扣叫 口= 采( v 岛) 一争影) 式中，三g ，表示曲线网格的面积b 2 ) ，乓，易分别为墨y 方向的滑动粘 滞系数。 通过j 下交曲线坐标变换，将复杂的实际不规则地形区域变换到规则的矩形 计算区域，将求解x - y 坐标系下变量z ，”，v 变换为求解f 一玎坐标系下的变量 z ，1 。 2 3 交错网格 对控制方程进行差分离散首先需要将求解变量布置到网格点上，常用的数 值离散网格有同位网格和交错网格两种。如图2 2 所示，同位网格把流速“、v 以及水位z 布置在同一个网格点，这样就有可能出现图2 2 中这种不合理的情 况：动量方程离散后的水位梯度在任何节点处均为0 ，和均匀的水位压力场一 致，不能正确反映重力项对水流的作用，但实际上在空间石，】，两个方向都存在 振荡的水位梯度。 河道= 维水流模拟高精度算法研究 r r t 5 s o 7 jt 0 7 5 r r 1 t 07 58 0t 5 s o p 7 5t 5& ot 5 一 8 07 5o z j ”、v ，z 节点 图2 2 同位网格示意图 为了解决上面提到离散后的动量方程不能反映水位梯度场的问题，研究者 提出使用交错网格来布置求解变量，把流速“、v 以及水位z 布置在不同的网格 点。如图2 3 所示，这样对动量方程的离散就能正确反映水位梯度场。交错网 格引入了三套网格系统，使得方程的求解变得更加复杂，增加了计算量，但由 于其可以成功地解决了水位梯度离散的问题，得到了广泛的应用。 iv 流速节点 _ “流速节点oz 水位节点 图2 3 交错网格示意图 2 4 控制方程求解方法 复杂的实际地形区域变换到曲线坐标系下规则的计算区域后，由于其计算 区域规则，原来在平面直角坐标系下规则区域的各种离散方法都可以采用，离 散形成的线性方程组求解也可以根据其特点采用直接求解或迭代求解。曲线坐 标系下的水流方程，由于做了坐标变换，采用的差分网格是曲线坐标网格，水 流控制方程与原方程相比多了几项，变得稍微复杂一些。水流方程求解计算量 1 4 第= 章河道= 维水流控制方程 有一点增加，但离散形成的矩阵规模较原来变小很多，增加的网格宽度和长度 计算是一次性的，占用内存较原来是节省的，总体计算量与原来相比也是节省 的。一般而言，直接对非守恒型方程进行差分离散求解，会出现水量不平衡现 象【埘。曲线网格下的连续方程采用守恒型，离散后守恒性很容易得到满足。动 量方程由于坐标是曲线坐标，采用非守恒型方程。一般对于河道水流数值模拟 而言，只要满足连续方程的守恒性就可以满足数值模拟精度要求。 对二维水流运动方程求解可以采用显式格式或隐式格式，但由于显式格式 稳定性差，时间步长不能取得很长，计算一个水流数值模拟方案需要很长的机 时，在工程实际计算中更多地是应用隐式格式或显隐混合格式。全隐非分步格 式的精度高，可以取很长的时间步长，但由于全隐非分步格式的数值离散求解 比较复杂，计算一个时间步长需要做大量的系数矩阵乘法及求逆运算，耗用的 计算机时比较长，工程界在实际应用中对数值求解精度和计算时间上做了一个 折衷，常采用分步法( 破开算子法) 来简化数值离散求解过程。 常用隐式求解方法有：a d i 法、分步法、全隐非分步算法。 2 4 ia d i 法( 交替方向法) p e a c e m a n ，r a c h f o r d 和d o u g l a s 于1 9 5 5 年提出了求解二维水流方程的a d i 法【3 9 】。将a d i 法应用到正交曲线坐标系下的浅水方程离散求解，就是把时间步 长f 分成2 个半步长等，前半时间步长在参，7 方向的偏导项分别用隐格式和 显格式离散，后半时自j 步长在量，7 方向的偏导项分别用隐格式和显格式离散。 控制方程采用( 2 9 ) 式，为了书写方便，以后出项的孓玎方向的流速。，v 都用扯，v 标记。离散求解过程如下： 在疗寸n + 去时间层： 河道= 维水瘴模撤离精a t 算法研究 l 出 2 西 l 锄 2 西 la v 2a + l 生( 。 j8 毕4 “加 + + g 8 “加 + 一一+ g a 砌) 一；品k 伸) 旦：型! ：! ：旦鱼：一三盟一竺坠 g f 鸳g 口却，却 + 【去毒睡一) 一i 1 瓦a 眩口) j g n 2 v 以2 + v 2 va v w 露“2a 洛g 。8 q j 粥 j 8 q 譬南一言善叫 + 导冬+ 夕 j8 e ( 2 1 0 ) 一羔喜一i f , g 口0 r l 对( 2 1 0 ) 式中三个方程等号左边采用隐式离散，等号右边采用显式离散。 先联立离散求解。动量方程和z 连续方程，采用追赶法可以解出z 4 气和。4 气，再 将z ”i 和”2 代入离散的v 动量方程求出y “气。 l 。二 l 在n + 三一n + 1 时间层：离散求解方法与前一步类似，先联立离散求解v 动 量方程和z 连续方程，用前半时问层的计算结果z ”j 、“”i 和y 4 气作为离散初始 值，采用追赶法可以解出z “和v ”，再将z ”和v ”1 代入离散的v 动量方程求 出u ”。 通过上面介绍的求解过程，对乒r l 两个方向的显隐交替求解，完成二维水 流方程一个整时间步长的隐式求解。 2 4 2 分步法 a d i 法及其演变出来的其它变形格式是隐式的，虽然理论上是证明是无条 件稳定的格式，但实际的计算表明，时间步长不允许取太大。当c o u r a n t 数超 过5 时由于相位误差太大等原因，计算可能不稳定，或者误差过大，从而使a d i 法精度不高，不能模拟水流的真实流态”。因此，b e n q u e 等提出了分型算法 4 0 l ， 这是一种物理概念上的分步格式，将水流控制方程按照物理过程特征分解为： 对流过程、扩散过程和波动过程。文科2 2 培合河道水流长宽比例较大的特性， 将水流控制方程分为：对流过程和传播过程，针对不同过程采用合适的求解方 法。文献旧j 介绍了分步求解的方法，过程如下： 1 6 第= 章河递= 维水漶彗e 制方程 ( 1 ) 对流过程，相应控制方程为 这是双曲型万程，- 3 - 以用矩阵娼魁珐求解。 ( 2 ) 扩散和波动过程，相应控制方程为 睁采脚南圳一o ia,+亟a,120tj 号薏+ 鱼弓等唼毒咖壶毒幢柚一品q 1 2 l ，鸳j j l ，； 岛鸳。岛鸳、。岛却” 睁篙玺or20t j j + 学+ 嚣= 咖+ 言南蚪言委l 鸳 岛却 。 岛却“岛筐” 先离散动量方程，分别建立、v 与z 的线性关系，再将线性关系代入连续 方程，采用矩阵追赶法求解出z ，代入线性关系求出“、v 。 2 4 3 全隐非分步算法 采用分步法的主要原因是由于直接全隐式离散浅水方程形成的方程组很复 杂，求解这样复杂的方程组需要做大量的矩阵运算，计算工作量很大，然而采 用分步法会引入分步误差。在一般流态中，分步法的分步误差较小，但在流态 变化较剧烈的情况下，当址取得较大时，其分步误差就很大，分步法计算的结 果就变得不合理。针对分步误差，文献 7 1 提出了全隐式的非分步算法。该方法 利用河宽与河长相比是小尺度的特性，对河道水流方程进行全隐式离散，构建 水位流速的递推关系，对动量方程和连续方程进行联立祸合求解。由于离散形 成系数矩阵都是小尺度的矩阵，推导出的矩阵相乘以及矩阵求逆运算也都是小 尺度矩阵运算，因此计算量相对直接全隐式离散求解较小。文献【n 提出了求解 河道二维水流方程的高离散精度全隐式矩阵追赶法，消除了分步误差，适用于 流态空间变化剧烈的河段。求解过程如下： 对控制方程( 2 1 ) 分别在”、v 、z 节点上隐式离散，将上游边界条件代入， l 互( o o = = 加一却c若一却 上岛上岛 加一西加一凿 。 旦旦 勿一a孤一a加一a 河琏= 维水流捐拟高精度算法研究 推导出一组上断面“、，、z 与下断面“、v 、z 的递推关系，最后把下游边界代入 递推关系，回代即可求解出所有断面的“、v 、z 。具体求解算法参考文献【1 4 】介 绍，这里不再详述。 2 5 动边界处理 采用贴体正交曲线坐标变换后，曲线网格与河道岸线以及河势或主流线之 间拟合较好，也能够较好的适应河道边界。在实际水流数值模拟中，由于水位 随着计算过程在不断变化，因而导致河道水流的计算边界也是在不断移动。处 理动边界前文提及已有多种方法，本文采用冻结法。在每步计算开始时，先比 较断面水位节点水位z 与河底高程耐之间的大小，然后根据比较结果改变相邻 计算流速网格点的糙率以及冻结或解冻计算水位节点水位。具体方法为：如果 z 耐，再判断相邻计算流速网格点是否为冻结点，是冻结点则将计算流速网 格点糙率设为正常值解冻该网格点，是正常节点则按正常计算过程计算；如果 z 0 ”2 0 n 2 _ = 0 v 墨2 。 0 其中：k = 0 5v 2 0 湘“= 0 10 0v 。o + 1 2 _ + l 0 ”2 0 h + 1 2 m l 20 “墨“2 m l 0 “2 0 h 2 _ = 0 “墨h o 将( 3 6 ) “动量方程按( 3 1 7 ) 离散可得： u a 2