第4讲-差分方法2

1、计算流体力学讲义计算流体力学讲义 第四讲第四讲 有限差分法（有限差分法（2）李新亮李新亮 ；力学所主楼；力学所主楼219； 82543801 知识点：知识点： 离散误差的离散误差的Fourier分析分析; 间断周围数值振荡的原因；间断周围数值振荡的原因； GVC格式；格式； 模型方程向模型方程向N-S方程的推广；方程的推广； 1讲义、课件上传至讲义、课件上传至 (流体中文网）流体中文网） - “流体论坛流体论坛” -“ CFD基础理论基础理论 ”Copyright by Li Xinliang 3.3 差分格式的进一步分析差分格式的进一步分析1. 耗散与色散误差耗散与色散误差2Copyrigh

2、t by Li Xinliang)sin()0 ,(0 xxuxutu精确解1阶迎风xu0123456-1-0.500.51Exact1st upwind2nd upwindCopyright by Li Xinliangt=solutions of 1D convection Eq.2阶迎风)2/()43(21xuuuujjjjxuuujjj/ )(1数值实验数值实验 boundaryperiodicxxxuxutu,2 , 0)sin()0 ,(0时间推进：时间推进： 3步步TVD型型Runge-Kutta, 且且时间步长足够小（误差忽略）时间步长足够小（误差忽略）空间离散：空间离散： 1

3、阶及阶及2阶迎风格式阶迎风格式 （20个网格点）个网格点）实验观察到的现象实验观察到的现象 两类误差：两类误差： 振幅误差振幅误差 相位误差相位误差 （波速误差）（波速误差） Copyright by Li Xinliang3对以上对以上“实验现象实验现象”进行理论分析进行理论分析半离散分析：半离散分析： 假设时间推进是精确的，仅假设时间推进是精确的，仅分析空间离散带来的误差（难度小、常用）分析空间离散带来的误差（难度小、常用）全离散分析：全离散分析： 同时分析时、空离散的误差同时分析时、空离散的误差 （难度大）（难度大）boundaryperiodicxexfcxfctfikx,2 , 0)

