高等计算流体力学-02

1、第二讲有限差分方法基本原理有限差分方法基本原理1有限差分方法概述有限差分方法概述2有限差分法概述有限差分法概述 22xutu1. 差分差分方法方法Ttx0 , 10离散点上利用离散点上利用Taylor展开，把展开，把微分微分转化成转化成差分差分! 2221xxuxuuxujjjj)(1xOxuuxujjj111112121njnjnjnjnjuuuxtuu j-2 j-1 j j+1 .)(! 31)(! 21)(3332221xxuxxuxxuuujjjjj)(2211xOxuuxujjj（等距网格）（等距网格）多维问题，各方向独自离散；（时间同样考虑）多维问题，各方向独自离散；（时间同样考

2、虑）)(1tOtuutunjnjnj比有限体积法计算量小；比有限体积法计算量小；便于构造高阶格式便于构造高阶格式；23311)(! 312xxuxuuxujjjj3基本概念：基本概念： a. 差分表达式及截断误差差分表达式及截断误差:截断误差截断误差差分表达式差分表达式（1阶）（2阶）b. 前差、后差、中心差前差、后差、中心差 j-2 j-1 j j+1 前前)(1xOxuuxujjj前差前差)(2211xOxuuxujjj中心差中心差)(1xOxuuxujjj后差后差其他：其他： 向前（后）偏心差分向前（后）偏心差分； 后后c.差分方程差分方程 经差分离散后的方程，称为差分方程经差分离散后的

3、方程，称为差分方程 精度精度如何确定精度？如何确定精度？ 1） 理论方法，理论方法， 给出误差表达式给出误差表达式 2）数值方法，）数值方法， 给出误差对给出误差对 的数值依赖关系的数值依赖关系x011xuuatuunjnjnjnj0 xuatu微分方程微分方程差分方程差分方程xxuxuuxujjjj221! 2123311)(! 312xxuxuuxujjjjxxuattuTEnjnj2222! 2! 21截断误差：截断误差：4d. 显格式及隐格式显格式及隐格式显格式：显格式： 无需解方程组就可直接计算无需解方程组就可直接计算n+1层的值；层的值；隐格式：隐格式： 必须求解方程组才能计算必须

4、求解方程组才能计算n+1层的值层的值011xuuatuunjnjnjnj01111xuuatuunjnjnjnje. 守恒型差分格式守恒型差分格式基本思想：基本思想： 保证（整个区域）积分守恒律严格满足保证（整个区域）积分守恒律严格满足 0 xuftu 定义：对于定义：对于上述上述守恒型方程守恒型方程，差分格式，差分格式njnjnjnjgghuu21211称为守恒型差分格式。称为守恒型差分格式。),(2121nljnljnljnjuuugg其中：nnNnjnjNjgggg2/12121211特点：特点： 消去了中间点上的值，只保留两端消去了中间点上的值，只保留两端物理含义：物理含义： 只要边界

5、上没有误差，只要边界上没有误差，总体积分总体积分方程方程不会有任何误差。不会有任何误差。1njju如果如果 是准确的，则是准确的，则 也是准确的也是准确的 （假设边界条件没有误差）（假设边界条件没有误差）njju守恒性的例子：守恒性的例子： 环形管道里的流动环形管道里的流动 总质量保持不变总质量保持不变 早期 极为强调守恒性 最近 重新认识5关于守恒性格式的一些注解关于守恒性格式的一些注解 xffxfjj2/12/1中的符号中的符号 2/1jf与函数与函数f 在在 点的值点的值无关无关！2/1j),.,(12/1ljljljjuuuff是是j点周围几个点上点周围几个点上 f (或者或者u)值的

6、函数，值的函数， 为一记号，请勿理解为为一记号，请勿理解为j+1/2点的值点的值 ！1）2） 常系数线性格式都是守恒的常系数线性格式都是守恒的)(126154132231jjjjjjjfafafafafafaxxf例如，差分格式：等价于xffxfjjj2/12/1其中2514312212/1jjjjjjfbfbfbfbfbf,.)3 , 2(;111kabbabkkk守恒方程守恒方程+ 守恒格式守恒格式= 守恒解守恒解67f. 传统型（非紧致）差分格式及紧致型差分格式传统型（非紧致）差分格式及紧致型差分格式传统型：传统型： 运用多个点函数值的组合逼近运用多个点函数值的组合逼近一点的导数一点的导

7、数 j-2 j-1 j j+1 123121.mkjmkjkjkjjfafafafaxfxfffjjj/ )(1xfffffffjjjjjjj60/ )3302060152(21123紧致型：紧致型： 多个点函数值的组合逼近多个点函数值的组合逼近多个点导数值的组合多个点导数值的组合11122134152jjjjjjjjFFFa fa fa fa fa fxffffFFFjjjjjjj36/)()(28(3/13/1221111例：例：例：例：xffFFFjjjjj2/ )( 34/14/11111xfffffFFjjjjjjj120/ )1236443(5/35/221111jjxfF联立求解

8、联立求解 ， 多对角方程多对角方程 追赶法求解（追赶法求解（LU分解法）分解法）jF 紧致格式：紧致格式： 同样的基架点，可构造更高阶格式同样的基架点，可构造更高阶格式 （因为自由参数更多）（因为自由参数更多） （最高）精度（最高）精度=自由参数个数自由参数个数-12. 构造差分格式的基本方法构造差分格式的基本方法 待定系数法待定系数法)(31431221Ouauauauaxujjjjj j-2 j-1 j j+1)()(! 31)(! 21)()()(! 31)(! 21)()()2(! 31)2(! 21)2(43)3(2143)3(2143)3(22 OuuuuuOuuuuuOuuuuu

9、jjjjjjjjj)(! 3/)10) 1()2(! 2/)10) 1()2()02()(4343332313)3(242322212432143211431221 Oaaaauaaaauaaaauaaaauuauauauajjjjjjjj0)1 (0)1()2(0)1 (0)1()2(1)1 (0)1()2(0)1 ()1()2(4333231342322212413121114032010aaaaaaaaaaaaaaaa3 , 2 , 1 , 0)3(41kbajjkjkotherkbk01/1解出解出akxffxfjjj2/12/1 （可选）化成守恒型（可选）化成守恒型小程序：小程序：

10、求系数求系数8) 1 . 3()(.123121mmkjmkjkjkjjOuauauauaxu更一般的情况：更一般的情况： m+1个基架点上构造的个基架点上构造的m阶差分格式：阶差分格式：要善于用数值计算的手段要善于用数值计算的手段研究研究CFD ， 不能仅限于不能仅限于用理论手段研究用理论手段研究CFD !基架点基架点 （stencil ）3. 复杂网格的处理方法复杂网格的处理方法1） 一维情况：一维情况： 非均匀网格非均匀网格 方法方法1 （常用）：（常用）： 网格（网格（Jacobian）变换）变换 j-2 j-1 j j+1 非均匀网格)(xx )1/()1(Nii0,1的均匀网格的均

11、匀网格)(iixxdxdfxf 将方程由物理空间变到计算空间将方程由物理空间变到计算空间 （以（以x 为自变量变为以为自变量变为以 为自变量）为自变量）dxd为已知函数为已知函数)(xx 常用的一维坐标变换函数：常用的一维坐标变换函数： 指数函数指数函数 双曲正切函数双曲正切函数9)tanh(/ )tanh(gjgjbbx物理坐标物理坐标 计算坐标计算坐标 要求：要求： 坐标变换必须足够坐标变换必须足够光滑，否则会降低精度光滑，否则会降低精度网格间距变化要缓慢，否则网格间距变化要缓慢，否则会带来较大误差会带来较大误差方法方法2） 在非等距网格上直接构造差分格式在非等距网格上直接构造差分格式 j

12、-2 j-1 j j+1 )(31431221Ouauauauaxujjjjj4131)3(21114131)3(21114232)3(2222)()(! 31)(! 21)()()(! 31)(! 21)()()(! 31)(! 21)(jjjjjjjjjjjjjjjjjjjjjjjjjjjjjjjjjxxOxxuxxuxxuuuxxOxxuxxuxxuuuxxOxxuxxuxxuuu 原理：原理： 直接进行直接进行Taylor展开，构造格式展开，构造格式 格式系数是坐标（或网格间距）的函数格式系数是坐标（或网格间距）的函数0)(0)()(0)(0)()(1)(0)()(0)()()(431

13、3323113242132221122411312111124013201102axxaaxxaxxaxxaaxxaxxaxxaaxxaxxaxxaaxxaxxjjjjjjjjjjjjjjjjjjjjjjjj解出系数解出系数jjjjaaaa4321,注：注： 系数随网格点系数随网格点(j)变化！变化！10 网格非光滑、间距剧烈变化不会降低精度；网格非光滑、间距剧烈变化不会降低精度； 随机网格都可保证精度随机网格都可保证精度2） 二维二维/三维情况三维情况坐标变换坐标变换 均匀的直角网格均匀的直角网格RAE2822翼型周翼型周围的网格围的网格),(),(),(zzyyxx123123UfffVV

14、Vtxyzxyz ,TUuvw E 21, ()Tfuupuvuw u Ep),(),(zyxxxxfffxf1111三个方向共需计算三个方向共需计算9次导数，次导数，计算量大计算量大yyyfffyf2222zzzfffzf3333对流项可组合，求对流项可组合，求3次导数即可次导数即可321321VVVffftUUJU1),(),(1zyxJ)(32111fffJfzyx)(32112fffJfzyx)(32113fffJfzyx11.0)()()(111xxxJJJ4. 时间项的离散时间项的离散1）直接离散法）直接离散法 把时间导数直接差分离散把时间导数直接差分离散 0 xuftu0)(1n

15、xnnuftuu0)(11nxnnuftuu0)()(2111nxnxnnufuftuu1阶阶Euler显格式显格式1阶阶Euler隐格式隐格式2阶阶Crank-Nicolson格式格式0)()1 ()(11nxnxnnufuftuu2） Runge-Kutta 格式格式)(ULtU目前最常使用的：目前最常使用的：3步步3阶阶TVD型型R-K)( 3/23/1)( 4/14/3)()2()2(1)1()1()2()1(UtLUUUUtLUUUUtLUUnnnnn推荐！推荐！半离散格式半离散格式!12 在某一点进行在某一点进行Taylor展开，构造格式展开，构造格式3） 时时-空耦合离散空耦合离

16、散0 xuctun+1nj-1 j j+1),(txu 蛙跳格式蛙跳格式 0 xuftu0221111xfftuunjnjnjnjn,jLax-Wandrof格式格式cucxtuuxtcuuxxnjnjnjnj2111212njnjnjnjnjffxtuuu112121212212121211njnjnjnjffxtuu半隐错点格式半隐错点格式022)(11111njnjnjnjnjnjuuuuxctuunjnjnjnjffxtuu11MacCormack格式格式1111121)(21njnjnjnjnjffxtuuu135 特征理论与差分格式特征理论与差分格式，CFL条件条件特征性质是双曲型

17、方程的重要特点。在构造差分格式是考虑到微特征性质是双曲型方程的重要特点。在构造差分格式是考虑到微分方程的数学物理性质，有助于得到性态较好的差分格式。分方程的数学物理性质，有助于得到性态较好的差分格式。 j-3 j-2 j-1 j j+1 n+1n000方程的特征线为，沿特征线。考虑的情形，参考上图，我们知道：uuaxatcuconstatxa1(,(1)(,)njjOPjuu xntuuu xa t n t 由于P点一般不在网格点上，P点的值必须由A、B、C等各点通过插值来获得。141）采用B，C两点线性插值，得令，上式为tcaxCBPBuuuuxa tx 111nnjjnnjjuuuuxa

18、tx 111nnnjjjuc ucu110nnnnjjjjuuuuatx. . .()()01(0)LT EOxtca 这就是所谓的迎风格式。其截断误差为，稳定性条件 j-3 j-2 j-1 j j+1 152DBPBuuuuxa tx 11111122nnnnnjjjjjcuuuuu111111()202nnnnnjjjjjuuuuuatx. . .()()1LT EOxtc Lax格式的截断误差为，稳定条件是2）采用）采用B,D两点进行线性插值两点进行线性插值该格式称为Lax格式（或LaxFriedrichs格式）。即或 j-3 j-2 j-1 j j+1 163)采用采用A.B.C.三点

19、进行二次插值。三点进行二次插值。 212(), (), ()取二次曲线，满足oAABBCCu xaa xa xu xu u xu u xu22222CAACBPBuuuuuuuxa txa txx 201212,oa a au xaa xa x定出系数。把P点坐标代入，得到：21122112122123422234222nnnnnnnnjjjjjjjjnnnnnjjjjjnnnjjjccuuuuuuuuuuuuuatauuutxx。. . .()() )02LT EOxtc 22Warming Beam格式的截断误差为稳定条件是（a0）该格式称为WarmingBeam格式。 j-3 j-2 j

