计算流体力学(宗智)_第03讲-差分方法

1、计算流体力学讲义计算流体力学讲义 第三讲第三讲 有限差分法（有限差分法（1）宗智宗智Tel: 84707694; 船池船池317房间房间知识点：知识点： 差分方法的理论基础差分方法的理论基础 （相容、收敛、稳定性；（相容、收敛、稳定性；Lax等价定理；精度、修正方程等价定理；精度、修正方程; 守恒性）守恒性） 差分格式的构造差分格式的构造 差分格式的差分格式的Fourier分析分析 1传统计算方法：传统计算方法： 有限差分法，有限差分法， 有限体积法有限体积法 ， 有限元法，有限元法， 谱方法（谱元法）等；谱方法（谱元法）等；最近发展的方法最近发展的方法: 基于粒子的算法（格子基于粒子的算法（

2、格子-Boltzmann, BGK），无网格），无网格 优点优点缺点缺点适用范围适用范围有限差分法有限差分法简单成熟，可构造高精简单成熟，可构造高精度格式度格式处理复杂网格不够灵处理复杂网格不够灵活活相对简单外形的相对简单外形的高精度计算高精度计算有限体积法有限体积法守恒性好，可处理复杂守恒性好，可处理复杂网格网格不易提高精度（二阶不易提高精度（二阶以上方法复杂）以上方法复杂）复杂外形的工程复杂外形的工程计算计算有限元法有限元法基于变分原理，守恒性基于变分原理，守恒性好好对于复杂方程处理困对于复杂方程处理困难难多用于固体力学多用于固体力学等等谱方法谱方法精度高精度高外形、边界条件简单外形、边界

3、条件简单简单外形的高精简单外形的高精度计算度计算粒子类方法粒子类方法算法简单，可处理复杂算法简单，可处理复杂外形外形精度不易提高精度不易提高复杂外形的工程复杂外形的工程计算计算 第三讲第三讲 有限差分法（有限差分法（1） （教材第（教材第3.1、3.2节及第节及第4章）章）23流体流体力学力学理论研究理论研究实验研究实验研究数值研究数值研究 计算流体力学计算流体力学（数值计算技术、（数值计算技术、计算方法研究）计算方法研究）理论研究：理论研究：格式推导、格式推导、 稳定性分析，稳定性分析，精度精度/误差分析，误差分析，实验研究：实验研究：数值实验，数值实验， 采用采用实际问题考核方法实际问题考

4、核方法的正确性的正确性数值研究：数值研究：采用数值计算推导格式、考察精度采用数值计算推导格式、考察精度/稳定性稳定性/分辨率分辨率“计算流体力学计算流体力学”作为一个学作为一个学科，科， 其研究手段依然包括理论、其研究手段依然包括理论、实验及数值模拟。实验及数值模拟。 与与 的依赖关系的依赖关系4举例说明：举例说明： 研究研究“计算流体力学计算流体力学” 学科的学科的理论手段理论手段、实验手段实验手段及及计算手段计算手段研究研究CFD的理论手段的理论手段例：例：Fourier分析分析线性系统：线性系统： 线性方程线性方程+ 线性格式线性格式 任意函数都可分解为三角函数的叠加任意函数都可分解为三

