(精品)第4讲-差分方法2.ppt

《(精品)第4讲-差分方法2.ppt》由会员分享，可在线阅读，更多相关《(精品)第4讲-差分方法2.ppt（32页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。

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

2、ang精确解1阶迎风2阶迎风数值实验数值实验 时间推进：时间推进：3步步TVD型型Runge-Kutta,且且时间步长足够小（误差忽略）时间步长足够小（误差忽略）空间离散：空间离散：1阶及阶及2阶迎风格式阶迎风格式（20个个网格点）网格点）实验观察到的现象实验观察到的现象 两类误差：两类误差：振幅误差振幅误差 相位误差相位误差（波速误差）（波速误差）Copyright by Li Xinliang3对以上对以上“实验现象实验现象”进行理论分析进行理论分析半离散分析：半离散分析：假设时间推进是精确的，仅假设时间推进是精确的，仅分析空间离散带来的误差（难度小、常用）分析空间离散带来的误差（难度小、

3、常用）全离散分析：全离散分析：同时分析时、空离散的误差同时分析时、空离散的误差（难度大）（难度大）考查问题：实际上就是普通三角函数，实际上就是普通三角函数，采用复采用复数形式仅仅是为了理论推导数形式仅仅是为了理论推导方便。方便。用实数形式用实数形式 sin(kx),cos(kx)推导形推导形式上略显繁琐。式上略显繁琐。精确解：差分格式：（1）其他格式 假设对于：假设对于：有有隐含假设隐含假设:线性差分格式线性差分格式 非线性系统作用于单波，会产生非线性系统作用于单波，会产生多个谐波多个谐波（2）差分没有误差差分没有误差Copyright by Li Xinliang4令：（1）式化为：“半离散

4、化半离散化”：空间导数差分计算，空间导数差分计算，时间方程（常微）精确计算时间方程（常微）精确计算如果如果 ，无误差无误差分析分析 （修正波数）与误差的关系（修正波数）与误差的关系 理想情况：理想情况：的误差导致解的幅值误差的误差导致解的幅值误差 耗散误差耗散误差 的误差导致解传播速度的误差的误差导致解传播速度的误差 色散误差色散误差假设对于：假设对于：有有 反映了一个波内的点数。反映了一个波内的点数。PPW （波内的点数）（波内的点数）=Copyright by Li Xinliang5耗散、色散误差分别由修正波数耗散、色散误差分别由修正波数的实部和虚部决定。的实部和虚部决定。关键参数：关键

5、参数：修正波数修正波数含义：含义：反应波数（谱）空间内差分的误差反应波数（谱）空间内差分的误差任意函数：定义：求导数，精确解求导数，精确解差分解差分解Fourier 分析的任务分析的任务计算出计算出 ，并考差其与，并考差其与 的逼近程度。的逼近程度。考察格式考察格式分辨率（分辨率（resolution）的重要指标的重要指标 精度：精度：反映反映 时的情况时的情况 分辨率：网格点数很少（例如波里面只有分辨率：网格点数很少（例如波里面只有6个点）时的性能个点）时的性能对于多尺度问题，分辨率更重要。对于多尺度问题，分辨率更重要。牺牲精度，提高分辨率牺牲精度，提高分辨率 优秀的差分格式，优秀的差分格式

6、，1个个波长里面波长里面6个点个点 即可即可精度精度 分辨率分辨率Copyright by Li Xinliang6如何计算修正波数？如何计算修正波数？定义：方法方法1.理论计算理论计算 根据差分具体表达式及定义计算根据差分具体表达式及定义计算例例1：令则：于是：1阶迎风阶迎风例例2：2阶迎风阶迎风Copyright by Li Xinliang7方法方法2：数值计算数值计算定义：Step 1）选取计算域）选取计算域0,2p p,计算网格（例如计算网格（例如64，128）Step 2）给定波数给定波数 k,生成函数值生成函数值Step 3)调用差分子程序，得到导数值调用差分子程序，得到导数值S

7、tep 4)通过通过Fourier反变换，得到谱：反变换，得到谱：假设已有求差分的子程序（黑箱，已知是线性的）假设已有求差分的子程序（黑箱，已知是线性的）线性黑箱强调：研究CFD本身，不能只使用理论手段，还要用数值手段根据修正波数的定义，根据修正波数的定义，有有Step 5)改变改变k的值，重复的值，重复2-5，得到得到 对于对于 的依赖关系。画图的依赖关系。画图非线性情况会产生高次谐波，造成非线性情况会产生高次谐波，造成 step 4 中隐含的假设无法成立中隐含的假设无法成立 将将Fourier分析手段拓展到非线性系统分析手段拓展到非线性系统 需要研究的课题需要研究的课题隐含条件：只有隐含条