20、-1 j j+1 174)采用采用B，C，D三点进行二次插值三点进行二次插值 212(), (), ()取二次曲线，满足oBBCCDDu xaa xa xu xu u xuu xu 201212,oa a au xaa xa x定出系数。把P点坐标代入，得到：22222BDCDBPBuuuuuuua ta txx 2111111211112222222nnnnnnnjjjjjjjnnnnjjjjnnnjjjccuuuuuuuuuuuatauuutxx。. . .()() )1LT EOxtc 22Lax Wendroff格式的截断误差为稳定条件是该格式称为LaxWendroff格式。 j-3

21、j-2 j-1 j j+1 180(通过上面的几个格式的构造，我们可以看出：(1) 当时，用及其左侧的点进行插值，得到的是迎风型格式 Warming-Beam可以看作二阶迎风格式)，其稳定性条件与a的正负有关；(2)当用及与之对称的点进行插值，得到的是中心型格式，此时，其稳定性条件与a正负无关。njnjauu0重要现：当差分格式稳定时，P点必然在插值点的中间， 即由插值点内插出P点的值，差分格式才可能是稳定的。可以验证，当P点的值由插值点外推得到时，格式必然不稳定，如，可以用C，D外推出P点的值，此时得到的所谓的顺风格式，该格式无条件不稳定。这种现象有一般意义,称为：Courant-Fried