4、0 ,()0(0考查问题：实际上就是普通三角函数，实际上就是普通三角函数，采用复采用复数形式仅仅是为了理论推导数形式仅仅是为了理论推导方便。方便。用实数形式用实数形式 sin(kx), cos(kx)推导形推导形式上略显繁琐。式上略显繁琐。 精确解：ikxikctctxikeeetxf)(),(差分格式：0jxjfctf（1）xfffjjjx/ )(1)2/()43(21xuuufjjjjx 其他格式 jikxjexu)(假设对于：假设对于：有有jikxjxexku隐含假设隐含假设: 线性差分格式线性差分格式 非线性系统作用于单波，会产生非线性系统作用于单波，会产生多个谐波多个谐波 （2） 差

5、分没有误差差分没有误差xikkCopyright by Li Xinliang4令：jikxjetutxf)( ),(jikxjxxektutxf/)( ),(（1）式化为：0/)( )( jjikxikxxektucedttud0/)( )( xktucdttud“半离散化半离散化”： 空间导数差分计算，空间导数差分计算，时间方程（常微）精确计算时间方程（常微）精确计算xctkeutu/)0( )( xctkikxikxjjjeuetutxf/)0( )( ),(ixikk如果如果 ， 无误差无误差分析分析 （修正波数）与误差的关系（修正波数）与误差的关系kirikkk)/(/)0( ),(

6、xkctkxikxctkjijreeutxf 理想情况：理想情况： 的误差导致解的幅值误差的误差导致解的幅值误差 耗散误差耗散误差 的误差导致解传播速度的误差的误差导致解传播速度的误差 色散误差色散误差xkkkir , 0rkikjikxjexu)(假设对于：假设对于：有有jikxjxexku 反映了一个波内的点数。反映了一个波内的点数。PPW （波内的点数）（波内的点数）= Copyright by Li Xinliang5jikxjexu)(jikxjxexku耗散、色散误差分别由修正波数耗散、色散误差分别由修正波数的实部和虚部决定。的实部和虚部决定。k关键参数：关键参数： 修正波数修正波

7、数含义：含义： 反应波数（谱）空间内差分的误差反应波数（谱）空间内差分的误差jikxkkjefxf)(任意函数：jikxkkjefikxf)(定义：kkfikf求导数，精确解求导数，精确解差分解差分解jikxkkjxjef xk ifF/kkfxkFFourier 分析的任务分析的任务计算出计算出 ，并考差其与，并考差其与 的逼近程度。的逼近程度。kixik 考察格式考察格式分辨率（分辨率（resolution）的重要指标的重要指标 精度：精度： 反映反映 时的情况时的情况 分辨率：网格点数很少（例如波里面只有分辨率：网格点数很少（例如波里面只有6个点）时的性能个点）时的性能对于多尺度问题，分

8、辨率更重要。对于多尺度问题，分辨率更重要。 牺牲精度，提高分辨率牺牲精度，提高分辨率0 xxk 优秀的差分格式，优秀的差分格式，1个个波长里面波长里面6个点个点 即可即可/2精度精度 分辨率分辨率Copyright by Li Xinliang6如何计算修正波数？如何计算修正波数？jikxjexu)(jikxjxexku定义： 方法方法1. 理论计算理论计算 根据差分具体表达式及定义计算根据差分具体表达式及定义计算例例1：xuuujjjx1令jikxjeu 则：jjjjikxiikxxxikikxjjjxexkexeeexxuuu)1 ()(1)(1xk于是：sin)cos1 (1ieki1阶

9、迎风阶迎风sin,cos1irkk例例2：)2/()43(21xuuuujjjjx2阶迎风阶迎风)43(2243221iiikxjjjjxeexexuuuuj2/ )43(2iieek2/ )2sinsin4(,2/ )2coscos43(irkkCopyright by Li Xinliang7方法方法2： 数值计算数值计算jikxjexu)(jikxjxexku定义：Step 1）选取计算域）选取计算域0,2 , 计算网格（例如计算网格（例如64，128）Step 2）给定波数给定波数 k, 生成函数值生成函数值Step 3) 调用差分子程序，得到导数值调用差分子程序，得到导数值Step

10、4) 通过通过Fourier反变换，得到谱：反变换，得到谱：假设已有求差分的子程序（黑箱，已知是线性的）假设已有求差分的子程序（黑箱，已知是线性的）x线性黑箱ju强调：研究CFD本身，不能只使用理论手段，还要用数值手段jxjuv).2 , 1(Njeujikxj).2 , 1(NjvjNjikxjkjevNv11根据修正波数的定义，根据修正波数的定义， 有有xvkk Step 5) 改变改变k的值，重复的值，重复2-5， 得到得到 对于对于 的依赖关系。画图的依赖关系。画图kxk非线性情况会产生高次谐波，造成非线性情况会产生高次谐波，造成 step 4 中隐含的假设无法成立中隐含的假设无法成立

11、 将将Fourier分析手段拓展到非线性系统分析手段拓展到非线性系统 需要研究的课题需要研究的课题隐含条件：只有隐含条件：只有波数为波数为k的那个的那个谱不为谱不为0 （线性（线性系统）系统）Copyright by Li Xinliang8中心差分格式的色散特性中心差分格式的色散特性0： 精确解；精确解； 1 ： 4阶普通阶普通2： 6阶普通；阶普通； 3：4阶紧致阶紧致4： 6阶紧致；阶紧致； 5：6阶超紧致阶超紧致迎风差分格式的色散特性迎风差分格式的色散特性0： 精确解，精确解， 1： 2阶迎风阶迎风2： 5阶迎风偏心阶迎风偏心 3： 3阶迎风紧致阶迎风紧致4： 5阶迎风紧致阶迎风紧致

