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

§10紧致格式一、紧致格式的定义及构造1.紧致格式的定义及待定系数的确定方法许多物理现象中包含有大范围的时间和空间尺度，湍流就是一个常见的例子，对于这些多尺度的物理现象的数值模拟要求适当地体现出所有相关的尺度，特别是小尺度的演化过程。这些要求促进了谱方法的研究，但是谱方法只能局限于简单计算区域和简单边界。为了克服谱方法的局限性，人们试图采用有限差分方法和谱元法，并在这方面取得了一定的进展。对于显式的高精度、高分辨率的有限差分法，必然需要较多的计算模板点，这样将带来边界点处理的困难。以均匀网格上的有限差分问题为例，在求某一点上物理量的导数（为简单起见，仅考虑一阶导数的情况）时，一般是利用该点附近的有限个网格点上的函数值构造导数的差分近似，即miiniiuuuLu,,,,'，(1)其中'iu代表对ixx处导数ux的差分近似。一般来说，差分算子L是线性的。此时，如果'()()iiuuOxx，(2)则mn。(3)也就是说，要构造阶精度的差分近似，至少需要1个模板点。目前，湍流的直接数值模拟、大涡模拟、计算气动声学领域对数值计算的精度提出了很高的要求。如果用常规的方法构造高阶格式，需要的模板点数也相应增加，这给边界附近数值方法的实施和边界处理带来了困难。紧致格式是一种可以用较少模板点逼近导数到高阶精度的差分格式，它除了具有边界处理相对简单的优点外，更重要的是具有类似谱方法的高分辨率；这一点对于计算具有复杂多尺度性质的流场（如湍流和声场）具有重要的意义。所谓紧致格式，是一种隐式的数值微分计算格式。即''',,,,,,,,iliikiniimQuuuLuuu(4)一般情况下，Q也是线性算子。于是，i点上物理量的导数值'iu不再能单独求出，而需要同网格线上的其它点（在附加的边界条件下）一起求解。下面讨论一般的线性紧致格式的构造过程。线性紧致格式可以写为：'''001()lilikikninimimuuuuuux。(5)如何确定(5)式的系数，使得'iu逼近ux到指定的阶数呢？我们把左右两侧在ixx处作Taylor展开，有'''2'''3'''''''2'''3''''1234012341iiiiiiiiicucxucxucxudududxudxudxux(6)对(6)式求导，有'''''2'''''''2'''3''''12301231iiiiiiicucxucxudududxudxux(7)''''''''''''''''120121iiiiicucxudududxux(8)……把(7)式乘以21cxc与(6)式合并，有2'2'''3''''22313411'''0212012112'''3''''22323411()()1()()()()iiiiiiiiccccucxucxuccdcdcdududxuxccdcdcdxudxucc(9)把(8)式乘以222311()xcccc与(9)式合并，有2'3''''22143112'''02120201231111222'''3''''221232223343111111((2))1()(())(())(())iiiiiiicccuccxuccdcdcdcdududcxuxccccdcdcdcdcdcxudcxucccccc(10)依次类推，可以消去(6)式中的高阶导数项。由(10)式，要使'iu逼近ux到一阶精度，需要满足：00211110()ddccddc要使要使'iu逼近ux到二阶精度，则应有00211112120223221110()(())0ddccddcdcdcdccdccc要使要使'iu逼近ux到三阶精度，则应有00211112120223221112221233331110()(())0(())0ddccddcdcdcdccdcccdcdcdccdccc可见，如果要构造阶精度的格式，只需(6)式两端前1个对应项的系数相等即可。在(5)式中有k+l+m+n+2个待定参数，按照上面的分析，该紧致格式最高可达k+l+m+n+1阶精度。当要求的精度低于k+l+m+n+1时，(5)中可以存在无穷多种满足精度要求的系数。特别的，当要求的精度为k+l+m+n阶时，(5)可以有一个系数可以任意给定；当要求的精度为k+l+m+n-1阶时，(5)可以有两个系数可以任意给定。其余可依此类推。在确定系数后，实际求解'iu的过程中，(5)式的右端是已知的。左端的未知量构成了一个有k+l+1个对角元素的方程组，需要结合边界条件联立求解。一般三对角线方程可以用高效的追赶法求解，所以常取k=1,l=1。2．中心型紧致格式如果在(5)式中,klmn且,jjjj，则称相应的紧致格式为中心型紧致格式。例如左侧三点、右侧5点的中心型紧致格式为'''11221124iiiiiiiuuuuuuuabhh(11)式中hx为网格间距，方程组的系数矩阵为三对角阵。将上式中的各项进行Taylor展开，并令各阶项（420,,hhh）的系数相等，得到!4!522!2!3222142bababa,,ba满足第一式时有二阶精度，满足前两式时有四阶精度，三式都满足时有六阶精度。四阶精度的格式满足)14(),2(3132ba，其中可以任取。其截断误差为4)5(!54)13(hu如果0，则格式退化为四阶中心差分格式，即huuuuuiiiii12882112'截断误差为4)5(!54hu如果31，那么第三式满足，格式具有六阶精度，此时截断误差为6)7(!74hu。从上面的例子可以看出，在相同模板点下，常规的差分方法精度最高为四阶，而紧致格式的精度可以达到六阶。2．迎风型紧致格式迎风型紧致格式是一种偏心（左右不对称）紧致格式。考虑方程0uuctx，0c时空间导数用下式计算''1121iiiiiiuuuuuuabhh(12)而0c时计算方法为：''1211iiiiiiuuuuuuabhh。(13)这种类型的格式即为迎风格式。以0c的情况为例分析其精度，利用Taylor展开并令各阶项（210,,hhh）的系数相等，得到21!38!31!23!211bababa,,ba满足前两式具有两阶精度，满足所有三式时具有三阶精度。由前两式得到21,23ba截断误差为2)3(121)53(hui如果0，则格式退化为二阶迎风格式，即huuuuiiii24321'截断误差为2)3(125hui如果35，那么第三式满足，格式具有三阶精度，此时截断误差为3)4(3623hui0c时，容易知道,,ba的关系与前述情况相同。二、紧致格式误差的Fourier分析1．分辨率分析设()fx是周期为L的周期性函数，把[0，L]上划分为N等分，网格间距h=L/N，对()fx作Fourier级数展开：/2/22ˆ()exp()kNkkNkxfxfLi(14)定义无量纲化波数与坐标分别为22khkLN和xsh，则解可以表示为：ˆ()exp()kfsfsi(15)对()fs求一次导数后，得到：ˆ()exp()kfsfsii(16)假定用紧致差分得到的导数可以写为：ˆ()exp()fdkfsfsii(17)于是可以用()来衡量差分后的导数值逼近真实导数值的程度。如果规定某个0值，则总可以找到一个区间2[,]fN，使得位于这个区间内的所有波数，逼近真实导数的误差都小于0。所以fr代表了格式在0条件下的分辨能力的有效范围，r越大，格式的分辨率越高。一般称r为分辨效率。下面推导'的公式。为了简便起见，以三点中心型紧致格式为例，即：112211'''24iiiiiiifffffffabhh(18)将方程当中各项作Fourier展开，只考虑一个频率成分，可得：112111(1)()exp()(1)()exp()(2)()exp(2)(2)()exp(2)''(1)'()exp()''(1)'()exp()iiiiiiiiiiiiiiiiiiffsfsffsfsffsfsffsfsffsfsffsfsiiiiiiii其中：ˆ()exp()kfsfis。把上面的关系代入(18)式，化简得到：sin()(/2)sin(2)(/3)sin(3)12cos()2cos(2)abc对于我们讨论的特定格式，为实数，但一般情况下，可能是复数。图1针对下面几种不同精度的格式，给出了曲线。另外，下面还给出了00.01时的分辨效率。显式格式（＝0）二阶：1,0,0,0.08abr（a）四阶：41,,0,0.2333abr（b）紧致格式（≠0）四阶：31,0,,0.3524abr（c）六阶：1411,,,0.5993abr（d）图1.中心型紧致格式的曲线2．色散性和耗散性考虑线性对流方程0uuctx，(19)其对应的半离散差分格式为0iiucut(20)其中的空间导数假定用中心型紧致格式'''11221124iiiiiiiuuuuuuuabhh(21)计算。首先我们讨论(19)式的精确解。设(19)式的解可以写为：(,)()exp()uxtvtsi(22)00.511.522.5300.511.522.53abcde代入(19)式，得：()()0cvtvthi。(23)积分上述方程并考虑初始条件(,0)(0)exp()uxvsi，有(,)(0)exp[()]ctuxtvshi(24)由(24)式，知(,)uxt的幅值为(0)v，相速度为c（与波数无关），因此方程(19)既无耗散，也无色散。下面讨论应用(21)式计算空间导数时，半离散方程(20)在时间方向的解析解。令(,)()exp()iiuxtvtsi，则(,)()exp()iiuxtvtsxii类似的，我们令'()exp()iiuvtsii，则利用上一小节的分析方法，易知sin()(/2)sin(2)(/3)sin(3)12cos()2cos(2)abc(25)由上面的关系，我们可以得到()()0cvtvthi(26)所以，(20)式的解为(,)(0)exp[()]iictuxtvshi。(27)对于我们讨论的紧致格式(21)式，是一个实数，但是一般情况下，可以是复数。为了使我们的讨论更具一般意义，我们把记为Re()Im()i，(28)则(27)式可以改写为Re()(,)(0)exp[Im()]exp[()]iictctuxtvshhi。(29)此时，(,)iuxt的幅值为(0)exp[Im()]ctvh，相速度为Re()c。当Im()0c时，差分格式有正耗散；对于我们讨论的中心型紧致格式Im()0c，无耗散；当Im()0c时，格式是不稳定的。当Re()10差分格式有正色散；R