《(精品)第8讲-差分方法4.ppt》由会员分享，可在线阅读，更多相关《(精品)第8讲-差分方法4.ppt（29页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、计算流体力学讲义计算流体力学讲义 第八讲第八讲 差分方法（差分方法（4）李新亮李新亮；力学所主楼；力学所主楼219；82543801 知识点：知识点：高精度激波捕捉格式高精度激波捕捉格式WENOWENO 激波捕捉方法激波捕捉方法 Godnov,MUSCL Godnov,MUSCL 常用的隐式处理方法常用的隐式处理方法LU-SGSLU-SGS 1Copyright by Li Xinliang讲义、课件上传至讲义、课件上传至 (流体中文网）流体中文网）-“流体论坛流体论坛”-“CFD基础理基础理”讲课录像及课件上传至网盘讲课录像及课件上传至网盘 http:/cid- by Li Xinliang
2、2知识回顾知识回顾1.Roe 格式格式 准确特征方向的守恒型格式准确特征方向的守恒型格式标量方程的情况：标量方程的情况：等价于等价于标准的标准的1阶迎风差分格式阶迎风差分格式f(u)的平均增长率平均平均斜率斜率 双曲守恒方程组的双曲守恒方程组的Roe 格式格式平均增长率矩阵平均增长率矩阵平均斜率平均斜率Roe 平均平均含义：含义：点点 处的瞬时增长率刚好为区间处的瞬时增长率刚好为区间 上的平均增长率上的平均增长率Copyright by Li Xinliang42.TVD 格式格式概念：概念：网格网格Reynolds数数 单调格式、保单调格式及单调格式、保单调格式及TVD格式格式 Harten
3、定理：定理：正系数原则正系数原则TVD保单调保单调单调单调TVD格式=1阶迎风阶迎风+j j*（LW格式格式-1阶迎风）阶迎风）二阶精度区二阶精度区TVD区区二阶精度二阶精度TVD区（二区（二者交集）者交集）Copyright by Li Xinliang5 8.1 WENO格式格式 高精度的激波捕捉法高精度的激波捕捉法 1.基本思路基本思路 j-3,j-2,j-1,j,j+1,j+2j-3,j-2,j-1,j；j-2,j-1,j,j+1；j-1,j,j+1,j+2五个基架点被分成三个组五个基架点被分成三个组1）若高精度逼近若高精度逼近 ，必然利用多个基架点必然利用多个基架点 2）如果该基架点
4、内函数有间断，会导致振荡如果该基架点内函数有间断，会导致振荡 3）间断不可能处处存在间断不可能处处存在 4）把基架点分成多个组（模板），把基架点分成多个组（模板），每个模板每个模板独立独立计算计算j点导数的逼近。点导数的逼近。得到得到多个差分多个差分 5）根据每个模板的光滑程度，设定权重）根据每个模板的光滑程度，设定权重 6）对多个差分结果进行加权平均对多个差分结果进行加权平均。光滑度越高，。光滑度越高，权重越大。如果某模板存在间断，则权重趋于权重越大。如果某模板存在间断，则权重趋于0；如果都光滑，则组合成更高阶格式。如果都光滑，则组合成更高阶格式。Copyright by Li Xinlia
5、ng62.WENO格式的原理描述格式的原理描述考虑线性单波方程：考虑线性单波方程：注：注：为了简便，以非守恒型形式为例讲授其思路，实际使用时，请采为了简便，以非守恒型形式为例讲授其思路，实际使用时，请采用下一节介绍的守恒形式用下一节介绍的守恒形式（1）确定网格基架点：确定网格基架点：6个点个点 j-3,j-2,j-1,j,j+1,j+2 构造出该基架点上的目标差分格式构造出该基架点上的目标差分格式计算这这6个点可构造个点可构造5阶迎风差分：阶迎风差分：该格式为该格式为WENO 的的“目标目标”格式，格式，即，即，光滑区光滑区WENO 逼近于该格式。逼近于该格式。利用利用Taylor展开，可唯一
6、确定系数展开，可唯一确定系数（可利用小程序（可利用小程序coeff-schemes.f)实际上，还可利用分辨率优化技术，可构造出新的目标实际上，还可利用分辨率优化技术，可构造出新的目标格式（降低精度、提高分辨率，见第格式（降低精度、提高分辨率，见第4讲）。目前大量讲）。目前大量WENO的优化版做这种工作。的优化版做这种工作。Copyright by Li Xinliang7(2)将这将这6个基架点分割成个基架点分割成3个组（称为模板）个组（称为模板）每个组独立计算每个组独立计算 的差分逼近的差分逼近 模板模板1模板模板2 模板模板3模板模板1：j-3,j-2,j-1,j 模板模板2：j-2,j
7、-1,j,j+1模板模板3：j-1,j,j+1,j+2利用这利用这三个三个模板的基架点，可构造出模板的基架点，可构造出逼近逼近 的的3阶精度差分格式阶精度差分格式计算计算j点的导数点的导数u,竟然算出了三个不同的值，怎么办竟然算出了三个不同的值，怎么办？ENO 方法：方法：选择最优（最光滑）的，舍弃其余两个选择最优（最光滑）的，舍弃其余两个WENO的处理方法：的处理方法：三个都要，加权平均它们。三个都要，加权平均它们。利用Taylor展开式，可唯一确定这些系数）（可利用小程序coeff-schemes.f)也可运用优化技术，降低精度、提高分辨率Copyright by Li Xinliang8
8、（3）对这对这3个差分值进行加权平均，得到总的差分值个差分值进行加权平均，得到总的差分值原则：原则：A.模板内函数越光滑，则权重越大；模板内函数越光滑，则权重越大；模板内有模板内有间断时，权重趋于间断时，权重趋于0 B.三个模拟内函数都光滑时，这三个模拟内函数都光滑时，这三个三阶精度三个三阶精度的逼近式可组合成的逼近式可组合成一个五阶精度一个五阶精度的逼近式。的逼近式。“理想权重”（3.1）确定理想权重确定理想权重令：5阶精度容易解出：Copyright by Li Xinliang9（3.2）度量每个模板内函数的光滑程度度量每个模板内函数的光滑程度 IS越大，表示越不光滑。越大，表示越不光滑
9、。光滑区，不同模板上的光滑区，不同模板上的IS趋近同一值。趋近同一值。具体形式见下一节。具体形式见下一节。（3.3）给出实际权重给出实际权重构造构造IS方法很多，方法很多，例如：例如：第k个模板光滑区逼近光滑区逼近 O（1）量级）量级间断区间断区 量级，很大量级，很大特点：特点：间断区权重很小间断区权重很小 光滑区，趋近于理想权重光滑区，趋近于理想权重（3.4）给出最终的差分逼近给出最终的差分逼近Copyright by Li Xinliang103.Jiang&Shu 的五阶的五阶WENO格式格式 守恒型；目守恒型；目 前使用的前使用的WENO格式均为守恒型格式均为守恒型针对方程：模板模板1
10、模板模板2 模板模板3构造差分格式如下：构造方法与前文相同构造方法与前文相同（但注意这里构造的是通量，（但注意这里构造的是通量，而前文是直接构造差分格式）而前文是直接构造差分格式）针对整个网格基，构造出针对整个网格基，构造出5阶精度的通量（理想情阶精度的通量（理想情况下的通量）况下的通量）并构造出每个模板上的通量，计算出理想权重。并构造出每个模板上的通量，计算出理想权重。仍利用程序仍利用程序coeff-schemes.f求系数求系数理想权重理想权重光滑度量因子实际权重Copyright by Li Xinliang11光滑度量因子的计算光滑度量因子的计算（Jiang&Shu）k=1 k=2 k
11、=3其中：其中：j-2 j-1 j j+1 j+2 是使用模板是使用模板k 得到的插值函数得到的插值函数 利用利用j-2,j-1,j点上的点上的值构造的插值函数值构造的插值函数 ，特点：特点：光滑区趋近同一个值光滑区趋近同一个值 非光滑区值远大于光滑区非光滑区值远大于光滑区 O(1)j 点一阶、二阶导数的差点一阶、二阶导数的差分逼近（用模板分逼近（用模板k计算）计算）代入Copyright by Li Xinliang12光滑因子光滑因子 ISk 的进一步讨论的进一步讨论（光滑区精度分析）（光滑区精度分析）同样：在光滑区，在光滑区，趋近于同一个值，相互差距是一个趋近于同一个值，相互差距是一个6
12、阶小量阶小量 是小量，可忽略5阶迎风格式的通量格式具有格式具有5阶精度阶精度Copyright by Li Xinliang13最终最终5阶阶 WENO 格式为格式为 正通量情况（正通量情况（a0）负通量情况负通量情况（a0）的方法计算）的方法计算采用针对负通量（采用针对负通量（a0）的方法计算）的方法计算可采用可采用Steger-Warming,L-F,Van Leer 等分裂，等分裂，见第见第4讲讲 1)计算计算，它们是，它们是 的函数（推荐使用的函数（推荐使用Roe平均计算平均计算 ）效果更好效果更好但计算量较大但计算量较大计算量小计算量小但效果略差但效果略差有轻微振荡有轻微振荡数值测试
13、：数值测试：1维维Sod激波管问题，网格点激波管问题，网格点128 差分方法差分方法:7阶精度经典阶精度经典WENO(Jiang&Shu)分裂方式：分裂方式：Steger-Warming FVS;特征投影分裂（算术平均特征投影分裂（算术平均/Roe平均）平均）t=0.1 时刻密度及速度分布时刻密度及速度分布Steger-Warming FVS+WENO 仍无法避免振荡仍无法避免振荡特征投影分裂特征投影分裂+WENO可避免振荡可避免振荡Copyright by Li Xinliang175.WENO格式的改进格式的改进 Martin P 的的WENO-SYMBO格式格式 （Martin P,JC
14、P 220,270-289,2006)思路：思路：原原WENO格式采用迎风型网格基，耗散偏大格式采用迎风型网格基，耗散偏大 新的新的WENO 格式采用格式采用对称型网格基对称型网格基，耗散小，耗散小 采用采用优化优化方法，提高格式的分辨率方法，提高格式的分辨率理想权重情况下差分格式的修正波数；理想权重情况下差分格式的修正波数；(SYMOO:未优化未优化 8阶中心差分；阶中心差分；SYMBO 分辨率优化后的分辨率优化后的4阶格式）阶格式）针对如下针对如下Sod 激波管问题激波管问题 用用5阶阶WENO格式格式计算其数值解，画出计算其数值解，画出t=0.14时刻密度、速度及压力的分布；并与精时刻密
15、度、速度及压力的分布；并与精确解进行比较（要求数值解与精确解画在同一张图上，便于比较）。确解进行比较（要求数值解与精确解画在同一张图上，便于比较）。要求：要求：1）空间网格数空间网格数100，时间推进格式选用时间推进格式选用3阶阶Runge-Kutta,时间步长自选。时间步长自选。2）结合使用结合使用Steger-Warming 流通矢量分裂，画出结果流通矢量分裂，画出结果 3）结合使用特征投影分裂（结合使用特征投影分裂（Roe平均计算平均计算Uj+1/2)，画出结果，画出结果18Copyright by Li Xinliang作业作业 8.1 针对如下针对如下Shu-Osher 激波激波-密
16、度扰动波干扰问题：密度扰动波干扰问题：用用5阶阶WENO格式格式计算其数值解，画出计算其数值解，画出t=0.1时刻密度、速度及压力的分布时刻密度、速度及压力的分布 要求：要求：1）空间网格数空间网格数200，时间推进格式选用时间推进格式选用3阶阶Runge-Kutta,时间步长自选。时间步长自选。2）结合使用结合使用Steger-Warming 流通矢量分裂，画出结果流通矢量分裂，画出结果 3）结合使用特征投影分裂（结合使用特征投影分裂（Roe平均计算平均计算Uj+1/2)，画出结果，画出结果 4）使用使用2000个网格点计算，其结果作为个网格点计算，其结果作为“精确解精确解”，与其它结果画在
17、一起，与其它结果画在一起，便于比较。便于比较。19Copyright by Li Xinliang作业作业 8.2初值为:20Copyright by Li Xinliang作业作业 8.3（选作题）（选作题）推导推导7阶精度的阶精度的WENO格式，给出详细的推导过程及格式的具体表达式格式，给出详细的推导过程及格式的具体表达式提示：提示：与与5阶阶WENO推导思路相同，但网格基架点扩大了，模板数目也增加了推导思路相同，但网格基架点扩大了，模板数目也增加了j-3 j-2 j-1 j j+1 j+2 j+3正通量正通量(a0)的的WENO通量通量 使用基架点使用基架点j-3,j-2,j-1,j,j
18、+1,j,j+2,j+3如上图示，将其分割为如上图示，将其分割为4个组（模板），每个模板上个组（模板），每个模板上4个基架点。个基架点。1）先构建整个网格基（先构建整个网格基（7个点）上个点）上7阶迎风格式的通量表达式阶迎风格式的通量表达式 2）对于每个模板，构造逼近对于每个模板，构造逼近j点导数的点导数的4阶差分格式的通量表达式阶差分格式的通量表达式 3）按照理想权重进行组合可得到理想格式（按照理想权重进行组合可得到理想格式（7阶迎风格式）阶迎风格式）对照具体表达式，求出理想权重对照具体表达式，求出理想权重Ck 4）构建构建WENO通量的加权表达式通量的加权表达式 5）仿照）仿照 本本PPT
19、第第35页的方法，给出光滑度量因子的具体表达式页的方法，给出光滑度量因子的具体表达式 可利用小程序求系数可利用小程序求系数coeff-schemes.f减小工作量减小工作量Copyright by Li Xinliang21 8.2 其他激波捕捉格式简介其他激波捕捉格式简介 1.Godnov格式格式思路：思路：利用精确利用精确Riemann解解难点：难点：精确精确Riemann解要求间断两侧物理量为解要求间断两侧物理量为常数常数 对策：对策：采用采用分片常数分片常数代替原先的函数代替原先的函数方法：方法：1）采用分片常数代替原先的函数采用分片常数代替原先的函数 即假设即假设xj-1/2,xj+
20、1/2 区间物理量为区间物理量为xj点上的值点上的值 2）在每个间断处，精确求解在每个间断处，精确求解Riemann问题，得到问题，得到D Dt时刻物理量的分布。时刻物理量的分布。（假设（假设D Dt足够小，各间断之足够小，各间断之间无相互影响，因此可独立求解）。间无相互影响，因此可独立求解）。3）在区域在区域xj-1/2,xj+1/2 上平均，得到上平均，得到D Dt时刻时刻xj点上物理量的值。点上物理量的值。4）重复重复1-3，直到指定时刻。，直到指定时刻。j-1 j j+1 间断1 间断2 间断3 k=1/3 三阶迎风三阶迎风 k=0 k=-1Copyright by Li Xinlia
21、ng222.MUSCL方法方法理解理解1：Godnov方法的改进方法的改进 分片线性函数代替分片常数（以后介绍）分片线性函数代替分片常数（以后介绍）理解理解2：一阶格式的基础上的改进一阶格式的基础上的改进二阶迎风二阶迎风Formm格式格式k=1 二阶中心二阶中心1阶迎风阶迎风修正部分修正部分间断处都会产生振荡间断处都会产生振荡思路：思路：振荡由修正项引起，需要对修正项进行限制振荡由修正项引起，需要对修正项进行限制Copyright by Li Xinliang23新格式新格式函数函数 minmod(a,b)a,b 符号相反，则返回符号相反，则返回0 符号相同，则返回绝对值小的符号相同，则返回绝
22、对值小的1阶迎风阶迎风修正部分修正部分j-1jj+1j-2j+1/2思路：思路：修正部分是由修正部分是由 j-1,j,j+1 三点计算而得，三点计算而得，是是 与与 的插值的插值旧格式旧格式二者的组合：j-1,j 计算；j,j+1计算 如果两个修正方案如果两个修正方案趋势趋势（符号）相反，干脆不修正（符号）相反，干脆不修正 如果如果趋势趋势相同，则取（绝对值）最小者相同，则取（绝对值）最小者 避免振荡避免振荡两个修正方案两个修正方案 选前者：二阶迎风选前者：二阶迎风 选后者：二阶中心选后者：二阶中心 组合之：组合之：Formm、3阶迎风阶迎风4 阶3 阶阶（TVD型）型）2阶Copyright
23、 by Li Xinliang24 8.3 时间推进方法时间推进方法 1.显格式显格式推荐方法：推荐方法：Runge-Kutta法法更高阶 Copyright by Li Xinliang252.隐格式算法简介隐格式算法简介（以二维（以二维Euler方程为例）方程为例）1）原理介绍原理介绍（1）时间离散后时间离散后方案方案 A：直接将（直接将（1）进行空间离散，得到）进行空间离散，得到 Un+1 的代数方程组的代数方程组困难：困难：大型大型非线性非线性方程组，求解困难方程组，求解困难方案方案B:设计一种迭代算法设计一种迭代算法令两端同时添加显式项两端同时添加显式项右端项，已知右端项，已知是已知
24、项，可采用某种差分方法是已知项，可采用某种差分方法显式显式计算得到计算得到（对算法无限制）（对算法无限制）Copyright by Li Xinliang26(2)已知，（已知，（2）为线性方程）为线性方程。（1）(2)是一个线性化过程是一个线性化过程 含义：含义：先用显格式计算，再先用显格式计算，再用隐格式计算修正量用隐格式计算修正量线性方程，离散求解线性方程，离散求解离散后为大型带状方程组，求解计算量大离散后为大型带状方程组，求解计算量大LU-SGS原理：原理：将矩阵分解为上、下三角阵，避免矩阵求逆运算将矩阵分解为上、下三角阵，避免矩阵求逆运算隐式问题显式化隐式问题显式化未知量：显式部分隐
25、式修正若隐式修正为若隐式修正为0，则为显格式，则为显格式Copyright by Li Xinliang272)LU-SGS 方法方法a)将矩阵将矩阵A,B分裂分裂为了简化计算，通常采用为了简化计算，通常采用L-F分裂分裂矩阵分裂矩阵分裂1阶迎风格式离散阶迎风格式离散整理整理迭代收敛迭代收敛后，不影响精度后，不影响精度Copyright by Li Xinliang28其中：其中：特点：特点：严格对角占优，收敛性好严格对角占优，收敛性好 稳定性好稳定性好可采用局部时间步长等加速收敛措施可采用局部时间步长等加速收敛措施近似近似LU分解分解 do j=1,ny do i=1,nxEnddoEndd
26、oCopyright by Li Xinliang29具体求解方法具体求解方法（1）计算右端项计算右端项显式，可采用前面构造的各种方法计算显式，可采用前面构造的各种方法计算（2）计算计算（3）计算计算 自变量采用自变量采用n时刻的值时刻的值由于由于L矩阵的下三角特性，可采用推进方法矩阵的下三角特性，可采用推进方法显式显式计算计算(4)计算计算显式计算显式计算（5）时间推进时间推进1阶精阶精度，空间精度由度，空间精度由RHS的算法决定。的算法决定。可用于定常及非可用于定常及非定常问题；定常问题；非定常问题通常非定常问题通常采用采用内迭代内迭代（保（保证精度）证精度）向上扫描向上扫描向下扫描向下扫描LU分解的对称分解的对称Gauss-Seidel迭代算法迭代算法