第4讲-差分方法2

计算流体力学讲义

第四讲有限差分法（2）李新亮lixl@

；力学所主楼219；82543801知识点：离散误差的Fourier分析;

间断周围数值振荡的原因；GVC格式；

模型方程向N-S方程的推广；

1讲义、课件上传至

(流体中文网）->“流体论坛”->“CFD基础理论

”CopyrightbyLiXinliang§3.3差分格式的进一步分析1.耗散与色散误差2CopyrightbyLiXinliang精确解1阶迎风2阶迎风数值实验

时间推进：3步TVD型Runge-Kutta,且时间步长足够小（误差忽略）空间离散：1阶及2阶迎风格式（20个网格点）实验观察到的现象——两类误差：

振幅误差

相位误差（波速误差）

CopyrightbyLiXinliang3对以上“实验现象”进行理论分析半离散分析：假设时间推进是精确的，仅分析空间离散带来的误差（难度小、常用）全离散分析：同时分析时、空离散的误差（难度大）考查问题：实际上就是普通三角函数，采用复数形式仅仅是为了理论推导方便。用实数形式sin(kx),cos(kx)推导形式上略显繁琐。

精确解：差分格式：（1）

其他格式……假设对于：有隐含假设:线性差分格式

非线性系统作用于单波，会产生多个谐波（2）差分没有误差CopyrightbyLiXinliang4令：（1）式化为：“半离散化”：空间导数差分计算，时间方程（常微）精确计算如果，无误差分析（修正波数）与误差的关系

理想情况：

的误差导致解的幅值误差——耗散误差

的误差导致解传播速度的误差——色散误差假设对于：有

反映了一个波内的点数。PPW（波内的点数）=CopyrightbyLiXinliang5耗散、色散误差分别由修正波数的实部和虚部决定。关键参数：修正波数含义：反应波数（谱）空间内差分的误差任意函数：定义：求导数，精确解差分解Fourier分析的任务计算出，并考差其与的逼近程度。考察格式分辨率（resolution）的重要指标

精度：反映时的情况

分辨率：网格点数很少（例如波里面只有6个点）时的性能对于多尺度问题，分辨率更重要。牺牲精度，提高分辨率

优秀的差分格式，1个波长里面6个点即可精度分辨率CopyrightbyLiXinliang6如何计算修正波数？定义：

方法1.理论计算

根据差分具体表达式及定义计算例1：令则：于是：1阶迎风例2：2阶迎风CopyrightbyLiXinliang7方法2：数值计算定义：Step1）选取计算域[0,2p],计算网格（例如64，128）Step2）给定波数k,生成函数值Step3)调用差分子程序，得到导数值Step4)通过Fourier反变换，得到谱：假设已有求差分的子程序（黑箱，已知是线性的）线性黑箱强调：研究CFD本身，不能只使用理论手段，还要用数值手段根据修正波数的定义，有Step5)改变k的值，重复2-5，得到对于的依赖关系。画图非线性情况会产生高次谐波，造成step4中隐含的假设无法成立

将Fourier分析手段拓展到非线性系统需要研究的课题隐含条件：只有波数为k的那个谱不为0（线性系统）CopyrightbyLiXinliang8中心差分格式的色散特性0：精确解；1：4阶普通2：6阶普通；3：4阶紧致4：6阶紧致；5：6阶超紧致迎风差分格式的色散特性0：精确解，1：2阶迎风2：5阶迎风偏心3：3阶迎风紧致4：5阶迎风紧致每个波长里面2个网格点，谱方法的分辨率，差分法分辨率的极限（只有无穷阶精度才能达到）20阶超紧致格式——接近谱方法CopyrightbyLiXinliang9不同差分格式的色散误差曲线结论：要求分辨率相同的情况下，