8、件：只有波数为波数为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阶迎风紧致阶迎风紧致 每个波长里面每个波长里面2个网格点，个网格点，谱方法谱方法的分辨率，的分辨率，差分法分辨率的极限（只差分法分辨率的极限（只有无穷阶精度才能达到）有

9、无穷阶精度才能达到）20阶超紧阶超紧致格式致格式 接近谱接近谱方法方法Copyright by Li Xinliang9不同差分格式的色散误差曲线不同差分格式的色散误差曲线结论：结论：要求分辨率相同的情况下，要求分辨率相同的情况下，采用高阶格式可放宽空间网格步采用高阶格式可放宽空间网格步长，从而减少计算量长，从而减少计算量重要方向：重要方向：高分辨率差分格式高分辨率差分格式0：精确解精确解1：2阶迎风阶迎风2:3阶迎风阶迎风3：3阶迎风紧致阶迎风紧致4：5阶迎风紧致阶迎风紧致 指定误差要求的情况下，不同差分格指定误差要求的情况下，不同差分格式能模拟的最大式能模拟的最大a (a a (a 越大，

10、所需网越大，所需网格越少）格越少）作业题作业题1：构造高分辨率差分格式，并进行理论分析及数值实验构造高分辨率差分格式，并进行理论分析及数值实验 针对单波方程针对单波方程：对于空间导数，构造出一种不超过对于空间导数，构造出一种不超过6点格式；并进行点格式；并进行Fourier误差分析，误差分析，画出画出kr,ki的曲线。的曲线。要求：精度不限；要求：精度不限；网格基架点数不超过网格基架点数不超过6个；个；能够分辨的波数范围尽量宽；能够分辨的波数范围尽量宽；（即（即kr,ki曲线近可能接近准确解）曲线近可能接近准确解）给出差分的具体表达式，给出差分的具体表达式，画出画出kr,ki的曲线；的曲线；说

11、明构造格式的阶数，并采用本说明构造格式的阶数，并采用本PPT第第5页的方法给出的精度验证；页的方法给出的精度验证；形如：另外，进行如下数值验证：另外，进行如下数值验证：空间采用空间采用20个网格点，采用新构造的差分格式离散；时间推进采用个网格点，采用新构造的差分格式离散；时间推进采用3步步Runge-Kutta方法，时间步长可足够小（例如方法，时间步长可足够小（例如0.01）。给出）。给出t=20,50两个时两个时刻的数值解，与精确解比较（画图），并给出数值解的刻的数值解，与精确解比较（画图），并给出数值解的L2模误差。模误差。10Copyright by Li Xinliang提示：1.如不

12、使用优化技术，则格式构造方法简单，Taylor展开后解代数方程组即可。2.建议尝试使用优化技术 例：假设格式形式如下如果要求其有5阶精度，则通过Taylor展开可得到6个方程，6个系数可直接解出。我们要求其有4阶精度（当然3阶，2阶也可），于是Taylor展开只能提供5个方程。6个未知数（a1-a6），5个方程；有1个自由参数。调整这个自由参数，使得kr,ki曲线最为理想。如何调整？1)可以人工调整，观察kr,ki曲线，选取满意的。2）可自动调整，设立一个优化目标函数。例如 调整自由参数，使得该目标函数取最大值。思路：牺牲精度，提高分辨率 11Copyright by Li Xinliang附

13、录：附录：部分差分格式部分差分格式 j-2 j-1 j j+1 表中的迎风差分格式均针对表中的迎风差分格式均针对 a0当a0 时，右侧为“前”）2）根据根据GVC的思想构造格式的思想构造格式间断前：快格式；间断前：快格式；间断后：慢格式；间断后：慢格式；格式 GVC23）改写成为守恒型改写成为守恒型非线性情况，通常守非线性情况，通常守恒型效果更好恒型效果更好NND格式 GVC2a17Copyright by Li Xinliang 4）a 0 时，时，同样思路构造同样思路构造 （利用对称性，仅需把下标（利用对称性，仅需把下标j+k换成换成j-k即可）即可）采用GVC2a(NND 2a)格式的计

14、算结果 消除振荡使用NND2a（守恒形式）；NND2（普通形式）及1阶迎风格式的计算结果将j换成j+118Copyright by Li Xinliang作业题作业题2：构造更高分辨率的构造更高分辨率的GVC格式格式 对于空间导数，构造出一种不超过对于空间导数，构造出一种不超过6点的点的GVC格式。要求：格式。要求：a.精度不限；精度不限；b.网格基架点数不超过网格基架点数不超过6个；个；c.求解模型方程求解模型方程 计算结果间断尽量保持计算结果间断尽量保持“锐利锐利”;计算结果振荡尽量小。计算结果振荡尽量小。振荡的定量判据：振荡的定量判据：总变差（总变差（Total Variation）：）

15、：间断间断“锐利锐利”的定量判据：的定量判据：间断区内的点数？间断区内的点数？（自行设计）（自行设计）给出差分格式的表达式、色散给出差分格式的表达式、色散/耗散分析耗散分析（ki,kr曲线）；曲线）；给出模型方程给出模型方程t=0.2的结果（空间的结果（空间100个网格点，计算域个网格点，计算域0,1，时间推进可采用，时间推进可采用3阶阶Runge-Kutta方法）；与精确解及方法）；与精确解及NND2a进行比较（画在同一张图上）进行比较（画在同一张图上）建议：利用优化方法19Copyright by Li Xinliang 3.5 从模型方程推广到从模型方程推广到N-S方程（方程（Euler

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

17、传播的波特征分解，找到独立传播的波常系数矩阵常系数矩阵A的情况的情况 完全解耦，独立求解完全解耦，独立求解变系数矩阵变系数矩阵A的情况的情况 局部讨论局部讨论20Copyright by Li Xinliang1.Jacobian 系数矩阵及其性质系数矩阵及其性质重要性质特点：特点：A 可以像常数一样，和求导运算交换可以像常数一样，和求导运算交换21Copyright by Li Xinliang2.对流项的分裂对流项的分裂 目的：目的：确定波传播方向，便于使用迎风差分确定波传播方向，便于使用迎风差分方法：方法：1）逐点分裂逐点分裂 2）严格特征分裂严格特征分裂 1）逐点分裂逐点分裂 利用性质

18、=+优点：耗散小缺点：导数间断方式 A:特点：特点：不必进行矩阵运算，计算量小不必进行矩阵运算，计算量小Steger-Warming 分裂分裂A:Steger-Warming 分裂分裂22Copyright by Li XinliangSteger-Warming 具体步骤具体步骤（以一维为例）（以一维为例）已知已知1）计算计算2）计算计算3）计算计算4）带入（带入（1）式得到）式得到5）利用不同的迎风格式，分别计算利用不同的迎风格式，分别计算 （1）（后差，前差）（后差，前差）6）计算计算7）时间推进时间推进23Copyright by Li Xinliang二维问题的二维问题的steger

19、-Warming 分裂分裂令：则：具体使用步骤，具体使用步骤，以计算以计算 为例为例1)令令 2)计算特征值计算特征值3)分裂特征值，计算分裂特征值，计算4)带入左式，计算正、负流通矢带入左式，计算正、负流通矢量量5)计算计算计算计算 设置设置 ，并注意，并注意对于曲线坐标系对于曲线坐标系仅需令 三维问题同样处理三维问题同样处理 二维、三维具体二维、三维具体 公式见傅德薰等公式见傅德薰等计算空气动力学计算空气动力学 4.7节节（158-162）书中公式有一定的排版错误，使用前务必书中公式有一定的排版错误，使用前务必重新仔细重新仔细推导！推导！24Copyright by Li Xinliang

