（计算数学专业论文）关于非线性双曲型偏微分方程数值解的研究.pdf

关于非线性双曲型偏微分方程数值解的研究 摘要 在航空、气象、海洋、石油勘探等方面的流体力学问题，在许多情况 下都归结成解非线性双曲型偏微分方程( 国外文献称为守恒律) 。这类方 程的基本困难是解出现了间断，当用高精度显式格式求解时，在问断处会 产生振荡。本文首先介绍了弱解、熵条件、物理解等概念和一些基本定 理，从理论上分析非线性双曲型微分方程产生间断解，以及高阶显式格式 在问断处产生振荡的原因；其次介绍目前关于这类方程数值解的一些结 果，包括显式人工粘性，隐式人工粘性和自动调节的混合格式；最后，介 绍在已有研究结果的基础上，我们给出的一种新的数值解法，即“二选一” 方法。与已有结果相比，“二选一”方法在保持精度高( 二阶精度) 的同 时在间断处不会产生振荡，并用实例说明所给解法的优越性。 关键词：双曲方程，间断，振荡，高阶显式格式，“二选一”方法 r e s e a r e ho nn u m e r i c a lm e t h o d sf o rn o n l i n e a r h y p e r b o l i cp a r t i a ld i f f e r e n t i a le q u a t i o n s a bs t r a e t m a n yf l u i dm e c h a n i c s p r o b l e m ss u c ha s a v i a t i o n ，s c e n e ，a n do i l r e c o v e r y ，e n di n n o n l i n e a rh y p er b o l i cp a r t i a ld i f f e r e n t i a le q u a t i o ns ( a ls o c a l lc o n s e r v a t i o n1 a w s ) p r e b l e ms t h eb a s i c d i f f i c u l t y w i t h t h o s ep r o b l e m sist h a tt h es o l u t i o n sd e v e l o pd i s c o n t i n u i t i e s ，i nt h i s p a p er ，c h a p t e r 1isd e v o t e dt om a t h e m a t i c a lt h e o r y w h i c hi n c l u d e w e e k l ys o l u t i o n ，e n t r o p yc o n d i t i o n ，p h y s i c a l s o l u t i o n ，a n d s o m e f u n d a m e n t a lt h e or e m t h e r e a s o nf ord i s c o n t i n u i t i e s a r i s i n g i n s o l u t i o no fh y p e r b o l i c e q u a t i o n sa n dw h yh i g h r e s o l u t i o ne x p l i c i t s c h e m e sg e n er a t eo s c i l l a t i o n sn e a rd is c o n t i n u i t i e s ，i st h e o r e t i c a l l y a n a l y z e d s o m e c u r r e n tm e t h o dsa r e g i v e n i n c h a p t e r i i t h es e i n c l u d ee x p l i c i ta r t i f i c i a lv i s e o s i t y ，i m p l i c i ta r t i f i e i a lv i s c o s i t y ， a n ds e l f - a d j u s t i n gh y b r i ds c h e m e s i nl a s tc h a p t e r ，w ec o n s i d e ra f l e w t e c h n i a u e “o n eo ft w o ”m e t h o dt h a tc a ns m o o t h0 s c i l l a t i o n sn e a r d i s c o n t i n u i t i e s 。a n dc o m p a r eo u rr e s u i t sw i t ht h e s eo ft h er e l a t e d m e t h o d s k e y w o r d ：h y p e r b o l i c e q u a t i o n s ，d i s c o n t i n u i t i es ，o s c i l l a t i o n s h i g h r e s o l u t i o ne x p l i c i ts c h e m e s ，“o n eo ft w o m e t h o d v 独创性声溺 本人声明所呈交的学饿谂文是本人在譬魉指导下进行的研究工 = 及淑褥的研究成柴。 据我所知，除了文中特别嬲以标注和致谢静地方外，论文中不包食其豫人已经发表或撰写 过的研究成果，也不包含为获得 盒腿王、业丕坐 溅其他教育机构的学位或证：瞎而使 用过的材料。与我一同工作的同志对本研究所做的任何赏献均已在论文中作了明确的说明 嚣表示溱意。 学接论文撵者签名：凳徭签字疆期2 庐秽年7 哆 强 学位论文版权使用授权书 本学位论文作者完全了解盒墅王壁太差有关保留、使用学位论文的规定，有权保留 莠离国家商关部f l 袋辊褥送交论文麓复秘箨和磁纛，允许论文搜查簇拳l 僚阕。车夫攫投垒 坦王、业本熬可以将学位论文的全部或部分内容编入有关数据库进行检索，可以采用影印、 缩印或扫描等复制手段保存、汇编学位论文。 ( 镰寮豹学谴论文在瓣密惹透蠲本授裁二器) 学位论文臻者签名：1 奚殛 签字日期：z o 够年7 月了日 学位论文作者毕业厉去向：笛等、掌赣 ： 作单位：亏舻互多承劳 通讯建戢：吾静；多式劳曩 噘疹黟；麓翕季 i i l 学9 雨签名： 签字日期 蹴：d 5 卜2 气e i 争s ， # 嚣编2 ，护矿护9 日 致谢 本文是在导师朱功勤教授、檀结庆教授指导下完成的。朱功勤教授是位学 识丰富，学术造诣颇深，他对问题的研究往往有着敏锐的洞察力，三年中他对 我孜孜不倦的教诲使我深深感动并且受益终身；同时，我很感谢檀结庆教授， 他治学严谨、平易近人，在学习上给了我很多帮助和启迪。我也衷心感谢苏化 明教授、邬弘毅教授、黄有度教授、周永务教授、汪泉副教授在学习方面对我 的指导以及给予我的关心和爱护。我要特别感谢林京博士给我的极大帮助，他 在学术上帮助我解决了许多问题，馒我顺利地进入了研究方向。 我还要感谢我的师兄李平和同学詹棠森、江平、芮义鹤，感谢他们在学习 和生活上给予我的帮助。 吴强 二零零三年五月 第一章绪论 1 1 问题的起源及介绍 在航空、气象、海洋，石油勘探等方面的流体力学问题，在很多情况 下都归结成双曲型偏微分方程( 国外文献称为守恒律) 的问题。由于这类 方程的解中常常出现间断和振荡( 通常称为激波或振荡问题) ，所以有关 它的数值解的研究直引起人们的关注。 我们以一维问题为例，简要睨明间断和振荡产生的原因。考虑如下的 双曲型方程 导+ 口罢= o x 佃o f t ( 1 1 ) a苏 。 一 ” u ( x ，o ) = 妒( x )f 1 2 ) 其中p ) 是给定的充分光滑函数。 容易验证，当日为常数时方程( 1 1 ) ，( 1 ，2 ) 的解为 “( x ，f ) = 妒( 工一a t ) 一0 0 x + c o0 0 是沿上轴正方向传播，当a 0 是沿x 轴负方向传播， 而波形保持不变。 在x i 上半平面上，方向为! ；= 口的直线为( 1 4 ) ，由于 d r d u ：丝+ 塑鱼：塑+ a 塑：o 一十一一+ 一= u d c8 l瓠d ta。x 所以“在x t 上半平面上任一斜率为口的直线( 1 4 ) 上其值为常数。当然c 不同，即直线不同，“的值也可以不一样。这簇直线覆盖了整个石，上半平 面，它可以刻画仞值问题( 1 1 ) 、( 1 2 ) 的某些特性，通常称此直线为特征 线。 在x i 上半平面上过任一点( x 。，“) 画特征线x a t = c ，即c = 一a t o 。 特征线 0 口一 0 r i r o 一 = x r 线卣 口 王个 的交点为( x 。一c l t 0 , 0 ) 。在这一点上“的值是已知的。即 u ( x o a t o ，0 ) = 妒( 工。一a t o ) 由于“沿着这条特征线的值为常数，所以 u ( x o ，t o ) = 妒( j o a t o ) 因为( 工。，气) 足任点，所以 u ( x ，t ) = 妒( 一a t ) ( 15 ) 就是( 1 1 ) 、( 1 2 ) 的精确解。 某些方程可以求得精确解的可以帮助我们验证数值解法的优越性，因 为方程( 11 ) 、( 1 、2 ) 的系数a 如果不是常数，则往往得不到类似( i5 ) 的精确 解，只能用数值解法。比如a 取未知函数“，即方程为 丝+ 甜塑；o 西a r u ( x ，o ) = e ( x ) x 忡0 ，r ( 16 ) ( 17 ) 在一f 上半平面上过任一点( x o ，。) 的特征方向是 象一( 砒) = t p ( 则有特征线 工= 妒( 工o ) r + 工o ( 18 ) 沿特征线有 a u ：0 础 即“是常数，为口( x o ) 。同样，对于另外任一点( ，( 而x - ) ，沿特征 线 工= 妒( 】) t + x 】 “也是常数，为妒( x ) 。如果有不等式 ! 堕) 二竺! 型 o 0 一1 x = 印( x o ) 七+ x o f ，_ ) 图1l两条相交的特征线 成立，那么分别通过x 。和工。的两条特征线相交( 见图1 1 ) ，交点为( 掌，玎) 。 这样来，u ( 4 ，7 7 ) 的值该如何确定。如果取p ( x o ) ，那么在过一的特征线上， “不再都是常数妒( x ，) ，“的值在( 告，叩) 处出现了间断，反之亦然。因此，在 这种情况下，无论初始函数妒( x ) 多么光滑，初值问题( 1 1 ) 、f 1 | 2 ) 的解不 再是连续可微的。这个事实特别在流体力学中有重要意义，它反映了自然 现象从连续到间断的转化。因此，如何计算间断解一直是人们关注的课题。 但是，文献 2 6 中指出，只能用一阶的显式格式( 截断误差对时间步长和 空间步长是一阶的) 计算间断解。如果采用高阶显式格式( 包括二阶) ， 那么数值解在间断处会出现振荡现象。因此，在间断处如何提高解的精确 度，引起了人们极大的兴趣。 探测这类问题始于二战时期，j o i n v o n n e u m a n n 1 8 ，2 0 ，2 1 试图用所谓 激波捕捉法去“捕捉”间断( 激波) 。虽然并不很成功，但他的激波捕捉 法的思想对这类问题的解决产生了重大影响。l a x 等人【1 2 ，1 3 ，1 4 ，1 7 提出 了弱解和物理解的理论并且对如何计算间断解进行了有益的探索。o r l enhh o a 2 7 提出了熵条件的理论，指出满足熵条件的弱解，即 物理解，是唯一的。虽然几十年来人们已经能够很好地计算间断解，但是 高阶显式格式在间断处产生振荡依然是无法克服的问题。 1 2 本文的主要内容及安排 本文主要讨论双曲方程数值解中如何计算间断解，以及在间断处高阶 ( 主要是二阶) 显式格式如何抹平振荡的问题。首先从理论上分析双曲方 程( 国外文献称为守恒律) 产生间断解，以及高阶显式格式在间断处产生 振荡的原因：其次介绍目前关于这类方程数值解的一些结果：最后，我们 给出的一种新的数值解法，即“二选一”方法，并用实例说明所给解法的 优越性。 在第二章中，给出了守恒律以及弱解、熵条件和物理解的基本概念， 以及一些基本定理。从理论上分析了双曲方程数值解中间断和振荡产生的 原因。 在第三章中，主要讨论了解决这类问题的常用方法一激波捕捉法，包 括显式人工粘性，隐式人工粘性和自动调节的混合格式。 在第四章中，给出了“二选一”方法的具体做法。“二选一”方法与已 有方法不同点就是一个二阶精度的格式能够很好地抹平在i q 断处的振荡 保持解的单调性，而我们目前见到的能够保持单调的方法都是一阶的。最 后，我们用具体例子来说明我们所给方法的优越性。 第二章关于双曲方程的基本理论 在上一章中，我们指出了双曲方程的解出现了问断a 但是对于间断解 来说，是不满足微分方程的。因此有必要拓广解必须满足的微分方程的古 典概念( 一般称为古典解) 。 2 1 守恒律与弱解 ( 臼一h 阶方程组称为守恒律( 有时也称为守恒形式或散度形式的方程) ，其 中向量= ( “”一“；) 7 是x = ( 。l 一，。) 和，的函数。而 则是x ，f ，h 的s 维 向量函数。这里”是空间变量的维数，s 是未知函数的个数。当五只依赖 于“时，方程还可以写成以下的形式： 塑+ 争4 丝：0 0 t 智0 x f 牲嘲。a 。是，关于未知向量函数“的j a c o b i 矩阵，如果对于任意的s 维实 函数u = ( 山一一，0 9 。) ，矩阵一，甜有s 个实特征值，则方程( 2 1 ) 或( 2 1 ) 。称 i t i 为双曲型方程组。如果这s 个实特征值彼此互不相同，则方程( 2 ，1 ) 或( 2 1 ) 称为严格双曲型方程组。 先研究一维( ，f 1 ) 单个( s = 1 ) 守恒律 a _ e u + 型：0 ，一 x 0( 2 2 ) 融a x 满足初始条件 u ( x ，o ) = p ( x ) ， 一 0 的上半平面中某个有界区域d 以外( 包括在区域d 的边界o d 上) 恒等于零的连续可微函数。区域d 称为的支集，有时用e 。表示。“( x ，) 满足关系式： ，蟥o 型f d x d t _ 0 v 掣口 犯6 如果实验函数集合甲，内除了p 中的函数外还包括这样的连续可微函 数p ( x ，f ) 它在f 0 的某个有界区域d 以及上除了和f - o 重合的一段以 外恒等于零，则当u ( x ，f ) 满足关系式： ，纱警“+ 芸f d x 出+ 胁，0 ) 贴舻。，v 斛口 就称( 工，f ) 为方程( 2 2 ) 满足初始条件( 2 3 ) 的c a u c h y 问题的弱解。 对于一维守恒律组 丝+ 旦业：o( 27 ) a a o 我们同样可以在分片光滑的向量函数类k 中给出三种弱解的定义来， 即如果 ( 1 ) 向量函数u ( x ，t ) e k 在其连续可微区满足微分方程( 2 7 ) ，而在u ( x ，) 的间断线x = f ( ) 上满足关系式 ( “l “一) 警_ ，_ 其中。+ 和h 一分别表示u ( x ，f ) 在间断线上两侧的右极限和左极限 f + = f ( u + ) f 一= f ( u 一) ( 2 ) 对于t 0 的半平面上与向量函数“ ，f ) 的间断线只相交于有限个 点的任意逐段光滑闭回路1 1 ，“( x ，) 满足以下关系式： ，d x f e t = 0( 29 ) ( 3 ) 设是具有紧致支集的实验函数集合，( z ，) 满足关系式： ，蟛“+ d x 删，v 。 ( 2 1 。) 上面给出了弱解的定义。但是，需要指出的是：这样定义的弱解是不 唯一的。例如对于任意常数d 1 ，定义在t 0 半平面上的单参数函数族 妒( 工，f ) = x 一三( 口一i ) f ， 1 ， 1 一口，一寺( 口一1 ) f x o ， 正 o x 妻。一i ) ， 中的每一个函数部是方程 坐+ “坐：0 3 t出 满足初始条件 贴，书拦 的c a u c h y 问题的弱解。在这个例子中弱解有无穷多个。 2 2 熵条件与物理解 上面指出了弱解是不唯一的，那么在这些弱解中，方程( 2 2 ) 、( 2 3 ) 是否 存在类似于古典解的解，并且是否唯一? o n en hh o a 2 7 提出了所谓的 熵条件，并论证了一维单个守恒律满足熵条件的弱解的唯一性。 设u ( x ，t ) 是定义在t 0 上半平面的，除了在有限条光滑曲线上间断外 的连续可微的函数。将这样的函数集合记为k 定理方程( 2 2 ) 满足初始条件( 2 3 ) 的弱解u ( x ，t ) ，如果在其间断线上满 足 f ( u - ) - f ( u * ) f ( u + ) - f ( u - ) 竺2 二! 生，v ，( 2 1 1 ) “一一c o“+ 一“一 “+ 一出 其中，= ( m i n u - , “+ ) ，m r x “- , “+ ) ) ，则这样的弱解在函数类k 。中是唯一的。 不等式( 2 1 1 ) 称为熵条件，以后用( o e ) 表示。我们以后将满足( o e ) 的 弱解称为物理解或广义解。 下面对熵条件( o e ) 作一些解释和分析。 在( “，y ) 平面上考察曲线y = ( “) ，当“+ 珊 “一时，( o e ) 可以写成 f ( c o ) 等m + ) + 芒芝m 一) ( 2 1 2 ) “一“一“ 这正表明曲线( “) 在区间，= u - , “+ ) 上位于联结 + ，f + ) 和0 - - 9 f 。) 两点的 直线的下面( 图2 1 ) ，而当“一 0 9 h + 时，( o e ) 可以写成 f ( c o ) 冬兰m 一) + 车兰他+ ) “一甜“一1 , 1 i 、 、 ： 、 、 jj f 寺古a 图2 1f ( u ) 在区间i = ( “一，u + ) 上 ( 2 1 3 ) 捌 一一 i ， 、 ，i ， ，7 7 i 1 ， l，7 l 卜j 7 ； 茸苜t 卜可 图2 2 f ( u ) 在区间，= 一，“+ ) 上 这正表明曲线f ( u ) 在区间，= ( u - , “+ ) 上位于联结 一，f 一) 和0 + ，f + ) 两点的 直线的上面( 图2 2 ) 。 熵条件中的堕型是间断传播的速度，用s 表示。在( 2 1 1 ) 左边 “。一“ 令- - + “一，同时在( 2 1 1 ) 右边令j “+ ，则得 a ( u 一) j a ( u + ) ( 2 1 4 ) 这可以理解为间断线两侧的特征线都走向间断线。不等式( 2 1 4 ) 是( o e ) 的 一个直接推论，有时也当作熵条件，记作( c e ) 。 当仁阶连续可微并厂”( “) 在区间肚不变号时，熵条件特别简单。如 果，”( “) 0 ，这时口( “) = f ( “) 是单调减的，从( o e ) 就推出“一 “+ 。又由于 厂 0 ，故对于任意的国，( 2 1 3 ) 成立，即( o e ) 成立。因此当f “ 0 时， 我们就推出 ( o e ) 寸( c g ) 斗“一 o 时，也可以顺序推出 “+ j f d e l 物理解不仅是唯一的，而且在厶范数意义下是连续依赖于初始条件 的，即是稳定的。这是因为如果h ( o ) o ( 2 15 ) o t积 出 定理如果当一0 时，方程( 2 15 ) 满足初始条件 “。( x ，o ) = “o ( x ) ( 2 t 1 6 ) 的解“，( x ，) 一致有界且几乎处处收敛到分片光滑的函数u ( x ，f ) ，则u ( x ，f ) 是 ( 2 2 ) 、( 2 3 ) 的物理解。 其中 对于多维双曲性方程 詈嘻詈脚州列，炉吣础”，e ( 0 ，r 】 口1 7 堑! 坐生：望! 坐生+ o f ( x , t , u ) 一o u 出良， 孔 缺， “( x ，0 ) ：妒( x ) 工r “ ( 2 18 ) 的c a u c h y 问题，c h k p y ) l c k 。b 2 8 】在有界可测函数类中定义物理解为 满足以下两个条件的函数“仁，砂： ( 1 ) 对于任意常数女和任意的非负的实验函数o ( x ，) 甲，都有 脚心，f ) 卅 ( u - k ) + 缸( x , i , u ) - - 脚，纠詈 一【主导，( x , t , b 1 ) + g ( 列】吖揪 0 ( 2 1 9 ) 其中积分区域仃为z ，) k r ”，os f 丁j ，曲的支集e n 。 ( 2 ) 在区间旧力上存在一个测度为零的集合r ，当，e o ，即s 时，函数 “肛在r ”上几乎处处有定义，并对任意的球彭，= 搁 ， c r 一 沏_ i ，r ) 妒( 刮出= 0( 2 2 0 ) ，1 0 、 按照c ，h ，kp y kob 的定义，当物理解“伍，口是分片光滑函数时，在 间断面上对任意的常数k ，有以下不等式： r、 s g n ( u + 一七) ( “+ 一k ) c o s ( 虬咖睡( x , t , l g * ) ，( 置f ，) 】c o sv , x i ) 、，2 _ s g n ( u 一女) ( “一k ) c o s ( y ，咖 f ( x , i , u - ) 一z ( z ) c o s ( ) l”i j ( 2 2 1 ) 其中“+ ，是f 1 1 f j 断两侧函数的极限值，v 是间断面上的法线方向。规定 从u - 指“+ 向的方向为正方向。不等式( 2 2 1 ) n t 直接从( 2 1 9 ) 得出。 对于一维情况( 即当n 取1 时) ，从( 22 1 ) ，可直接推1 t i ( 2 1 1 ) ，即( o e ) 来。 最后我们讨论维双曲型方程组的熵条件。假定一维守恒律组 等+ 掣：o ( 2 2 2 ) a fa ￥ p 是严格双曲型方程组，则矩阵 4 f h l ：塑- 盟 a ( h f ，以) 有j 个互不相同的实特征值，按大小排序为 1 ( “) 2 ( “) ，t ( n ) 由微分方程 妄吲口) 拄l ，2 ，s 定义的曲线称为方程( 2 2 2 ) 1 拘特征线。 l a x 9 称方程组( 2 2 2 ) 的弱解中的一个强阃断为激波，如果对于间断 线z = 孝( ，) ，存在某一个t o s 七s ) ，使得下列不等式成立 州“- ) 警 拍一) 引“) 警“) ( 2 2 3 ) 这时间断线x = f ( f ) 称为激波。从不等式( 2 2 3 ) n - 以得到 枷+ ) 警 o r 2 3 3 ) r 2 3 4 ) 232 ( 一2 丑 对 互 甜 1 x + z 一 一2 + 一2 ” o 的，差分方程的解也应该保持单调的性质。r ony hob c ，k f 2 6 就称这种保持单调性的格式为单调差分格式。 定理差分格式 “川= 口f “川月 ( 2 3 5 ) ， 是单调差分格式的充要条件是所有的口，0 。 此外，ro y h 0b c k 在 2 6 】中特别证明了常系数单调差分格式的 截断误差只能是一阶的。 r 、o 且y h oh ( g o d u n o v ) 定理如果( 2 3 5 ) 是- - 个和( 2 3 3 ) n 容的，并且 极大范数有界以及至少有两个系数是非零的差分格式，那么格式( 2 3 5 ) 是一阶精确度的。 o j e n n i n g s 10 将r oa y hob 的常系数差分格式的概念推广到非线型 差分格式上。他称单个守恒律( 2 2 ) 的守恒型差分格式( 2 2 8 ) 是单调的， 如果 l j n + i ( “川，“川，。) 对其每一个分量的导数都大于等于零，即 盟o m = - t , - t + 1 ，f o u m ， ( 2 3 6 ) 和常系数单调差分格式一样，h a r t e n ，h y m a n ，l a x 5 等x i i e 明了非线型单 调差分格式的截断误差也只能是一阶的。 定理设守恒律( 2 2 ) 的守恒型差分格式( 2 2 8 ) 是单调的，n ( 2 2 8 ) 的截 断误差是一阶的。 在【5 】中h a r t e n 等人还证明了一维当单个守恒律的单调守恒型差分 格式的解，如果几乎处处有界收敛，则极限函数是物理解。 c r a w l 等人 2 证明了多维单个守恒律的单调守恒型差分格式的解一 定收敛于物理解。 正是由于保持单调的守恒型差分格式的截断误差只能是一阶的，所以 当用二阶的l a x - w e n d r o f f 格式解方程( 2 2 9 ) 、( 2 3 0 ) 时，得到的是( 2 3 2 ) ， 而不是物理解( 2 3 1 ) 。订三式山于高阶显试格式不能保持在间断解的单调性， 因此它们必定会在间断处产生振荡，从( 2 3 1 ) 可以很明显地看出这一点。 第三章数值解法概述 上一章从理论上解释了无论双曲方程的初始条件多么光滑，它的解可 能会产生问断，以及高阶显式格式在间断处产生振荡的原因。那么如何用 数值方法计算双曲方程间断解，以及在间断处如何提高差分格式的精确度 的同时又能够保持解的单调性，是人们面j | 刍的首要任务。目前大多数采用 的是激波捕捉法。 激波问题创始人v o n n e u m a n n 2 1 1 提出了“捕捉”激波的思想，即： 无论双曲型方程是否有间断解，都不加区别地按统一的格式进行计算。后 来人们把这类方法统称为激波捕捉法。采用激波捕捉法，间断是作为解的 一部分进行计算的。问断解被表示为具有一定过渡层的连续解，这个连续 解基本上能反映出本来的间断特征，这类方法的最大优点是简单。由于计 算间断不作特别处理，所以程序编制容易，机时也省。其次，守恒型差分 格式收敛性已有保证。当然这类方法也有缺点，主要是有时把间断的过渡 区拉得过宽，有的虽然抹平了振荡，但精确度都是一阶的，而高阶的方法 在间断的前或后都会产生不应有的振荡。 激波捕捉法主要包括显式人工粘性、隐式人工粘性和自动调节的混 合格式。 3 1 显式人工粘性 显式人工粘性方法是在要求解的方程中另外加上一定的扩散项，比如 a 2 t ， 等。在数值平滑计算中，这种人为引进的项起着粘性作用，所以称作 “工 人工粘性项。其效果是把间断解化为具有过渡层的连续解。上述r 是一个 小参数，它随步长r ，h 趋于零面趋于零。 这种方法是 c o nn e u m a n n 和r i c h t m y e r 2 i 为求解气体动力学方程组 最先使用的。在l a g r a n g e 坐标下气体动力学方程为 堡一塑：o 西苏 塑+ 望：o 甜缸 堡+ 。丝：0 魂。瓠 其中v 为气体比容，“为速度，p 为压力，e 为比内能。 v o n n e u m a n n 和r i c h t m y e r 引入人工粘性项g 后，方程组为 生一塑：o e t瓠 丝+ ! ! 旦盟：o 西ao 宴+ ( p + g ) 娑= o 对于人工粘性q ，司选择如f 形式 ：一1 2l , _ c a u ，, 2 l , 0 = u 一, q j2 一l _ ，l = 一j 其中 ( 参= 1 ，如果罢 o 肌瓦j - 川禾瓦 u l 是具有长度量纲的正常数。 下面我们给出( 3 1 ) 的一个差分格式 ! 比：l 二坚。! 型! ：! 二! ! 二坐 f h 堑塑! 堑! 型一! 业：! 二坠： f h = 0 ( 3 1 ) = 0 ( 3 2 ) 生字啊忆。竿2 o “+ 1 4 一“月0 + 1 ，一“m 0 r 为网格时间步长，h 为空间步长。 引进人工粘性项是重要而技巧性很强的方法，它不应使间断强度过于 减弱，也不应使过渡区拉得太宽。如何构造出好的k - r a 性，自v o n n e u m a n n 和r i c h t m y e r 的开创性工作以来很受人们的重视。 3 2 隐式人工粘性 隐式人工粘性不是在微分方程中加一个另外的项来使间断具有光滑 的过渡，而是对微分方程进行离散时，差分格式本身就带有类似于粘性这 样的项。因此也称为隐式人工粘性为格式粘性。有了这样的项之后，就可 去 以使间断有一个光滑的过渡。下面我们通过一些例子来说明某些差分格式 确实存在这样的项。 首先我们考察逆风格式。为了简单，我们取常系数双曲方程 竺+ 口竺：o ，口 of 3 3 1 a孤 ? 其逆风差分格式是 坠型二堑+ d ! ! 二堑：! ：o fh ( 3 4 ) 这个格式的截断误差是阶的。为了进行考察，我们把改写成( 3 4 ) 下面的 形式 坠： 坚+ 口堕：! 二垡：! f2 丝堑：! 二塾! 堑：1 2h 2 上式右端项可以看作( ) 器2 的离散，因此逆风差分格式可以看作微分方 程 a “a “ a “ + a 一= s ( 厅) 了 a苏氖2 显示格式。占( ) 窑称为隐式人工粘性项或格式粘性项。我们注意到，当 0 x 。 h 斗o 时( ) 一0 。但在具体计算中，h 0 。因此格式粘性项是会起作用 的。 现在我们考察逼近守恒律 塑+ 型：0( 3 5 ) 甜出 的i ，a x f r j e d r i c h s 格式 垫竺：竺! ! + 丑! 二红! r 2 h = 0( 3 6 ) 这个格式的截断误差对时问步长f 是- g r n ，对空间步长h 是二阶的。上 式可以写成 因e v a ( 3 6 ) 和( 3 7 ) 可阻看作微分方程 堑! 二! 弩! 竺坐( 3 7 ) h 2 詈+ 掣叫n ，磐0 xa反 、。 2( 3 8 ) 的简单显式格式，其中s ( ) ：笔。s ( ) 鲁就是格式粘性项。正因为有了 z r 0 霄。 它，可以用l a x f r i e d r i c h s 格式计算间断解。需要注意的是，由于 l a x f r i e d r i c h s 格式对时间步长r 的截断误差是一阶的，精度不高，从一 些实际例子可以看出，在解的间断处，连续的过渡区比较大。 最后我们简单考察一下l a x w e n d r o f f 格式。对于常系数对流方程组 o u + a o u ：0 西融 ( 3 9 ) 的l a x w e n d r o f f 格式是 1 厂爿去( 。叶o + 妄以一2 m + h j - l , n ) 这个格式的截断误差是二阶的，它相当于用高于二阶的精度逼近方程组 坐+ 4 坐：d “ 出揪一 其中 如一i ia ( h 2 _ r 2 a 2 ) 窘一扣1 1 2 - r 2 a 2 ) 窘 由此可知，利用l a x w e n d r o f f 格式进行计算相当于逼近一个具有小参数 项的微分方程。从而l a x w e n d r o f f 格式也包含了隐式人工粘性项。 由于l a x w e n d r o f f 格式的截断误差是二阶的，因此它在间断解的计 算中经常使用。需要注意的是，它的计算结果具有明显的振荡，比如2 3 中的例子。 3 3自动调节的混合格式 我们从前面的讨论知道，利用阶精度格式( 如l a x f r i e d r i c h s 格式) 计算间断解，间断的过渡区比较宽。而在解的光滑区，计算精度不高。但 是计算得到的解不会出现振荡。用高精度格式( 如l a x w e n d r o f f 格式) 进行计算，在间断前后产生振荡，但间断的过渡区比较窄并且在光滑区精 度高。 a h a r t e n 6 提出了一种自动调节的混合格式，这种格式是在解的光滑 区用高精度格式以提高计算精度，而在解的间断处则用一阶精度格式以消 除振荡。它们之f 司的转换是使用丌关函数来自动完成的。 设三和k 是逼近方程组 o u 研+ 掣：o ( 3 1 0 ) 品缸 7 两个差分格式，并且上，是一阶精度的，。是r ( r 2 ) 阶精度的。令 ( 口) ，= 一一五( 辟：，：一碍：，：)( 3 1 1 ) ( 工。“) = h 一五( ，l j 2 ：一厅j 2 ：)( 3 1 2 ) 其中兄= 是时问步长与空间步长之比，即网格比。自动调节的混合格式 定义如下 ( 工h ) = “一a ( j + l ，2 一h 一i ，2 ) ( 3 13 ) 其中 h 卜，：= o + 。，2 矗j ? ，2 + ( 1 一嘭+ ，2 ) ，l j 嚣2 0 臼1 。口就是自动开关，其构造方法如下：在间断的过渡区取0 ；】， 而在解的光滑区取0 z0 。这样一来，在光滑区和间断的过渡区就分别用 高精度格式和一阶精度格式进行计算。 可以把混合格式( 3 1 3 ) 改写成 ( 工) ，= ( 工。) 厂五咿j + l 2 ( 螂，：一 2 ) - o j 。：( ，j ”- i ，2 盘：) ( 3 1 4 ) 我们看到，在解的光滑区取0 = o ( a x ) ，p r 一1 ，那么有 ( 上“) ，= ( 月) j + 0 ( 船) 于是，在光滑区内l 和。的截断误差是同阶的。 从( 3 ，13 ) 我们还可以看出，混合格式可以看作l 。加上了一个人工粘性 项。作为例子我们取。为l a x w e n d r o f f 格式，此时r = 2 ， ( l 2 u ) ，= 旷兰( “一几) + 等 4 二。( t t j + 。- t t ) ) 一4 m u ) - t 1 ) _ t ) ( 3 15 ) 其中 一( ) =a ( ， o ( u i ， ，正) ，) 爿( 毕) 取厶为l a x - f r i e d r i c h s 格式 ( “) ，= 百1 ( “川+ “川) 一鲁( + 一 一，) ( 3 1 6 ) 为了便于比较，可以把上式改写为 ( 厶) ，= u 1 一鲁( + i 一 一1 ) + i 1 ( “川一2 “j + “川) ( 3 1 7 ) 这样有 = 去( + 乃+ 。) 一害爿二。( “川一h j ) ，l ：= 丢( + 厶+ ，) 一鲁( “，) 因此混合格式可以写作 ( 上h ) = ( 工2 h ) ，+ 去【目p i ，2 ( j 一 2 q ，2 “，2 ) ( h p l 一1 , 1 1 ) 一o j - t l2 ( ，一 2 爿卜2l ，2 ) ( “一u j - i ) 】 其中，表示单位矩阵。于是，构造混合格式就相当于在l a x w e n d r o f f 格式 中 加上了如下的人工粘性项 垡墅呈一t 2 az ) 当 2d rc ! r 可以看出，如果0 = d ( 缸) ，那么和l 2 是有同样精度。 下面我们考虑混合格式的稳定性。取0 为常数，0 0 1 ，那么混合 格式( 3 13 ) 就化为 上= o l l + ( 1 一o ) l w 由此可以得到 = o l l l 。i i + ( 1 一o ) l l l 。 因此，如果厶和l 。是稳定的，那么l 也是稳定的。而且l 的稳定性条件是 和。中较严的一个。作为例子，考虑常系数线性方程组，即方程( 3 9 ) ， 爿为p x p 常数矩阵。k 取l a x w e n d r o f f 格式( r = 2 ) ( l 2 u ) ，：“一要( r 一，) + 冬 爿( 一“，) 一4 ( j - u i _ , ) 1 我们知道，其稳定性条件是 五p ( 4 ) 1 l ，取一阶精度格式 ( 工1 h ) j = ( 2 “) j + i ( “卜l 一2 “+ h 一l 可以求得其稳定性条件是 如( 爿) 塑2 ( 3 18 ) 对应于上述和厶的混合格式三是 1 ( l u ) j2 ( 三2 “) + 吉 臼卜1 ，2 ( “j + l u j ) 一8 j - i 2 ( “，- - i , i 卜1 ) 】 显然，上在条件( 3 1 8 ) 之下是稳定的。 最后，我们考虑怎样来选取自动开关函数0 ，方便起见仅对标量情况 进行讨论。设a ( u ) 是“的一个具有间断的函数。令盯，= g ( u ) ， a ，+ l ，2 盯= 盯，+ 1 一o - ， 只= 卜。刮一卜。刮1 晤。：0 1 + k 抛。1 l o f ：0 1 知：和占 令w 爿叫知2 c f 占 其中s 是一个适当选取的常数。容易看出，当解光滑时，有色= o ( ( 缸) ) 。 由百，的定义有o 西j 1 ，而且当a j + l 2 盯= o 和i 川，：仃l s 或者当i a ) + i 2 盯l f 和a l _ f 2 t y = o 时有幺= 1 。现在定义自动开关函数如下 0 ，2 = m a x ( 舀，亩川) 由此得出，0 0 n ，2 蔓1 ，满足我们的要求。 注：1 9 7 2 年r i c h t m y e r 2 4 1 认为激波拟合法是解决一维问题的最好方法。 m o r e t t i 1 9 】做出了在二维问题的网格上如何拟合一个单激波。 激波拟合法的思想是在解的光滑区采用高阶格式，而在间断处采用熵 条件( 2 8 ) 或r a n k i n e h u g o n i o t 关系式【1 4 来连接间断两边的间断点。这类 方法的困难在于首先要找到间断点，其次是该如何连接光滑区和间断区， 而且得到的不是一个统一的格式，分柝该格式的稳定性是非常困难的。尽 管如此，激波拟合法毕竟是对如何提高间断处的精确度、抹平振荡进行了 可贵的探索，由于资料缺乏，这罨不作详细讨论了。 9 第四章 一种新的数值解法一一“二选一”方法 前一章给出了计算双曲方程间断解的些常用方法。然而，这些方法 都无法解决高阶显式格式在间断处的振荡问题。为了解决这个问题，我们 提出了一个“二选一”方法，这个方法是在每一网格点上用两个高阶格式 分别进行计算，得到两个数值解然后根据给出的选取准则，选取其中之 一作为该点的最后数值觯。这个方法可以在间断处不会产生振荡，保持了 解的单调性。最后用些实例说明了它的优越性。 4 1 一维闯题 考虑如下的一维对流方程 塑+ 口竺= 0一 工 懈0 fs7 1 8 ca x “( 工，0 ) = 缈( ) 其中“= 撕，“，百g 3 1 , g ，璺o x ， 对于区域：= f f x 选取网格如下： ，= ， ： 【，。= r 是初始条件，u ( x ，) 是待求函数。 ，t ) i 一一 0 一4 、。巩，+ 竿( 1 1 j + l , n - - l d j , a ) 一等( 1 一巾厂2 ”) 口j ，。 0 ( 45 ) 格式i ，i i 的截断误差都是二阶的，稳定性条件是i 竿卜1 。如果( 4 ，1 ) 中的a 含有 ，则格式( 4 3 ) 一( 4 5 ) 的口。中作下面的代换 一。丝 缸 。万1 ( “巾一“巾) 二选一 在每一个网格点( x ，t 。) ，分别用上面的两个格式进行计算，得到“，+ 的两个数值解。将它们分别代入 取使该值最小的一个作为“。+ 的最后解。 数值例子 我们选择如下的对流方程： 竺+ 0 0 0 5 工竺= 0- - o o 工 栅0 f t a苏 初始条件选作 “( x ，。) = e l ”i 。_ e 净- l o o t 时间步长r ：0 0 9 ，空问步长h = 0 i 。r 。= 4 5 对，用l a