（数学专业论文）二维浅水波方程的高阶有限体积法.pdf

国防科学技术大学研究生院学位论文 摘要 浅水波方程的数值求解越来越受到人们的重视，已成为计算数学的一个热门课题非 结构网格的迅速发展推动了有限体积法的发展，非结构网格上的高阶有限体积法的研究变 得越来越重要，本文主要研究高阶有限体积法在浅水波方程数值求解中的运用，本文的主 要工作有： 考虑二维浅水波方程，在非结构网格上，根据复合型有限体积法的构造思想，对不同 的重构函数，分别选择l a x - f r i e d r i c h 数值流函数和r o e 的r i e m a n n 解算子，构造了几种 高精度的复合型有限体积格式。以一阶单调格式为基础，在每一个三角形单元上对各变量 作单调线性重构函数，构造了具有= 阶精度的t v d 型有限体积法。通过对每一个单元构 造线性插值多项式，构造了具有二阶精度的e n o 型有限体积法。在每一个单元上构造加 权的二次插值多项式，并用两点高斯积分来计算交界面的流通量，得到了具有三阶精度的 w e n o 型有限体积法。我们采用的方法是时间和空间分开处理，时间离散采用多步d r u n g e k u t t a 方法，空问离散是对重构函数，数值流函数，和积分的离散格式分别进行考 虑。 对溃坝问题进行了分析，包括边界的处理，虚拟网格的生成，以及不可容许基点出现 的可能和回避的方法；运用本文构造的几种方法，在规则和非规则两种网格上，对平底溃 坝问题进行了数值模拟，并对结果作了详细的分析和比较，结果表明本文构造的格式确实 具有高精度，高分辨率的性质，特别是w e n o 型有限体积法能很好地解决间断问题；最 后对非平底部溃坝问题作了数值模拟，表明本文构造的方法能解决具有任意几何形状的问 题。 关键词：浅水波方程，非结构网格，t v d 型有限体积法，e n o 型有限体积法，w e n o 型有限 体积法 第i 页 璺堕型兰垫查奎堂翌壅竺堕兰堡笙茎 a b s t r a c t i nt h ep r e s e n tw o r k ，t h en u m e r i c a ls i m u l a t i o no ns h a l l o ww a t e re q u a t i o n sh a s a t t r a c t e dm o r ea t t e n t i o n ，a n dh a sb e c o m eah o tt o p i ci nc o m p u t a t i o n a lm a t h e m a t i c s w i t ht h ed e v e l o p m e n to ft h et e c h n i q u e so fm a k i n gu n s t r u c t u r e dm e s h ，t h er e s e a r c h o ff i n i t ev o l u m em e t h o dw i t hh i g ho r d e ro fa c c u r a c yb e c o m e sm o r ea n dm o r ei m p o r t a n t i nt h i st h e s i s 。w em a i n l yd i s c u s st h eu s i n go ff i n i t ev o l u m em e t h o d w i t hh i g ho r d e r a c c u r a c yi ns h a l l o ww a t e re q u a t i o n s t h ef o l l o w i n gw o r kh a sb e e nd o n e ： o nu n s t r u c t u r e dt r i a n g u l a rm e s h ，b a s i n go nt h ec o n s t r u c t i o no fc o m p o s i t e f i n i t ev o l u m em e t h o d ，b yc h o o s i n gl a x f r i e d r i c hf l u xf u n c t i o na n dr o e sr i e m a n n s o l v e rf o rd i f f e r e n tr e c o n s t r u c t i o nf u n c t i o n s ，w eg i v es e v e r a lf i n i t ev o l u m e s c h e m e sw i t hh i g ho r d e ro fa c c u r a c yf o r2 ds h a l l o ww a t e re q u a t i o n s f i r s t ，b a s e d o nt h em o n o t o n es c h e m eo fo n eo r d e ra c c u r a c y al i m i t e d1 i n e a rr e c o n s t r u c t i o nf o r e a c hv a r i a b l ei sm a d eo ne v e r yt r i a n g u l a rm e s h 。t h u sw eo b t a i nat 卜t y p ef i n i t e v o l u m em e t h o dw i t hs e c o n do r d e ro fa c c u r a c y 。s e c o n d l y 。t h r o u g hc o n s t r u c t i n gl i n e a r i n t e r p o l a t i o nf o rv a r i a b l e so ne v e r ym e s h ，w eb r i n gf o r w a r de n o - t y p ef i n i t ev o l u m e m e t h o dw i t hs e c o n do r d e ro fa c c u r a c y t h i r d l y ，i no r d e rt og e ts c h e m e sw i t hh i g h e r o r d e ro fa c c u r a c y ，as i m p l i f i e dw e n o - t y p es c h e m eisa d v a n c e db yg i v i n gaw e i g h t e d q u a d r a t i ci n t e r p o l a t i o no ne a c hm e s h i na d d i t i o n ，w es e r v es p a c ed i s c r e t i z a t i o n b yc o n s i d e r i n gr e c o n s t r u c t i o n ，f l u xf u n c t i o na n di n t e g r a lf o r m u l a t i o nd i v i d u a l l y a n dt i m ed i s c r e t i z a t i o nb yu s i n gt v dr u n g e k u t t am e t h o d w ea n a l y z ep a r t i a ld a mb r e a kp r o b l e mi nt h ef o l l o w i n ga s p e c t ：t r e a t i n go f b o u n d a r yc o n d i t i o n ，c o n s t r u c t i o no fd u m m ym e s h ，a n dt h ep o s s i b i l i t yo fs t e n c i l t ob en o ta d m i s s i b l ea n dw h ot op r e v e n t w i t ht h es c h e m e sp r o p o s e dh e r e w em a k e n u m e r i p a ls i m u l a t i o na b o u tp a r t i a ld a mb r e a ko nt w ok i n d so fu n s t r u c t u r e dg r i d ( r e g u l a rg r i da n di r r e g u l a rg r i d ) 。a n dc o m p a r et h er e s u l t so fd i f f e r e n ts c h e m e s t h en u m e r i c a lr e s u l t ss h o wt h a tt h ep r o p o s e ds c h e m e sa r ea c c u r a t ea n do fh i g h r e s o l u t i o n ，e s p e c i a l l y ，w e n o - t y p ef i n i t ev o l u m em e t h o dw h i c hh a sb e t t e ra b i l i t y t oc a p t u r ed is c o n t i n u i t i e s a tl a s t ，w em a k en u m e r i c a ls i m u l a t i o na b o u tp a r t i a l d a mb r e a kw i t ht h ev a r y i n gb o t t o mo ni r r e g u l a rg r i d 。w h i c hm e a n st h a tt h ep r o p o s e d s c h e m e sc a ns o l v ep r o b l e m sw i t ha r b i t r a r yg e o m e t r y k e yw o r d s ：s h a ii o ww a t e re q u a t i o n s ，u n s t r u c t u r e dm e s h t v d - t y p ef i n i t ev o i u m e m e t h o d 。e n o - t y p ef i n i t ev o l u m em e t h o d w e n o t y p ef i n i t ev o i u m em e t h o d 第i i 页 国防科学技术大学研究生院学位论文 图目录 图1 1 浅水波的自由面的描述。2 图1 2 格子中心型控制体积示意图：。5 图1 3 格子顶点型控制体积示意图一5 图2 1 局部三角形单元7 图2 2 二维溃坝问题的平面图1 4 图2 3 固壁边界流通量处理的示意图i 1 4 图2 a 非结构三角形网格1 5 图2 5 特殊单元的虚拟网格1 6 图2 7t v d 型f v m 在规则网格上求解二维溃坝问题的等高线图( 仁7 2s ) 1 8 图2 8 分片常数型f v m 在非规则网格上求解二维溃坝问题的水位图( t - 7 2s ) 1 8 图2 9t v d 型f v m 在非规则网格上求解二维溃坝问题的等高线图( 卢7 2s ) 1 9 图3 1 不可容许基点图例2 1 图3 2 正的v o n n e u m a n n 邻域( 左图) 和m o o r e 邻域的一部分( 右图) 2 2 图3 3 局部网格图2 4 图3 4e n o 型f v m 在规则网格上求解二维溃坝问题的数值结果( t - 7 2s ) 2 6 图3 5e n o 型f v m 在非规则网格上求解二维溃坝问题的数值结果( 卢7 2s ) 2 6 图4 i 典型的基点2 9 图4 2 二次加权多项式的基点3 0 图4 3w e n o 格式空间离散的示意图3 1 图4 4w e n o 型f v m 在规则网格上求解二维溃坝问题的数值结果( 户7 2s ) 3 2 图5 1 几种格式在规则网格上求解二维溃坝问题的水位的侧面图( t = 7 2s ) 3 5 图5 2 特殊区域的示意图3 6 图5 3 规则网格上的分片常数型f v m 的含局部极值点的局部图3 8 图5 4 规则网格上的t v d 型f v m 的含局部极值点的局部图3 9 图5 5 规则网格上的e n o 型f v m 的含局部极值点的局部图4 0 图5 6 规则网格上的w e n o 型f v m 的含局部极值点的局部图4 1 图5 7 非规则网格上的分片常数型f v m 的含局部极值点的局部图4 2 图5 8 非规则网格上的t v d 型f v m 的含局部极值点的局部图4 3 图5 9 非规则网格上的e n o 型f v m 的含局部极值点的局部图4 4 图5 1 0 非规则网格上的w e n o 型f v m 的含局部极值点的局部图4 5 图5 1 l 非平底部二维溃坝问题的流场水下底部曲面4 7 图5 1 2 高阶有限体积法在非规则网格上求解非平底部二维溃坝问题舻7 2s ) 4 8 第1 i 页 独创性声明 本人声明所呈交的学位论文是我本人在导师指导下进行的研究工作及取得的研 究成果尽我所知，除了文中特别加以标注和致谢的地方外，论文中不包含其他人已 经发表和撰写过的研究成果，也不包含为获得国防科学技术大学或其它教育机构的学 位或证书而使用过的材料与我一同工作的同志对本研究所做的任何贡献均已在论文 中作了明确的说明并表示谢意 学位论文题目：三丝遗壅遮左焦盟商睑煎匮旦显鎏 学位论文作者签名：塞玺霾 日期：矽衫年厉月加日 学位论文版权使用授权书 本人完全了解国防科学技术大学有关保留，使用学位论文的规定本人授权国 防科学技术大学可以保留并向国家有关部门或机构送交论文的复印件和电子文档，允 许论文被查阅和借阕；可以将学位论文的全部或部分内容编入有关数据库进行检索， 可以采用影印、缩印或扫描等复制手段保存，汇编学位论文 ( 保密学位论文在解密后适用本授权书) 学位论文题日：三丝鸯盔选友焦盟苞睑查匦佳塑鎏 学位论文作者签名：丞聋丕 作者指导教师签名：社 日期：矽口z 年陟月，矿甘 日期：舶5 年i 月f 口日 里堕型兰垫查查堂堕垄竺堕兰垡笙奎 第一章绪论 浅水波方程是各种浅水环境流体运动的数学描述，浅水波方程的数值求解可用来研究 诸如潮汐、涌浪、溃坝、水环境污染扩散等人们关心的实际问题。运用数值方法对浅水波 方程进行研究已成为计算数学，计算流体力学的一个热门话题。近年来，非结构网格有限 体积法的发展非常活跃，它能处理具有任何几何形状的问题，这给浅水波方程求解注入了 强大的动力 1 1 浅水波方程的数学模型及应用 当水域沿特征方向的水平尺度比水的深度大得多时，可称为浅水环境，其流体运动的 规律可用二维浅水波方程来描述。二维浅水波方程的形式有很多种，这里只介绍两种典型 的形式：一种是由连续方程和动量方程组成的形式，另一种是守恒律方程的形式。 设石= “，恐) r 2 ，善“，恐，) 为基准水平面以上的自由面的高度，h b ( x l ，屯) 为基准水 平面以下的深度，h = 善+ 吃，设脚= ( “，屯，) ，v ( 五，而，f ) ) 1 为深度平均水平速度。那么由 连续性方程和动量方程组成的二维浅水波寺程1 3 6 1 ： 等+ v ( 以) = 0 祟+ ( v ) + g v f z a u + r , u + f = o 这里 d 是粘性系数常数，g 是重力加速度，( 掌，) ：c ，= f 衲是底部摩擦函 数，f = ( 正七“- h + 跑一g v ) 是作用力函数，由c o r i o l i s 力，表面风切应力，表面 大气压力和潮汐势等体积力和面积力组成。基于特征方向的g a l e r k i n 方法【3 6 1 是在此形式的 浅水波方程上构造的方法。不过，有限体积法用得比较多的是下述的守恒律形式1 1 0 1 1 3 z l ： 等+ 掣+ 掣：s ( 劬，) ( l d m蕊鲫 积分形式：杀眨删q + 唾f 材r = i l ( s ( z ，乃f ) p q ，f ( = g ) fh1 u = hi ，( c ，) = i h v j u h 鼎+ 丢劝2 一( 帆) m , h - ( h v v ，) ，g ( u ) = v h u r h - ( h v u ，) v 2 i i + 丢砂2 一( 鸭) ( 1 2 ) 第1 页 璺堕型兰垫查奎兰翌茎尘堕堂堡笙奎 s ( x , y ，) = o g h s s , + q v g h s 口+ h c s u ( 1 3 ) 其中，q 是控制体，r 是q 的边界，j i = ( 以，b ) 是r 的单位外法向量，v 分别是流速沿工 m y 方向的分量，g 为重力加速度，f = 譬) 是通过边界r 的通量函数。h = z - z # ，磊是 流扬水下底部曲面，z 是自由面水位，一砂( 鲁) ，一咖( 参 分别是水下底壁作用力沿棚 y 的分力；一荫啄，一曲分别是水下底壁的摩擦力沿x ，j ，方向的分力，是工，y 方 向的摩阻啉：毒警，妨：毒掣，k 是摩擦瓢c ，胁蒯i 。 参数，p 是( e d d yv i s c o s i t y ) 。 图1 i 浅水波的自由面的描述 浅水波问题模型最常见的有五种：一维溃坝模型，二维不可压流模型、后台阶河题、 二维溃坝模型和倾斜水跃问题。当前许多水动力学数值模拟方法的研究中，常常采用二维 溃坝模型，在结构网格和非结构网格的应用，边界条件的数值实验研究中，这个模型也很 有意义。本文主要是对s = 0 ，无黏性，无摩擦，河底无倾斜的二维浅水波方程构造有限 体积格式，对二维溃坝模型进行了数值模拟。 浅水波方程可以结合实际的问题建立模型，从而来模拟各种各样的问题，比如对长江 流域的部分江段进行水流模拟；为了研究工程，工业废水及生活污水的水环境的问题，对 水流过程和相应的污染物输运扩散过程的模拟；运用泥沙数学模型，对水库，河道，湖泊 及河口的河床冲淤变形的模拟等等，浅水波方程的应用范围很广，应用价值很大。 第2 页 、， 锄瓦锄一砂 ，。l，_=_ 舶 鼢 一 一 孓笾一赵 国防科学技术大学研究生院学位论文 1 2 浅水波方程求解方法的研究与进展 浅水波方程也是双曲型守恒律方程；所以，研究浅水波方程求解方法时，有必要回顾 一下双曲型守恒律的有限体积法。 双曲型守恒律的积分形式仅当解是充分光滑时，才能约化成等价的微分形式。流体力 学的基本规律是积分守恒律，它们表示在表面限定的体积内的质量、动量和能量的守恒， 所以，从某种意义上讲积分型的守恒律是更根本的。 有限体积方法就是在物理空间中选定的控制体积上把积分型守恒律方程直接离散的 类数值方法。它综合了有限元法在几何上的灵活性和有限差分法在规定离散流动变量的 灵活性的优点。由于它可以在结构网格和非结构网格上处理问题，从而非常适合处理具有 复杂区域的问题，目前已成为计算流体力学的一种重要方法。 有限体积法从5 0 年代末就开始发展，如g o d u n o v 格式【3 】就是一种从非定常流体力学 的积分守恒律出发进行离散的数值方法。但一般作者认为m c d o n a l d 2 0 1 和m a c c o r m a r k 及 p a u l l y l 2 1 】关于二维非定常e u l e r 方程的数值解的论文是f v m 的先驱。李德元【3 3 1 等的非定常 流体力学数值方法的著作中也给出了有限体积法的数值方法。j a m e s o n 和m a v r u p l i s t 2 3 l 关于 二维e u l e r 方程的非结构三角形网格的有限体积的工作使非结构网格有限体积法引起了广 泛的注意，目前已成为高效数值模拟复杂的高速流动问题的重要方法。国内外目前对非结 构网格有限体积法的研究工作发表的有： 2 0 世纪8 0 年代初期，t v d ( t o t a lv a i l a t i o nd i m i n i s h i n g ) 格式 4 1 提出之后，s a n t l l m 等【l 】在三角形网格上对双曲型守恒律方程提出了迎风型有限体积法；宋松和，李荫藩，陈 矛章等【2 s 1 1 3 4 1 对二维标量双曲型方程发展了满足极值原理的无结构网格有限体积法，将t v d 型有限体积法运用到平面激波反射和空穴流动等问题的计算中。 t v d 格式在光滑鳃的局部极值点处降阶，为了克服这个弱点，h a r t e n 提出了e n o 格 式1 7 】，r a b g r a l l t 8 】在非结构网格上提出e n o 格式和讨论了基点的选择，并对激波管和激波 反射问题作了模拟。这方面的工作另见文献 2 9 1 1 6 1 5 。w e n o 格式是基于e n o 格式发 展起来的，与e n o 格式有相同的可选模版的情况下，在光滑区域，w e n o 格式具有更高 的精度，而在间断附近，仍保有e n o 格式的性质。l i u ，o s h e r 和c h a n t l 7 1 首先提出了一维 空间上的三阶有限体积格式，然后j i a n g ，s h u t z s l 的在多维空间上构造了高阶w e n o 差分 格式。并提出了光滑因子和非线性权的构造。f f i e d r i c h 在文献【8 】的基础上，构造了二维三 角形网格上的w e n o 型有限体积格式。 有限体积格式在解决实际问题的过程中也得到了发展，其中，用f v m 的数值方法求 解浅水波方程推动了f v m 的发展。 浅水波方程也是双曲型守恒律方程，所以双曲型守恒律方程的数值方法都可以用来求 解浅水波方程。目前，已有许多数值方法可以对浅水波问题进行数值模拟，如有限差分法， 特征线法，有限元法，谱分析法和有限体积法g l a i s t e r t l 2 墚用了一种通量分裂技术：汪继 文，刘儒勋提出了一种基于特征方向的g a l e r k i n 方法【3 6 1 ；张理论等给出了谱元素方法离散 的具体过程 3 5 1 。近些年来，有限体积法引起了人们的广泛关注，在模拟二维浅水波问题上 已取得一些成果。 a l c r u d o f 和g a r c i a - n a v a r r o p 运用r o e 的r i e m a n n 解算子对二维浅水波方程提出了 高精度的g o d u n o v 型m u s c l 有限体积格式【l l 】；艮a n a t a s i o u 和c t c h a r t 提出了二阶迎 第3 页 国防科学技术大学研究生院学位论文 风型有限体积法p j ；w a n gj i w e n 和l i ur u x u n 提出了l a x - f f i e d f i c h 格式与l a x w e n d r o f f 格式相结合的复合型有限体积格式【6 1 ；d a m b r o s i 采用了r o e 型m u s c l 方法【1o 】；目前， 处理任意地形的水下地形仍然是没解决的困难问题，这方面的成果并不多，v m c e n tg u i n o t 讨论了对一维r i e m a n n 问题加权的带角度的二维r i e m a n n 问题，并对含源项的浅水波方程 进行了离散洲；e m m a n u e la u d u s s e 和m a r i e - o d i l eb r i s t e a u 在非结构网格上对具有任意地 形的含源项的浅水波方程进行了分析和离散1 2 5 1 ；蔡启福等【3 2 】提出了局部平底化概念，构造 了一种适用于有任意水下地形的高精度，无数值振荡的有限体积法。 目前，将t v d 型有限体积法，e n o 型有限体积方法和w e n o 型有限体积法运用到二 维浅水波方程求解中的文章并不多，本文主要讨论了这些有限体积法的构造和在浅水波方 程中的运用，并模拟了二维溃坝问题，测试了这些格式处理激波问题的能力，证明这些格 式的确具有高精度，高分辨率的性质。 1 3 二维标量双曲型守恒律的非结构网格有限体积法 这一节，我们结合二维标量双曲型方程，简单地介绍非结构网格有限体积法的构造。 有限体积法的构造过程是：首先，对计算区域进行网格剖分，然后对积分方程进行空间离 散和时间离散，其中，空间离散包括构造重构函数，构造数值流函数和选择积分离散公式 三卜方面。 网格生成是有限体积法的重要组成部分， 模拟和数值结果的好坏，是非常关键的因素， 工作量的相当大的部分。 网格生成的质量和效果，能直接关系到数值 而且，网格生成的任务和工作量常常占据总 一种网格称作结构的，假如它的顶点间的连结是有限差分型的，否则就是非结构的。 对于结构网格，结点问的连通性是( i j ) ( 二维问题) ，( i j 均( 三维问题) 。在二维情况下， 若计算区域被剖分成非结构的( 三角形) 网格，我们必须以列表方式保存每个格子的顶点 的编号，每格的边的编号，以及相邻单元等等。 结构网格对简单几何区域能直接生成可利用结构网格特点发展算法。但对复杂的计算 区域生成适当的结构网格是困难的。如果对网格放弃任何内在的结构要求，采用完全非结 构网格，对二维问题用三角形网格，对三维问题用四面体网格，原则上在处理几何上的奇 异性和不规则性上不存在困难。 由于非结构网格舍去了网格节点的结构性限制，易于控制网格单元的大小、形状及网 格点的位置，因此比结构网格具有更大的灵活性，对复杂外形的适应能力非常强。非结构 网格实际上是有限元型的，在解决复杂几何外形问题时，结构网格型方法所遇到的困难， 促使人们尝试发展非结构网格型的方法，使之成为综合有限差分法和有限元法优点的一种 有力的方法。但是非结构网格算法的精度和效率一般不如结构的，改善的途径有多重网格 技术，自适应网格方法和移动网格方法。 守恒型问题的有限体积法是典型的有限体积法。考虑二维标量双曲型守恒律： 珥+ 1 ，【“jj + l g l “”= 0 ( 1 ，4 ) 这里虬厶g 是标量，u 是待求变量。积分型： 第4 页 国防科学技术大学研究生院学位论文 昙( 矩“出砂) + l ( 厂 ) a y - g ( “) 出) = o ( 1 5 ) 将计算区域剖分成网格后，可取该网格的单元本身作为控制体，守恒量的平均值储存 在三角形格子的形心，这样得到的有限体积法叫做格子中心型格式；如图1 2 ，a 以单元本 身为控制体，相邻单元为b ，c ，d ；另一类控制体是由网格周围三角形单元的整体或按统 一的规则或取一部分合并而成，守恒量储存在三角形格子的顶点，这样得到的有限体积格 式叫做格子顶点型，其控制体常用的有如下三种：( a ) 经过每个顶点的所有三角形格子的 并集；( b ) 由围绕每个顶点的周围三角形单元的形心连接构成控制体的边界线；( c ) 由围绕 每个顶点的三角形格子的形心和从该点发出的边的中点相连接而构成控制体的边界线。 b 图1 2 格子中心型控制体积示意图 (a)(” ( c ) 、 图1 3 格子顶点型控制体积示意图 禹散方式：记f = 少+ ，x = “力，考虑格子中心型格式，对( 1 1 ) 在三角形单元爿上积分， 并由g r e e n 公式得 昙( 瓦) 一砉。聃西 ( 1 6 ) 其中h 表示三角形单元爿的面积，k 表示三角形单元4 的第k 条边，瓦是爿的第k 条边 的单位外法向，瓦= 扯出| 爿l 是“在鲋b c 上的网格平均值。设 ( q ，呸) 是满足l i p s c h i t z 连续的单调流函数，即 ( q ，吐) 对q 非减，对吧非增，则对( 1 3 ) 式右边我们进行如下近似： 一塞。f 瓦凼* 一塞h 粪厅( 以。，瓯，。) ；z 饬) ( t 7 ) 其中，女挑。表示在4 中构造的加权高次插值在4 的第k 条边上的第q 个高斯点处的值， 第5 页 彭 致幻表示与一共边的三角形单元上的加权高次插值在该第q 个高斯点处的值，而且， 哆o ，q = l n = i ，即为取一个高斯点，q ( - l ；n = 2 ，即每边取两个高斯点，q = i 1 ； 1 吐2 互。 在构造有限体积法时，人们通常采用时空分离法，时间方向的离散一般采用s h u ”1 所 提出的t v d 型r u n g e - k u t t a 时间离散方法，对于i e l 折t 晖- ，它和稳定的空间离散相结合是非 线性稳定的。例如，一个二阶的t v d 型r u n g e k u t t a 格式是 瓦o k 乃”+ 出z ( 死8 ) ， 瓦一：吾瓦一+ 吾( 瓦m + 出z ( 乃m ) ) o 名 另一个三阶t v d 型r u n g e - k u t t a 格式是 乃o = 瓦4 + 出z ( 瓦“) 。 瓦2 ) = ；瓦”+ ( 死1 ) + 出z ( 瓦1 ) ) ( 1 9 ) 瓦”l = 瓦4 + 寻( 瓦但，+ 血z ( 瓦让，) ) 其中z 是空问离散算子。 在空间离散过程中，关键是重构步，即从单元平均值重构解的表达式，进而得到高斯 积分点的值。对于向量方程组，重构步是对各物理量分别进行的，在下面两章我们将结合 二维浅水波方程组，构造重构函数来得到各种复合型有限体积格式。 本文做了以下工作： 第二章我们对二维浅水波方程组，在非结构三角形网格上，构造了具有二阶精度的 t v d 型有限体积格式，并对二维溃坝问题进行了数值模拟。我们采用的方法是时间和空间 分开处理，时间离散采用t v dr u n g e - k u t t a 方法，空间离散对重构函数，数值流函数，和 积分的离散格式分别进行考虑。 第三章讨论了e n o 型有限体积法的构造，并给出了具有二阶精度的简化e n o 型有限 体积格式，并对二维溃坝问题进行了数值模拟 第四章在e n o 型有限体积法的基础上构造了具有三阶精度的w e n o 型有限体积格 式，并对二维溃坝问题进行了数值模拟。 第五章对结果进行了分析和比较，给出了高阶有限体积法的应用，并对二维非平底溃 坝问题进行了数值模拟。 第6 页 里堕型堂垫查奎兰坚茎竺堕竺垡笙奎 第二章浅水波方程的t v d 型有限体积法 自1 9 8 3 年h a r t e n 提出了t v d 4 1 格式后，出现了一些高分辨率的流体力学数值方法。 特别是，近年来，随着网格技术的发展，尤其是非结构网格的发展，使f v m 进入了一个 快速发展时期。有限体积法结合其他一些数值方法，如u p w i n d ( 迎风) 格式，t v d 格式， e n o 格式，产生了一系列新颖有效的方法。 t v d 格式保持了数值解的单调性，能有效抑制间断附近的振荡。近来，文献 2 8 提出 了一种求解二维标量方程的。t v d 型有限体积法。这一章，我们主要把该方法运用到求解 二维浅水波方程组中，并根据文献 2 7 q h 提出的复合型有限体积法的构造思想，得到了两 种具有高精度的复合型有限体积格式。这些格式基于非结构三角形网格，能够处理具有任 意复杂计算区域的问题。用它们对溃坝问题进行了数值模拟，结果表明，这些方法具有高 精度且稳定。 2 1 浅水波方程及其离散 二维浅水波方程的守恒形式为 q + ( f ( u ) l + ( g ( u ) ) ，= o 其中，u 是守恒物理量，f ( u ) 是x 向通量，g ( u ) 是y 向通量 h u = ih u l 训 f = p 弓曲2 i z f ， g = h v h u v h v 2 + l g h 2 ( 2 1 ) 式中，h 是水的深度，z f 和v 分别是x 和) ，向垂线的平均流速分量，g 为重力加速度。 本文采用格子中心型的有限体积格式，即将三角形单元本身作为控制体，如图2 1 。 d 图2 1 局部三角形单元 在控制体a a b c 上对( 2 1 ) 式积分，得 昙( l 嬲呦) + j 椰c ( f ( u ) ) ，+ ( g ( u ) ) ，p 谚= o 利用g r e e n 公式并化简，得 第7 页 里堕型堂垫查奎兰塑窒生堕堂垡丝塞 ，西o - - u 一+ 网1 l 。( e g ) 砌= 。 ( 2 2 ) 其中，玩盯2 匹去可k u d q 是u 在崩占c 上的网格平均值，面是边界- a a b c 的单位外法 向量。 设k ，k ，乙是削b c 的边，死。，哆。，再。是各边的单位外法向量，则有 昙珏= 一函南( l ( e g ) 讲+ l ( e g ) 刃+ ( 只g ) 酝讲) 从而得到其半离散形式 争= z ( 玩鲇) - 一南 日( u ) ( 如) ) k 。 日( 二脚( k ) ，三脚( ) ) k + 日( 三脚( 瓦。) ，厶。( 。) ) 乙】 ( 2 3 ) ( 2 4 ) 其中，日( ，) 为数值流函数，j ，j 乞为a b ，b c , c a 的中点，。( x ) 表示物 理守恒量u 在a , 4 b c 中的重构函数。 重构函数的精度是由重构多项式逼近各物理量时产生的截断误差的阶来决定的。 例如，分片常数逼近型是把单元, s a b c 的重构函数取为： “( x ) = 死孵 一 这是一阶精度方法，稳定性好，但精度不高。 ， 文献 2 6 1 中关于有限体积法空间离散的精度分析，有如下定理： 定理1 空间离散所用的积分法是q 阶的；对每个固定时刻t 和单元z ，重构函数 是，阶的；u 和数值流函数日是充分光滑的，则空间离散就是m i n ( q ，) 阶的j 从该定理可知高阶有限体积法的空间离散关键在于高阶重构函数的构造。下面两节将 给出构造具有高精度的重构函数的几种方法。 由定理1 可知，空间离散所用的积分法的精度，需要等于或高于重构函数的精度。一 般情况下，( 2 2 ) 式中的线积分，对于线性重构函数，采用一点高斯积分公式，高斯点为各 边的中点，比如对边k ，g = 去( j 乙+ ) ，对二次重构函数，采用两点高斯积分公式，高 二 ，石 斯点为g l = c 磁+ ( 卜c ) x 8 ，g 2 = ( 1 _ c ) k + c k ，c = 毒+ 半，高斯积分权为 l q 2 哆2 j 。 数值流函数日( ，吸) 是( e g ) 引拘近似，这里给出三种形式 第8 页 国防科学技术大学研究生院学位论文 第一种形式，算术平均形式： 日( 巩，) ：f ( u r ) 2 f ( u l ) ，型掣1 菇 一 j 第二种形式，r o e 的r i e m a n n 解算子。 在推导二维浅水波方程的r o e 的r i e m a r m 解算子之前，我们先介绍一下r o e 方法的基 本思想。 考虑r i e m a n n 间断分解问题： u + ( ( u ) ) 。= 0 训加舷：；： u ：( 妒，) r ，( u ) ：( 石( 吩，厶( ) r ，彳( u ) ：亘坚婴 c t u 此间断问题数值求解的难度，主要是因为非线性项f ( u ) 引起的，也就是a ( u ) 的非线性， r o e 方法的根本点在于，利用左右函数的状态u ，u 去构造a ( u ) 的替代矩阵，从而将非线 性问题线性化。r o e 提出的构造线性化替代矩阵的原则有三点： ( 1 ) 相容性：j ( 珥，q ) 一4 ( u ) ，当奶，以- * u ； ( 2 ) 相似性：j ( 奶，u ) ，( u ，一q ) = z 一石； ( 3 ) j ( q ，q ) 的特征向量反，七= 1 ，2 ，m 是线性无关组。 r o e 给出了一种构造j ( 奶，珥) 的途径： ( 1 ) 引入了参向量历作为过渡的量； ( 2 ) 对间断跳跃量【u 】和【八u ) 】用参向量表示： 【，】= 以一q = 秀( q q ) ，( u ) = 厂( u ，) 一，( u ) = e ( 缉一q ) 则夤，q ) = 茚一； ( 3 ) 确定j 的特征值和特征向量，互，五，k = l ，2 ，掰； ( 4 ) 将未知函数跳跃量p 】写成上述的特征向量的线性组合： m = u r u = 瓯五 ( 5 ) 确定数值流函数的特征表示 第9 页 国防科学技术大学研究生院学位论文 兀= 吾( ，) + 厂( u ) 一薹瓯f 互【瓦) 下面我们进行二维浅水波方程的近似r i c m a n n 解算子的推导； 二维浅水波方程：q + ( ，( u ) ) ，+ ( g ( u ) ) ，= o m u = ih u l 【加j i 甜 - 1 f = h u + 2 9 h 2 则j a c o b i a n 矩阵为 。 印，= 掣：i p j = ( b ( u ) ，c ( u ) ) 亓= 掣： 圳 。g = h v h u v 厅v 2 + 昙劝： o ( c 2 一“2 ) 以一哪 哪+ ( c 2 一v 2 ) 哆 其中，f = 厨，再= k ，b ) 。孑的特征值和特征向量为： 以n 2 u n x + v n y“n r 1 u n x + 2 = , 1 x + w 0 一c ，也= “以+ w 0 ，乃= u f x + w 0 + 册 q 怯卜卧恻 下一步，我们寻找p ( u ) ，c ( 【，) ) 露的线性化矩阵i = ( 雪( u ，u ) ，e ( q ，以) ) 疗。 q = 再+ 压 = i “何+ 蚱历) ( 何+ 历) 一历+ 咋压) ( 压+ 厄) 适当地定义“当地声速”孑，可以使线性化矩阵i 与j 的形式相同，这里取 。：坪， 第l o 页 r 伊。叫：一 ， v 警 国防科学技术大学研究生院学位论文 则b ( u ) 和c ( 【，) 的平均线性化替代矩阵雪( u ，q ) 和0 ( u ，以) 分别为 。云c 奶，珥，= 2 弱，。c 奶，以，= 每乏0 ：“vv “u 一“ 身= ( 秀，e ) 拈 o n ln v ( j 2 一存2 ) 心一嘲2 讥+ 鸭鸭 一磊溉+ ( 0 2 一f 2 ) b 饥 纸+ 2 鸭 j 的特征值为互= 讥+ 一声，互= 纸+ ，互= 讥+ 弧+ 5 磊= ；刍 ，五= 高 ，磊= ；刍 可以验证取；= 丛掣得到的复合型线性化替代矩阵i ：( 雪，。) 菇满足r 0 e 的近 似m e m 锄解算子的原则，即 ， ( 1 ) i 满足相容性。因为当q ，u ，一u 时，露哼，哥一v ，o _ c ，则 j = ( 雪，e ) 疗_ j = ( b ，c ) 露 ( 2 ) i 满足相似性。可以验证下式成立 ( ( 豆o ) j i ) ( u u ) = ( 以) 一f ( u ) ，g ( u r ) 一g ( 弘) ) 矗 所以，跳跃量【u 】和 ( f ( u ) ，g ( u ) ) 露 满足下式； i 【u 1 = ( ( 云，e ) 蔬) 【，】 ( f ( u ) ，g ( u ) ) 厅 ( 3 ) i 的特征向量如乞，毛是线性无关的。 3 一以= 幺磊 k = l 其中展开系数为： 第1 1 页 里堕型堂垫查奎兰塑塑竺堕堂垡丝奎 磊= 等一去( ( 砌) 以+ ( 自v ) 嘭一( 讥+ 鸭) 幽) 幺= 吉( ( ( v ) 一幽) 吃一( ( 粤) 一幽) b ) 岛= 等+ 去( ( 砌) 以+ ( v ) b 一( 讥+ ) 幽) 从而得到r o e 的r i e m a n n 解算子的逼近形式5 日( 叽，) = 兰( ( f ( ) + f ( ) ，g ( ) + g ( 巩) ) 疗一喜瓯| 互k ) 第三种形式，l a x f f i e d f i c h 数值流函数。 日( 眈，) = 寺( p ( ) + f c u 。) ，g ( ) + g ( ) ) 亓一口( 一眈) ) 其中a 是j a c o b i a n 矩阵! f ( ! ：! ! ：塑的复合型线性化替代矩阵i 的谱半径。 采用不同形式的流函数和g r 不r 同的重构方法，就可以得至q 各种各样的复合型有限体积格 式，如分片常数逼近型l a x f r i e d r i c h 有限体积格式，分片常数逼近型r o e 有限体积格式， t v d 型l a x f f i e d f i c h 有限体积格式，t v d 型