12、每个波长里面每个波长里面2个网格点，个网格点， 谱方法谱方法的分辨率，的分辨率， 差分法分辨率的极限（只差分法分辨率的极限（只有无穷阶精度才能达到）有无穷阶精度才能达到）2)int(WavePersPoPPW20阶超紧阶超紧致格式致格式 接近谱接近谱方法方法Copyright by Li Xinliang9不同差分格式的色散误差曲线不同差分格式的色散误差曲线结论：结论： 要求分辨率相同的情况下，要求分辨率相同的情况下， 采用高阶格式可放宽空间网格步采用高阶格式可放宽空间网格步长，从而减少计算量长，从而减少计算量重要方向：重要方向： 高分辨率差分格式高分辨率差分格式0： 精确解精确解1： 2阶迎

13、风阶迎风2: 3阶迎风阶迎风3： 3阶迎风紧致阶迎风紧致4： 5阶迎风紧致阶迎风紧致 指定误差要求的情况下，不同差分格指定误差要求的情况下，不同差分格式能模拟的最大式能模拟的最大 ( ( 越大，所需网越大，所需网格越少）格越少） /2PPW作业题作业题1： 构造高分辨率差分格式，并进行理论分析及数值实验构造高分辨率差分格式，并进行理论分析及数值实验 针对单波方程针对单波方程：0 xutu 对于空间导数，构造出一种不超过对于空间导数，构造出一种不超过6点格式；并进行点格式；并进行Fourier误差分析，误差分析，画出画出kr,ki的曲线。的曲线。 要求：精度不限；要求：精度不限； 网格基架点数不

14、超过网格基架点数不超过6个；个； 能够分辨的波数范围尽量宽；能够分辨的波数范围尽量宽； （即（即kr,ki曲线近可能接近准确解）曲线近可能接近准确解） 给出差分的具体表达式，给出差分的具体表达式， 画出画出kr,ki的曲线；的曲线； 说明构造格式的阶数，并采用本说明构造格式的阶数，并采用本PPT第第5页的方法给出的精度验证；页的方法给出的精度验证； 26154131231jjjjjjjjuauauauauauaxuu16514233241jjjjjjjjuauauauauauaxuu形如：另外，进行如下数值验证：另外，进行如下数值验证：)sin()0 ,(2 , 0, 0 xxuxxutu空间

15、采用空间采用20个网格点，采用新构造的差分格式离散；时间推进采用个网格点，采用新构造的差分格式离散；时间推进采用3步步Runge-Kutta方法，时间步长可足够小（例如方法，时间步长可足够小（例如0.01）。给出）。给出t=20,50两个时两个时刻的数值解，与精确解比较（画图），并给出数值解的刻的数值解，与精确解比较（画图），并给出数值解的L2模误差。模误差。10Copyright by Li Xinliang提示： 1. 如不使用优化技术，则格式构造方法简单， Taylor展开后解代数方程组即可。 2. 建议尝试使用优化技术 26154131231jjjjjjjjuauauauauauaxu

16、u例： 假设格式形式如下如果要求其有5阶精度，则通过Taylor展开可得到6个方程，6个系数可直接解出。我们要求其有4阶精度（当然3阶，2阶也可），于是Taylor展开只能提供5个方程。6个未知数（a1-a6）， 5个方程； 有1个自由参数。 调整这个自由参数，使得kr,ki曲线最为理想。 如何调整？ 1) 可以人工调整，观察kr,ki曲线，选取满意的。 2）可自动调整，设立一个优化目标函数。 例如 调整自由参数，使得该目标函数取最大值。思路：牺牲精度，提高分辨率 05. 0)(:*ik11Copyright by Li Xinliang附录：附录： 部分差分格式部分差分格式 j-2 j-1

17、j j+1 0 xuatu表中的迎风差分格式均针对表中的迎风差分格式均针对 a0当a0 时， 右侧为“前”） 2） 根据根据GVC的思想构造格式的思想构造格式jjjjjjjjjjjjjjxuuuuwhenxuuuuuuwhenxuuuu11111121)2/()()2/()43(间断前：快格式；间断前：快格式；间断后：慢格式；间断后：慢格式；格式 GVC23） 改写成为守恒型改写成为守恒型非线性情况，通常守非线性情况，通常守恒型效果更好恒型效果更好jjjjjjjjjjjjjuuuuwhenuuuuuuwhenuuF1111112/12/ )(2/ )3(xFFujjjx/ )(2/12/1NN

18、D格式 GVC2a17Copyright by Li Xinliang 4） a 0 时，时， 同样思路构造同样思路构造 （利用对称性，仅需把下标（利用对称性，仅需把下标j+k换成换成j-k即可）即可） 00axuatujjjjjjjjjjjjjuuuuwhenuuuuuuwhenuuF1111112/12/ )(2/ )3(采用GVC2a (NND 2a)格式的计算结果 消除振荡xu00.8100.51NND 2aNND 2Exact solution1st upwindDeveloping of initial discontinuity2nd NND(C) Li Xin

19、liangt=0.2使用NND2a （守恒形式）；NND2 （普通形式）及1阶迎风格式的计算结果1211121212/12/ )(2/ )3(jjjjjjjjjjjjjuuuuwhenuuuuuuwhenuuF将j换成j+118Copyright by Li Xinliang作业题作业题2： 构造更高分辨率的构造更高分辨率的GVC格式格式 对于空间导数，构造出一种不超过对于空间导数，构造出一种不超过6点的点的GVC格式。要求：格式。要求： a. 精度不限；精度不限； b. 网格基架点数不超过网格基架点数不超过6个；个； c. 求解模型方程求解模型方程 1,0axuatu5 . 015 . 00

20、)0 ,(xxxu计算结果间断尽量保持计算结果间断尽量保持“锐利锐利”;计算结果振荡尽量小。计算结果振荡尽量小。振荡的定量判据：振荡的定量判据： 总变差（总变差（Total Variation）：）： 间断间断“锐利锐利”的定量判据：的定量判据： 间断区内的点数？间断区内的点数？ （自行设计）（自行设计）111NjjjuuTV给出差分格式的表达式、色散给出差分格式的表达式、色散/耗散分析耗散分析 （ki,kr曲线）；曲线）；给出模型方程给出模型方程t=0.2的结果（空间的结果（空间100个网格点，计算域个网格点，计算域0,1，时间推进可采用，时间推进可采用3阶阶Runge-Kutta方法）；与

21、精确解及方法）；与精确解及NND2a进行比较（画在同一张图上）进行比较（画在同一张图上）xu00.8100.51NND 2aNND 2Exact solution1st upwindDeveloping of initial discontinuity2nd NND(C) Li Xinliangt=0.2建议： 利用优化方法19Copyright by Li Xinliang 3.5 从模型方程推广到从模型方程推广到N-S方程（方程（Euler方程）方程）0 xuatu?0?,0aa格式F+格式F-visf(U)UxtTTpEupuufEuU)(,(,),(2xf(U)xua

22、（教科书第6章）对流项：信息（波）从上游传至下游对流项：信息（波）从上游传至下游 上游信息更重要上游信息更重要 迎风差分迎风差分扩散项：扩散项： 信息从中心向周围扩散信息从中心向周围扩散 不区分上、下游不区分上、下游 中心差分中心差分迎风差分优点：迎风差分优点： 有效利用信息传播的方向，增强稳定性有效利用信息传播的方向，增强稳定性微分与差分方程的影响域N-S方程：单波方程：单波方程单波方程 一个波，容易判断波传播方向一个波，容易判断波传播方向N-S对流项（对流项（Euler） 方程组：多波问题，方程组：多波问题， 复杂复杂双曲方程组的原则双曲方程组的原则 特征分解，找到独立传播的波特征分解，找

23、到独立传播的波0000 xvtvxtxtxtkkkVVUSSUUAU1常系数矩阵常系数矩阵A的情况的情况 完全解耦，独立求解完全解耦，独立求解变系数矩阵变系数矩阵A的情况的情况 局部讨论局部讨论20Copyright by Li Xinliang1. Jacobian 系数矩阵及其性质系数矩阵及其性质xxUAf(U)Uf(u)Af(U)Uf)(f(U)UUU)f()()(1AUUUf(U)f(U)(xxxUAAUf(U)(重要性质uucucuuu222322311221)3(23010AtttUAAUf(U)(特点：特点： A 可以像常数一样，和求导运算交换可以像常数一样，和求导运算交换vis

24、f(U)UxtTTpEupuufEuU)(,(,),(221Copyright by Li Xinliang2. 对流项的分裂对流项的分裂 目的：目的： 确定波传播方向，便于使用迎风差分确定波传播方向，便于使用迎风差分方法：方法： 1） 逐点分裂逐点分裂 2） 严格特征分裂严格特征分裂 1） 逐点分裂逐点分裂 利用性质)(AUUAf(U)xxxSSA1kdiagkdiag2,2kk=+优点：耗散小缺点：导数间断SSA1SUSUAf1方式 A:2)(2/122kkk) 1(2)(3(232cw特点：特点： 不必进行矩阵运算，计算量小不必进行矩阵运算，计算量小Steger-Warming 分裂分裂

25、wcucuucucuu232221321321)(2)(2) 1()()() 1(2) 1(22)(f)(fffffA: Steger-Warming 分裂分裂22Copyright by Li Xinliang Steger-Warming 具体步骤具体步骤 （以一维为例）（以一维为例）已知已知TEuU),(1） 计算计算2） 计算计算3） 计算计算4） 带入（带入（1）式得到）式得到5） 利用不同的迎风格式，分别计算利用不同的迎风格式，分别计算 pu,cucuu321,2)(2/122kkk)3 , 2 , 1(,kkkwcucuucucuu232221321321)(2)(2) 1()(

26、)() 1(2) 1(22)(f) 1(2)(3(232cw)(),(ffff（1）ff ,fffxxff,xxff0a0 xuatu0a（后差，前差）（后差，前差） 6） xxxfff计算计算7） 时间推进时间推进0 xtfU23Copyright by Li Xinliang二维问题的二维问题的steger-Warming 分裂分裂021yxtffUTpEuuvpuu)(,21fTpEvpvuvv)(,22fTEvu,U21fff令：则：WvuvuVvvvuuu2222)(2222421213202413024130430f) 1(2)(3(243cW10) 1(211kcuu12kcuu

27、/1k/2k21kcvv22kcvv22具体使用步骤，具体使用步骤， 以计算以计算 为例为例x1f1) 令令 2) 计算特征值计算特征值3) 分裂特征值，计算分裂特征值，计算4) 带入左式，计算正、负流通矢带入左式，计算正、负流通矢量量5) 计算计算0, 1cucuu4321,2)(2/122kkk)4 , 3 , 2 , 1(,kkk)(1ffxxx111fffy2f计算计算 设置设置 ，并注意，并注意1, 0cvcvv4321,对于曲线坐标系对于曲线坐标系021ffUt21111fffyxJJ仅需令yxJJ11, 三维问题同样处理三维问题同样处理 二维、三维具体二维、三维具体 公式见傅德薰

28、等公式见傅德薰等计算空气动力学计算空气动力学 4.7节节 （158-162）书中公式有一定的排版错误，使用前务必书中公式有一定的排版错误，使用前务必重新仔细重新仔细推导！推导！24Copyright by Li XinliangB: L-F分裂分裂AUf(U) 2/ )(, 2/ )(*IAAIAA特点： A正特征值 负特征值A2/ )(*UffUA=+缺点：耗散偏大局部局部L-F分裂分裂，每个点上计算，每个点上计算 全局全局L-F分裂分裂，全局（一维）上计算，全局（一维）上计算 cu *足够大cu *数学性质（光滑性）数学性质（光滑性）最好，但耗散偏大最好，但耗散偏大)(max*cux常数与

29、迎风格式结合，等价于人工粘性与迎风格式结合，等价于人工粘性UfUfUfUfUfUffffxxxxxxxxxxxxxx*0*)(2)(21)()(21)()(21例如，可取例如，可取0 xtf(U)U22*xxtUf(U)U25Copyright by Li Xinliang2,2kkkkkkkkk方式很多=+S-W:L-F:=+2,2*kkkk0, 0kkVan Leer:=+221112) 1() 1(22) 1(cufcufff2121Macfff0f1Ma1MaffcuMa/0f1Ma26Copyright by Li Xinliang 分裂后分裂后 失去了失去了A的性质（可以的性质（可

30、以像常数一样与求导交换）像常数一样与求导交换）逐点分裂：逐点分裂： 优点：优点： 无需矩阵运算，计算量小无需矩阵运算，计算量小 缺点：缺点： 分裂后改变了特征方向，分裂后改变了特征方向， 耗散大耗散大xxx(U)f(U)ff(U)SUSUAf1利用了性质)(AUUAf(U)xxx一般情况下：一般情况下：xxxUAUAf变系数，变系数， 不能与导数交换不能与导数交换AAxBxUf AB实质：实质： 没有做到解耦；没有做到解耦； 只是把原变量重新组合，组合后波的传播方向的保证只是把原变量重新组合，组合后波的传播方向的保证 f+ 向正向正向传播，向传播，f-向负向传播向负向传播 缺点：缺点： 由于未

31、解耦，各变量的误差会相互传递由于未解耦，各变量的误差会相互传递 27Copyright by Li Xinliang概念澄清：概念澄清： 流通流通矢量分裂本身不带来耗散矢量分裂本身不带来耗散， 但其会但其会影响到差分的耗散影响到差分的耗散;举例：0)(xUftUfff分裂过程),(21),(21UffUfffffxxxUaftUxxx00耗散如果差分格式无耗散（例如都用中心差分），则通量分裂不带来耗散。ffffffxxxxx0000)(=+向上平移向下平移分裂差分格式耗散分裂后的流场越偏离原先流场，则分裂后的流场越偏离原先流场，则总体耗散越大总体耗散越大fff精确满足，不引入误差！如使用低精差

32、分度格式，如使用低精差分度格式， 则对分裂形式敏感则对分裂形式敏感 （推荐使用特征分裂）（推荐使用特征分裂）如使用高精度格式（低耗散），则对分裂形式不敏感如使用高精度格式（低耗散），则对分裂形式不敏感 （可使用逐点分裂）（可使用逐点分裂）28Copyright by Li Xinliang2. 严格特征分裂严格特征分裂 基架点上冻结系数基架点上冻结系数常系数常系数方程组：方程组：00 xtxtUAUf(U)U0000 xvtvxtxtxtkkkVVUSSUUAU1完全解耦变系数情况变系数情况 局部冻结系数局部冻结系数0jjjxtUAU j-2 j-1 j j+1 在基架点上系数 不变jjxUA

33、jjxUAjA计算：在差分基架点上在差分基架点上Aj 不变，不变， 可按常矩阵处理可按常矩阵处理jj1jSSAj)()(jjjjjjjjxxxxxxVVSVSUSSUSSUA1j1j1jj1j局部冻结系数分别采用后差和前差USVj优点： 严格保证（局部）特征方向，数值解质量好；缺点： 大量矩阵运算，计算量大。29Copyright by Li XinliangUSVjUSVj通常写成守恒型差分通常写成守恒型差分，计算计算x)(Ufxxjjj/ )()(2/12/1ffUf j-2 j-1 j j+1 在基架点上系数 不变jjxUAjA)(2/12/112/12/1jjjjSVVfUSV2/1j

34、具体步骤：具体步骤： 假设已知假设已知 U， 且针对模型方程（线性单波方程）且针对模型方程（线性单波方程） 已构造出差分格式已构造出差分格式xvvvxvvvjjjjjj/ )(;/ )(2/12/12/12/10 xvatv（1）1） 计算出计算出2/12/112/1,jjjSS教材教材130页的公式页的公式(6.1.11-6.1.13), 式中用到各变量在式中用到各变量在j+1/2的值（例如的值（例如 ） 可使用可使用j, j+1 点值的算术平均点值的算术平均 （如（如 ） 或或Roe平均平均 （教材（教材6.4节）；节）； 由由 计算；方法很多，例如前面介绍的计算；方法很多，例如前面介绍的 或或 2/1ju2/ )(12/1jjjuuu30Copyright by Li