参考文献CFD-trainingm-UCD5-CN

1 保单调迎风紧致差分格式保单调迎风紧致差分格式 何志伟 1 李新亮1* 傅德薰2 马延文2 1 中国科学院力学研究所 LHD 实验室 2 中国科学院力学研究所 LNM 实验室 * 联系人，Email: lixl 摘要：摘要： 针对迎风紧致格式，结合保单调思想，构造了一种 5 阶精度的非线性保单 调迎风紧致格式。该差分格式保留了原迎风紧致格式波数分辨率高的优 点，并大幅提高了激波分辨能力，有效克服了过激波产生的非物理振荡。 与经典的 5 阶 WENO 格式相比，新的格式的耗散小且计算量低。并且本文 通过逐元素以及特征分解两种形式，将该方法推广应用到求解 NS/Euler 方程组。通过 Sod 激波管问题，ShuOsher 激波密度波干扰问题及二维双 Mach 反射问题，展示了该方法良好的波数分辨率及激波捕捉能力，可用 来模拟含激波的复杂流动。 关键字：迎风紧致格式； 保单调； 可压缩流； 激波捕捉格式 1 引言引言 航空航天领域的应用为可压缩复杂流动数值模拟技术提出了强劲需求。 而这 些流动中包含了湍流、激波、分离流等多种复杂因素，流动具有多尺度、强间断 等复杂特征，各类波系在其中充当着不同的角色。为了真实刻画流动细节，人们 纷纷要求采用高精度格式。因为高精度格式的耗散误差和色散误差更小，在同样 条件下能分辨出更加精细的流场，捕捉到更细微的结构。 一个好的格式应该能得到原物理问题锐利的、逼真的物理图像，它习惯地被 称为具有高分辨率（high resolution）3。根据上述定义可知，一个高分辨率的格 式必须具有如下的特征4： （1）在光滑区至少有 2 阶及以上的精度； （2）格式不 会有或者有小的数值振荡； （3）格式对间断要有高的分辨率，穿越间断所需的网 格点数与一阶精度的单调格式相比要小的多。所以，高（阶）精度未必意味着高 分辨率。 Harten3在 1983 年首先提出了总变差不增（total variation diminishing, TVD）的概念，并通过对一阶精度的迎风格式的数值通量进行反扩散 （antidiffusion） 处理 （本质上亦为流通量矫正迁移 （flux correction transport, FCT） 方法）使的其过大的数值耗散反掉，从而得到了一个 2 阶精度的 TVD 格式。在 Harten 给出的 TVD 概念及其充分条件3以及运用限制器（limiter）技术的思想3 基础上，Sweby5 提出了采用流通量限制器（flux limiter）的方法来构造一系列 的 2 阶 TVD 格式。此外，Van Leer6提出构造具有 TVD 性质格式的 MUSCL 2 （monotonic upstreamcentered scheme for conservation law）方法，即首先将网 格节点上的物理量外插至网格单元边界面处，并进行适当的限制，最后再在边界 面处采用流通量矢量分裂（FVS）方法，这样就可以获得具有良好的色散和耗散 特性的一种激波分辨率较高的格式。从此，即 20 世纪七八十年代之后，国内外 研究者提出了各种具有 TVD 性质的格式，如 NND7。对于激波捕捉，TVD 格式已 经被证实是可靠的和有效的。但是，为了格式满足总变差不增的数学性质，所有 的 TVD 格式在极值点附近的精度退化为 1 阶。为了解决这个问题，人们提出了 各种方法。如，Harten8提出的本质无波动（essentially nonoscillatory, ENO）的 概念，其放宽了 TVD 中强调消除波动的限制，允许格式产生高阶虚假波动。从 此，人们构造了各种具体形式的高阶 ENO8,9格式及 WENO10格式。此外,常见还 有总变差有界（TVB）11，保单调（MP）12等方法。 作为高精度格式的一种，紧致格式由于使用了网格上的导数信息，因而采用 相同的网格基架点构造出的紧致格式比传统格式有更高阶的精度、 更高的尺度分 辨率以及更小的色散/耗散误差。 紧致格式可以分为两大类： 对称型紧致格式 （如 1） ，及迎风型紧致格式（如2） 。紧致格式同样面临过间断的非物理振荡，且 精度越高，振荡越难以克服。Halt13认为紧致格式的进一步发展在于对无波动激 波捕捉技术的研究。为了能使紧致格式较好的捕捉激波，需要对紧致格式的非线 性处理。目前，在这方面的工作较少，主要有以下研究： 1）Cockburn14通过 TVD 与总变差限制 TVB 概念建立了 3 阶和 4 阶精度的非 线性紧致格式。但是其 4 阶格式在激波附近会产生明显的数值振荡。 2) Ravichandran15将 minmod 通量限制器应用到高阶紧致格式，构造对于一 维标量方程具有 TVD 性质的高精度紧致格式。但是，这种方式造成在极值处降 为 1 阶。 3）Adams 和 Shariff 16、Pirozzoli17、Ren 等18分别构造了高精度紧致格式和 ENO、WENO 的杂交格式来数值模拟含激波的流动。这种格式依赖于开关函数的 设计，而开关函数里往往含有经验参数。 3）Deng 等 19针对一类特殊的单元中心型紧致格式，采用自适应非线性差 值技术，构造了一类非线性紧致格式。但是该格式的计算量较大。 4）Ma 与 Fu20采用群速度控制方法构造了一种群速度控制的紧致格式。该 格式计算量很小但鲁棒性有待提高。 本文在 Fu 与 Ma21的迎风紧致格式基础上，结合保单调思想12，构造了一 种 5 阶精度的保单调迎风紧致格式， 并通过算例进行了验证比较。 本文结构如下： 第二节我们将具体地介绍 5 阶保单调迎风紧致格式的构造方法。在第三节中，通 过采用逐元素和特征分解两种形式， 我们将新格式运用到求解 N-S/Euler 方程上。 第五节分析了一些经典的测试算例和数值结果。 第四节和第六节则分别讨论了边 界格式的处理和结论。 3 2 基本原理基本原理 2.1 基本方程基本方程 我们考虑如下的标量守恒律方程： ( ) 0 uf u tx (2.1) 离散点 j x为均匀网格点，间距为 h. 方程(2.1)的半离散的守恒型差分格式可以写为： 1/21/2 1 ()0 jj u ff th (2.2) 其中 1/2 j f 为数值通量。如果 1/21/2 1 ()() k jj j f ffO h hx (2.3) 那么这个格式是k阶精度。在本文中，时间方向我们使用 3 阶 TVD RungeKutta 法22。令 1/21/2 1 ( )() jjj L uff h (2.4) 则 3 阶 TVD RungeKutta 法为： (1) (2)(1)(1) 1(2)(2) () 311 () 444 122 () 333 nn jjj n jjjj nn jjjj uut Lu uuut Lu uuut Lu (2.5) 2.2. Fu和和Ma的五阶精度迎风紧致格式的五阶精度迎风紧致格式 紧致格式一般可以在较小的基架上获得相对较高的精度。 Fu 和 Ma21构造了 一种 5 阶迎风紧致有限差分格式(UCD5)。这类格式的一般构造过程可以总结为： 首先进行通量分裂 fff ，其中0 df du 、0 df du 。然后对 f 与 f ，分别 构造如下形式的差分表达式 2112 1 3443612 23 5560 jjjjj jj fffff FF h (2.6a) 2112 1 3443612 23 5560 jjjjj jj fffff FF h (2.6b) 4 其中 j j f F x 。 下面以如下线性单波方程的初值问题为模型，讨论数值方法的波数分辨率。 00 uu aaconst tx 初始条件：( ,0) ikx u xe 该方程的精确解为： () ( , ) ik x at u x te 考察 1 阶迎风格式、5 阶迎风偏斜格式和 5 阶迎风紧致格式的色散和耗散误差。 如果记 j F为逼近 j u x 的各种差分表达式，则可得出 j ikx je Fke，其中 eri kkik。则我们可以求出0 j j du aF dt 的数值解为 1 () ( , ) ri at kik x akt hh u x tee 。 不同的格式的色散误差和耗散误差可以通过 i k和 r k来反映23。 对于精确解有 0 r k ， r kkh。1 阶迎风格式、5 阶迎风偏斜格式和 5 阶迎风紧致格式的色 散和耗散误差曲线如 Fig.1 所示。可以看出：5 阶迎风紧致差分格式(UCD5)比相 应的普通 5 阶迎风偏斜差分格式(UD5)相比，具有更高的分辨力和更小的耗散。 (a) (b) (a)修正波数的虚部（色散） (b)修正波数的实部（耗散） Fig. 1. 各种差分格式的修正波数 表 1 则给出了 1 阶迎风格式、5 阶迎风偏斜格式和 5 阶应风紧致格式的 i k和 5 r k与精确解不超过 5%和 2%时所对应的最大值()kh。可以看出：5 阶迎风 紧致格式能够模拟的波数范围要明显地高于普通 5 阶迎风偏斜有限差分格式。 5% r k 2% r k (1/)5% i k(1/)5% i k 1 阶迎风阶迎风 0.318 0.201 0.552 0.348 5 阶迎风偏斜阶迎风偏斜 1.29 1.08 1.49 1.25 5 阶迎风紧致阶迎风紧致 1.57 1.35 2.01 1.71 表 1 i k和 r k与精确解不超过 5%和 2%时所对应的最大值 不仅如此，该类迎风紧致格式还有一个特别之处，即可以将其数值通量写成 “拟显式”的形式。即，计算出上游点上的通量后，通过如下递推关系就可以得 到本网格点上的通量。 112 11 22 34711 52 3605 jjjj jj ffff ff (2.7a) 21 13 22 34711 52 3605 jjjj jj ffff ff (2.7b) 而通常的对称紧致格式需要求解多对角方程组获得数值通量， 虽仍可以利用追赶 法获得递推关系式，但通常为双向的递推关系（ “追”过程和“赶”过程） ，求解 效率，尤其是并行效率不如 Fu 和 Ma 的迎风紧致格式。 Fu 和 Ma21迎风紧致格式的单向递推特点为采用限制器提供了便利。 2.3 保单调限制保单调限制 针对线性方程，文献12提出了一种“1 阶精度保单调”思路。文献24证明 了此方法与 TVD 的等价性。不同于12，24，本文将该方法进行了推广及修正， 给出了通用的保单调限制方法。并且将上述方法运用到有限差分形式的 Fu 和 Ma21迎风紧致格式中。 假设需求解的方程为（2.1） ，且 0 df aconst du 即讨论正通量 1 2 j f 的情况。对于0a （即 1 2 j f ）的格式，依据格式的对称性容 易获得。 传统的 2 阶 TVD 格式可以看作是对 1 阶迎风的修正。所以在迎风基架 6 211 , jjjj xxxx 上，我们可以写出如下的数值通量的一般表达式： 1111 222 1 (, ) () jjj jjj T ffrff (2.8) 其中a h 是 CFL 数， 1 1 12 nn upwjj nn j locjj uu r uu . 其中的 T 是传统的 TVD 限制器的上限，即： 1 2 0 BT j (2.9) 利用 Harten 定理3，得出一步格式 TVD 的一般限制条件为: 11 122 2 11 0()1 jj TT j r (2.10) 则有： 1 2 0 BT j (2.11a) 11 22 1 0 ST jj r v (2.11b) 上述两个式子 （2.11a、 2.11b） 决定了一个区域， 即为 TVD 区域， 见 Fig.2.(a)。 定义 11 22 1 (, ) jj T r (2.12a) 11 22 1 (, ) jj S r (2.12b) (2.11a)和(2.11b)可以化简为： 01 (2.13a) 01 (2.13b) 我们一般地重新定义： 11 2 (, ) () ul jjjj j ffrff (2.14) 其中 1 12 2 (, ) S j T j r r 7 则(2.8)经过运算可以变成为： 11 2 () jjj j ffff (2.15a) 1 2 () ul jjj j ffff (2.15b) 利用(2.13a)和(2.13b)，有： 11 2 , jj j fff (2.16a) 1 2 , ul j j fff (2.16b) 其中 12 ,., k z zz表示区间 1212 min( ,.,),max( ,.,) kk z zzz zz。 也就是说，如果数值通量在上述两个区间的交集内，则这个格式是 TVD 的。而 文献12中的4正好即为 1 2 (,0.2) j r 时的取值。 (a) (b) Fig. 2 (a) TVD 区域 (b)TVD 区域决定的开关函数 最后，我们分析保单调的思路，并对其进行改进。 (1) 保单调的过程，精度不是考虑的因素。其本质是在单调区利用 TVD 性质来提 供一个区间，即 1 2 j f 必须在某个区间内。而这个区间的上界/下界是： 1j f 和 ul f，下界/上界是： j f。下面我们分别讨论这两个界限值的起到的作用。 上界/下界： 23 1 2 1111 1(1)(1)(1)() 8 ul j TTT ffh fhfO h 8 23 1111 222 11 () 28 j jjj ffh fhfO h 对应的下界/上界是： 23 111 222 11 () 28 j jjj ffh fhfO h 所以可以得到满足 TVD 性质的区间：, mpmp ab ff ，而文献12取得是： mp aj ff 1 minmod(,) mpul bjjjj ffffff 进一步我们可以用 11 22 ()() mpmp ab jj ffff 来作为一个开关函数来区分单调区域和极值与间断区域， 见 Fig.2(b)。 所以在 单调区可以用各类高（阶）精度格式，而其他区域则可以用 ENO/WENO 这 种处理极值和间断较好的格式。通过这种杂交技术，可以显著地降低 ENO/WENO 格式的计算量。 (2) 我们新定义其中的不再是常数，而更一般地定义为( , )r的函数。那么 各种传统的 TVD 限制器均可以求出对应唯一的( , )r函数。而那么上述 过程的只相当于用了其中一种特殊情况ULTRABEE 限制器。Roe25给出了 一种针对半离散形式的 TVD 条件：( )max(0,2 )rr。显而易见，当0r 和 r 时应该有同样的处理机制，即对称性。所以我们根据对称性，可知， ( )min max(0,2 ),2rr， 这样， 我们的得到了一种半离散的 1 阶保单调方法， 而时间方向没有固定是 RK 方法。对于这种情况，那么根据本文的定义和讨 论，上述函数简单的取常数 1。数值实验表明，本文的修正运用到上述 5 阶迎风紧致格式（得到的格式记为 mUCD5）略微增加了耗散，但是对应的 逐元素形式和和特征分解形式的差分格式在0.2CFL 下都可以获得较好的 分辨率和鲁棒性。而采用原来的方法24运用到上述 5 阶迎风紧致格式（得 到的格式记为 MPUCD5）则不适合逐元素形式的有限差分格式，会产生剧烈 的数值振荡。而且即使运用到特征分解形式的有限差分格式上时，仍然会产 生较为明显的数值振荡。 2.4.精度保持精度保持 9 同传统的 TVD 类型格式，上述方法在极值处降为 1 阶精度。为了解决此问 题，文献12提出了一种方法：扩大上述限制区间的范围，使得在极值处原格式 算出的界面值不受影响。从而达到精度保持。文献24讨论并修正了这种方法。 文献12最先根据对抛物插值的几何分析，给出了如下的思路： 将上述两个区间分别扩大为： 1 , md jj fff 、 , ullc j fff 。 其中： 11 2 11 () 22 mdMM jj j fffd (2.17a) 11 2 18 () 23 lcMM jjj j ffffd (2.17b) 文献24修正了上述两个参量的定义： 11 2 11 () 22 mdMM jj j fffd (2.18a) 1 2 1 1 (1) () 2 MM j lcul jj jj d ffff ff (2.18b) 并在数学上证明了上述两个修正后的参量具有如下的性质：在单调区，这两个用 来扩大的区间的两个参量 md f和 lc f满足： 1 , md jj fff (2.19a) , lcul j fff (2.19b) 这里，根据我们的一般性地讨论，本文新定义如下的两个参量： 11 2 11 () 22 mdMM jj j fffd (2.20a) 11 22 11 ()(, ) 22 lculMM jj jj fffrd (2.20b) 其中： 11 2 minmod(,) MM jj j ddd ， 11 2 jjjj dfff 。在实际的应用过程中，可 以采用文献12一种更为苛刻的 4 1 2 M j d ，其定义为： 4 1111 2 minmod(4,4,) M jjjjjj j ddddddd 新定义的(2.19b)中 1 2 (, ) j r 函数为： 10 ( , )0 ( , ) 0 rr r r (2.21) 上述新的定义：在极值处，参数是一个与极值有关的经验参数。但对于极值情 况各异，难以一般的统一地描述；在单调区，不难证明：本文新定义的两个参量 md f和 lc f仍满足上述性质。 根据我们的一般性的讨论和上述对半离散形式进行的相应地修订， 本文新定 义的 lc f最终简化为： 1 2 11 () 22 lculMM jj j fffd (2.20b) 从而上述一阶精度保单调算法中的两个区间分别扩大为： 11 2 , md jj j ffff (2.22a) 1 2 , ullc j j ffff (2.22b) 所以上述扩大区间的过程在单调区间其实并没有改变上述保单调区间。 下面我们讨论一下，上述扩大区间的过程在极值点处的精度问题。文献12 利用几何方法针对抛物插值曲线论述了在极值处可以实现精度的保持， 但对于一 般的极值情况迥异，难以统一地描述。所以很难给出一般性的数学证明。但数值 试验25表明上述方法在极值处的精度虽有所下降，但总体而言，是可行的。 最终，利用本方法计算通量最终，利用本方法计算通量 1/2 j f 的具体算法为：的具体算法为： （1） 依据(2.7a)算出数值通量 1/2 j f （2） 如果 11 22 ()() mp jb jj ffff 成立，则计算结束， 1/2 j f 为最终的 数值通量；否则进入第三步 （3） 计算(2.14), (2.18a), (2.20b), 然后计算： max 1 max(min(,),min(,) mdullc jjjjjjj fffffff min 1 min(max(,),max(,) mdullc jjjjjjj fffffff minmax 1111 2222 minmod(,) jj jjjj ffffff minmod 函数采用文献12的定义，即： 11 123123 minmod( ,)min(,) kk z zzxszzzx 而其中的s为： 12131 111 ( )()( )().( )() 222 k ssign zsign zsign zsign zsign zsign z 3 推广到推广到 Euler 方程组方程组 3.1 Euler 方程组及其半离散格式方程组及其半离散格式 一维 EULER 方程组为： ( ) 0 tx UF U (3.1) 其中( ,)TuE U是守恒量矢量， 2 ( )(,() )TuupEp uF U是相应的通 量矢量。而、u、E、p则分别表示密度、速度、单位体积的能量和压力。 方程(3.1)的半离散的守恒型差分格式可以写为： 1/21/2 1 ()0 jj th U FF (3.2) 其中 1/2 j F为数值通量。 如果 1/21/2 1 ()() k jj j O h hx F FF (3.3) 那么这个格式是k阶精度。 3.2 推广到推广到 Euler 方程组方程组 本文使用 Steger-Warming 分裂。先将流通量分为正负流通量 FFF。接 下来有两种方式可以将上述格式应用到求解 Euler 方程组上来。 我们分别讨论之： 第一种方法：对于正负流通量的每一个分量分别直接使用上述格式来计算。 此方法得到的结果即为逐元素形式的五阶非线性迎风紧致格式。记为 （m-UCD5,component-wise） 第二种方法：即采用特征分解形式。这种方法不但符合从标量方程推广到方 程组的思路，而且往往会获得比逐元素形式更好的计算结果。但是由于我们用的 是紧致格式，紧致格式严格说来是依赖于所有离散点的。所以与传统格式的特征 分解形式有所不同。下面给出特征分解形式的五阶非线性迎风紧致格式 （m-UCD5,characteristic-wise）的具体算法: 12 (1). 在每一个固定点 1/2j x ，计算平均状态 1/2j U：可以用 j 点及 j+1 点值的 算术平均或者 Roe 平均。 (2). 计算相应的特征值 ( ) 1/2( 1,2,3) i j i 、 左特征向量 ( ) 1/2( 1,2,3) i j li 和右特征向 量 ( ) 1/2( 1,2,3) i j ri 。 (3). 流通量函数进行局部特征分解： ( )( ) 1/2 ,1,2,3;1,.,2 ii mjm wlimjj F ( )( ) 1/2 ,1,2,3;2,.,1 ii mjm wlimjj F 此外，不同于传统的差分格式，五阶非线性迎风紧致格式的迎风点处数 值通量需重新进行局部特征分解：对于 1 2 j F和 3 2 j F在1/ 2j处重新进行局部 特征分解，得 ( ) 1 2 i j w 和 ( ) 3 2 i j w 。 (4).在局部特征变量 ( )( 1,2,3) i m wi 上分别计算： ( )( )( )( ) 112( )( ) 11 22 34711 52 3605 iiii jjjjii jj wwww ww ( )( )( )( ) 211( )( ) 13 22 34711 52 3605 iiii jjjjii jj wwww ww (5). 进行限制过程。 (6). 变换回物理空间计算出 1/21/21/2 jjj FFF。 4 边界格式 4 边界格式 对于边界格式，我们采用文献27里的处理方式。即使用六点模板 21123 (,) jjjjjj xxxxxx 来近似计算 j F ，其中 j j f F x 。 211235 126512060203 () 60 jjjjjj j ffffff FO h h 可以得出： 121012 2 1 ( 32747132) 60 ffffff 同理可得： 13 12112 2 1 (21347273) 60 NNNNN N ffffff 通过采用虚拟网格即可获得上述的数值通量。 5 数值结果数值结果 本文中，我们将使用上述两种形式的 5 阶保单调迎风紧致差分格式计算一 维、二维算例。并且和原来的 5 阶迎风紧致差分格式及经典的 WENO5 格式进行 了比较。 5.1 算例算例 1: 一维单波方程的初始间断问题一维单波方程的初始间断问题 为了更好的与 WENO5 比较间断的分辨率和格式的计算量，我们考察一维单 波方程的初始间断问题。问题具体如下： 控制方程：0 uu tx 初始条件为：在区间 1,1上，初始条件为： 2 0 0 0 0 ( )exp( log(2)(0.7) /0.0009)0.80.6 ( )10.40.2 ( )110(0.1)00.2 ( )1 1 uxxifx uxifx uxxifx ux 2 1/2 0 00(0.5) 00.2 ( )0 xifx uxotherwise 空间离散分别采用本文提出的保单调 5 阶迎风紧致格式（下文记为 mUCD5） 及 5 解精度的 WENO 格式（下文记为 WENO5） ，时间推进采用 3 阶精度 TVD 型 RungeKutta 方法。空间网格为 200 个。 Fig.3. mUCD5 和 WENO5 的计算结果 14 表 2 mUCD5 和 WENO5 求解算例 1 所需的 CPU 时间（推进 2000 时间步） 方法 CPU 时间 （毫秒） mUCD5221.755 WENO5 353.396 图 3 给出了算例 1 的数值解及精确解。 从图 3 可以明显地看出，5 阶保单 调迎风紧致格式 （mUCD5） 的计算效果要好于同阶的WENO差分格式 （WENO5） 。 mUCD5 对于方波和尖波的分辨率更为出色，其耗散明显小于 WENO5。 表 2 显示了采用两种不同方法求解算例 1 所需的时间。 该测试采用的网格点 为 200 个，推进的时间步为 2000。 可以看出，mUCD5 差分格式的计算时间则 要小于同阶的 WENO 差分格式。采用 mUCD5 所需的 CPU 时间仅为采用 WENO5 的 62%。 5.2 算例算例 2: Sod 激波管问题激波管问题 这是一个典型的算例。其控制方程为一维的 Euler 方程组，初始条件为： (1,0,1)0.5 ( , , ) (0.125,0,0.1)0.5 if x u p if x 最终的结果是在计算区间0,1x上计算结束时间为0.14t 。 时间步由以下 公式求出： max () jj j h tCFL ua 。 其中的 CFL 取 0.2。 15 Fig.4. Sod 问题，网格点 N=100，0.14t 时刻密度、速度、压力分布。逐元素形式。 图4显示逐元素形式的5阶保单调迎风紧致格式(mUCD5)与原来的5阶迎风 紧致格式(UCD5)的数值结果。从图 4 中可以看出，UCD5 的数值结果产生了明显 而又剧烈的数值振荡，而 mUCD5 的数值结果则避免了绝大部分的数值振荡，取 得了很好的数值结果。 图 57 分别比较了逐元素形式的 5 阶保单调迎风紧致格式 (mUCD5)与经典 5 阶 WENO 格式(WENO5) 计算的密度、速度和压力分布。其中 图 57(a)为全场分布，图 57(b，c)为局部放大图。从图 47 可以得出结论，新的 格式（mUCD5）很好地解决了原格式（UCD5）的振荡问题。而且 mUCD5 在计 算量小于 WENO5 的情况下，取得了与 WENO5 相当甚至更好的数值结果。 (a) (b) (c) Fig.5. Sod 问题，N=100，0.14t 时刻密度分布。逐元素形式。 16 (a) (b) (c) Fig. 6. Sod 问题，N=100，0.14t 时刻速度分布。逐元素形式 (a) 17 (b) (c) Fig. 7. Sod 问题，N=100，0.14t 时刻压力分布。逐元素形式 最后， 我们给出计算量略大但是计算效果更好的特征分解形式的 5 阶保单调 迎风紧致格式(mUCD5)的数值结果，见图 8。可以看出，基于局部特征分解理论 的特征分解形式的 mUCD5 格式确实可以有效地抑制数值振荡。 Fig.8. Sod 问题，网格点 N=100，0.14t 时刻密度、速度、压力分布。特征分解形式。 5.3 算例算例 3：ShuOsher 激波激波密度波干扰问题密度波干扰问题 这个问题的控制方程为一维的 EULER 方程组，初始条件为： (3.857143,2.629369,10.333333)4 ( , , ) (1 0.2sin(55),0,1)4 if x u p xif x 最终的结果是在计算区间 5,5x 上计算结束时间为1.8t 。 时间步由以下 公式求出： 18 max () jj j h tCFL ua 。 其中的 CFL 取 0.2。 原来的 5 阶迎风紧致格式(UCD5)由于剧烈的数值振荡造成计算发散而无法计 算这个问题。引用保单调限制后，新格式（mUCD5）可很好地计算该问题。图 Fig.9(a)给出了逐元素形式的 5 阶保单调迎风紧致格式(mUCD5)与 WENO5 的计算 结果的比较。其中图 9(b)为 Fig.9(a)的局部放大。从图图 9 可以看出，经过修正 后的逐元素形式的 5 阶保单调迎风紧致格式(mUCD5)不但可以算， 而且还取得了 与 WENO5 相当甚至是好于 WENO5 的数值结果。 (a) (b) Fig.9. ShuOsher 问题，1.8t 时刻的密度分布。N=200。其中的精确解是由 N=4000 时 WENO5 算出。逐元素形式。 (a) (b) Fig.10. ShuOsher 问题，1.8t 时刻的密度分布。N=200。逐元素形式。其中的精确解 是由 N=4000 时 WENO5 算出。 图 10 给出了网格点 N=400 时逐元素形式的 5 阶保单调迎风紧致格式 (mUCD5)与 WENO5 的计算结果的比较。 19 (a) (b) Fig.11. ShuOsher 问题，1.8t 时刻的密度分布。(a)N=200， (b) N=400。特征分解形 式。其中的精确解是由 N=4000 时 WENO5 算出。 图 11 给出计算量略大但是计算效果会更好的特征分解形式的 5 阶保单调迎 风紧致格式(mUCD5)的数值结果， 可以看出本方法具有较高的分辨率和激波捕捉 能力。 原格式（UCD5）无法计算这个问题。新的格式不但可以算，而且从上述的 数值结果中可以看出新的格式的计算效果要好于 WENO5。 5.4 算例算例 4：双马赫反射问题：双马赫反射问题 这是一个典型的测试格式的二维算例，来自于文献27。其控制方程为二维 的 EULER 方程组。 计算区域为0,4 0,1。 一平板位于计算区域的底部6/1x处。 初始时刻有一道位于平板左侧、与x轴夹角为60o、马赫数10M的激波。计算 时间为0.2t 。示意图如下： Fig.12 双马赫反射问题示意图 时间步由以下公式求出： xy xy tt tCFL tt 其中 20 , , max() x j kj k j k x t ua CFL 取 0.2。 原格式（UCD5）计算该问题会造成计算发散（出现负的压力及温度） ，从而 无法计算；而进行本文构造的 mUCD5 可以很好地计算该问题。 (a) mUCD5 (b) WENO5 ( c) mUCD5 局部放大 (d) WENO5 局部放大 Fig.13. 双马赫反射问题，0.2t 时刻的密度分布（值从 1.731 到 20.92 共 30 条等间距 分布的等值线） 。网格点 960240。逐元素形式。 图 13 给出了 mUCD5 与 WENO5 的比较结果。 从图中可以看出 mUCD5 对于 21 激波与接触间断具有更高的分辨率。且数值耗散小于同阶的 WENO 格式。 逐元素形式 mUCD5 的的缺点就是流场中含有数值振荡， 采用特征分解形式 可有效解决该问题。 图 14 给出计算量略大但是计算效果更好的特征分解形