22、richsLewy(CF件象L)条：a 19定理：双曲型方程的差分格式收敛（在定理：双曲型方程的差分格式收敛（在Lax等价性定理满足时等价性定理满足时亦即稳定）的必要条件是差分格式的依赖域包含微分方程的依亦即稳定）的必要条件是差分格式的依赖域包含微分方程的依赖域。赖域。一阶迎风(a0)差分格式依赖域差分格式依赖域CFL条件条件01c稳定条件稳定条件01c差分格式差分格式Warming-Beam（a0）02c02cEuler显式1c 无条件不稳定：CFL数a tcx01/tax20差分方法理论基础差分方法理论基础21 1） 相容性：相容性： 当差分方程中当差分方程中 ，时间与空间步长均趋近于，时

23、间与空间步长均趋近于0 时，差分方程的时，差分方程的截断误差截断误差也趋近于也趋近于0，则称差分方程与原微分方程是，则称差分方程与原微分方程是相容相容的。的。差分差分方法理论基础方法理论基础2）收敛性：）收敛性：0lim0,uuhtxL2 模：212)()(dxxuxu 模：)(max)(xuxuxLxuxujjh2)(jjhuxumax)(22 当时间与空间步长均趋近于当时间与空间步长均趋近于0 时，差分方程的时，差分方程的解解趋近于微分方程的解，趋近于微分方程的解，则称差分方程的解则称差分方程的解收敛收敛于原微分方程的解。于原微分方程的解。注意！注意！ 方程互相趋近方程互相趋近 解互相趋近

24、解互相趋近（多值性、奇异性多值性、奇异性 ）)()(lim00 xfxfxx不一定等于不一定等于只有连续函数才满足只有连续函数才满足 （根据（根据Lax等价定理，只有稳定性条件满足的等价定理，只有稳定性条件满足的情况下，方程趋近才能保证解趋近）情况下，方程趋近才能保证解趋近）含义：含义： 方程趋近方程趋近含义：含义： 解趋近（更强）解趋近（更强）)(),(xuxuh分别为差分方程和微分方程的解1. 相容、收敛、稳定性与相容、收敛、稳定性与Lax等价定理等价定理相似的例子：3） 稳定性稳定性：xt 和xt 和0021,exp,txuttcctxuhnnh定义：称差分方程的初值问题定义：称差分方程

25、的初值问题是稳定的，如果当是稳定的，如果当 做够小时，存在于做够小时，存在于 无关的常数无关的常数C1和和C2使得使得:含义：含义： 在差分方程的求解过程中，如果引入的误差随时间的增长有界，在差分方程的求解过程中，如果引入的误差随时间的增长有界，则称差分方程是稳定的。则称差分方程是稳定的。234） Lax 等价定理等价定理 如果线性微分方程如果线性微分方程的的初值问题初值问题是适定的，差分方程是相容的，则是适定的，差分方程是相容的，则差分方程解的差分方程解的收敛性收敛性与与稳定性稳定性是等价的。是等价的。含义：含义： 如果微分方程不出问题（适定），差分方程性质好（稳定），如果微分方程不出问题（

26、适定），差分方程性质好（稳定），则则方程逼近方程逼近就可保证就可保证解逼近解逼近。 如果方程逼近就可以导致解逼近，则差分方程的性质肯定好（稳定）如果方程逼近就可以导致解逼近，则差分方程的性质肯定好（稳定）2. 差分格式稳定性分析方法差分格式稳定性分析方法 njnjnjnjuuxcuut1111Fourier分析法：分析法： 基本思想：基本思想： 初始时刻引入单波扰动，考虑其随时间的变化初始时刻引入单波扰动，考虑其随时间的变化 原理：原理： 任何扰动都可认为是单波扰动的叠加；任何扰动都可认为是单波扰动的叠加； 线性情况下不同波之间独立发展线性情况下不同波之间独立发展。 引入单波扰动，带入差分方程

27、，如果其振幅放大，则不稳定；否则稳定引入单波扰动，带入差分方程，如果其振幅放大，则不稳定；否则稳定00cxuctujikxnnjeG引入单波扰动引入单波扰动：inneGGG11解出放大因子：解出放大因子：)2/(sin)1 (4122G1024xeecGetGGxikikxnikxnnjj)1 (11G稳定性条件：对所有稳定性条件：对所有 带入差分方程例例1： 考察右式给出差分格式的稳定性考察右式给出差分格式的稳定性一些注解一些注解：1,nnGGG 通常为复数；通常为复数； 可反映可反映振幅及相位；振幅及相位；xtc/ 称为库朗数称为库朗数稳定条件：稳定条件：xk 有效波数有效波数/2PPW

28、一个波里面的网格点数一个波里面的网格点数3. 差分方程的修正方程差分方程的修正方程修正方程修正方程 差分方程准确逼近（无误差逼近）的方程差分方程准确逼近（无误差逼近）的方程0 xuatu011xuuatuunjnjnjnjtttttnjnjnjututtuutu6221差分方程截断误差微分方程微分方程=差分方程差分方程+截断误差截断误差通常要求：通常要求： 修正方程中不出现时间的修正方程中不出现时间的高阶导数高阶导数项项 （便于进行空间分析）（便于进行空间分析）25.6221xxxxxnjnjnjuxuxxuuxu2202266txttxxtttxxxta xta xuauuuuu与差分方程等

29、价的微分方程称为修正方程!表表 1 自循环消去过程自循环消去过程 tu xu ttu txu xxu tttu ttxu txxu xxxu （*） 1 a 2t 0 2xa 62t 0 0 62xa *2tt 2t 2a t 0 42t 0 4xta 0 xta2 2a t 22at 0 24a t 0 42xta 22212tt 2121t 122ta 0 0 221(*)3a tt x 231ta 2213at 0 22221(*)34a t xatx 43122xtata 431223xtata 按列求和按列求和 1 a 0 0 12cxa 0 0 0 132622ccxa 22(1)

30、( 231)26txxxxxxa xa xuauc uccua tcx110nnnnjjjjuuuuatx2202266txttxxtttxxxta xta xuauuuuu264 差分格式的耗散和频散 考虑线性波动方程 的某一差分格式，其修正方程为：0uuatx2lllluuuatxx研究差分格式解的性质，相当于研究修正方程解的性质。 的初值问题存在解析解。设其初始条件为：2lllluuuatxx0( ,0)( )u xux假定 是周期为 的周期性函数，则当长度为 的区间划分为M等分时， 可以展开为离散的Fourier级数：0( ,0)( )u xuxLL0( ,0)( )u xux/20/

31、2 1( )(0)mMk xmmMuxAei27由由Fourier 级数级数的正交性，我们可以只考虑波数为的正交性，我们可以只考虑波数为 的分量的分量，此时，上式可以简写为mk0( )(0)mk xmuxAei2lllluuuatxx记 的解的波数为 的分量为：mk( , )( )mk xmu x tA t ei带入方程 得2( )()( )() ( )mmmk xk xk xlmmmlmmldA tea kA t ekA t edtiiiii2( )()() ( )lmmlmlmdAta kkdtAt ii22() () ln( )()() ( )lmlmllmmlmla kktmdA ta

32、 kkdtA tCe iiii28由初始条件，可知 ，所以(0)mCA2() () ( )(0)lmlmla kktmmA tAeii即2() () ( , )(0)lmmmlmakktk xmu x tAeeiii上式也可以写为：2222111221( 1) ( 1) ( 1) 2211( , )(0)(0),( 1)lllll mmlmllmlll mlktkxaktkx atmmmktllmlmlu x tAeeAg egeaak其中，ii在相同的初始条件下， 的解为0uuatx( , )(0)mkx atmu x tAei显然，数值解的幅值可能随时间变化，而解析解的振幅不随时间变化，二

33、者的比值是 。利用 可以分析差分格式的耗散（dissipation）特性。221( 1) lll mlktmgemg29 当 ，即 时，数值解的振幅是随时间衰减的，我们称差分格式具有正的耗散； 当 ， 即 时，数值解的幅值不变，这时称差分格式具有零耗散。 当 ，即 时，数值解的振幅随时间增长的，我们称差分格式具有负耗散。 具有负耗散的差分格式是不稳定的，具有正耗散的差分格式是稳定的，具有零耗散的差分格式是中性稳定的。利用这个性质，我们也可以通过修正方程进行稳定性分析。1mg 221( 1)0lllmlk1mg 221( 1)0lll mlk1mg 221( 1)0lll mlk注意到：只有修正

34、方程中的偶数阶导数项影响差分格式的耗散特性，我们称这些项为修正方程中的耗散项。对于一阶精度的差分格式，如迎风格式，其修正方程耗散项的主项是 ，所以 的符号主要由 决定，即如果差分格式是耗散（正耗散或者中性耗散）的，必须有 。222ux221( 1)lll mlk22030通过考察一阶迎风格式的修正方程式，差分格式有非负耗散则要求 。这和由Von Neumann方法得到的稳定性条件是一致的。22(1)( 231)26txxxxxxa xa xuauc uccu1c 对于二阶精度的格式，如LaxWendroff格式， ，所以耗散项的主项是 ， 即如果差分格式是耗散的，必须有 。由LaxWendro

35、ff格式的修正方程可得耗散非负的条件为 ，仍然与Von Neumann方法的结论相同。20444ux402322(1)(1)68txxxxxxxxa xa xuauc ucc u 1c 一般的，如果修正方程中所有 阶项的系数均为非负，所有 阶项的系数均为非正，则数值解的振幅不会随时间增长，差分格式是耗散的。在实际应用中，我们常常只考虑主项的影响，如在上面所作的那样。42 (1,2,)nn4 (1,2,)n n 31差分格式的耗散性质，造成数值解和解析解之间幅值的差别。而数值解和解析解之间相位的差别，则取决于差分格式的频散（又称色散）(dispersion)特性。 为了解释所谓差分格式的频散性质，我们首先引入相速度的概念。所谓相速度，是指波数为 波动的传播速度。显然，由修正方程的解析解，知道相速度为 。如果相速度与 有关，称差分格式有色散（或者频散）。如果 ，则称差分格式有正的色散；如果 ，则称差分格式有负的色散 。 的精确解中，各个波数的波的传播速度都是 ，因此，精确解的波形在传播的过程中保持不变。但是，差分方程的数值解由于色散的作用，不同波数的波的传播速度是不同的，从而造成波形的畸变。mk2211( 1)ll