高等计算流体力学讲义(6)

高等计算流体力学讲义（6）§9ENO格式一、ENO方法概述二阶或二阶以上精度的格式可能在间断附近产生振荡（Gibbs现象），这些振荡并不会随着网格间距的变小而减小，并经常导致包含间断的非线性问题的数值不稳定。早期消除或减小这种间断附近振荡的尝试基本上是基于两种方法：人工粘性和限制器。前一种方法是增加人工粘性项，通过调节使其在间断附近足够大以便消除或是减小这种振荡，而在光滑的地方则很小以便保持较高精度。这种方法的一个缺点是调节过程中包含自由参数，计算结果的好坏与自由参数的取值有关，而自由参数的“最佳值”往往是与具体问题有关的。后一种方法就是利用限制器来减小振荡。从效果上看，它降低了在间断附近的精度。使用经过精心设计的限制器，一维非线性问题或者线性系统可以获得TVD性质，这类格式对于多维线性问题在最大值范数（范数）的意义下是稳定的。而且，在限制器中一般没有多余的参数供调节。这种方法的缺点是在光滑的极值点附近会降为一阶精度。ENO（EssentiallyNon-Oscillatory）格式是1987年由Harten，Engquist，Osher和Chakravarthy等人最先提出的。他们提出的ENO思想似乎是第一个获得具有自相似性（与网格尺度无关）、一致高阶精度而且本质上无振荡的分片光滑函数插值方法的尝试。ENO格式特别适合于既包含复杂激波又包含光滑区域的流场。自从Harten等人发表他们的第一篇关于ENO的文章以来，他们和许多其他研究者们继续在研究ENO的道路上前进，改进算法和扩大该方法的应用领域。Shu和Osher发展了基于点的ENO重构，时间上用TVDRunge-Kutta离散的方法，对于多维问题该方法可以大大节省计算成本；Liu，Osher和Jiang，Shu等人则发展了一维和高维的加权ENO（WENO，WeightedENO）格式，可以提高计算精度。在应用方面，ENO和WENO方法已经成功的用于模拟激波和湍流的相互作用以及直接模拟可压缩湍流运动等方面。人们已经看到，ENO和WENO已经相当广泛的应用于不同的领域。用ENO和WENO解决的问题的共同特征是含有很强的激波和复杂的光滑区域，低阶精度格式在这方面有着难以解决的困难，而像ENO和WENO这样的高精度格式就显示出了它们的优势。二、一维重构基于有限体积的ENO格式主要解决数据重构的问题，通过本质上无波动的高阶重构方法，构造高阶的有限体积格式。通常，一维重构问题的提法是：给定网格13112222...NNaxxxxb(1)我们定义单元，单元中心和单元大小为1111112222221[,](),1,2,...,2iiiiiiiiiIxxxxxxxxiN，，(2)并且把最大单元尺寸表示为1maxiiNxx(3)给定函数()vx在单元内的平均值12121()1,2,...,iixixivvdiNx，(4)要求在单元iI内求得一个至多k-1次的多项式()ipx，在iI内逼近函数()vx到k阶精度，即()()()1,2,...,kiipxvxOxxIiN，，(5)特别的，这种重构给出了单元边界上的()vx函数值的近似11112222()()1,2,...,iiiiiivpxvpxiN，，(6)它们都是k阶精度的11112222()(()(1,2,...,kkiiiivvxOxvvxOxiN)，)，(7)这里的多项式()ipx也可以换成其他的简单函数，例如三角函数。单元边界上的近似函数值作为Riemann问题或通量分裂方法的左右状态，可以构造高阶数值通量，从而得到高阶的有限体积格式。下面讨论这个重构问题的解法。给定单元位置iI和精度k，我们先选取一族模板点，从iI左边选取r个单元，右边选取s个单元，连同iI一起有100rskrs，，，即(){,...,}irisSiII(8)为了方便，我们称()Si为一个模板点集。可以证明，存在一个唯一的至多1krs次多项式()ipx（后面在不致引起误会的地方就省去下标i，记为()px）。它在()Si的每个单元内的平均值都等于()vx在相应单元内的平均值，即12121(),...,，jjxjxjpdvjirisx(9)[我们令1011()kkpxaaxax，则有1krs个待定系数，而(9)式正好提供了1krs个方程，所以通过(9)式可以唯一确定待定系数011,,,kaaa。但是，实际应用中，我们有更直接的方法确定这些系数，这个方法将在后面介绍]容易证明，只要()vx在()Si中足够光滑，这个()px就正是我们要寻找的k阶精度近似。对于()px为多项式的情况，(7)式中从()Si的平均值jv到12iv和12iv的映射是线性的，所以存在与r，精度k和单元大小jx有关但与函数()vx无关的常数rjc和rjc，使得11110022kkirjrjirjrjiijjvcvvcv，，(10)系数之间满足下述关系：1,rjrjcc。为了证明这个关系，我们注意到:1122()iiivpx，11122()iiivpx，1,2,...,iN。如果1()()iipxpx，重构的模板点集不同,则一般1122iivv。但是,对于确定的模板点集()Si，()px是唯一确定的，而()Si可由i,r,k完全确定。所以，(10)式实际上隐含了插值的模板点集()Si。注意到如果(){,...,}irisSiII，(1)(1)11(1){,...,}irisSiII，则()(1)SiSi，1()()iipxpx。此时111/2021111/21,(1)(1)02()()krjirjiiijkiirjirjijvpxcvvpxcv。所以1,rjrjcc。既然(10)式的右端项隐含了插值的模板点集()Si，所以我们可以简单的记i+1/2处的函数值为1102kirjrjijvcv，而不必区分从左侧还是从右侧的插值多项式出发进行计算。在后面的WENO格式中，我们将用这种写法。rjc和rjc的值与()Si有关。对于均匀网格情形，k=2时，()Si有以下两种可能的取法：1(){,}iiSiII和1(){,}iiSiII，即固定k和i,r分别取0r和1r。容易得到，在这两种情况下，分别有11112211312222iiiiiivvvvvv，；11112213112222iiiiiivvvvvv，。一般的，对于给定的i，k，()Si有k种不同的取法。可以证明（请自己证一下）：对于任意选取的r和k，rjc（j＝0,1,…,k-1）可以写成下面的一般形式：1/21/200,11/21/20()()kkiirqlqklmqmlrjirjkmjirmirlllmxxcxxx对于均匀网格，ixx，此时rjc与i及x无关：00,10(1)()kklqklmqmlrjkmjllmrqcml通过引入原函数确定()px比直接求解（9）式更为方便。首先我们定义函数()vx的原函数为1/2()()xxVxvd(11)则：1/21()iijjjVxvx所以通过单元内()vx的平均值我们可以得到其原函数在半节点上的精确值。在给定模板点集的情况下，我们可以构造单元iI中1/2()()xxVxvd的Largrange插值多项式()iPx（也记做()Px），则()()iiidPxpxxIdx要达到k阶精度，()iPx的阶数为k。很显然，插值多项式()iPx的模板()Si与()Si不同。但二者之间存在一一对应的关系，如果(){,...,}irisSiII则1/21/2(){,...,}irisSixx。三、一维ENO重构一维ENO重构的基本思想就是采用自适应模板点，尽可能的避免将不连续的单元包含到模板点集里面，以便消除或减小在间断附近的振荡。为了自适应调节模板，采用Newton插值多项式计算()iPx往往更加方便。通过考察()Vx的Newton插值多项式，我们可以发现Newton均差是函数在模板点集内是否光滑的量度。均差的定义为：1/21/21/21/21/23/21/21/21/21/2[]()(0)[,,][,,][,,]()iiiijiijiijijiVxVxthorderVxxVxxVxxjthorderxx这样我们就可以得到ENO重构的基本过程为：假设我们的目的是寻找一个由包含12ix和12ix的k+1个连续点构成的模板点集，使得()Vx在这个模板点集内相对其它可能的模板点集是最光滑的。我们将这项工作分几步来完成，每一步都只添加一个点到模板点集中。首先从包含两个点的模板点集21122(){}iiSixx，(12)开始，这个模板点集内的Newton插值多项式可以写成111112222()[][]()iiiiPxVxVxxxx，(13)下一步，我们有两种选择来向2()Si中添加一个点，如果添加左边的点32ix，那么得到的插值多项式为13111122222()()[]()()iiiiiRxPxVxxxxxxx，，(14)如果添加右边的点32ix，那么得到的插值多项式为11131122222()()[]()()iiiiiSxPxVxxxxxxx，，(15)我们注意到(14)、(15)式中都是以1122()()iixxxx分别乘上311113222222[][]iiiiiiVxxxVxxx，，，，，(16)它们分别是不同模板点集内()Vx的二阶Newton均差。从前面我们可以知道，较小的均差表示函数在该模板点集内比在其它模板点集内光滑。所以，要决定添加哪个点，我们可以通过比较这两个均差，并选取其绝对值较小的一个。如果311113222222[][]iiiiiiVxxxVxxx，，，，(17)那么就选取包含3个点的模板点集为3311222(){}iiiSixxx，，此时，2()()PxRx。否则，就选取3113222(){}iiiSixxx，，此时，2()()PxSx。这个过程可以继续进行下去，直到包含预定数目的点的模板点完全确定为止。一旦找到了()Si，也就是找到了(8)式中的()Si，我们就可以利用常数rjc来计算在单元边界上的重构数值了。对于分片光滑函数()Vx，以上重构过程具有如下几条性质：1、精度条件1()()()kiPxVxOx(18)对于任何不包含间断的单元iI都是成立的2、在包含()Vx的间断的任何单元里，()iPx都是单调的。3、该重构是总变差有界的（TVB，TotalVariationbounded），也就是说，存在函数()zx，如果对任何单元（包括包含间断的单元）都满足1()()()kiizxPxOxxI，则有()()TVzTVV成立。这几个性质对于理解ENO格式十分重要。性质1显然成立，性质2则不容易从直观上把握。性质2的完整证明比较复杂，我们仅就一种简单的情况给出说明。假定()Vx是一个阶越函数0,0()1,0xVxx，()Px是在1k个点1/23/21/2kxxx上对()Vx的插值。在001/21/20jjxx间存在间断。对于不包含间断的单元1/21/20[,]jjxxjj有1/21/21/21/2()()()()jjjjPxVxVxPx。因此，至少有一个点1/21/2(,)jjjxx，成立()0jP。由于有k