（计算数学专业论文）level+set方法在双曲守恒律中界面追踪的运用研究.pdf

摘要 在一维双曲守恒律方程的的数值计算中，众多学者研究并建立了各种各样的 差分格式，这些格式不断趋于成熟，分辨率更强，精度更高。e n o 格式和w e n o 格式是目前应用非常广的高分辨率格式。但是，我们知道，即使用这两种格式来 计算双曲守恒律方程的初值问题，在间断面的附近仍会发生数值耗散，也就是说 在间断处的分辨率降低了。 19 8 8 年，o s h e r 和s e t t i a n 提出了l e v e ls e t 方法，在许多复杂的界面追踪问 题中得到成功应用。他们把运动界面定义成l e v e ls e t 函数的零等值面，然后始 终保持它是零等值面。l e v e ls e t 方法能够准确地追踪运动界面。 本文将高分辨率激波捕捉格式与l e v e ls e t 方法结合起来计算双曲守恒律方 程( 组) 。对标量守恒律方程、守恒律方程组分别构造了一种虚拟区域，将一个 方程( 组) 转化成两个方程( 组) ，对这两个方程( 组) ，我们仍然使用高分辨率 激波捕捉格式，而l e v e ls e t 方程用来追踪间断的位置，原方程( 组) 的解最后 由l e v e ls e t 函数决定；这样做弥补了高分辨率激波捕捉方法在间断附近发生数 值耗散的缺陷，提高间断处的分辨率。通过实验可以看出，对标量双曲守恒律方 程应用l e v e ls e t 方法后，我们得到了锐利的图像，界面得到了准确追踪，而不 像高分辨率激波捕捉方法，在界面附近有一个过渡区域。我们把l e v e ls e t 方法 运用到守恒律方程组中，分别对接触间断和激波进行追踪，均取得了显著的效果。 关键词：双曲守恒律e n o 格式w e n o 格式l e v e l s e t 方法虚拟区域 a b s t r a c t f o rt h en u m e r i c a lc o m p u t a t i o no fo n e d i m e n s i o n a lh y p e r b o l i cc o n s e r v m i o nl a w s ， v a r i o u sd i f f e r e n c es c h e m e sw e r ee s t a b l i s b e d t h e s es c h e m e sa r eb e c o m i n gm o r ea n d m o r ep e r f e c tw h i c hh a v eh i g h e rr e s o l u t i o na n da c c u r a c y e n oa n dw e n o a r et w o t y p e so fr e p r e s e n t a t i v e so fh i g h r e s o l u t i o ns c h e m e sn o w h o w e v e r w h e nw eu s et h e s e s c h e m e st oc o m p u t et h ei n i t i a lp r o b l e mo fh y p e r b o l i cc o n s e r v a t i o nl a w s ，t h e r ei ss t i l l n u m e r i c a ld i s s i p a t i o nn e a rt h ei n t e r f a c e ，t h a ti st os a y , t h er e s o l u t i o ni sd e c r e a s e dn e a r t h ei n t e r f a c e i n1 9 8 8 ，o s h e ra n ds e t h i a np u t e df o r w a r dl e v e ls e tm e t h o da n d g a i n e ds u c c e s s i nm a n yc o m p l i c a t e di n t e r f a c e - t r a c k i n gp r o b l e m s t h e yd i f i n e da n dk e p tt h em o v i n g i n t e r f a c ea sz e r ol e v e lo fl e v e ls e tf u n c t i o n l e v e ls e tm e t h o dc a nm i c kt h ei n t e r f a c e a c c u r a t e l y i nt h i st h e s i s ，w ec o u p l eh i g hr e s o l u t i o ns h o c k c a p t u r i n gm e t h o da n dl e v e ls e t m e t h o dt oc o m p u t eh y p e r b o l i cc o n s e r v a t i o nl a w s f o rs c a l a re q u a t i o na n ds y s t e mo f e q u a t i o n s ，w eb u i l dd i f f e r e n tg h o s tf i e l d s ，t r a n s l a t eo n ee q u a t i o n ( s y s t e m ) i n t ot w o e q u a t i o n s ( s y s t e m ) w es t i l lu s eh i g hr e s o l u t i o ns h o c kc a p t u r i n gm e t h o dt oc o m p u t e t h et w oe q u a t i o n s ( s y s t e m ) ；l e v e ls e te q u a t i o ni su s e dt ot r a c kt h ei n t e r f a c e ，a n dt h e r e s u ao fo r i g i n a le q u a t i o n ( s y s t e m ) i sd e t e r m i n e db yt h el e v e ls e t f u n c t i o n t h u s w e e l i m i n a t et h en u m e r i c a ld i s s i p a t i o nw h i c hh i g hr e s o l u t i o ns h o c kc a p t u r i n gm e t h o d c a n n o ta v o i dn e a rt h ei n t e r f a c e ，a n dt h er e s o l u t i o ni se n h a n c e d f r o m t h ee x p e r i m e n t s i nh y p e r b o l i cc o n s e r v a t i o nl a w s ，w ec a l ls e et h a ti fw e o n l yu s eh i g hr e s o l u t i o ns h o c k c a p t u r i n gm e t h o di ns c a l a re q u a t i o n ，t h e r ew o u l dh a v eat r a n s i t i o nn e a rt h ei n t e r f a c e ： b u tw h e n c o u p l e d w i t hl e v e ls e t m e t h o d ，w ec a ng e tas h a r pi n t e r f a c e ，t h ei n t e r f a c e i st r a c k e d e x a c t l y ，a n dw h e nw eu s el e v e ls e tm e t h o dt o t r a c kt h ec o n t a c t d i s c o n t i n u i t y , t h es h o c k r e s p e c t i v e l yi ns y s t e mo f e q u a t i o n ，w ea l lh a v e g o o dr e s u l t s k e yw o r d s ：h y p e r b o l i cc o n s e r v a t i o nl a w s ，e n os c h e m e ，w e n os c h e m e l e v e ls e t m e t h o d ，g h o s t f l u i d 两北工业大学硕士毕业论文 第一章绪论 1 9 1 7 年，英国科学家l f r i c h a r d s o n 作了一个预报天气的报告，他是第一 个通过数值求解偏微分方程来预报天气的科学家。他揭开了计算流体力学、科 学计算的序幕。随着计算机科学和技术的日新月异的迅速发展，计算机已经成 为探索流体力学中的各种物理现象以及应用计算流体力学方法解决各类工程实 际问题的强有力的物质基础和重要手段。随着计算方法及其相关技术的不断改 进，计算流体力学已发展成了一门重要的学科。 计算机问世以前，虽然有些先驱者曾探讨过用数值计算方法求解流动问题， 但只能在计算方法上做出有意义的研究成果。在这里应该特别提到c o u r a n t ， f r i e d r i c h s 和l e w y 的工作，他们证明了连续的椭圆型、抛物型和双曲型方程 组解的存在性和唯一性定理；针对线性方程的初值问题，首先将偏微分方程离 散化，然后证明了离散系统收敛到连续系统，最后用代数方法确定了差分解的 存在。他们还讨论了双曲型方程的特征性，提出了特征线方法，给出了著名的 稳定性判别条件c f l 条件。这些工作，结合其他一些数学家研究的偏微分 方程数学理论，构成了有限差分方法的数学理论基础。随后，v o nn e u m a n n 、 r i c h t m y e r 、l a x 等研究并建立了非线性双曲型方程守恒律的数值方法理论，特 别是弱解的理论，为含有激波以及其他间断的气体流动数值模拟打下了理论基 础。y o nn e u m a n n 还提出了线性稳定性的分析方法，提供了分析线性稳定性问题 的较简单的实用的方法，直到目前在计算流体力学中分析线性稳定性时仍广泛 地使用着这些方法。随后l a x ，k r e i s s 等给出了非定常偏微分方程差分逼近的 稳定性理论，进一步促进了双曲型方程或双曲一抛物型方程的时间相关方法的发 展。时间相关法的基本思想是从非线性e u l e r 或n - s 方程出发，利用双曲型方 程或双曲一抛物型方程的数学特性，沿着时间方向推进求解，由此得到对于时间 t 趋近于无穷大的渐近解即为所要求的定常解。此方法既能取得流动定常解， 又能模拟流体运动的非定常过程；若流场中存在间断面激波，它能自动捕 捉到激波，因而是应用范围极广的一般性方法。关于激波捕捉，l a x 在1 9 5 4 年 1 西北工业人学硕士毕业论史 首先提出了一种采用守恒形式的主控方程而不对激波作特殊处理的计算带有激 波流动的方法，v o nn e u m a n n 和r i c h t m y e r 则明确提出用人工粘性的方法计算 带激波的流动，他们这些工作是“激波捕捉”法的起源。随后g o d o n o v 提出了 迎风格式，l a x 和w e n d r o f f 提出了二阶精度的差分格式，改进了带激波流动的 数值计算的精度。与激波捕捉技术发展的同时，g r a y 提出了激波拟合技术的想 法，可避免激波捕捉技术中激波被抹平的缺点。m o r e t t i 等应用和发展了这种 技术，计算了多维超声速流动，但相比于激波捕捉法来说，其计算复杂的多，因 而目前主要还是应用前者计算带有激波的流动。 经过这些年的努力，计算流体工作者发展了相当数量的高精度、高分辨率 的差分格式，如总变差减少( t v d ) 格式，本质无震荡( e n o ) 格式，无波动、 无自由参数的耗散( n n d ) 格式，耗散比拟方法以及矢通量分裂格式等，形成了 第二代差分格式。这些格式可以模拟包含激波、漩涡等现象的非光滑流场。9 0 年代以来，为进一步提高对粘性流场的高分辨率，人们正努力研究和发展更高 精度( 二阶以上) 的计算格式和方法。h a r t e n ，e n g q u i s t ，o s h e r ，c h a k r a v a r t h y ， s h u c - w ，l i u ，x 一d ，e t a d m o r 等在这方面做了大量的工作 1 1 5 。尽管 如此，间断面附近的数值耗散现象仍没有消除，分辨率尚待提高。 1 9 8 8 年o s h e r 和s e t h i a n 首先提出了l e v e ls e t 方法 1 6 。其基本思想是： 把一个隐式表面定义成一个符号距离函数o ( 哥) 的零等值线。l e v e ls e t 方法在 最近1 5 年中得到了人们的关注，被广泛地应用于流体力学、图象处理、计算可 视化等方面 1 6 儿1 7 1 8 。其主要优点是运动前沿或问断面不用显式追踪，而 是用一函数来定位，通过此函数，我们可以知道间断面的准确位置。假定r ( t ) 是 一运动界面，n ( t ) 是界面包围的区域，我们构造一函数中，即l e v e ls e t 函数， 巾是l i p s c h i t z 连续的且满足下列方程： f 驴( x ，t ) 06 2 ( 0 区域外 这样我们就可以直接追踪到包括激波在内的间断。相对于激波捕捉方法来说， 两北e 业大学硕士毕业论文 结合l e v e ls e t 方法后，激波附近的过渡区域变得更窄，从而提高了分辨率。 以上各种方法我们将在后面作详细的叙述。 本文主要研究l e v e ls e t 方法在双曲守恒律中的界面追踪。对标量双曲守 恒律方程来说，运用l e v e ls e t 方法后，界面得到了准确追踪，而不像高分辨 率方法，在界面附近有一个过渡区域；把l e v e ls e t 方法运用到守恒律方程组 后，我们分别对接触问断和激波进行了追踪，均取得了显著的效果。 本文安排如下：第二章介绍拟线性偏微分方程的基本理论；第三章主要讲 各种高分辨率激波捕捉格式及l e v e ls e t 方法；第四章介绍本文的数值离散格式； 第五章研究l e v e ls e t 方法在标量双曲守恒律中的应用；第六章研究l e v e ls e t 方法在双曲守恒律方程组中的应用。 西北工业人学硕士毕业论文 第二章双曲型偏微分方程的介绍 本章主要介绍一下双曲型偏微分方程的基本性质。这些特性对于分析流体 以及构造数值方法具有非常重要的作用。有关偏微分方程的基本理论，读者可 参阅以下文献 1 9 、 2 0 以及 2 1 。在这里我们主要研究双曲型偏微分方程以及 双曲型偏微分方程组，有下述几个主要原因：( i ) 在不考虑粘性和热传导的情 况下，可压缩流体的方程可简化成双曲型方程；( i i ) 从数值的角度来看，流体 偏微分方程中双曲项的离散要求要比其它项严格的多：( i i i ) 有关双曲型偏微分 方程及双曲型偏微分方程组的基础理论要比其他数学模型( 像n a v i e r s t o k e s 方 程) 完善的多。另外，近些年中，众多科研工作者们在双曲型问题方面作了大 量的研究，取得了丰硕的成果。 第一节拟线性偏微分方程的基本概念 我们考虑如下形式的方程 堕a t + 静埘儿) 等哪 以) - o ( 叫，( 2 1 ) 其中“，是x 、，的函数，即“，= 。( x ，r ) ，婴，婴分别是珥关于，、x 的偏导数。 o fo x 如果令 u = 睦 b = 则方程( 2 2 ) 可写成如下形式 u ，+ a u ，+ b = 0 b l 6 2 ，a = 口1 口2 m ： 口m ( 2 2 ) ( 2 3 ) 西北工业大学硕士毕业论文 如果矩阵a 中的每个元素口。和向量b 中的每个元素6 ，均为常值，则我们称方 程( 2 3 ) 为线性常系数方程：如果口，= 口，( x ，) ，b ，= b ，( z ，r ) ，则方程( 2 3 ) 为线性 变系数方程；如果矩阵a 是向量( ，的函数，即a = a ( u ) ，那么方程( 2 3 ) 为拟 线性方程( 一般来洗，拟线性方程是非线性方程) ：如果b = 0 ，则方程( 2 3 ) 称之为齐次方程。 线性对流方程 罢+ 口_ o u ：0( 2 4 ) 西反 。 和无粘b u r g e r s 方程 孚+ “娑：0 ( 2 _ 5 ) 研融 。 是方程( 2 1 ) 的两个简单例子，也是数值方法中经常谈论的两个方程，是验证 我们数值格式是否有效的经典方程。 定义( 守恒律)守恒律是指可以写成以下形式的方程或方程组 u ，+ f ( u ) ，= 0 ( 2 6 ) 此处， 阽臣 f i f ( ：l 参 | = i s ( 2 7 ) 我们把u 称为守恒变量，f = f ( 称为矢通量，( 的每一个分量，均是守 恒变量各元素“，的函数，即，= ，( “，“：，“。) 。 西北工业人学硕士毕业论文 定义( 雅可比阵)方程( 2 6 ) 中通量函数f ( u ) 的雅可比阵是指 4 r u l 一o f 、。 a u 萌奶甄 抛l动2a “， 识蔹荫 乩t加2孔。 可m 耐。茸m a “1a “20 “。 矩阵4 ( 【，) 中每个元素口。是矢通量f ( 中元素，关于相应的守恒变量u 中的 元素“，的偏导数，也就是= 暑。 我们可以看到，如果b = 0 ，通过 o f ( u ) o f o u _ 一 良a u 反 方程( 2 6 ) 可以写成如下拟线性形式 u + a ( u ) u ：= 0 ( 2 8 ) 这是拟线性方程( 2 3 ) 的一个特例。相反，线性对流方程( 2 4 ) 和无粘b u r g e r s 方程( 2 5 ) 可写成守恒律形式 o u 。+ + 掣_ 0 ，m ) 喇 ( 29 ) o ro x 詈+ 掣地m ) = 丢“2 ( 2 1 0 ， 定义( 双曲型方程组) 如果矩阵彳( z ，f ，“) 有m 个实特征值五。，五2 ，五。和 相应的卅个线性无关的右特征向量k ”，足”，彤【，我们说方程( 2 3 ) 是双 曲型的。如果m 个特征值互异，则称方程( 2 3 ) 为严格双曲型的。 西北工业大学倾士毕业论文 定义( 特征值)矩阵4 的特征值旯，是指特征多项式 l 以一“ = d e t ( a 一加) = o ( 2 1 1 ) 的解。这里，是指单位矩阵。方程( 2 3 ) 中的系数矩阵4 的特征值也常称为 方程( 2 3 ) 的特征值。 定义( 特征向量) 如果向量膏= ( 七：，) 7 满足彳髟= 五k ( “，那 么k 。就是一的属于特征值 的右特征向量：类似的，如果三2 a = 五三，则掣) 称之为a 的属于特征值 的左特征向量。 对于标量方程，我们易知，线性对流方程( 2 4 ) 的特征值五：“：无粘 b u r g e r s 方程( 2 5 ) 的特征值五= “。下面我们看一下一维欧拉方程组 一维可压缩e u l e r 方程的守恒形式是 u + f ( ，= 0 ( 2 1 2 ) u 措小= 暖 其中，p ，甜，p 和e 分别表示密度，速度，压力和单位体积的总能，。，日表示 比内能和总焓a e = p ( e + u 2 2 )，h = ( e + p ) p 。对理想气体， e2 么0 一1 ) p ) ，其中，y = ，c p q 分别表示等压和等容比热，本文中y = 1 4 。 方程( 2 1 2 ) 的雅可比矩阵a 为 a ( t o = o 妒s l 掣+ ( y p p 2 + 倒 ，。 o _ “列 y r飞莉惜 。玎卜 卜l 尹 丝q 丫_ 、心百 也一啦 d 西北工业大学硕一l 毕业论丈 其特征值为“一a ，“，u + a ，相应的左右特征向量为 三( = 口“2 2 4 日 k ( 1 j = 三( 2 ) = 在燧蜥憨d = 居。 2 a u2 2 a h ( i u 一2 d h 一2 a u 口“2 2 a l l 2 a 口“2 2 a l l 1 “ 1 2 一“ 2 ( 3 ) = 第二节r i e m a n n 问题 一! “3 + h h 一土。“2 至至 口“2 2 a l l 土“2 一h + 口“ 2 为什么人们都喜欢讲“非线性”? 其实线性的问题，如线性波、线性表示 等等，也都是现实非线性的近似，是人们处理问题的一种简化方式。越接近非 线性，就越走向现实世界。非线性现象，非线性问题是实际，那么要处理非线 性，就会遇到许多困难，如梯度的变化、奇性、激波、涌浪、介质面、运动界 面等等间断问题。 数值离散的思想似乎有两种观点，一是把间断看成梯度很大的连续性情况， 通常的许多差分方法、有限元方法实际就是这样，它利用格式的数值耗散效应 把问断“光滑”了些许，从而使得数值解得以实现。特别是过去的“人工粘性 法”，更是显而易见。另一种观点是把问题处处当作间断处理，这就是由g 1 i m m 等人于1 9 6 5 年提出的所谓r i e r n a n n 问题思想 2 2 ，即把离散后的问题，对于每 一个网格的交界面都看作是间断点( 面) ，从而处处需要进行间断的分解和装配。 r 口一 一 之丽竺柑害等 一 一“ 2 2 一一k l “ 一2 一 一 一 西北工业大学颂士毕业论文 特别是c h o r i n 的随机选择法 2 3 ，r o e 提出的参向量解算子法 2 4 2 5 ，v a n l e e r 的m u s c l 方法 2 6 和c o l l e l a 等人的p p m 方法 2 7 等等。 下面我们以后一种观点介绍这种r ie m a n n 间断分解的思想和方法，以及如 何求解，这对于我们理解和构造各种高分辨率方法有很大的帮助。 以守恒律方程的柯西问题为例 “，+ 厂( “) ，= 0 如果给它左、右不同的常数态初值条件 出，o ) - 心x ? 【u r x 0 ( 2 1 3 ) ( 2 1 4 ) 一般对于这种左、右不同常数态的守恒律初值问题，我们称为r i e m a n 。问 题。这显然是一种间断分解问题，因为它理应满足r a n k i n e h u g o n i o t 间断关系 s u l = ( “) 】，或j ( “l 一“月) = 厂( “ ) 一f ( u r ) ( 2 1 5 ) 其中，j 表示该间断线的斜率或激波的传播速度；或者必须是一种中，巴, n n n 的结构，符合中心稀疏波的相似解特征。 o 图2 1 r i e m a n n 问题 这种r i e m a n n 间断分解问题，本来可以利用特征线方法及特征关系，并结 合r a n k i n e h u g o n i 。t 关系，进行特征装配方法求局部解析解的。如图2 1 所示， 9 一 西北工业大学硕士毕业论史 首先判断间断分解的结构，如果是激波结构，就找到激波位置及发展方向，而 后确定左右状态；如果是中心稀疏波结构，就在中心稀疏波内采用相似的解形 式给出解析解 2 5 ，左、右状态而后确定。而左、右状态的确定又有不同的形 式。 ( 1 ) 简单地以吡，u 。给出，或者以某种平均形式给出。 ( 2 ) 进行单侧的线性化外推，这就是b e na r t z i 的g r p 方法的思想。 ( 3 ) 进行随机选择定解，这就是c h o r i n 的r c m 方法。 ( 4 ) 节点模板选择与插值类方法，无振荡类高分辨率方法，如e n o 和 w e n o 方法。 我们以无粘b u r g e r s 方程为例作简单说明 x o z 0 ( 2 1 6 ) ( 2 1 7 ) 我们可以根据特征线的分布进行解的分析。如果在( x ，f ) 平面上的原点邻域引入 特征线 d x = u d t ( 2 1 8 ) 由于左、右的常数态，那么显然特征线分布可以分为两种情况( 如图2 2 所示) ， 即激波和中心稀疏波。 ( 波结构) 如果 ，间断分解为激波结构。 如图2 2 a 所示，左、右特征线必然相交，即形成激波线 ，i n o i l l r r 、lllrl，f；，l 一2 p 二一、 斗 # “ “ 西北：r 业大学硕 ：毕业论文 ) 图2 2 激波和中心稀疏波 d x = s 出，s = 三( “。+ “。) ( 2 1 9 ) 那么在间断分解点x = 0 的左、右两侧状态将以激波线为界，向上和左右 开拓，从而我们有 工s f 工 s t ( 中心稀疏波) 如果“。 u 。，间断分解为中心稀疏波结构。 ( 2 2 0 ) 如图2 2 b 所示，由于在左、右两侧临界特征线o a ，o b 之间有一个稀疏波 结构才是物理解，也就是唯一的真解，所以在间断解的存在唯一性问题上，将 这种情况确立为中心稀疏波，就一定满足间断解或弱解的存在唯一性的熵条件。 在这种情况下，解有唯一的相似解 l “l “( 州) = 彳 【“o x s u ，r “l 甜a r ( 2 2 1 ) 【 r “ 甜 r，、【 j j ) fx ( “ 两北工业大学硕上毕业论文 第三章高分辨率激波捕捉方法 我们知道，拟线性方程组有间断解，如气体动力学中的激波。即使初始条 件是充分光滑的，也有可能在之后的时间内自发地形成间断解。因此拟线性方 程组的计算都要考虑这个事实。当然，计算间断解的方法很多，在这里我们主 要采用激波捕捉方法。激波捕捉方法的基本思想是不论拟线性方程组的解有无 间断解，不加区别地按统一的格式进行处理。此时，间断是作为解的一部分进 行计算的。间断激波被表示为具有一定过渡层的连续解，这个连续解基本上能 反映出本来的间断特征。这个方法的最大优点是简单。由于计算间断不作特殊 处理，所以程序编制简单，机时也省。其次，守恒型差分格式收敛性已有保证。 当然这个方法也有缺点，主要是有时把间断的过渡区拉地过宽，有时在间断的 前、后会产生不应有的振荡。见图3 1 ，图3 2 。 |辛 i i 8 嗥， o5 一j o f + 1 05 。 j 。 l 。! 卜斌1 百5 。一也 图31i _ - f 格式的计算结果图3 2l 0 格式的计算结果 以上两个例图是根据线性对流方程的r i e m a n n 问题 丝+ 塑：o 西融 1 1 0 x 0 “= l 1 0 x 0 在t = 0 5 时的结果。两个例图分别采用一阶精度的l a x f r i e d r i c h s 格式 西北1 二业大学硕士毕业论文 “? + - 一；( 甜二。+ “三) + 丛二亟：o f2 力 和二阶精度的l a x w e n d r o f f 格式 “? “= “? 一云( “三。一“二) + 嘉( “晶砌? + “：- ) 计算而得。 高分辨率是由h a r t e n 3 5 等人首先提出的概念。所谓高分辨率，就是说在 计算与模拟激波等大梯度、间断问题时，能够得到比较锐利的图像，而不像 “人工粘性方法”那样使激波受到过分的光滑。其实，许多学者构造的方法， 如b o r i s 和b o o k 的f c t 方法，v a nl e e r 的m u s c l 方法，c o l l e l a 等的p p m 方法， 以及其他二阶g o d u n o v 型方法，也属于这种高分辨率方法。本文就是基于 r i e m a n n 问题，构造一种“虚拟网格点”，将高分辨率格式和l e v e l s e t 方法结 合起来，对间断进行捕捉，力求使得间断附近的过渡区域尽可能的小，从而提 高间断处的分辨率，得到相当锐利的图像。 在解决包括间断( 激波，接触间断以及其他间断) 的数值问题中，一种典 型的情况就是双曲守恒律的解。它或是单个标量方程，或是一个双曲守恒方程 组；它可以是一维空间情形 “，+ ( “) ：= 0 ( 3 1 ) 也可以是多维( 例如三维) 空间情形 “。+ 厂( “) ，+ g ( “) ，+ ( “) ：= 0 ( 3 2 ) 激波捕捉方法的一个主要因素是格式的守恒性，即一个逼近( 3 1 ) 的守恒格式为 百d u j + 去( z + 必一乒一) = 。 ( 3 3 ) 其中“，是( 3 1 ) 中精确解在x ，处f 时刻的点值u ( x ，r ) ，或单元平均值 玎( x j , t ) = 石1 睡舳 ( 34 ) 西北工业大学硕卜毕业论文 夕是依赖于相邻点上的数值通量 o + i fl = 厂( “j 一女，“，一k + 【，u j + m ) j + ： 。 ( 3 5 ) 并且夕满足两个条件：数值通量，与物理通量厂 ) 相容；数值通量夕对 ，+ j 于它的所有变量均为l i p s c h i t z 连续。我们注意到( 3 3 ) 是一个半离散格式，然 而在实际计算时，时间变量也必须进行离散。由著名的l a x w e n d r o f f 定理知， 如果一个守恒格式收敛，则必然收敛于( 3 1 ) 的个弱解。故守恒格式( 3 3 ) 是 能够处理具有激波的解的。当然在一定的条件下，对严格守恒放松限制将是非 常有益的，并且不会影响关于弱解的收敛。有兴趣的读者可参阅o s h e r 和 c h a k r a v a r t h y 2 8 的“弱守恒型”格式，以及f e d k i w 的“虚拟流体”方法 2 9 。 第一节一阶单调格式 在上世纪7 0 年代末和8 0 年代初，关于方程( 3 1 ) 和( 3 2 ) 的一阶单调格式 的构造是一个非常活跃的主题。这些格式具有单调的激波跳跃，并且收敛于物 理上的弱解 3 0 。g o d u n o v 格式是一阶单调格式中数值耗散最小的单调格式。 但是，对于复杂的通量函数f ( u ) ，其计算量较大，并且其通量函数不光滑，仅 仅是l i p s c h i t z 连续。l a x f r i e d r i c h s 格式的通量函数f ( u ) 可以很容易地计 算，但数值耗散较大。 在 3 1 和 3 2 中，e n g q u i s t 和o s h e r 对跨音速位势方程和一般的标量守恒 律构造了一种单调格式。这种单调格式相对于g o d u n o v 格式来说，对复杂的通 量函数f ( u ) 容易计算且通量函数c 1 光滑，并且具有几乎和6 0 d u n o v 同样的数值 耗散。 后来，o s h e r 3 3 ， 3 4 1 把这些格式推广到双曲守恒律方程组，得到了o s h e r 格式。 o s h e r 格式比o o d u n o v 格式更光滑并且比最简单的l a x f r i e d r i c h s 格式 西北工业大学顾卜毕业论止 有更小的数值耗散。关于o s h e r 格式的应用可见 3 3 。 第二节高分辨率 r v d 格式 一阶单调格式具有非常好的稳定性并且收敛于丁f 确的满足熵条件的物理 解，但在许多应用过程中，它们耗散太大。因此，我们必须采用更多的网格点 才能得到一个比较合理的解，这就严重地限制了在多维模拟中的应用。 1 9 8 3 年，h a r t e n 在他的具有历史意义的论文中 3 5 ，第一次提出了高分辨 率( h i g hr e s o l u t i o n ) 和总变差不增( t o t a lv a r i a t i o nd i m i n i s h i n g ) 的概念， 这在数值方法的发展历程中具有十分重要的意义，在其后相当长的时期内，成 为人们关注、发展和应用的课题 7 3 6 3 7 。t v d 、e n o 、和w e n o 方法都是从 守恒形式出发，通过构造矢通量函数的数值矢通量的方式，来设计高分辨率格 式。 这些“高分辨率”格式具有至少二阶精度；当激波出现时是稳定的。t v d 格式是一个无伪振荡格式，在解的光滑区域里具有至少二阶的精度，在间断处 有较高的分辨率。 考虑简单的标量守恒律方程( 3 1 ) ，一般显式t v d 格式的形式为 “? “= “? 一尝( 也一，：) ( 3 6 ) 其中，：必为格式的数值矢通量，一般取节点模板形式 ：必2 厂( “，“，。，“；+ t + - ) ， 女，12 0 ( 3 7 ) 要求满足相容性条件 ：魂( “) = ，( ”一，“，“) = ，( “) ( 3 8 ) 通常取节点模板很小，如三节点模板，f _ o ，k = 1 ，即 ：k = 厂，”“) ( 3 9 ) 这时格式为 “? “= “? 一- - 兰 f ( “? ，群二- ) 一厂( “：，h ? ) ( 3 1 0 ) 西北工业大学硕上毕业论文 若定义总变差为 t v ( u ”) = 血h ( 3 1 1 ) 所谓t v d ，即要求 w ( u ) t v ( u ”) ( 3 1 2 ) 人们又常把格式写成变差形式 “? = “? 一c 卜z “卜+ d ，+ 烂“，+ ( 3 1 3 ) 可以证明，格式为t v d 的一个充分条件是 卜o ，d f + 必洲 ( 31 4 ) 0 “，+ + d l + 1 t v d 性质能保证格式是单调的，其解在光滑区域能够至少达到二阶精度， 在间断附近能够抑制过分的振荡和保持比较好的锐利形态。前面我们提到的 l a x f r i e d r i c h s 就是t v d 格式的。 8 0 年代中期，人们认识到，尽管t v d 格式具有非常好的稳定性和高分辨率 特征，但也有不足，这是因为在解的极值处，格式降低为一阶精度 5 。在 2 】 中，h a r t e n 和o s h e r 放松了t v d 限制，采用了u n o 限制，这样构造的u n o 格式在解的极值处也是二阶精度的，因此是一个一致二阶精度格式。然而人们 很快就认识到u n o 的限制条件太强，并且不能构造出高于二阶的一致精度格 式。在2 0 世纪8 0 年代和9 0 年代初，大量的研究论文，各种类型的改进，优化 形式的t v d ，t v b 和大量的应用文章出现在各类数值方法和科学计算的杂志上 【3 6 3 7 1 1 2 5 7 。 第三节e n o 格式 前面我们介绍的高分辨率格式的精度不会很高，其关键在于，被用来构造 格式的节点，或者单元的模板是很小的。因为它们所依赖的用以构建的网格点 个数，或者积分控制元个数为2 至3 个，能够构造的逼近多项式不会超过2 阶。 如果考虑到适当的自适应调节，譬如引入限制函数等等，总体精度也只能低于 - 1 6 一 西北工业大学硕上毕业论文 2 阶。因此，即使我们使出浑身解数，也不可能提高格式的精度。 为了构造高精度格式，唯有扩展节点模板，增加节点的个数，或者说增加 控制元的范围，即单元模板的个数。但是，采用这样的方式扩展模板来构造的 更高阶多项式，常常又会引入新的极值，而加强了数值解的振荡机制，反而得 不偿失。 既可以扩展节点模板，又能够尽可能地减少引入数值解的振荡机制，这是 可能的。在v a nl e e r 构造的m u s c l 格式 2 6 和c o l l e l a 等提出的p p m 格式 2 7 中，我们已经看到他们的有效措施，即通过减小n e w t o n 插值逼近多项式的各阶 差商绝对值的幅度，来达到减弱振荡的效果。2 0 世纪9 0 年代发展起来的e n o 和w e n o 方法，更加圆满地解决了这个困难。 1 9 8 7 年，h a r t e n ，e n g q u is t ，o s h e r 和c h a k r a v a r t h y 1 提出了e n o 格式， 即本质无振荡格式。其主要思想就是根据解的局部光滑性，自适应地选择光滑 度最高的节点模板构造插值多项式。 仍然以简单的标量方程( 3 1 ) 为例，假定己知f 。时间层上有一个阶梯函数 ( 网格积分平均值) 分布，即 “( x ，。) = 玩“，石，巧 0 时，u 2 是蚝的激波状态， 反之亦然。 我们举个例子来说，假定给定初始条件是 f 3 s i n 3 x ， 羔0 舻1 2 2 2 8 ， x 0 则我们可以构造虚拟状态“。，“：为 西北工业大学硕e 毕业论文 “= 3 s i n 3 x “，= 2 x 3 8 l e v e ls e t 函数巾的初始条件可定为 巾o2 一x r 0 工o 图50 a 初始条件 图5o b 虚拟状态 此时，方程( 5 i ) 就成了 ( “1 ) ，+ ，( “1 ) ，= 0 ( “2 ) ，+ f l u 2 ) ，= 0 3 0 ， ( 5 6 ) ( 5 7 ) 西l t i 业大学硕士毕业论文 我们分别求解这两个方程，最后以( 5 4 ) 式确定真解。在此，我们需要注意的 是：我们给定的虚拟状态是否满足熵条件 f 。( “，) 厂+ ( ，) ( 5 8 ) 其中 旷仨三黑 、 驴恢三： 在中 0 ( 由0 ) 的地方，如“：( ) 不满足熵条件，n + u ：= ( 确= 甜：) 。 一旦熵条件满足，我们就可以知道激波的速度： c f ( u ，) 一f ( u f ) l ) 一一 “r 一玑 ( 5 9 ) 从而可以求解l e v e ls e t 函数。 中，+ 舯。= 0 ( 5 1 0 ) 5 3 1 线性对流方程 第三节数值试验 首先我们看一下本章的算法在线性对流方程中的应用。线性对流方程是简 单的双曲守恒律方程，它具有解析解，通过与解析解进行对比，我们可以充分 地检验我们的算法。 甜f + a u ，= 0 ( 5 1 1 ) 其初始条件为“。( x 。) ，那么其解析解为 d ( x ，) = d o ( x o ) = “o ( x a t ) 本节中我们考虑的是具有周期性边界条件的线性对流方程 “，+ ”，= 0 ( 5 1 2 ) 两北工业大学顾士毕业论文 初始条件为 2 ( x + 1 ) 一言s i n ( 3 石( x + 三) )一1 x 一圭 七一扣n 互3 砸一尹12 ) 一j 1 z 丢 2 石( ： s i n ( 2 7 r ( x z ( x 一一 一x ) ) 1 i i 一= 川 土s i n r 3 万h 一 6 、 、 = 1 x i 1 ( 5 1 3 ) oz 1 5 一 工 一 26 二 z 1 6 显然此方程具有周期性边界条件“( 一1 ) = “( 1 ) ，初始条件在上= 一i 1 ，x = 吉，x 处不连续，见图5 l - 口i 06 7 吐 图5 1 初始条件 通过l e v e ls e t 函数，函数u 可以表示为 2 ( x + 1 ) 一吉s i n ( 3 万( j + 兰) ) s i n ( 2 石( = 1 一x ) ) 2 ( z 一1 ) 一i 1s i n ( 3 石( z 一言) ) 一3 2 1 x 一二 6 一土 x 三 ( 5 1 4 ) 63 = x 1 两北工业人学硕士毕业论丈 “2 = 其中l e v e ls e t 函数中为 0 - o2 _ o4 o = s i n ( 2 厅( x 一言) )一1 x 一詈 七一争s i n ( i 3 m 一扣一i 5 x ；( 5 1 5 ) 如( 2 砸一扣； 川 x + 一 o l x 一一 l x 一一 6 1 一x + 一 x 一一 6 一l x 一一 5l 一 x 一一 66 11 一一 x 一 63 12 一 z 一一 33 = x 1 3 | i n山 ： j p j 专 争 。1a - i r 图52 aw e n 0 5n x = 5 0t i m e o u t = 20 ( 5 t 6 ) 图5 2 、5 3 是我们在，= 2 0 时的解。实线代表解析解，圆圈是数值解。 图5 2 a 、图5 2 b 、图5 2 c 是我们用5 0 、1 0 0 、1 0 0 0 个网格点，用五阶w e n o 格式得到的各种图像。我们可以看出随着网格的加密，数值解越来越接迓解析 - 3 3 i1叫f叶 两北l 业大学硕士毕业论文 解。但即使如此，在间断处的数值解仍发生很强的数值耗散。当我们采用l e v e l s e t 方法对间断追踪后，我们发现采用5 0 个网格点得到的数值解就比1 0 0 0 个 网格点的五阶w e n o 格式