5、角函数的叠加差分系统差分系统(解差分解差分方程方程)初始值初始值数值解（特数值解（特定时刻离散定时刻离散的函数值）的函数值）jjuv记为： 是差分算子，把离散函数是差分算子，把离散函数（有限点列）（有限点列） 映射为另映射为另一个离散函数一个离散函数jujvkikxkjjeuukikxkjjevv kvkuvi 与与ui的依赖关系的依赖关系线性系统，可大为简化线性系统，可大为简化波数空间单一的依赖关系：波数空间单一的依赖关系：)(kkufv ),.)(),( (2211ufvufv原理：原理： 线性系统，输入一个波，只能输线性系统，输入一个波，只能输出一个波（且波数不变）。出一个波（且波数不变

6、）。 非线性系统会产生多个谐波非线性系统会产生多个谐波 线性差分系统： 针对一个单波一个单波， 研究经过差分系统后的变化就可以了解该系统。 Fourier误差分析； Fourier稳定性分析理论分析的局限性：理论分析的局限性： 对于复杂系统（非线性方程、非线性格式）非常困难对于复杂系统（非线性方程、非线性格式）非常困难5研究研究CFD的实验手段的实验手段例： 精度分析思想：思想： 通过具体算例来研究（考核，分析通过具体算例来研究（考核，分析）差分方法）差分方法典型的文章： 提出方法+理论分析 + 算例验证dxdu差分离散jxjuF理论方法，理论方法，Taylor展开，求余项展开，求余项 。 对

7、对于复杂（如非线性）格式，难度大。于复杂（如非线性）格式，难度大。实验方法，通过算例考核精度实验方法，通过算例考核精度)cos()sin(xdxduxu精确解)cos(jjxFerrj 为该离散函数的模为该离散函数的模 常用的模：常用的模： 1 模：模： 2 模：模： 无穷模无穷模：jjj1jjjN221jjjmax计算误差 分析误差对网格步长的依赖关系xlgerrlg 斜率为精度的阶数（通常用最小二乘法计算）lglglgxnerrxerrn斜率为精度的阶数nCopyright by Li Xinliang6常用的验证算例（常用的验证算例（“实验验证实验验证”） 考核方法通常找一些难度大的（条

8、件苛刻、极端）的算例。否则，考核方法通常找一些难度大的（条件苛刻、极端）的算例。否则，无法突出方法的优越性。无法突出方法的优越性。1维算例：维算例： Sod 激波管，激波管， Shu-Osher, 方波方波/尖波尖波 2维算例：维算例： 前前/后台阶、双马赫反射、二维后台阶、双马赫反射、二维Riemann问题、漩涡问题、漩涡-激波干扰、翼激波干扰、翼型扰流、圆柱绕流型扰流、圆柱绕流3维复杂算例：维复杂算例： 各向同性湍流的各向同性湍流的DNS, 槽道湍流的槽道湍流的DNS, 激波激波-边界层干扰的边界层干扰的DNS Shu-Osher问题的计算结果问题的计算结果 （Li et al. Init

9、. J. Num. Fluid. 2005)航空领域权威的考核算例航空领域权威的考核算例 DPW标准计算模型标准计算模型常用一、二维算例整理后已发到流体中文网常用一、二维算例整理后已发到流体中文网 Copyright by Li Xinliang7研究研究CFD的计算手段的计算手段例：例： 差分格式构造差分格式构造理论方法：理论方法： 手工推导系数（工作量大）手工推导系数（工作量大） 数值方法：数值方法： 通过数值手段推导系数通过数值手段推导系数3 , 2 , 1 , 0)3(41kbajjkjk)(31431221Ouauauauaxujjjjj数值求解，获得系数数值求解，获得系数 格式优化

10、；格式优化； 通过数值计算手段进行通过数值计算手段进行 Fourier分析分析; 3.1 有限差分法基本原理有限差分法基本原理22xutu1. 差分方法的基本原理差分方法的基本原理Ttx0 , 10离散点上利用离散点上利用Taylor展开，把展开，把微分微分转化成转化成差分差分! 2221xxuxuuxujjjj)(1xOxuuxujjj111112121njnjnjnjnjuuuxtuu j-2 j-1 j j+1 .)(! 31)(! 21)(3332221xxuxxuxxuuujjjjj)(2211xOxuuxujjj（等距网格）（等距网格）多维问题，各方向独自离散；（时间同样考虑）多维

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

12、程，称为差分方程经差分离散后的方程，称为差分方程 精度精度如何确定精度？如何确定精度？ 1） 理论方法，理论方法， 给出误差表达式给出误差表达式 2）数值方法，）数值方法， 给出误差对给出误差对 的数值依赖关系的数值依赖关系x011xuuatuunjnjnjnj0 xuatu微分方程微分方程差分方程差分方程xxuxuuxujjjj221! 2123311)(! 312xxuxuuxujjjjxxuattuTEnjnj2222! 2! 21截断误差：截断误差：9Copyright by Li Xinliangd. 差分方程的修正方程差分方程的修正方程修正方程修正方程 差分方程准确逼近（无误差逼近

13、）的方程差分方程准确逼近（无误差逼近）的方程0 xuatu011xuuatuunjnjnjnjtttttnjnjnjututtuutu6221.62622211tttttxxxxxnjnjnjnjututuxuxxuuatuuxuatu差分方程截断误差微分方程微分方程=差分方程差分方程+截断误差截断误差.62622211tttttxxxxxxtnjnjnjnjututuxauxuauxuuatuu 差分方程差分方程=微分方程微分方程-截断误差截断误差 新的微分方程（修正方程）新的微分方程（修正方程）通常要求：通常要求： 修正方程中不出现时间的高价导数项修正方程中不出现时间的高价导数项 （便于进

14、行空间分析）（便于进行空间分析）011xuuatuunjnjnjnj等价于0.626222tttttxxxxxxtututuxuxuau修正方程0.662222xxxtttxxttxtuxutuxcutuau),(),(32xtOucuxtOucuxxxtttxxtt,) 132(612322322ttxtxxOuxcuxccuuxxxxxxtxtc修正方程主导项：主导项： 1阶；阶； 耗散型耗散型10.6221xxxxxnjnjnjuxuxxuuxud. 显格式及隐格式显格式及隐格式显格式：显格式： 无需解方程组就可直接计算无需解方程组就可直接计算n+1层的值；层的值；隐格式：隐格式： 必须

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

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

17、 f (或者或者u)值的函数，值的函数， 为一记号，请勿理解为为一记号，请勿理解为j+1/2点的值点的值 ！1）2） 常系数线性格式都是守恒的常系数线性格式都是守恒的)(126154132231jjjjjjjfafafafafafaxxf例如，差分格式：等价于xffxfjjj2/12/1其中2514312212/1jjjjjjfbfbfbfbfbf,.)3 , 2(;111kabbabkkk3） 关于关于2/1jf得到得到 后，将后，将j替换成替换成j-1即可得到即可得到 2/1jf2/1jf无需单独计算无需单独计算 ！ （白白增加计算量）（白白增加计算量） 2/1jf守恒方程守恒方程+ 守恒

18、格式守恒格式= 守恒解守恒解12Copyrigh by Li XinliangCopyright by Li Xinliang13f. 传统型（非紧致）差分格式及紧致型差分格式传统型（非紧致）差分格式及紧致型差分格式传统型：传统型： 运用多个点函数值的组合逼近运用多个点函数值的组合逼近一点的导数一点的导数 j-2 j-1 j j+1 123121.mkjmkjkjkjjfafafafaxfxfffjjj/ )(1xfffffffjjjjjjj60/ )3302060152(21123紧致型：紧致型： 多个点函数值的组合逼近多个点函数值的组合逼近多个点导数值的组合多个点导数值的组合2514312

19、2111jjjjjjjjuauauauauaFFFxffffFFFjjjjjjj36/)()(28(3/13/1221111例：例：例：例：xffFFFjjjjj2/ )( 34/14/11111xfffffFFjjjjjjj120/ )1236443(5/35/221111jjxfF联立求解联立求解 ， 多对角方程多对角方程 追赶法求解（追赶法求解（LU分解法）分解法）jF 紧致格式：紧致格式： 同样的基架点，可构造更高阶格式同样的基架点，可构造更高阶格式 （因为自由参数更多）（因为自由参数更多） （最高）精度（最高）精度=自由参数个数自由参数个数-1Copyright by Li Xinl

20、iang14一些记号一些记号 约定：约定： jjxxff为一阶偏导数的差分算子为一阶偏导数的差分算子jjxxxff22为二阶偏导数的差分算子为二阶偏导数的差分算子xfffjjjx1xfffjjjx1分别为分别为一阶精度一阶精度前、后差的前、后差的差分算子差分算子（本讲义中，上面两个算子表示的差分格式形式可以任（本讲义中，上面两个算子表示的差分格式形式可以任意，意， 包括线性包括线性/非线性、低阶非线性、低阶/高阶、普通高阶、普通/紧致紧致）jjjjxfxfff2)(2110为二阶中心差分算子为二阶中心差分算子上面三个算子有固定含义上面三个算子有固定含义。2. 构造差分格式的基本方法构造差分格式

21、的基本方法 待定系数法待定系数法)(31431221Ouauauauaxujjjjj j-2 j-1 j j+1)()(! 31)(! 21)()()(! 31)(! 21)()()2(! 31)2(! 21)2(43)3(2143)3(2143)3(22 OuuuuuOuuuuuOuuuuujjjjjjjjj)(! 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(

22、)2(0)1 ()1()2(4333231342322212413121114032010aaaaaaaaaaaaaaaa3 , 2 , 1 , 0)3(41kbajjkjkotherkbk01/1解出解出akxffxfjjj2/12/1 （可选）化成守恒型（可选）化成守恒型小程序：小程序： 求系数求系数15Copyright by Li Xinliang) 1 . 3()(.123121mmkjmkjkjkjjOuauauauaxu更一般的情况：更一般的情况： m+1个基架点上构造的个基架点上构造的m阶差分格式：阶差分格式：要善于用数值计算的手段要善于用数值计算的手段研究研究CFD ， 不能

23、仅限于不能仅限于用理论手段研究用理论手段研究CFD !基架点基架点 （stencil ）3. 复杂网格的处理方法复杂网格的处理方法1） 一维情况：一维情况： 非均匀网格非均匀网格 方法方法1 （常用）：（常用）： 网格（网格（Jacobian）变换）变换 j-2 j-1 j j+1 非均匀网格)(xx )1/()1(Nii0,1的均匀网格的均匀网格)(iixxdxdfxf 将方程由物理空间变到计算空间将方程由物理空间变到计算空间 （以（以x 为自变量变为以为自变量变为以 为自变量）为自变量）dxd为已知函数为已知函数)(xx 常用的一维坐标变换函数：常用的一维坐标变换函数： 指数函数指数函数

24、双曲正切函数双曲正切函数16Copyright by Li Xinliang)tanh(/ )tanh(gjgjbbx物理坐标物理坐标 计算坐标计算坐标 要求：要求： 坐标变换必须足够坐标变换必须足够光滑，否则会降低精度光滑，否则会降低精度网格间距变化要缓慢，否则网格间距变化要缓慢，否则会带来较大误差会带来较大误差方法方法2） 在非等距网格上直接构造差分格式在非等距网格上直接构造差分格式 j-2 j-1 j j+1 )(31431221Ouauauauaxujjjjj4131)3(21114131)3(21114232)3(2222)()(! 31)(! 21)()()(! 31)(! 21)

25、()()(! 31)(! 21)(jjjjjjjjjjjjjjjjjjjjjjjjjjjjjjjjjxxOxxuxxuxxuuuxxOxxuxxuxxuuuxxOxxuxxuxxuuu 原理：原理： 直接进行直接进行Taylor展开，构造格式展开，构造格式 格式系数是坐标（或网格间距）的函数格式系数是坐标（或网格间距）的函数0)(0)()(0)(0)()(1)(0)()(0)()()(4313323113242132221122411312111124013201102axxaaxxaxxaxxaaxxaxxaxxaaxxaxxaxxaaxxaxxjjjjjjjjjjjjjjjjjjjjjjj

26、j解出系数解出系数jjjjaaaa4321,注：注： 系数随网格点系数随网格点(j)变化！变化！17Copyright by Li Xinliang 网格非光滑、间距剧烈变化不会降低精度；网格非光滑、间距剧烈变化不会降低精度； 随机网格都可保证精度随机网格都可保证精度2） 二维二维/三维情况三维情况坐标变换坐标变换 均匀的直角网格均匀的直角网格RAE2822翼型周翼型周围的网格围的网格),(),(),(zzyyxx123123UfffVVVtxyzxyz ,TUuvw E 21, ()Tfuupuvuw u Ep),(),(zyxxxxfffxf1111三个方向共需计算三个方向共需计算9次导数

27、，次导数，计算量大计算量大yyyfffyf2222zzzfffzf3333对流项可组合，求对流项可组合，求3次导数即可次导数即可321321VVVffftUUJU1),(),(1zyxJ)(32111fffJfzyx)(32112fffJfzyx)(32113fffJfzyx18Copyright by Li Xinliang.0)()()(111xxxJJJ4. 时间项的离散时间项的离散1）直接离散法）直接离散法 把时间导数直接差分离散把时间导数直接差分离散 0 xuftu0)(1nxnnuftuu0)(11nxnnuftuu0)()(2111nxnxnnufuftuu1阶阶Euler显格式

28、显格式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推荐！推荐！ 在某一点进行在某一点进行Taylor展开，构造格式展开，构造格式3） 时时-空耦合离散空耦合离散0 xuctun+1nj-1 j j+1),(txu 蛙跳格式蛙跳格式 0 xuftu0221111xfftuunjnjnjnjn,j

29、Lax-Wandrof格式格式cucxtuuxtcuuxxnjnjnjnj2111212njnjnjnjnjffxtuuu112121212212121211njnjnjnjffxtuu半隐错点格式半隐错点格式022)(11111njnjnjnjnjnjuuuuxctuunjnjnjnjffxtuu11MacCormack格式格式1111121)(21njnjnjnjnjffxtuuu 1） 相容性：相容性： 当差分方程中当差分方程中 ，时间与空间步长均趋近于，时间与空间步长均趋近于0 时，差分方程的时，差分方程的截断误差截断误差也趋近于也趋近于0，则称差分方程与原微分方程是，则称差分方程与原

30、微分方程是相容相容的。的。 3.2 差分方法理论基础差分方法理论基础2）收敛性：）收敛性：0lim0,uuhtxL2 模：212)()(dxxuxu 模：)(max)(xuxuxLxuxujjh2)(jjhuxumax)(21 当时间与空间步长均趋近于当时间与空间步长均趋近于0 时，差分方程的时，差分方程的解解趋近于微分方程的解，趋近于微分方程的解，则称差分方程的解则称差分方程的解收敛收敛于原微分方程的解。于原微分方程的解。注意！注意！ 方程互相趋近方程互相趋近 解互相趋近解互相趋近（多值性、奇异性多值性、奇异性 ）)()(lim00 xfxfxx不一定等于不一定等于只有连续函数才满足只有连续函数才满足 （根据（根据Lax等价定理，只有稳定性条件满足的等价定理，只有稳定性条件满足的情况下，方程趋近才能保证解趋近）情况下，方程趋近才能保证解趋近）含义：含义： 方程趋近方程趋近含义：含义： 解趋近（更强）解趋近（更强）)(),(xuxuh分别为差分方程和微分方程的解1. 相容、收敛、稳定性与相容、收敛、稳定性与Lax等价定理等价定理相似的例子：3） 稳定性