20、B:L-F分裂分裂特点：正特征值 负特征值=+缺点：耗散偏大局部局部L-F分裂分裂，每个点上计算，每个点上计算 全局全局L-F分裂分裂，全局（一维）上计算，全局（一维）上计算 足够大数学性质（光滑性）数学性质（光滑性）最好，但耗散偏大最好，但耗散偏大常数与迎风格式结合，等价于人工粘性与迎风格式结合，等价于人工粘性例如，可取例如，可取25Copyright by Li Xinliang方式很多=+S-W:L-F:=+Van Leer:=+26Copyright by Li Xinliang 分裂后分裂后 失去了失去了A的性质（可以的性质（可以像常数一样与求导交换）像常数一样与求导交换）逐点分裂：

21、逐点分裂：优点：优点：无需矩阵运算，计算量小无需矩阵运算，计算量小 缺点：缺点：分裂后改变了特征方向，分裂后改变了特征方向，耗散大耗散大利用了性质一般情况下：一般情况下：变系数，变系数，不能与导数交换不能与导数交换实质：实质：没有做到解耦；没有做到解耦；只是把原变量重新组合，组合后波的传播方向的保证只是把原变量重新组合，组合后波的传播方向的保证 f+向正向正向传播，向传播，f-向负向传播向负向传播 缺点：缺点：由于未解耦，各变量的误差会相互传递由于未解耦，各变量的误差会相互传递 27Copyright by Li Xinliang概念澄清：概念澄清：流通流通矢量分裂本身不带来耗散矢量分裂本身不