采用高阶格式可放宽空间网格步长，从而减少计算量重要方向：高分辨率差分格式0：精确解1：2阶迎风2:3阶迎风3：3阶迎风紧致4：5阶迎风紧致指定误差要求的情况下，不同差分格式能模拟的最大a(a越大，所需网格越少）

作业题1：构造高分辨率差分格式，并进行理论分析及数值实验针对单波方程：对于空间导数，构造出一种不超过6点格式；并进行Fourier误差分析，画出kr,ki的曲线。

要求：精度不限；

网格基架点数不超过6个；

能够分辨的波数范围尽量宽；

（即kr,ki曲线近可能接近准确解）

给出差分的具体表达式，画出kr,ki的曲线；

说明构造格式的阶数，并采用本PPT第5页的方法给出的精度验证；

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

例：假设格式形式如下如果要求其有5阶精度，则通过Taylor展开可得到6个方程，6个系数可直接解出。我们要求其有4阶精度（当然3阶，2阶也可），于是Taylor展开只能提供5个方程。6个未知数（a1-a6），5个方程；有1个自由参数。调整这个自由参数，使得kr,ki曲线最为理想。

如何调整？1)可以人工调整，观察kr,ki曲线，选取满意的。2）可自动调整，设立一个优化目标函数。例如

调整自由参数，使得该目标函数取最大值。思路：牺牲精度，提高分辨率11CopyrightbyLiXinliang附录：部分差分格式…j-2j-1jj+1…表中的迎风差分格式均针对a>0当a<0时，需把下标的“j+k”换成“j-k”（例如把j+2换成j-2,把j-1换成j+1);并在表达式前加上“-”号。例：迎风偏斜格式：上游的基架点更多些（或上游权重更大）12CopyrightbyLiXinliang§3.4数值解的群速度及间断处数值振荡来源对于：有修正波数数值解色散误差：数值解传播的速度与精确解不一致数值解传播偏快数值解传播偏慢0：精确解；1：2阶迎风；2：5阶迎风偏心3：3阶迎风紧致；4：5阶迎风紧致

快格式(FST)：

慢格式(SLW):

混合格式(MXD):

特点：波数越高，误差越严重1.色散误差与群速度13CopyrightbyLiXinliangt=0.5时刻的精确解及数值解空间离散：五种不同格式；时间推进：3阶Runge-Kutta【数值实验】波的传播问题观察现象：1）高波数成分误差严重，低波数成分误差不明显；2）二阶Pade格式的解传播速度快于精确解，其余格式偏慢；3）迎风型格式有耗散，尤其是二阶迎风格式；

概念：群速度——波包传播的速度14CopyrightbyLiXinliang2.间断附近数值振荡的来源【数值实验】间断的传播计算域[0,1];计算网格点100

时间推进：3阶Runge-Kutta

空间离散：1）二阶中心差分

2阶迎风及2阶中心格式的色散特性2）二阶迎风差分15CopyrightbyLiXinliang过激波数值振荡的根源——色散误差导致群速度不一致快格式慢格式波前振荡波后振荡=+++…群速度控制的基本思路（群速度控制GVC:Fu&Ma）：

间断前、后分别采用快格式和慢格式，可有效抑制振荡Zhuang&Zhang:抑制波动原则

示意图：间断的Fourier分解好思路16CopyrightbyLiXinliang利用GVC的思想构造可计算间断的差分方法1）间断的前后判据

简易方法：则j点在间断左侧j-1jj+1则j点在间断右侧原理：越靠近间断，振荡越剧烈（a>0时，右侧为“前”）2）根据GVC的思想构造格式间断前：快格式；间断后：慢格式；格式GVC23）改写成为守恒型非线性情况，通常守恒型效果更好NND格式GVC2a17CopyrightbyLiXinliang

4）

a<0时，同样思路构造

（利用对称性，仅需把下标j+k换成j-k即可）

采用GVC2a(NND2a)格式的计算结果——消除振荡使用NND2a（守恒形式）；NND2（普通形式）及1阶迎风格式的计算结果将j换成j+118CopyrightbyLiXinliang作业题2：构造更高分辨率的GVC格式对于空间导数，构造出一种不超过6点的GVC格式。要求：a.精度不限；b.网格基架点数不超过6个；c.求解模型方程

计算结果间断尽量保持“锐利”;计算结果振荡尽量小。振荡的定量判据：总变差（TotalVariation）：

