（计算数学专业论文）浅水波方程有限体积方法的数值比较.pdf

摘要 摘要 本文讨论有限体积方法求解一维浅水波方程组时各种数值流通量和限制器 的数值表现差异。具体的处理方法有两种：其一是状态外插值方法，即对于单元 界面处的变量值进行某种捅值重构使得重构后的点值( 外捅值) 有比界面左右单 元平均值更高的精度，比如v a i ll e e r 的的m u s c l 方法即属于此类；其二是流通量 限制器方法，即对低阶耗散格式引入反扩散机制以达到高阶。这两种做法都能达 n 2 阶精度，亦能有效控制解的大梯度区域的数值振荡，t 开在不同的流通量和限 制器的组合下会有不同程度的数值表现。 关键词：有限体积方法浅水波方程组限制器数值流通量 a b s t r a c t a b s t r a e t i nt h i st h e s i s ，w ei n v e s t i g a t et h ed i f f e r e n c e su s i n gv a r i o u sc o m b i n a t i o n so ft h en u m e f i c a lf l u x e sa n dl i m i t e r sw h e nn u m e r i c a l l ys o l v i n gt h eo n e d i m e n s i o n a ls y s t e mo f s h a l l o ww r a t e re q u a t i o n s w ew i l lu s et w od i f f e r e n ta p p r o a c h e st oo b t a i no u rf i n i t e v o l u m es c h e m e ：t h ef i r s tw a yi st oa p p l yt h em e t h o do fe x t r a p o l a t i o no fs t a t ev a r i a b l e s ， t h a ti st os a y , w ew i l le x t r a p o l a t et h es t a t ev a r i a b l e sa c c o r d i n gt ot h ec e l la v e r a g e ss o a st og e th i g ho r d e ra c c u r a c y 鲥dv a l u e sa tt h ec e l li n t e r f a c e sw h i c hw ew i l lu s ei n t h en u m e r i c a lf l u x e s t h ef a m o u ss c h e m eo fv a nl e c r sm u s c l ( m o n o t o t i cu p s t r e a m c e n t e r e ds c h e m ef o rc o n s e r v a t i o nl a w s ) i so ft h i sc l a s s t h e0 t h e ri st ou s et h ef l u x l i m i t e rm e t h o dw h o s e g u i d e l i n ei st oa d da na n t i - d i f f u s i o nt e r mt ot h el o wo r d e rf l u xt o o b t a i nh i g ho r d e ra c c u r a c yi nt h es m o o t hr e g i o no ft h es o l u t i o n w | eb o t hg e ts e c o n d o r d e ra c c u r a c yu n d e rt h ec o n d i t i o no fs m o o t ht e s tp r o b l e m sa n dt h et w om e t h o d sc a n a l s os u p p r e s sn o n p h y s i c a lo s c i l l a t i o n si nt h er e g i o no fh u g eg r a d i e n t ，b u tt h ep e r f o r - m a n t e si nt h ev a r i o u sc o m b i n a t i o n sa r ed i f f e r e n t k e yw o r d s f i n i t ev o l u m e m e t h o d s h a l l o ww a t e re q u a t i o n sl i m i t e rn u m e r i c a lf l u x 南开大学学位论文版权使用授权书 本人完全了解南开大学关于收集、保存、使用学位论文的规定， 同意如下各项内容：按照学校要求提交学位论文的印刷本和电子版 本；学校有权保存学位论文的印刷本和电子版，并采用影印、缩印、 扫描、数字化或其它手段保存论文；学校有权提供目录检索以及提供 本学位论文全文或者部分的阅览服务；学校有权按有关规定向国家有 关部门或者机构送交论文的复印件和电子版；在不以赢利为目的的前 提下，学校可以适当复制论文的部分或全部内容用于学术活动。 学位论文作者签名：京稚彪 2 驴拜夕月勿日 经指导教师同意，本学位论文属于保密，在年解密后适用 本授权书。 各密级的最长保密年限及书写格式规定如下： 南开大学学位论文原创性声明 本人郑重声明：所呈交的学位论文，是本人在导师指导下，进行研究工作 所取得的成果。除文中已经注明引用的内容外，本学位论文的研究成果不包含 任何他人创作的、已公开发表或者没有公开发表的作品的内容。对本论文所涉 及的研究工作做出贡献的其他个人和集体，均已在文中以明确方式标明。本学 位论文原创性声明的法律责任由本人承担。 学位论文作者签名： 农聪彪 2 d 湃夕月多护日 第章引言 第一章引言 本文将针对几种简单实用的数值流通量以及几种常用限制器之间的组合搭 配进行讨论，寻求一种最佳搭配。本文采用一维浅水波方程做为模型方程。这 方面的工作还是很有明显的使用价值的，如文献【1 8 】比较了r k d g 方法结合不同 数值流通量以及不同限制器函数之间的搭配的一些数值结果，文献【8 】就各种求 解r i e m a n np r o b l e m 的方法所得到的数值流通量进行了比较研究。 针对目前高分辨牢有限体积方法的两种主要做法：一种做法是对于单元界 面处的状态变量值进行某种重构使得重构后的点值( 外捅值) 有比界面左右单元 状态平均值更高的精度，这个精度的提高就体现在我们将这些重构点值代入到 通常的数值流通量中可以达到高阶( 当然一般我们讨论的是发展问题，这个高阶 的得到还需要有相应的时间离散方法) ，这个思想首先由v a nl e e r 2 3 提出，即他 的m u s c l 格式。另外一种做法就是流通量限制器做法 1 5 ，2 ，5 】。这种方法的思想 在于在解的光滑区域我们要追求高阶精度，t 开是在大梯度区域我们不再一味追求 这个目标( 众所周知：l a x w e n d r o f f 格式是高阶线性格式，它在大梯度区域会产生 数值振荡) 而是着眼于解的单调，变差不增。为了这两个目标同时能实现，我们 需要利用解的光滑性来自适应地选择高阶还是低阶方法( 这样最后得到的格式是 非线性格式，即使对于线性方程) 。本文中的高阶格式采用由l w 格式，记相应数 值流通量为f 日；低阶格式由相应的低阶数值流通量得到，记为f l 。最后我们得 到的格式的数值流通量为，y d = f z + 咖( 口) ( f 日一f l ) ，其中0 为解的光滑性因 子或者叫做光滑性指标，西为限制器函数。 文章内容安排如下：第一章我们对浅水波方程的数学表达式及物理背景做出 介绍，同时给出有限体积方法的简单介绍及其发展，最后介绍本文主要工作；接下 来的第二章给出本文所采用的有限体积方法的“b u i l d i n gb l o c k s ”即数值流通量以 及限制函数的表达式；第三章给出采用文中所说方法的数值精度验证以及针对几 个浅水方程不同初值问题的一些数值效果的比较。在最后一章给出本文的结论。 1 1 一维浅水波方程数学模型 本节简要给出不带源项的一维浅水波方程的守恒形式的推导【l6 4 ，1 5 。 1 第一。章引言 假定水流在一横截面积为矩形的均匀水渠中流动，并且还假定在垂直于水渠 底平面的截面上水体速度分量为0 ，在任意垂直于水渠底平面的截面水流只在平 行于水渠底平面的方向速度分1 i t t # 零且为常数u ( x ，) 。以上假设在浅水波动问题 中是成立的：即水波振幅与水波波长相比起来较小时我们的假设是合理的。在上 述假设之下，再加上水的不可压缩假设，也就是p 为常数，我们用h ( x ，) 表示水高， 设水渠宽为常数枷。在每一时刻t ，区f 3 x m ，z 2 】的水的总质量为：ep w h ( x ，t ) 如， 在z 点t 时刻质量流通量为p w h ( x ，t ) u ( x ，t ) 全f ( x ，) 。根据质量变化率= x l 处质量流 通量一z 2 处质量流通量，于是有 o，啦 豢p w h ( x ，t ) d x = ，( t ) 一，( t ) “。j x l 上式对于任何区问k l ，z 2 】以及在任何时刻t 都成立，约掉w ，厕得微分形式的描述 质量守恒的方程： + ( 舰k = 0 ( 1 1 ) 在每。一时刻t ，位于陋1 ，z 2 】区问的动量为譬p w h ( x ，t ) n ( z ，t ) 如，动量流通量 为p w h ( x ，t ) u 2 ( z ，t ) + p 全r ( x ，t ) ，其中p = p g h 2 表示静水压力。根据动量变化 率1 处动量流通量一z 2 处动量流通量得： 耗2 删叫m 州肛巾小) - r 一 上式对于任何区问p l ，z 2 】以及在任何时刻t 都成立，约掉钳，p 我们得到微分形式的 描述动量守恒的方程： ( 珏) t + ( u 2 + 互1 ，h 2 ) 霉= 0 ( 1 2 ) ( 1 1 ) 与( 1 2 ) n j 以写为以下的方程组形式： u t + f ( u ) 户训。+ k 枷u h2 础2 = 0 ( 1 3 ) 其中h 代表水体高度，让表示水流速度。上述双曲微分方程组中第一个方程表示质 量守恒，第二个方程表示动量守恒。对于光滑解，我们也可将浅水波方程组写为 以下拟线性形式或叫做非守恒形式( 1 4 ) 。但是对于间断解( 即初值间断或者即使 初值光滑解也会在某时刻会出现不光滑) 的情况，我们必须采用守恒形式来进行 数值计算，因为此时采用非守恒形式的话计算出来的问断解是不正确的，不满足 跳跃条件，因为此时和控制体的物理守恒律不一致【4 ，1 6 】。 2 第。章引言 + b ：9 九三雠卜 m 4 ， 令 a c u ，= 一u 2 ：9 那么a 的左右特征向量可取为： r = 乱三c 乱三c 全t r t 您，l = 磊1 一u 。让+ 一c c ，1 1 全 。1 。1 在上述特征向量下a - - y 以对角化为： a = 苫0 = 牡一0c 钍+ 0c = 三a r 其中c = 、而为当地声速。 1 2 有限体积方法简介 有限体积法( f i n i t ev o l u m em e t h o d ，简记为f v m ) 又称为控制体积法，是八十 年代以来发展起来的一种新型的偏微分方程的空间离散方法，具有独特的优点， 目前已成为偏微分方程问题和计算流体力学问题的数值计算、数值模拟中一个重 要的方法。有限体积法设计原理很容易于理解。由它得到的离散方程，就是状态 变量在控制体积中的守恒原理，这就如同微分方程表示状态变量在无限小的控制 体积中的守恒原理一样。众所周知，对于求解双曲型守恒方程而言，格式的守恒 性是一个很重要的标准。基于有限体积法得出的离散方程，要求因变量的积分守 恒对任意一组控制体积都得到满足，因而对于整个计算区域自然守恒性也得到满 足。这是有限体积法吸引人的优点。有限体积法即使在较粗网格情况下，也能显 示出准确的积分守恒。有限体积方法构造守恒数值格式，很简便直观，即有限体 积在求解守恒方程方面有先天的优势。 有限体积方法的基本设计原理是：将计算区域划分为一系列不重复的控制 体积单元( 规则或者不规则) ，并使每个网格节点周围有一个控制体积，然后认为 状态变量在每一个控制体积内为常数( 这个常数即为变量在单元的平均值) 。将 待求解的微分方程在空问的每一个控制体积积分，便得出一组关于单元平均值 的半离散常微分方程组，进而再在时间上离散，比如显式离散( 在实际中应用较 3 第章引言 多) ，譬如简单的e u l e r 向前，显式r u n g e k u t t a 方法等或者隐式时间离散方法比如 说c r a n k n i c o l s o n 时问离散格式。 有限体积法可视作有限元法和有限差分法的杂交产物，它融合了两者的一些 优点。有限元法必须假定状态变量近似解在网格点之间的变化规律。有限差分法 只考虑网格点上的点值而不考虑整体函数表达式在网格点之间如何变化。因而有 限体积方法具有有限元方法在对求解区域的适应性强的优点，即对于无结构网格 也可以方便进行离散化计算求解。 比起最近兴起的间断有限元方法【7 1 而言，有限体积方法仍具有一定优势即： 需要求解的未知量或者自由度数量少要很多，因而工作量要小。特别是近几十年 来发展起来的高分辨牢有限体积方法，有效提高了早期有限体积方法的精度方面 的不足。这些高分辨率方法巧妙融合了迎风思想与使用限制器的做法，使得高分 辨率格式对于求解间断问题可以有所谓总变差不增的性质，而且在解光滑的区域 有高的精度，从而使数值解很好的收敛于熵解。 f v m 的原始雏形从6 0 年代开始出现，如h a r l o w 等人提出的p i c il 】，m a c 1 2 】 和f l i c i o i 方法，7 0 年代的m a c c o r m a r k 1 7 和p a t a n k a r 3 1 等人的f v m 思想8 0 年 代以来，网格生成技术、特别是无结构网格生成技术的发展，给f v m 注入了强 大的生命力，使f v m 进入了一个快速发展时期同时，f v m 结合其他一些数值 方法，如有限元方法，u p w i n d ( 迎风) 格式，r u n g e k u t t a 方法，还有新发展起来 的m u s c l 2 3 1 、r o e 1 9 、t v d 1 3 1 、e n o 1 4 等方法，产生了一系列新颖有效的方 法。从上述文献中我们可发现，有限体积方法的发展也是伴随人们对于数值格式 精度的不断追求而不断完善进步的，如何处理高精度与数值振荡之间的矛盾就成 为人们努力的方向。 在这个发展的过程中，人们发现在解间断的地方仍要求高精度是没有什么太 大好处，因为此时在间断附近高阶格式所赖以存在的数学基础：t a y l o r , 展开不再成 立。相反在这些区域使用低阶精度格式则更有利于数值格式的稳定发展。即在解 光滑的区域我们应该追求高阶收敛率，而在解不光滑的地方我们更需要追求的是 数值解的不振荡。于是h a r t e n 提出d 和h i g hr e s o l u t i o n 等概念 1 3 l ，在具体的实 现这些提法的手段就是采片j 限制器 2 1 1 。一种实现手段是在分片常数基础上进行 ( 高阶) 线性重构( 本文所讨论仅限于线性重构) 得到单元界面处的高阶近似。重 构的过程我们由斜率限制器( 所产生的格式为非线性格式【l ，2 5 】) 保证其t v d 性 质。另一种做法就是将低阶有限体积格式与高阶有限体积格式作适当组合，可以 4 第，章引言 认为是在低阶高耗散格式基础上添加适当的反扩散项，进而最后得到形式上的高 分辨率数值格式。 本文将沿着这个思路进行一些数值比较的工作，试图给出求解浅水波方程的 最佳匹配方案。 5 第二章有限体积方法以及数值实现 第二章有限体积方法以及数值实现 本章将简要介绍如何用有限体积方法的具体实现过程，并给出常见的一些数 值流通量和限制器。我们讨论如下5 种数值流通量和4 种限制策略，即u ：数值流 通量 1 5 。1 6 1 ，l l f 数值流通量【1 5 ，1 6 1 ，r o e 数值流通量【1 9 ，2 2 ，1 6 ，7 】，h l l e 数值 流通量【9 ，1 6 】和f o r c e 数值流通量【2 2 ，1 8 ，以及i i l i n m o d 限制器，s u p e r b e e 限制 器，m c 限制器和v a nl e e l 限制器，这方面许多文献都有相应介绍，比如 1 5 ，2 2 ，4 ， 1 6 ，5 】。 2 1 有限体积法的基本框架 下而给出求解守恒方程的有限体积方法基本实现框架。设双曲守恒律方程为 阢+ f ( u ) 霉= 0 ， 其中是欲求解的标量( 或向量) 函数。 首先将求解区域剖分为一系列单元厶= ( 墨一1 2 ，x i + l 2 ) 的并。在每个单元 ： 秋分微分方程，并引进数值流通量可得 掣= - ( 2 ( 沪) ( 2 i ) 其中 嘴) = 击z 哪) 如 表示阢在单元厶的平均值， f + l 2 c t ) = f ( u ( z i l 2 ，t ) ，u ( z 0 1 2 ，t ) ) 是所谓的数值流通量，它仅仅依赖边界两侧的取值。u ( z i l 2 ，孟) 表示( 戤+ 1 2 ，t ) 的 左极限，u 扛善l ，2 t ) 表示u ( z 件l 2 ，t ) 的右极限。 在格式的具体设计中，空间精度紧密依赖于数值流通量e + l 2 ( ) 的设计选取。 回顾历史上数值流通量的设计我们发现提高精度并有效抑制非物理振荡的方法 主要有两种：一是提高f ( ，) 的两个变元的精度，即采取所谓斜率重构做法或叫 6 第二章有限体积方法以及数值实现 状态值外捅方法。二是采用对于单元甲均值有高空间精度的数值流通量( 如l a x - w e n d r o f f ) 加上反扩散修正( 可使在解的大梯度区域退化为低阶精度从而抑制非 物理振荡) 最终的流通量形式为：，v d = f l + ( p ) ( f 日一f l ) ，其中f h ，f l 为 高阶低阶数值流通量，0 为解的光滑性网子或者叫做光滑性指标，为流通量限制 器函数。无论哪种做法都涉及两个重要概念：数值流通量以及限制器。下面我们 将分别介绍之。 2 2 数值流通量 数值流通量在有限体积方法中占据重要的地位，它保证格式具有良好的数值 效果和局部的守恒性。以下介绍的各种数值流通量都是一阶的单调型数值流通量， 使得有限体积格式可以很好地跨越激波区，不产生虚假的数值振荡。 1 l a x f r i d r i c h s ( l f ) l o c a ll a x f r i d r i c h s ( l l f ) 数值流通量 对于标量问题，l f 数值流通量的数值粘性是所有常用流通量中最人的。l l f 数 值流通量中引入的的粘性比l f 数值流通量要小，但是这种局部的特征值上界已经 能够提供足够的粘性以保持数值解的不振荡。 l f 数值流通量和l l f 数值流通量分别定义为： 1 f 仆( 魄，) = 去【r ( u l ) + f ( u r ) 一q ( 一魄) 1 二 1 f 乩f ( 吮，u r ) = 去【r ( u l ) + f ( ) 一q ( 一吮) 】 二 其中的q 是j a c o b i 矩阵a ( u ) 的特征值的最大绝对值。对于l f 数值流通量，u 取遍整 个值域；对于l l f 数值流通量，u 取遍观与界定的区间。 2 l l f 熵修正的r o e 数值流通量 这里采用了r o e 线性化方法对非线性r i e m a n n 问题作线性化以便近似求解。这 种近似导致的近似线性问题只考虑了间断，没有对实际的r i e m a n n 问题之解可能 出现的稀疏波情况进行处理，这里我们采用u f 熵修正来增加粘性，使得线性 化矩阵的特征值与0 的距离不至于太近，从而使数值解收敛至粘性解( 熵解) 。具 体做法：首先对a ( c 厂) 作r o e 线性化处理：如果我们已经有 阢，l = 1 ) 以及 令厶= ( 戤一v 2 ，x i + t 2 ) ，阢= ，t l n 在每一个单元界面z h 2 ，令： 1 h = 言( 氐一l + 乜) 厶 7 第二章有限体积方法以及数值实现 以及 厕i _ l + x 冗i t t i 也一l + 乜 那么a ( u ) 在x i 一1 2 线性化为： 五一t ，2 = 一铲：而刍 易知上述矩阵的特征值以及相应的特征向量： 娃1 2 = 岔一磊礓l 2 = 鑫+ 8 吐。，z = 矗二爸 ，象。，2 = 鑫三 】 l l l 2 。麦【鑫+ 爸一l1 ，l 乙l 2 。麦f 一( 舀一刁1 】 其中a = 折 最后我们得到带熵修正的流通量为： f f ( e l ) ， i f 0 a i _ l 2 ， f r 。e l l f ( u l ，) = 一工f ( 吮，) ，i f 礓l 2 垂0 礓l 2 ， 【f ( u r ) ，i f 疑l 2 0 3 h l l e 数值流通量 这坐采用另外一种近似求解r i e m a n n 问题的方法，即将以相应速度传播的所 有的波近似为以最小与最大波速( 这两个速度需要作适当估计) 传播的两个波。 从而对r i e m a n n 问题求解带来方便。采用这种方法最后得到的数值流通量为t 一 黯等一， 其中的s 一和s + 定义如下： 8 一= m i n ( m i n ( a l l ，, x t l v 2 ) ) p 矿= m a x ( m a x ( a ，鼍l 2 ) ) p 一 8 i f 0 s 一， i f s 一05s + ， i fs + 0 第二章有限体积方法以及数值实现 这里的骘仍然为j a c o b i 锄f ，( ) 的第p 个特征值，礓l 2 为r 0 e 线性化矩阵五一1 2 的 第p 个特征值 4 一阶中, b ( f o r c e ) 数值流通量 它是l f 数值流通量与r i c h t m y e r 数值流通量的算术平均，即 f f o n c e ( 观，u n ) = 去( f 卯( 观，) + 胪( 吮，) ) 其qj r i c h t m y e r 数值流通量定义为 产( 吮，) - - - f ( 矿) ，旷= 三( 吮+ 一是( f ( ) ) 一( f ( 观) ) ) 换言之，一阶中，i 二, ( f o r c e ) 数值流通量是对l f 数值流通量的某种修正，降低了数 值粘性。 2 3 限制器函数 为提高有限体积方法的数值解精度，我们试图找到高阶的数值流通量，或者 将分片常数解提升为高阶元( 如线性函数) 。但是这样的处理均不可回避地遇到高 阶与数值振荡，单调却低阶的矛盾，为此计算工作者提出一种自动调整的策略，将 数值格式在不同的光滑区域采用不同的数值格式。总体上有两种主要实现途径； 对于线性守恒律问题，两者是等价的。下面分别简要介绍之。 ( 一) 数值流通量限制 传统做法是简单地取流通最中的吮，为网格界面左右单元的。、f 均值以一l ，阢。 为提高数值解精度，我们对上述数值流通量进行某种合适的非线性修正，其实质 是自适应地增加一个修正项【6 ，1 5 ，1 6 ，5 】，使数值流通量在光滑区域中变为某个 高阶数值流通量，而在激波区域是单调数值流通量。本文的高阶数值流通量取 为l a x w e n d r o f f 数值流通量，它在光滑区域具有二阶精度。 这样处理的机制是在两种数值流通量中获取某种平衡，以避免高阶格式在激 波区域可能产生强烈的数值振荡的，最终目的是使数值解满足所谓的t v d 性质， 即： t y ( 矿+ 1 ) st v ( q n ) 对于离散网函数q 而言，其总变差( t o t a lv a r i a t i o n ) 定义为 = t v ( q ) = i 酝一q t 1 1 = 一 9 第二章有限体积方法以及数值实现 它在某种程度上刻画了q 的振荡程度。 这个实现过程要借用于t v d 流通量限制器【6 ，1 5 ，1 6 ，5 】，它们的图像均落在所 谓的s w e b y 区域( 2 玢t v d ，图2 1 中折线i l l i n m o d 和折线s u p e r b e e 之间的区域) 中， 使得通过限制得到的数值结果可以达到形式_ k 水j 2 阶精度并且保证数值解的总变 差不增。 图2 1 l e f t ：m i n m o d ，s u p e r b e ea n dm cl i m i t e r s ；r i g h t ：m i n m o d ，s u p e r b e ea n dv a nl e e r l i m i t e r s 其中 1 m i n m o d ( o rm i n b e e ) ： 该限制器图形( 分段线性函数) 落在s w e b y 区域的下边缘。具体表达式： 砂( p ) = m i n m o d ( 1 ，口) m i n m o d ( a , b ) = 襄 i f 川a 0 2 s u p e r b e e ： 该限制器图形( 分段线性函数) 落在s w e b y l x 域的上边缘。具体表达式： ( 秽) = m a x ( o ，r a i n ( 1 ，2 0 ) ，r a i n ( 2 ，口) ) 3 m c ( m o n o t o n i z e dc e n t r a l - d i f f e r e n c e5 m i t e r ) ： 该限制器图形( 分段线性函数) 落在s w e b y 区域。具体表达式： 4 v a l l l e 宅r ： ( p ) = m a x ( o ，m i n ( ( 1 + o ) 2 ，2 ，2 p ) ) 1 0 第二章有限体积方法以及数值实现 该限制器图形( 分段光滑函数) 也落在s w e b y 区域。具体表达式： 础) = 端 以上表达式中的0 为数值解光滑性指标。本文我们不取逐分量即h ，h u 分别对 其计算光滑性。而是采用将u = 阢一阢一1 分解到特征方向吐i 2 , 吃1 2 上去，对 得到的展开系数在界面z i - 1 2 处与该点的迎风方向相应展开系数信息作比较来确 定解的光滑性。 ( 二) 斜率修正限制器 第二种常见的处理方式是采用状态外插值方法 2 3 ，1 4 & p 提高流通量作用到 的两个变元吮，的精度。为达t u 2 阶，我们需要根据当前已知的各单元平均值进 行满足t v d 性质的线性重构： 观：以+ 等吼一l ，：阢一等吼 其中以为单元厶之重构斜率。这个重构过程目标是一方而在光滑区域提高精度，另 一个方面在间断( 大梯度) 区域抑制数值振荡。对应于上述的流通量限制函数，满 足这两个要求的斜率限制线性重构有： 1 m i n m o d ( o rm i n b c e ) ： 其中 2 s u p e r b e e ： 。r = m i n m 硼c 竽， 吼= m a x i 孙) d ( 蠢，蠢2 ) = 曲删( 攀捌訾，) 蠢2 ) = 商n r o o d ( 2 c 警，訾) 3 m c ( m o n o t o n i z c dc e n t r a l d i f f e r e n c el i m i t e r ) ： 以= 面姗嘲( 訾叛訾犯c 警，) 1 1 第二章有限体积方法以及数值实现 更一般的斜率限制器可以由流通量限制器得到：【1 5 ，2 】 r rr 厂 u i = ( 型羔半) ( 一i l ， 注意到我们现在是在数值求解方程组，所以需要对未知函数的两个分量都做出斜 率限制。对应于流通量限制做法这里同样有两种途径：一种叫做逐分量限制，这 种做法得到的数值解会出现某种“w r i g g l e s ”【7 1 ；另一种则是在特征方向上进行投 影【7 ，2 2 ，1 6 ，5 】然后限制，最后再返回到原来的数值解，本文将采取这种做法。 2 4 全离散数值格式 下而给出本文最终用到的数值格式，即给出时间的离散方式。基于不同的处 理机制，分别陈述如下。 ( 一) 基于斜率限制器的有限体积格式 基于斜率限制的有限体积数值格式如下： w = 叼一石a t ( 研t 2 ( 叨，v a ) 一f t - 1 2 ( 睨，) ) w = w 一塞( 只+ l 2 ( 咙，) 一f - 1 2 ( 咙，) ) 叼+ 1 = 去( 叼+ 嵋) 其中 观：阢+ 等以一l ，u r ：阢一等吼 咙：w + 百a x 吒_ 1 ：叼一百a x 吼, 以，西分别是由 阢 墨1 ，【w 】鉴1 是按照上述斜率修正限制而重构得到的斜率而数 值流通量f 则取以上提到的5 种流通量中的任意一种。以上全离散格式事实上采 用的是显式2 阶s s pr u n g e k u t t a 时间离散方法 2 0 1 。这是为达到整体误差的高阶精 度，我们同样需要对时间变量做出相应于空间精度的时间离散。这是因为对于一 阶双曲守恒问题，为满足c f l 稳定性条件，我们需要使得与z 满足：c a t a z c f l ，其中c 为j a c o b i a n 的绝对值最大特征值。 ( 二) 基于流通量限制的有限体积格式 首先我们给f l j l w ( l a x w e n d r o f o 格式所对应的数值流通量： e 2 = i ( f ( u i 1 ) + f ( 阢) ) 一互1 石a tl a 一1 2 1 2 ( 阢一阢一1 ) 第二章有限体积方法以及数值实现 其中a l 2 为上面提到的r o e 线性化矩阵，i a i 一1 2 1 定义如下： i a x 2 i = r i a i l 其中l a l = d i a g ( i a l l 2 i ，l a i l 2 1 ) 最后给出本文最终的数值流通量 只一1 2 = 硭l 2 + ( p ) ( e l w l 2 一硭l 2 ) 其中硭l 2 取上面提到的5 种数值流通量当中的任何一种，而可以取上述4 个流通 量限制器中任意一个。从而最终的全离散格式为 叼+ l = 叼一a z t ( f t + 1 1 2 - - e l 2 ) 采取这种做法，时问离散采用e u l e r 向前即可满足时空方向的整体2 阶精度。事实 上在解光滑区域，这种做法退化为著名的l w 格式。 1 3 第三章数值模拟结果比较 第三章数值模拟结果比较 我们主要针对四组不同初值( 解的结构不同) 进行刁同数值流通量与限制 器之间的组合搭配，进而观察其对于问断解的数值效果即分辨率的高低。为此， 我们将考虑不同初值情况进行相关的数值比较。所有算例均取空间区域为i = ( - 5 0 ，5 o ) ，重力加速度g 取为1 0 。相应的边界条件) 白d i f i c h l e t 边界条件，即 h ( x = - 5 0 ，t 0 ) 兰h ( x = - 5 0 ，t = 0 ) h u ( x = 5 0 ， 0 ) 兰h u ( x = 5 0 ，t = 0 ) 本文将考虑如卜网种初值，其定义如下： ( 1 ) 光滑初值 ( 2 ) 溃坝初值 fh ( x ，江0 ) = 1 0 + 0 4 e - 8 - 眙2 i - h u ( x ，= 0 ) = 0 0 k 芝：) ( 3 ) - - 相向水波初值 0 ， i f 0 i f z 0 i fz 0 ( 4 ) 二相离水波初值 三篆兰：) ： i ：s i ；f fx z o 。 相应的问题简称为模型( i ) 一模型( i 。 1 4 & l d o o l l d r l “卜 ) 0 划 江 0 扛畎 危 九 ，l_jll 第三章数值模拟结果比较 3 1 数值格式精度验证 我们利用问题i 进行精度验证，此时的问题真解是光滑的。对于精确解我们 采用3 阶的r k d g 方法数值解作为近似真解考察我们的格式的精度。网格数我们 选取为n = 2 0 0 ，4 0 0 ，8 0 0 ，1 6 0 0 计算终止时间为t = o 5 我们选取l 1 ，三2 ，l 范数去比 较 矿= ( 莩z 肼一) 。肛m 由表3 2 可知，数值解在厶l ，l 2 ，l 范数度量下均达n - - 阶精度，南表格3 4 我们同 样可以看到3 种误差度量下均有2 阶精度。但注意表格3 4 中我们只比较了3 种误差 度量下l w 格式的误差精度阶，这是因为在光滑区域所有流通量均变为l w 数值流 通量。在表3 1 和表3 3 我们还附带给出了程序的运行时间( 以秒为单位) 。数值试 验程序是在一台c p u 主频1 8 6 g ，内存l g 的机器上进行计算的，使用的足c o m p a q v i s u a lf o t r a n6 5 编译器。 i l f l l f r o e l l f h l l e f o r c e ic p u 时间 0 1 4 0 6 2 57 8 1 2 5 e 0 0 27 8 1 2 5 e 撇1 1 0 9 3 7 5o 1 2 5 表3 in = 1 6 0 0 时的程序运行时问：斜率限制器修i 3 2 基于斜率限制器处理的数值效果比较 本文更关心非光滑解时的数值表现，下面将对后面三个模型问题进行数值计 算模拟。此处讨论采用斜率限制器的处理方法，下一节讨论流通限制的处理方法。 为此本文采取先固定一种数值流通量，进而比较各种不同限制器的数值表现 差异；然后我们固定一种限制器，进而比较不同数值流通量的数值表现。最后，我 们发现在上述限制器中，s u p e r b e e 限制器得到的数值解耗散最弱，而m i n m o d 限制 器的数值耗散最强。m c 限制器则表现适中，也可视为一种整体效果不错的限制 器。对于算例n 和算例，当固定限制器后，在所有5 种数值流通量中，f o r c e 流 通量则表现地更好一些，而对于算例i ，在所有5 种数值流通量中，h l l e 流通量 则表现地更好一些。 1 5 第三章数值模拟结果比较 nf l u xl l e r r o ro r d e rl oe l t o ro r d e rl 2 e r r o ro r d e r 2 0 0 l f 0 0 0 4 40 0 0 5 2 o 0 0 2 8 2 0 0l l f0 0 0 4 40 0 0 5 2 0 0 0 2 8 2 0 0r o e0 0 0 4 40 0 0 5 20 0 0 2 8 2 0 0h l l e0 0 0 4 40 0 0 5 3o 0 0 2 8 2 0 0f o r c e0 0 0 4 40 0 0 5 3o 0 0 2 8 4 0 0l f0 0 0 l l1 9 6 0 9o 0 0 1 4i 8 7 9 47 3 2 1 5 e 0 0 4 1 9 2 8 5 4 0 0l l f0 0 0 l l1 9 6 2 60 0 0 1 41 8 8 7 17 3 3 5 7 e 0 0 4 1 9 3 0 6 4 0 0i 的e0 0 0 l l1 9 6 2 60 0 0 1 41 8 8 7 l7 3 3 5 7 e 0 0 41 9 3 0 6 4 0 0陇j eo 0 0 l l1 9 5 9 70 0 0 1 41 8 9 3 67 3 3 5 5 e 0 0 41 9 3 1 2 4 0 0f o r c eo 0 0 l l 1 9 6 0 0o 0 0 1 41 8 9 8 77 3 6 2 5 e - 0 0 41 9 3 5 9 8 0 0l f2 8 5 8 2 c 0 0 41 9 7 3 03 6 i ) 9 9 e - 0 0 41 9 5 8 71 8 6 2 8 e - 0 0 41 9 7 4 6 8 0 0l l f2 8 6 0 4 e - 0 0 41 9 7 3 l3 6 l8 7 e - 0 0 41 9 6 2 41 8 6 4 9 e - 0 0 41 9 7 5 9 8 0 0r o e2 8 6 0 4 e - 0 0 4l 。9 7 3 l3 6l8 7 e - 0 0 41 9 6 2 41 8 瞅4 ) 0 41 9 7 5 9 8 0 0h l l e2 8 5 8 5 e 0 0 4 1 9 7 2 2 3 6 2 2 6 c 0 0 41 9 6 4 61 8 6 4 7 e - 0 0 41 9 7 5 9 8 0 0f o r c e2 8 6l3 e - 0 0 41 9 7 1 93 6 3 3le - 0 0 41 9 6 8 31 8 6 8 4 e 0 ( m1 9 7 8 4 1 6 0 0l f7 4 7 2 3 e 0 0 51 9 3 5 59 1 6 8 8 e 0 0 51 9 7 7 24 7 9 1 2 e 伽1 9 5 9 1 1 6 0 0 l l f 7 4 8 5 0 e 0 0 5 1 9 3 4 l9 1 8 0 6 c 加1 0 51 9 7 8 84 7 9 4 c 0 0 51 9 5 8 4 1 6 0 0r o e7 4 8 5 0 e 0 0 51 9 3 4 l9 1 8 0 6 e 0 0 51 9 7 黯 4 7 9 8 4 e 0 0 5 1 9 5 4 8 1 6 0 0h u 。e7 4 8 2 9 e 0 0 51 9 3 3 69 1 8 5 2 e 舢1 9 7 9 6 4 7 9 8 2 e 0 0 51 9 5 8 4 1 6 0 0f o r c e7 5 0 3 6 e 0 0 51 9 3 l o9 1 9 7 2 e 0 0 51 9 8 1 94 8 0 9 2 e 删1 9 5 7 9 表3 2 斜率限制器仃限体秘方法的精度验证：时f h j 为t = 0 5 ，n 是等分单元个数 i l fl i j pr o e l l fh l l ef o r c e ic p uw , t n 0 2 3 4 3 7 5o 2 1 8 7 5o 2 1 8 7 50 2 5o 2 5 表3 3n = 1 6 0 0 时的程序运行时间：流通最限制 3 2 1 算例 下面给出数值结果的一些图形，计算终止时刻均取为t = 1 0 。 1 6 第三章数值模拟结果比较 nl 1 ，l l e r r o ro r d e rl 2 e r r o rl 2 0 r d e rl o ce r r o ro r d e r 2 0 0l 1 1 7 l ，o 0 0 2 9o 0 0 1 90 0 0 3 0 4 0 du 7 6 2 0 3 e - 0 0 41 9 l c r 74 9 5l3 e - 0 0 41 9 4 5 97 8 蛾- 0 0 4 1 9 2 3 2 8 0 0u 1 9 8 7 l e 删1 9 3 9 21 2 5 9 2 e - 0 0 41 9 7 5 32 0 1 3 3 e - 0 0 4 1 9 6 1 5 1 6 0 0l i 5 51 0 0 e 0 0 51 8 5 0 63 3 3 8 8 e 加1 0 51 9 1 5 l5 4 6 l l e - 0 0 51 8 8 2 3 表3 4 流通量限制的数值格式的精度验证：此时的数值流通量均为l w 图3 1固定流通量为l f ，四利限制策略的数值比较 图3 1 是固定流通量l f ，分别采取四种限制策略所得到的数值结果比较，其中 网格数是5 0 。其中左边第一个图是比较精确解( 高阶数值解) 、m i n m o d 限制策略数 值解和v a nl e e r 限制策略数值解的比较，图中可以看出在稀疏波和激波部分，v a n l e e r 限制策略数值解明显比m i n m o d 限制策略数值解更加接近于精确解，也就足说 在算例i i 中，当固定流通量l f 时，采用v a nl e c f 限制策略数值效果更优。然后我们 在看右边的两个图，也可以看出m c 限制策略优于v a nl e e r 限制策略，s u p e r b e e 限 制策略优于m c 限制策略。最后得出结论是，对于算例l i ，采用斜率限制方法、固 定流通量是l f 时，s u p e r b e e 限制策略最好。 图3 2 固定流通量为l l f ，四种限制策略的数值比较 图3 2 是固定流通量u j ，分别采取四种限制策略所得到的数值结果比较，其 中网格数是5 0 。其中左边第一个图是比较精确解( 高阶数值解) 、m i n m o d 限制策略 数值解和v a nl 胱r 限制策略数值解的比较，图中可以看出在稀疏波和激波部分，v a n l e e r 限制策略数值解明显比m i n m o