（计算数学专业论文）求解双曲守恒律方程的群速度修正法.pdf

求解双曲守恒律方程的群速度修正法 专业：计算数学 硕士生：周超红 指导教师：朱庆勇教授 中文摘要 本文研究了利用群速度控制法的思想构造双曲守恒律方程高分辨 率格式的一般方法。在前人研究的基础上，本文研究了群速度控制项 的性质，进而指出，普通差分算子均可以利用群速度修正算子进行修 正而得到相应的群速度控制格式利用群速度修正的思想，本文得到 了三阶迎风群速度控制格式、三阶迎风紧致群速度控制格式等多种高 分辨率格式。数值实验结果表明，这些格式均可以很好地控制非物理 振荡并能保证精度。通过引入色散度量以及各向异性度量的概念，本 文还为评估差分算子对微分算子的逼近程度提供了客观标准。 关键词：高分辨率格式，群速度控制法，双曲守恒律方程，群速度修 正，各向异性度量 g r o u pv e l o c i t ym o d i f i e ds c h e m e f o rs o l v i n gh y p e r b o f i c c o n s e r v a t i o nl a w s m a jo r ：c o m p u t a t i o n a lm a t h e m a t i c s n a m e ：z h o uc h a o h o n g s u p e r v i s o r ：p r o f e s s o r z h uq i n g y o n g a b s t r a c t b a s e do nt h eg r o u pv e l o c i t yc o n t r o lm e t h o d ，t h et h e s i si sd e v o t e dt os t u d yt h eg e n e r a lm e t h o d o f c o n s t r u c t i n gh i g h - r e s o l u t i o ns c h e m e st os o l v eh y p e r b o l i cc o n s e r v a t i o nl a w s o nt h eb a s i so f o t h e r s s t u d i e s ，t h et h e s i sr e s e a r c h e st h ep r o p e r t yo fg r o u pc o n t r o li t e m ，a n dt h e np o i n t so u tt h a t c o m m o nd i f f e r e n c eo p e r a t o rc a nb em o d i f i e db yt h eg r o u p v e l o c i t ym o d i f i e di t e mt oc o n s t r u c tt h e c o r r e s p o n d i n gg r o u pc o n t r o ls c h e m e b ym a k i n gu s eo ft h et h e o r yo f g r o u pv e l o c i t ym o d i f i c a t i o n ， t h et h e s i so b t a i n ss o m eh i g h - r e s o l u t i o ns c h e m e ss u c ha st h r e e - o r d e ru p w i n dg r o u pc o n t r o ls c h e m e , t h r e e - o r d e ru p w i n dc o m p a c tg r o u pc o n t r o ls c h e m e ，a n ds oo n t h er e s u l to f n u m e r i c a le x p e r i m e n t d e m o n s t r a t e st h a tt h e s es c h e m e sc a nc o n t r o ln o n p h y s i c a lo s c i l l a t i o na n de l l b x l r ea c c u r a c y b y m e a n so fi n t r o d u c i n gc o n c e p t i o n so fd i s p e r s i o nm e a s u r ea n da n i s o t r o p ym e a s u r e , t h et h e s i s p r o p o s e sa no b j e c t i v ec r i t e r i o nt oe v a l u a t ea p p r o x i m a t i o nd e g r e eo fd i f f e r e n c eo p e r a t o rt o d i f f e r e n t i a lo p e r a t o r k e yw o r d s ：h i g h - r e s o l u t i o ns c h e m e ，g r o u pv e l o c i t yc o n t r o lm e t h o d ，h y p e r b o l i c c o n s e r v a t i o nl a w s ，g r o u pv e l o c i t ym o d i f i e d ，a n i s o t r o p ym e a s u r e 原创性声明 本人郑重声明：所呈交的学位论文，是本人在导师的指导下，独立进行研究 工作所取得的成果。除文中已经注明引用的内容外，本论文不包含任何其他个人 或集体已经发表或撰写过的作品成果。对本文的研究作出重要贡献的个人和集 体，均已在文中以明确方式标明。本人完全意识到本声明的法律结果由本人承担。 日期：洲。年r 月1 日 学位论文使用授权声明 1， 一躲触渊嗽：骞易、锣 日期：列睥f 月工7 日 日期：工叶年f 月上7 日 1 1 选题背景和意义 第1 章引言 许多物理过程的数学模型都可以写成偏微分方程的形式，针对这些偏微分 方程的不同性质，人们把偏微分方程分为椭圆型、抛物型、双曲型。守恒律方 程则是指出物理过程满足的守恒性质直接得到的方程。流体力学中的方程大多 为双曲守恒律方程。 求解双曲守恒律方程的数值方法主要有有限差分法、有限元法、有限体积 法三大类。有限差分法以经典的数学逼近理论为基础，直接将微分问题变为代 数问题，数学概念直观，表达简单，方法灵活简便，易于实现且具有高度的通 用性，是发展最为成熟和应用最为广泛的数值方法。有限元方法则是借助于变 分原理，将微分方程离散求解。有限元方法最早应用于结构力学，后来随着计 算机的发展慢慢用于流体力学的数值模拟。有限体积法可视作有限单元法和有 限差分法的中间物。但有限体积法中因变量的积分守恒对任意一组控制体积都 得到满足，与有限差分法相比，有限体积法对守恒型双曲方程的求解有着很大 的优越性，越来越受到人们的关注。 早期针对双曲型方程的研究主要偏重于理论方面，h a d a m a r d 、f r i e d r i c h s 和c o u r a n t 等人研究了偏微分方程的基本特征、数学提法的适定性、物理波的传 播特性和解的唯一性问题。c o u r a n t 、l e w y 和f r i e d r i c h s 1 证明了连续的椭圆型、 抛物型和双曲型方程组解的存在性和唯一性：他们还进一步研究了双曲型方程 的特征性时，进而提出了特征线方法，给出了著名的稳定性判别条件- c f l 条 件。这些工作结合其他一些数学家研究的偏微分方程的数学理论，构成了有限 差分方法的数学理论基础。随后v o nn e u m a n n 、r i c h t m y e r , h o p f 和l a x 等人建立 了非线性双曲守恒律方程的数值方法理论，为含有涌波和其他间断的流动数 值模拟打下了理论基础。v o nn e u m a n n 还提出了一种时间推进问题数值方法线 性稳定性的分析方法，该方法至今仍被广泛使用【2 】。之后，g o d o n o v 提出了迎 风格式，l a x $ 1 ：i w e n d r o 躜出了二阶精度的差分格式，改进了计算激波流动的精 度。 对于非线性双曲型方程，不管初始值如何光滑，其解都可能出现间断。早 1 2 国内外研究现状 期提出的各种求解非线性双曲型方程的差分格式，对于没有大梯度的定常光滑 流动的求解，计算结果尚令人满意，但当出现强间断时，数值解就会产生虚假 振荡。而各种低阶格式虽然可以有效消除数值振荡，但往往耗散过大，造成格 式精度下降及间断抹平的现象，与实际物理现象不符。间断的处理问题是非线 性双曲型方程数值模拟的关键所在，也是难点所在。间断问题的处理，需要一 种既能有效消除非物理振荡，又能具有比较高的精度的高分辨率方法。 1 2 国内外研究现状 自二十世纪八十年代h a r t e n 第一次提出高分辨率方法和总变差递减( t o t a l v a r i a t i o nd i m i n i s h i n g ，t v d ) 的概念以来，在计算流体力学领域，掀起了研究有 限差分方法的新高潮。t v d 格式能够保持数值解的单调性，因而可以有效地抑 制间断附近产生的振荡。然而t v d 格式的精度至多只能达到二阶，而且多维问 题不存在高于一阶的t v d 格式。为了获得更高精度的格式，h a r t e n 、e n g q u i s t 3 】 等又提出本质无振荡( e s s e n t i a l l yn o n o s c i l l a t o r y ，e n o ) 格式。与t v d 格式相比， e n o 格式并不要求总变差严格递减，允许总变差有微小的增加。为达到高阶精 度并消除虚假振荡的目的，e n o 格式在流通矢量重构的过程中，通过比较差商 绝对值的大小来自适应地选择插值模板。由于传统e n o 格式求解多维问题时计 算量较大，之后s h u 和o s h e r 4 又提出了基于点值重构和t v d 龙格库塔时间离 散的通量形式e n o 格式。为解决e n o 格式模板选择以及算法并行计算的问题， l i u 、o s h e r 等人提出了加权e n o ( w e n o ) 格式 5 1 。w e n o 的基本思想是将e n o 格 式只选择最光滑模板改进为选择每一个可能模板的加权平均，其中权值为模板 的光滑程度的度量。目前，e n o 和w e n o 格式已被广泛应用于双曲守恒律方程 的求解中，并取得了良好的效果 6 】。 傅德薰和马延文等人【7 】通过对差分格式非物理振荡产生原因的分析，提出 了群速度控制方法，该方法通过在激波前后使用不同的算子控制振荡波的群速 度，从而抑制激波附近的非物理振荡。文【8 】采用具有三阶精度的迎风紧致格式 求解驱动方腔流动问题。文【9 】提出了一种用于空气动力学方程组高精度的群速 度数值方法。文【1 0 】采用了群速度控制方法求解二维r i e m a n n 问题。文【1 1 】从六 阶精度的紧致格式和五阶精度的迎风紧致格式出发，引入一个函数用以控制数值 解的群速度，构造出了具有六阶精度的紧致型差分格式。文 1 2 1 提出了八阶群速 度控制方法，用于对可压缩衰减湍流的直接数值模拟。 此外，国内外学者还提出了c s c m 方法【1 3 】、n n d 格式【1 4 】、m m b 格式 1 5 1 、 精致差分格式以及耗散比拟方法【1 6 】等方法。 2 第1 章引言 1 3 本论文研究的主要内容 本文将研究使用群速度控制的思想构造双瞳守恒律方程高分辨率格式的一 般方法，并给出几种高分辨率格式的具体参数。 论文的具体章节安排如下； 第一章介绍了论文研究的背景和意义，双曲守恒律方程的国内外研究现 状，同时简要介绍了本文的研究内容。 第二章我们介绍本文研究的一些基本预备知识，包括曲线坐标的网格变 换的方法、群速度控制法的一些基础知识以及一阶导数的高阶逼近格式具体表 达形式和基本性质。 第三章为本论文的重点，第一节研究了群速度控制项的性质，第二节我 们提出了算子修正的思想，给出了构造高精度群速度控制格式的一般方法。第 三节，我们首先引入了色散度量、各项异性度量的定义，进而提出了确定群速 度修正参数的方法，然后利用给出的方法求解出了几种常用算子对应的最佳的 群速度修正系数的值。 第四章介绍了无粘性、无摩擦、河底无倾斜的理想情况下二维浅水方程 的形式，曲线坐标网格变换下二维浅水方程的形式及其流通矢量的一种分裂格 式。最后，给出了利用本文提出的方法求解二维浅水方程的几个数值算例。 第五章总结分析了本文提出的群速度修正法的特点及意义。 3 第2 章预备知识 2 1 曲线坐标网格的变换 在流体力学的数值模拟中，由于流体的物理区域的几何形状比较复杂，物 理量在区域内的变化较大，边界条件复杂，这样，很多在规则区域能够较好使 用的计算方法无法使用。为了解决这样的问题，我们希能通过坐标变换，将几 何形状复杂的区域变成规则区域，将复杂的物理区域上的定解问题转化为规 则的变换区域的定解问题，然后再选择有限差分法，或有限元法等数值算法求 解变换区域的问题。美国密西西比大学的j f t h o m p s o n ，f c t h a m e s 和c w m a s t i n 提出了边界拟合法。边界拟合法给出了从复杂不规则的物理区域与其对应的规 则且易于计算的区域的变换。利用复变函数的思想，t o m p s o n 他们得出，把物理 平面的求解区域变换成变换平面的矩形区域的交换f ( z ，秒) ，叩( z ，可) ，是下列拉普 拉斯方程组边界问题的解【1 7 】： fv 2 = 0 1 2 栌o ： ( 2 1 ) 蝴 = = 手 u j 【，7 i 边界t ：= 亓 实际上，我们希望在变换平面上得到一个均匀的网格，也就是说需要在变 换平面网格已知的情况下求解物理平面上与之相对应的网格节点。将式( 2 1 ) 自 变量与因变量调换，就得到了如下方程： lq z 毖一2 p z e 卵+ ，y z 7 聊= 0 q 珧f 一2 p 欢叩+ y y , m = 0 ( 2 2 ) l z ，y i a 界 = = 佩雪) 其中： a = x ；+ 砺2 ， p = 歌十耿蜥，7 = x ；+ 鳕 式( 2 2 ) 是一个非线性方程组，利用如下差分逼近： z e2 三生皆，珧2 丝学，z 叼2 兰尘铲，5 丝皆， 5 2 2 群速度控制法简介 z 翱2 耿叩2 z 5 珧2 ( x i + l ，j + l x i 十1 ，j 一1 )一( x i 一1 ，j + 1 一x i 一1 ，j 一1 ) 2 叩 ( y i + l ，j + l y i + l ，卜1 )一( y i 一1 , j + l y i 一1 , j 一1 ) x i + l ，j 一2 x i ，j + 既一1 ，j 2 ( 越) 2 y i + l ，j 一2 y i ，j + y i 一1 ，j 2 ( 世) 2 2 a a r l ，z 聊2 x i ，j + l 一2 x i ，j + 舰，j 一1 2 ( 叩) 2 y i ，j + l 一2 y i ，j + 犰，j 一1 2 ( z x , 7 ) 2 将式( 2 2 ) 化为非线性差分方程组，可利用超松弛迭代法求解。图( 2 1 ) 显示了倾 斜水跃问题的曲线网格变换。 x 图2 1 倾斜水跃问题的曲线网格变换 2 2 群速度控制法简介 考虑单波方程及其半离散方程 ( 2 3 ) u ，j = ， 仉 = 6 ，一za 一拶 + 丝珧 第2 章预备知识 尝+ 曼：0 (24)a 挑z 、 7 这罩乃z 是一阶导数o f a x 的差分逼近，定义差分算子如下： 砖办= 土( 办士一办) 鲤= ( 砝+ ) 1 2 = ( 乃+ 1 一办一1 ) 2 ( 2 5 ) 磋= 砝e = ( 办+ 1 2 办+ 办一1 ) 2 弓可以表示为乃= 风砝力+ m = 风( 乃+ m + 一办+ m ) 。 mm 下面我们利用傅里叶分析法对差分方法进行理论分析。 考虑y ( u ) = n u ，其c a y 9 常数，初值条件u ( x ，0 ) = ( z ) ，考虑以2 丌为周期 o 。 的解，初始条件可进行傅里叶级数展开，u o ( z ) = e x p ( i k x ) ，这里忽表示 七= 一0 0 波数。不难得到方程( 2 3 ) 的精确解 0 0 u ( x ，亡) = c ke x p i k ( x 一酬】_ 呦( z a t ) ( 2 6 ) 七= 一o o 由此可知，具有不同波数的波，全都以相同的速度a 传播。 对于半离散方程( 2 4 ) ，我们设z 方向网格步长为：a x = 努，= j a x ，j = 0 ，1 ，n 设其初始条件为 n2,v2 = e x p ( i k x j ) ，马= n c k k e e x p ( i k x j ) ，k e = + i 觑 ( 2 7 ) k = - n | 2k = - n | 2 其中： b = 6 , c o s k ( m + 1 ) a m 一c o s k m a x = f i r s i n 陋( m + 1 ) a x 一s i n k m a x 则半离散方程( 2 4 ) 的解是 ( 2 8 ) 不失一般性，我们只考虑单波解= e x p ( i k x j ) ，那么此时半离散方程的解为： = e x p ( 一口惫加x p i k ( x ：i n 函k i 圳 ( 2 9 ) 7 d 旦心 一 z p x 幻 旦如 一x 舭叫 = 嘶 2 3 一阶导数的高阶逼近格式 对比( 2 6 ) 和( 2 8 ) ，可知对于任意的k ，k a x 0 说明格式是正耗散的，n 0 说明格式是负耗散的，而a k ，= o n 说 明格式是非耗散的。负耗散的格式是不稳定的，不能用于数值求解。数值色散 造成不同波的传播速度不一致，从而产生非物理振荡。而数值耗散对数值解有 光滑效果，可以抹平色散造成的振荡，但可能会降低数值解的精度。高分辨率 格式要求尽量消除数值色散引起的非物理振荡，同时也要求尽量减弱数值耗散 对数值解的抹平作用。 波的群速度，是指波振幅外形上的变化在空间中传播的速度。根据以上的 讨论，我们可定义群速度为s 口= d k t d n ，其中q = k a x 。对于( 2 3 ) 的精确解( 2 6 ) ， k i = o l ，s 口= 1 。对( 2 8 ) 的分析我们不难发现，由于差分离散的误差导致对不同 的波数k 对应的波的传播速度不同，即群速度不一致，从而引起数值解产生非物 理振荡。 定义2 1 【7 】如果8 9 1 ，0 o l 1 , 0 q 7 r ，则格式称为快格式，相应的差 分算子称为快算子；若存在使得当0 k a x k o a x 1 ，而 当k o a x k a x 7 r 时8 9 1 ，那么格式称为混合格式，相应的差分算子称为混 合型算子。 傅德熏、马研文在分析了数值振荡产生的原因之后，提出了群速度控制的 思想，简单说就是在激波前使用慢算子，在激波后使用快算子，使振荡向激波 靠拢，从而达到消除振荡的目的。实际上，单纯的高阶快算子并不存在，由于 耗敖的存在，混合算子的振荡基本都出现的激波前，因而一般可使用混合算子 替代快算子。 2 3 一阶导数的高阶逼近格式 为了构造双曲型方程的高精度求解算法，需要了解一些逼近一阶导数项的 常用算子的群速度性质。 以下f j a z 是- - 阶导数o y o x 的差分逼近。 ( 1 ) 四阶中心差分逼近格式：乃= 丧【f j + 2 + 8 ( 办+ 1 一办一1 ) + 乃一2 】 由此得到 = 0 = ( 一s i n 2 a + 8 s i n a ) 6 8 9 = 一4 ( c o sa 一1 ) 2 + 6 6 l 当o a 1 ，当詈 0 觑= ( 8 s i n a s i n 2 a ) 6 8 9 = 一吾2 ( c o s o ! 一1 ) 2 + 1 0 ) ： 地乃= 喜乃+ 言乃一， 元j c f = ( 丢篮+ 吾砝) 乃 鳄= ) - 1 瓦 9 2 3 一阶导数的高阶逼近格式 由差分格式可知 图2 2 一阶导数差分算子的色散性质 妊鼎 。 k i = 岛2 8 s i n q + s i n q c o s q 4c o sq + 5 ( 1 0c o s ( n ) 2 + 4 0c o s ( i x ) + 4c o s ( q ) 3 + 2 7 ) ( 1 6c o s ( q ) 2 + 4 0c o s ( a ) + 2 5 ) l , a e ( 0 警) 可见，格式是耗散的，对低波数波，算子为快算子，实际应用中可视为快算 子。 1 0 第3 章群速度修正格式 傅德薰、马延文提出的群速度控制法，基本思想是在激波前使用慢算子， 激波后使用快算子。然而，根本不存在真正的快算子，因而直接构造出的群速 度控制格式往往仍然存在非物理振荡。马研文【l8 】、朱庆勇【1 9 ，2 0 等人在构造 具体的群速度控制格式时均是在右端引入修j 下项，以使得修正后的差分算子在 激波前表现为慢算子，激波后表现为快算子。 群速度控制法的思想，其实也可以看做是引入了群速度修j 下算子将原有的 差分进行修正。群速度修j 下算子在激波前后的不同性质，使得原有的差分格式 所引起的激波前后的振荡向激波靠拢，从而达到消除激波振荡的目的。 本章将着重研究群速度控制法中经常引入的群速度修正项的性质，进而提 出一种构造群速度控制格式的一般方法。同时本章还将针对利用普通差分算子 构造群速度控制格式的一维、二维问题，提供寻找最合适的群速度修正参数的 方法。为了行文方便，我们把离散格式中未经过群速度控制时使用的用于逼近 一阶导数的算子称为基本算子。 3 1 群速度控制项的性质 为了说明引入新的格式，我们还要回到方程( 2 3 ) 及其半离散格式( 2 4 ) ： 窑+ 筹_ 0 ，训u ) ( 2 3 ) 丝+ 墨：o( 2 4 ) o t a z 。 、- 。7 在对群速度控制法的研究中我们发现，朱庆勇在文【1 9 】中构造三阶迎风紧致群 速度控制格式 等+ 乏：石o 硝0s 鲥一圳。i 鹾u j 一丝坠卺掣型( 3 1 ) 其实就是针对三阶迎风紧致格式引入了如下的群速度控制项 酽t y0 m 训h 黝一型学 1 1 3 1 群速度控制项的性质 对半离散方程( 2 4 ) 进行修j 下。参考上式，我们考虑群速度修正项 l g = 堡掣一鹾 s 勺咖h 鹾】 其中s 彤= s i a n ( u j + 1 一) ( + 2 一u j + i u j + + 2 ) 】，称为激波结构函数。咖满 足在振荡区域奶= 1 ，平滑区域咖= 0 ，称为开关函数1 。我们利用第二章介绍 的知识来分析一下岛的性质。 考虑，( 让) = a u ，其中n 为常数，并先假设n 0 。通过分析s 勺的不同取值有 l f j + l + 3 乃一3 乃一1 + 办一2 ， s s j + l = s s j 一1 = 1 ； r j 一办+ 1 + 2 乃一力一1 ，8 8 j + 1 = 1 ，s s j 一1 = 一1 ； l g 一1 办+ 2 3 办+ 1 + 3 乃一3 乃一1 + 力一2 ，s s j + l = 一1 ，8 s j 一1 = 1 ； 【f j + 2 3 力+ 1 + 3 办一厶一1 ，8 s j + 1 = s s j 一1 = - 1 ( 1 ) 在激波后光滑区域，s s j + 1 = s s j 一1 = 1 。 b ( l 9 ) = 2 ( c 0 8 1 ) 2 0 ( l g ) = 2s i n c e s i n2 a s 9 ( l 夕) = - 4 c o s 2n + 2c o s c e + 2 0 ，a ( 0 ，等) 从而可知，岛是币耗散的，且会加快基本算子。换句话说，岛一方面可以 耗散掉振荡，另一方面又可以使得振荡向激波靠拢。 ( 2 ) 在极值点，8 s j + 1 = 1 ，s s j 一1 = - 1 。 辟( l 9 ) = 2 2c o s a 0 ( l 9 ) = 0 s g ( l g ) = 0 可得出岛为j 下耗散算子，极值点一般为振荡波位置，岛的耗散性质正有利于振 荡波的消除。 ( 3 ) 在拐点，s 勺+ 1 = 一1 ，s s j 一1 = 1 。 ( l 9 ) = 4 c o s 2q 一6c o s + 1 k i ( l 9 ) = 0 s 9 ( l g ) = 0 f1 ， 1 可取咖= l 0 ， ( 1 一u ) ( i u j + 1 一 t j + 1 一j 。- l u 2 u 3 + 1 l u j + 。1 l + l + 2 “j + ”j l i ( 1 一u ) ( i t j + 1 一u j+ l u j t 一1j ) + u i u j + 1 + 2 u j + u f 一1 l 1 2 ， 0 w 1 k = 0 0 5 一 第3 章群速度修正格式 c o s 0 ( 等乎，1 ) 时，b ( 岛) 0 k i ( l 9 ) = s i n 2 ( 1 2 s i n a s g ( = 4 c o s 20 1 - - 2 c o s o ! - - 2 岫( 0 ，警) 这表明l 。是j 下耗散的，且会拖慢基本算子。也就是说，l 。即可以耗散掉振 荡，又可以使振荡波指向激波。 根据以上分析我们得知，在平滑区域，l g = 0 ；在激波前后，l 。能使得振 荡波向激波靠拢；在极值点，己。只有耗散而无色散，有利于振荡的消除；而在 拐点，l 口呈现负耗散性质，则有利于激波形状的保持。根据l 。的这些性质，我 们得知，己。其实可以对任意差分格式进行修正，而不仅仅限于三阶迎风紧致 1 3 3 2 群速度修正格式的构造 格式。同时，我们注意到，如在激波前后的性质具有明显的对称性，因而可推 知厶对一些在激波前后出现对称性振荡的算子应当有相对比较好的修正效果。 3 2 群速度修正格式的构造 我们将方程( 2 3 ) 的半离散格式重新写成 鲁+ 警= o 2 ， 在对方程( 2 3 ) ( 厂= u ) 进行数值实验中，我们发现群速度控制项岛表现出 的是很明显的耗散性质，岛可看作一个光滑算子。 ( i ) 四阶紧致差分算子( i i ) 四阶紧致著分+ 群速度修正 ( 1 1 1 ) 三阶迎风筹分算子( i 三阶迎风差分+ 群速度修正 注：蓝色线为t = 5 s u , t 结粜红色线为t = 5 0 s 结粜 图3 2 群速度修正项的耗散性 由图3 2 可知，无论基本算子性质如何，加入群速度修正之后都能有效消除 振荡。但根据( i i ) ( i v ) ，我们发现经过修j 下后的算子的表现出来的耗散性，使得 我们基本上无法区分所使用的基本算子。因此，直接使用群速度控制项固然能 消除振荡，但大大降低了基本算子原有的精度。 1 4 第3 章群速度修正格式 为了兼顾振荡控制和分辨率，我们采用在群速度修正项前加修正系数的方 法。即采用半离散格式： 等+ f j f + c r l 9 = 0 ，其中弓为基本算子，0 仃 1 ( 3 3 ) 我们改称式( 3 3 ) 中的l 。为群速度修正算子，称盯为群速度修正系数。式( 3 3 ) 禾l j 用 算子修正的观点来构造群速度控制格式，有利于对构造出的格式进行分析。 图3 3 显示了仃= 0 2 时，各基本算子修正后的计算结果。 ( i ) 阴阶紧致筹分算f( i i ) 三阶迎风紧致群速度控制 ( i i i ) 三阶迎风群速度控制( i v ) 四阶迎风群速度控制 注：蓝色线为t = 5 s 时结粜红色线为t = 5 0 s 时结果 图3 3 基本算子的群速度修正 从图3 3 可知，基本算子加修正系数盯进行群速度修j 下后，四阶紧致、四阶 中心差分在开始时都是有振荡的，但n t = 5 0 s 时，振荡都基本消失了，计算 精度也得到了很好的保持。我们注意到，作为更高阶逼近的( i ) ( i v ) 的逼近精度 与( i i ) ( i i i ) 相比，并不见得更好。因为四阶中心格式和四阶中心紧致格式的振荡 几乎都在激波后，群速度修萨项没能够抵消振荡。( i i i ) 的计算结果精度最好，这 跟我们之前的理论分析相符。 对于半离散格式( 3 3 ) ，使用不同的基本算子和群速度修正系数o r ，我们可以 构造出不同的群速度控制格式。我们将使用格式( 3 3 ) 构造群速度控制格式的方 1 5 3 3 群速度修正系数的确定 法，称为群速度修正法。 3 3 群速度修正系数的确定 为了确定不同基本算子的修币格式中群速度修正参数的大小，我们引入 色散度量函数，用来度量逼近一阶微分的差分算子的色散误差的大小。定义 色散度量d m f = pf 下( s 口一1 ) 2 d a ，其中s 。为差分算子的群速度，是q 的函数。 b 。 丁= e c o s ( a ) 一e - 1 称为波成分比例因子，表示数值解中与q 对应的波的比例。肛称 为色散显现因子，用于度量耗散对色散的压制力度。对于基本算子与其群速度 控制格式之间色散程度的比较，“可取为1 v e 。显然d i n 为盯的函数。 波成分比例冈子色散显现闪子 图3 4 波成分比例冈子与色散显现冈子 下面以三阶迎j x l 群速度控制格式为例，利用色散度量函数，寻找最佳的群 速度修f 参数。由于三阶迎风格式的振荡主要在激波后，我们考虑，激波后三 阶迎风群速度控制格式的差分算子 b = 丢( 2 办+ 1 + 3 办一6 办一。+ 办一2 ) + 口( 一办+ + 3 f j 一3 乃一+ 乃一2 ) 1 6 第3 章群速度修正格式 s 9 ( 岛) = ( 一百a ( c o sq 一1 ) 2 + 1 ) + 盯( 一4 c o s 2a + 2c o s + 2 ) u 给定不同的口，将8 口代入d m f 右端，对d m f 的右端求数值积分( 可取q n = n t r l o ，礼= 0 ，1 ，1 0 ) ，再画出盯- d m f 的曲线，d m f 最低点对应的仃，即为最 优的修j 下参数。由于群速度控制格式中的数值耗散随着仃的增大而增大，因而若 最低点偏左侧点的纵坐标与最低点的纵坐标相差不大，最优的修正参数则应取 偏左侧点对应的盯。 图3 5 群速度控制算子的色散度量 对于二维问题，我们还需要考虑算子的各项异性效应【9 】。考虑二维- 0 - 波方 程 尝+ 。不o u + 6 豢：o ， n ，6 为常量( 3 4 ) 瓦+ o 瓦+ 6 瓦2 o ，n ，6 为常量( 3 4 ) 初值u ( x ，y ，0 ) = e x p ( i k x ) 其中k = k l ，k 2 r , x = p ，y t , k 1 ，k 2 分别为z ，y 方向 上的波数。定义 z = 丽a ，丽br 1 7 3 3 群速度修正系数的确定 并将它重写为f = c o s ，s i n 0 t ，这罩的秽就是方位角。该问题的精确解可以表示 为 u ( x ，t ) = e x p ( i k x 一 0 2 + b 2 k 冽) ( 3 5 ) 方程( 3 4 ) 的半离散格式为 鲁+ n 警舶哿= 。 6 ， 其中毛乱巧a x 、毛u 巧夕分别为挑如、d 乱a 秽的差分逼近。；b n o 6 ) 的精确解 可以表示为 u ( x ，亡) = e x p ( 一、石雨k 1 # ) e x p ( i k x 一 孑了万k 2 卅) ( 3 7 ) 其中 扣【譬c o s 伊，譬咖旷 讧 譬c o s 口，譬咖旷 q = k l a x ，p = k 2 a y 这罩的姆) ，七：，舻) ，庇：2 分别为差分算予如，如傅立叶分析中的实部和虚部。比较 式( 3 5 ) 、式( 3 7 ) 我们得知，考察算子的各项异性效应只需考察 丽l l i = ：1 k , ) c o so + 2 ) s i n0 对川l 的逼近程度。 为度量差分算子的各项异性，仿照前面色散度量的定义，我们定义 埘= 丁( | l z 卜赫) 2 批 为各项异性度量。其中丁= e e o s ( 。) 一e 一1 称为波成分比例因子，表示数值解中u 与 对应的波的比例。z ，= 1 盯4 3 称为各向异性显现因子，用于度量耗散对各向异性 的压制力度。 由于方向角的不确定性，对于群速度控制法，可假设在z 方向与! ，方向的群 速度修j 下参数相等，均为盯，则各项异性度量a e f 就是仃的函数。同样，对于不同 的盯，我们可以利用数值积分求解a e ，得到盯与a ，之间的对应关系，进而找 出最佳的群速度修正参数盯。图3 6 显示了三阶迎风群速度的不同盯值对应的榆与 第3 章 群速度修正格式 圈3 6 三阶迎风群速度控制格式的箨向异性效应 方位角之间的关系，可看出盯= 0 1 时，三阶迎风群速度控制算子的各项异性效 应最小。 根据本节提供的方法，我们可以得到各种高阶群速度控制法的最佳群速度 修正参数，具体为表( 3 1 ) 。 由于开关函数的存在，平直区域将不进行群速度控制，因而只需考虑主要 振荡区域对应的最佳群速度修正参数的值。对于一维问题只需考虑色散度量下 的值。例如，三阶迎风算子的振荡主要出现在激波后，因而推荐值取为激波后 色散度量下的最佳值0 2 0 。对于二维问题，要同时考虑色散度量和各向异性度 量。如四阶迎风算子的振荡主要出现在激波前色散度量最佳值为0 2 0 ，各向异 性度量下的最佳值为0 1 0 ，推荐值取两者的平均值0 1 5 。对于三阶迎风紧致算 子，实验发现盯取0 0 5 时色散比较明显，故推荐值取为0 1 0 。 从表( 3 1 ) 可得出，三阶迎风紧致群速度控制格式具有最好的性质。而四阶 中心紧致群速度控制格式，由于各向异性效应大，不太适合模拟二维问题。令 人意外的是，四阶中心群速度控制格式也将是不错的格式，这个格式将是大量 流通矢量不易分裂的方程构造差分格式的首选。 1 9 3 3 群速度修正系数的确定 面墨i 刨 趟 琳 星 谑 。 图3 7 群速度控制算子的各向异性度量 表3 1 基本算子最佳群速度修正参数对应表 第4 章二维浅水方程 浅水方程是计算流体力学中的一类重要方程，是一种特殊的非线性双曲守 恒律方程。浅水方程的应用范围十分广泛，可用于河道流量、溃坝决堤、污染 物输运扩散、河道输沙、洪水预报、河口潮汐、涌浪等大量问题。由于浅水方 程的应用背景十分广泛，其数值方法的研究自然成为了是计算流体力学和计算 水动力学的重要课题。由于二维浅水方程的一般形式比较复杂，本文中使用在 无粘性、无摩擦、河底无倾斜的理想情况下的二维浅水方程作为求解模型。 4 1 浅水方程模型 二维浅水方程在无粘性、无摩擦、河底无倾斜的理想情况下， 的齐次守恒形式如下： 掣+ 娑+ 榘：o 8 td z d n 其中， 其数学模型 丫= ( 兰：) ，f = ( 三兰+ ；夕舻) ，g = ( 三+ ；夕舻) ( 4 1 ) 为水深，牡为z 方向的速度大小，秽为秒方向的速度大小，夕为重力加速度。方 程( 4 1 ) 的雅可比矩阵为 这里c = “_ o 矩阵a 的特征矩阵、右特征向量矩阵、左特征向量矩阵分别是 、l，-、 1 u 孙 o 仃0 u 秒 o ，j-il_lilil、 i i g 一扩丝彬 = b 、llllj， 0 o 珏 1 轧 秽 u ” 0 一 wo扛叫 ，。 = f u竺彬 l | a 、liiij， c 0 o +u 0 缸0 c 扒 一 0 0 珏 ，i。一 i | ld 4 2 浅水方程的曲线坐标变换 耻彳1 卜文兰 。2 = ( 吾c 三移呈c ) 耻1 鼍1 卜文羔 4 2 浅水方程的曲线坐标变换 在对实际流场进行数值模拟时，往往会碰到区域边界形状复杂和物理变量 在区域内变化大的情况。这时，为了较好的进行模拟，需要先建立曲线坐标网 格，即把不规则物理域先映射到规则计算域，然后才能在计算域上进行差分计 算。当然，相应的流体力学方程组也需要变换到计算域。 假设物理域上的( t ，z ，y ) 到方形计算域上的( 丁，舌，7 7 ) 之间的全局同胚映射关系 为： 7 - = t ，荨= ( 亡，z ，可) ，叩= 叩( 亡，z ，y ) 则由方程( 4 1 ) 可导出二维曲线坐标系下的守恒型浅水波方程如下 2 1 1 ： 口 1 ，2 尝+ 豢+ 娑：o ( 4 2 ) 百+ 页+ 百2 o ( 4 2 ) j = f 绷l - 1o0 z 7 z z ，7 翟f 营乏 譬n r 、l ， 地也叽 ，f =u ， = y 中其 、，i， 口， ， 地嵋 恐 9 口卜12 z + 一 h d 。：_ht犰。喵“ 9 f扫口 他 一 研 1 ”旭 一衫地睨 观 跏跏跏 + + + 们 仍 仍 蜘蜥蜥 z z z 一 一 一 z z z 玑坍坍 ，il、 要求 的流 分裂 第4 章二维浅水方程 对于双曲型方程，为了保证计算我们常常需要使用迎风方法，但迎风方法 对流通矢量f 4 、口进行分裂【2 2 】。但浅水方程与欧拉方程不同，浅水方程 通矢量f + 、g 4 都不是基本矢量y 的一次函数，因而不能通过将雅可比矩阵 的方法来分裂流通矢量。文 2 3 1 提出了针对二维浅水方程的有效的分裂方 参照该文中的方法，我们给出曲线坐标下二维浅水方程流通矢量的分裂格 令 则满足f + = a v 。矩阵a 的特征矩阵和右特征向量矩阵分别为： 。，= ( 口喜矿弓u ，兰矿) ，冗，= ( 耋1 j 1 ：1 霎量 其中 u = ( ( 蜘一z r 蜘) + 蜘耽v 1 - - x ，i v 3 u 1 ) j ，c 宰= 、( x 2 + y 2 ) g v l 2 j j 同样，令 b ：；f 珧c 器y r x f 学- - ，x 二r y 取 警扩二篡酱妥1 珧警+ 欢( 一番+ 努)一耿嚣d + 欢嚣 则b + 满足g + = b + y 。矩阵b 的特征矩阵和右特征向量矩阵分别为： 。z = ( 彬i 矿寻叫，呈扩) ，r 2 = 1 0 | 硼4 = ( ( 可下z 一z r 珧) 一y e v 2 口，+ z 秒t ) j ，d = 、厂i 乏- i i 丽j 2 3 、l， 们了 ， 嵋 0 1 2 托什 d f_=l机旧舭 地 + f ， 嘞仇加 + 衫抛也 线纵佻 一 一 一 仉 化 协 必必纵 z z z 一 一 一 z z z 可可可 ，iiilllii l j l lg ， o 法式 叩嚣晴 一 a l蜥协晴 弩哟 十 喝哟f + 哪 舭卦一 一碧9警 蜥 咯 、liiij， 1 _ 1 - ”一c v c 9 丐g i 翌j 旦j黔筹 + 一 s l m 监 、lliij， 1 _ 1 一 一，o”一d g 虿9 i ，j i ，鲈一互一互l 一一 一 + 盟砚塑m 4 3 数值算例 针对上边我们提供曲线坐标系下的二维浅水方程流通矢量的一种表示形式， 我们给出流通矢量的分裂格式。设入七( 七= 1 ，2 ，3 ) 是a + 的特征值，可令 牡生学坐 其中是一个小参量 2 4 】。令 导曼j ) ，a 士= 兄。手r f l ，f + 士= a + 士y 有a + = a 4 + + a “，f = f + + + p 一，这样就得到了流通矢量f + 的一种分裂格 式。同理可得到g ”、g ”。 利用第三章构造的方法和本节给出的二维浅水方程流通矢量的分裂格式， 我们就可以得到求解曲线坐标网格下的浅水方程的具体方法。下一节，我们将 给出几个具体的数值算例。 4 3 数值算例 算例l ：全溃坝问题 - - 5 0 m x 2 0 m 的矩形断面规则河道，在x = 2 5 m 处发生瞬问全溃，初始上游水 深2 m ，下游水深1 5 m ，不计底面坡度。z 、y ：z 向空间步长均取为1 m ，时间步长 取为0 0 5 s ，计算