间断“锐利”的定量判据：间断区内的点数？（自行设计）给出差分格式的表达式、色散/耗散分析（ki,kr曲线）；给出模型方程t=0.2的结果（空间100个网格点，计算域[0,1]，时间推进可采用3阶Runge-Kutta方法）；与精确解及NND2a进行比较（画在同一张图上）建议：利用优化方法19CopyrightbyLiXinliang§3.5从模型方程推广到N-S方程（Euler方程）格式F+格式F-（教科书第6章）对流项：信息（波）从上游传至下游——上游信息更重要——迎风差分扩散项：信息从中心向周围扩散——不区分上、下游——中心差分迎风差分优点：有效利用信息传播的方向，增强稳定性微分与差分方程的影响域N-S方程：单波方程：单波方程——一个波，容易判断波传播方向N-S对流项（Euler）——方程组：多波问题，复杂双曲方程组的原则——特征分解，找到独立传播的波常系数矩阵A的情况——完全解耦，独立求解变系数矩阵A的情况——局部讨论20CopyrightbyLiXinliang1.Jacobian系数矩阵及其性质重要性质特点：A可以像常数一样，和求导运算交换21CopyrightbyLiXinliang2.对流项的分裂

目的：确定波传播方向，便于使用迎风差分方法：1）逐点分裂2）严格特征分裂

1）逐点分裂利用性质=+优点：耗散小缺点：导数间断方式A:特点：不必进行矩阵运算，计算量小Steger-Warming分裂A:Steger-Warming分裂22CopyrightbyLiXinliangSteger-Warming具体步骤（以一维为例）已知1）计算2）计算3）计算4）带入（1）式得到5）利用不同的迎风格式，分别计算

（1）（后差，前差）6）计算7）时间推进23CopyrightbyLiXinliang二维问题的steger-Warming分裂令：则：具体使用步骤，以计算为例令

计算特征值

分裂特征值，计算

带入左式，计算正、负流通矢量

计算计算设置，并注意对于曲线坐标系仅需令三维问题同样处理

二维、三维具体公式见傅德薰等《计算空气动力学》4.7节（158-162）书中公式有一定的排版错误，使用前务必重新仔细推导！24CopyrightbyLiXinliangB:L-F分裂特点：正特征值负特征值=+缺点：耗散偏大局部L-F分裂，每个点上计算全局L-F分裂，全局（一维）上计算足够大数学性质（光滑性）最好，但耗散偏大常数与迎风格式结合，等价于人工粘性例如，可取25CopyrightbyLiXinliang方式很多=+S-W:L-F:=+VanLeer:=+26CopyrightbyLiXinliang

分裂后失去了A的性质（可以像常数一样与求导交换）逐点分裂：

优点：无需矩阵运算，计算量小

缺点：分裂后改变了特征方向，耗散大利用了性质一般情况下：变系数，

不能与导数交换实质：没有做到解耦；只是把原变量重新组合，组合后波的传播方向的保证f+向正向传播，f-向负向传播

缺点：由于未解耦，各变量的误差会相互传递27CopyrightbyLiXinliang概念澄清：

流通矢量分裂本身不带来耗散，但其会影响到差分的耗散;举例：分裂过程耗散如果差分格式无耗散（例如都用中心差分），则通量分裂不带来耗散。=+向上平移向下平移分裂差分格式耗散分裂后的流场越偏离原先流场，则总体耗散越大精确满足，不引入误差！如使用低精差分度格式，则对分裂形式敏感（推荐使用特征分裂）如使用高精度格式（低耗散），则对分裂形式不敏感（可使用逐点分裂）28CopyrightbyLiXinliang2.严格特征分裂——基架点上冻结系数常系数方程组：完全解耦变系数情况——局部冻结系数…j-2j-1jj+1…在基架点上系数不变计算：在差分基架点上Aj

不变，可按常矩阵处理局部冻结系数分别采用后差和前差优点：严格保证（局部）特征方向，数值解质量好；缺点：大量矩阵运算，计算量大。29CopyrightbyLiXinliang通常写成守恒型差分，计算…j-2j-1jj+1…在基架点上系数不变具体步骤：

假设已知U，且针对模型方程（线性单波方程）

已构造出差分格式（1）1）计算出教材130页的公式(6.1.11-6.1.13),式中用到各变量在j+1/2的值（例如）

可使用j,j+1点值的算术平均（如）或Roe平均（教材6.4节）；由计算；方法很多，例如前面介绍的或30CopyrightbyLiXinliang

均可本人