22、带来耗散，但其会但其会影响到差分的耗散影响到差分的耗散;举例：分裂过程耗散如果差分格式无耗散（例如都用中心差分），则通量分裂不带来耗散。=+向上平移向下平移分裂差分格式耗散分裂后的流场越偏离原先流场，则分裂后的流场越偏离原先流场，则总体耗散越大总体耗散越大精确满足，不引入误差！如使用低精差分度格式，如使用低精差分度格式，则对分裂形式敏感则对分裂形式敏感 （推荐使用特征分裂）（推荐使用特征分裂）如使用高精度格式（低耗散），则对分裂形式不敏感如使用高精度格式（低耗散），则对分裂形式不敏感（可使用逐点分裂）（可使用逐点分裂）28Copyright by Li Xinliang2.严格特征分裂严格特征

23、分裂 基架点上冻结系数基架点上冻结系数常系数常系数方程组：方程组：完全解耦变系数情况变系数情况 局部冻结系数局部冻结系数 j-2 j-1 j j+1 在基架点上系数 不变计算：在差分基架点上在差分基架点上Aj 不变，不变，可按常矩阵处理可按常矩阵处理局部冻结系数分别采用后差和前差优点：严格保证（局部）特征方向，数值解质量好；缺点：大量矩阵运算，计算量大。29Copyright by Li Xinliang通常写成守恒型差分通常写成守恒型差分，计算计算 j-2 j-1 j j+1 在基架点上系数 不变具体步骤：具体步骤：假设已知假设已知 U，且针对模型方程（线性单波方程）且针对模型方程（线性单波

24、方程）已构造出差分格式已构造出差分格式（1）1）计算出计算出教材教材130页的公式页的公式(6.1.11-6.1.13),式中用到各变量在式中用到各变量在j+1/2的值（例如的值（例如 ）可使用可使用j,j+1 点值的算术平均点值的算术平均（如（如 ）或或Roe平均平均（教材（教材6.4节）；节）；由由 计算；方法很多，例如前面介绍的计算；方法很多，例如前面介绍的 或或 30Copyright by Li Xinliang 均可本人感觉会限制精度（例如本人感觉会限制精度（例如2阶）阶）；但；但数值实验没发现问题。数值实验没发现问题。张树海张树海 JCP 2009 对其进行了深入探讨对其进行了深

25、入探讨；欢迎；欢迎对该问题进行研究对该问题进行研究。2）在网格基上计算在网格基上计算 j-2 j-1 j j+1 计算fj+1/2用到的点注意，在该网格基上（例如注意，在该网格基上（例如k=j-1,j,j+1）保持不变保持不变例如：3)利用已构造好的差分格式，计利用已构造好的差分格式，计算通量算通量4）得到总通量得到总通量 5）计算差分计算差分（j点处）点处）步骤的算法描述步骤的算法描述 （注意：（注意：实际上是两重循环）实际上是两重循环）do j=1,N do k=j-1,j+1 (网格基，可以是更多或更少点)enddo enddo do j=1,N enddo 需要多次矩阵运算，需要多次矩

26、阵运算，计算量大计算量大 守恒性好，耗散小，守恒性好，耗散小，数值解质量好数值解质量好31Copyright by Li Xinliang作业题作业题3：针对如下针对如下Sod 激波管问题激波管问题 计算其数值解，画出t=0.14时刻密度、速度及压力的分布；并与精确解进行比较（要求画在一张图上）。要求：1）空间网格数100，时间推进格式选用3阶Runge-Kutta,时间步长自选。2）可选用逐点分裂，也可选用特征分裂。3）建议采用本讲作业题2（或作业题1）自行构造的差分格式计算。（作业题2是激波捕捉格式，效果应当会好些）。如果作业题1和作业题2遇到困难，也可采用现有的差分格式。32Copyright by Li Xinliang