CFDfudan可压缩方程离散李新亮老师课件

2013复旦大学暑期课程：计算流体力学李新亮lixl@Tel:82543801力学所主楼219

1CopyrightbyLiXinliang讲义、课件上传至

（流体中文网）->“流体论坛”->“CFD基础理论”下载地址：/?cid=1cc0dcbff560c149&id=1CC0DCBFF560C149%21128&sc=documents第四讲可压缩流体力学方程组的离散方法（上）

1.通量分裂技术:流通矢量分裂FVS2.Roe格式，通量差分分裂3.隐式时间推进：LU-SGS（下）

有限体积法初步2CopyrightbyLiXinliang§1通量分裂技术格式F+格式F-对流项：信息（波）从上游传至下游——上游信息更重要——迎风差分扩散项：信息从中心向周围扩散——不区分上、下游——中心差分迎风差分优点：有效利用信息传播的方向，增强稳定性微分与差分方程的影响域N-S方程：单波方程：单波方程——一个波，容易判断波传播方向N-S对流项（Euler）——方程组：多波问题，复杂双曲方程组的原则——特征分解，找到独立传播的波常系数矩阵A的情况——完全解耦，独立求解变系数矩阵A的情况——局部讨论3CopyrightbyLiXinliangCopyrightbyLiXinliang4Step1针对模型方程构造差分格式Step2将格式推广到Euler方程方法1：流通矢量分裂(FVS)

特点：对流通矢量f(U)本身进行分裂

Steger-Warming,vanLeer,Lax-Friedrichs利用通量分裂技术，将模型方程推广到Euler(或N-S)方程格式1格式2格式1格式2CopyrightbyLiXinliang5方法2：通量差分分裂（FDS）特点：对流通矢量f的导数进行分裂

Godnov,Roe,HLL,HLLC

方法3：AUSM类方法vanLeer分裂法+压力项单独处理

利用Riemann解利用差分表达式，计算求解Riemann问题，获得通量1.Jacobian系数矩阵及其性质重要性质特点：A可以像常数一样，和求导运算交换6CopyrightbyLiXinliang2流通矢量分裂(FVS)利用性质=+优点：耗散小缺点：导数间断方式A:特点：不必进行矩阵运算，计算量小Steger-Warming分裂A:Steger-Warming分裂7CopyrightbyLiXinliangSteger-Warming具体步骤（以一维为例）已知1）计算2）计算3）计算4）带入（1）式得到5）利用不同的迎风格式，分别计算

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

计算特征值

分裂特征值，计算

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

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

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

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

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

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

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

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

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

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

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

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

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

均可2）

在网格基上计算…j-2j-1jj+1…计算fj+1/2用到的点注意，在该网格基上（例如k=j-1,j,j+1）保持不变例如：3)利用已构造好的差分格式，计算通量4）得到总通量

5）计算差分（j点处）步骤的算法描述（注意：实际上是两重循环）doj=1,Ndok=j-1,j+1(网格基，可以是更多或更少点)

enddoenddodoj=1,N

enddo需要多次矩阵运算，计算量大

守恒性好，耗散小，数值解质量好16CopyrightbyLiXinliangCopyrightbyLiXinliang171.单方程的Roe格式线性化，用平均变化率代替(j,j+1)之间的变化率a(u)“平均斜率”，不等于“斜率的平均值”，也不等于中点处的斜率§2Roe格式——FVS分裂技术非线性情况根据Langrage中值定理，[uL,uR]之间必有一点uRoe,该点处的斜率为平均斜率;二次函数f(u)=u2中点处的斜率=平均斜率CopyrightbyLiXinliang182.方程组的情况平均斜率线性化，以平均增长率代替瞬时增长率[j,j+1]区间内连续，且可通过相似变换对角化

应当具有的性质常系数方程的Riemann解CopyrightbyLiXinliang19平均斜率xj+1/2常系数单波方程的Riemann解Roe格式：微分型近似Riemann解CopyrightbyLiXinliang203.矩阵的构造关键：“向量除以向量”？直接求平均增长率：uf(u)uLuRuRoeRoe点的斜率为平均斜率（根据拉格朗日中值定理，[UL,UR]区间内肯定存在Roe点）思路1：在UL与UR之间寻找一个点URoe,该点处的增长率为平均增长率f(u)=u2u二次函数——Roe点与中点重合标量函数的启示：Roe点肯定存在（Langrage中值定理）

二次函数的中点即为Roe点思路2：进行坐标变换，得到一个二次（齐）函数引入如果

是二次（齐）函数，则其中点即为Roe点重要启示更准确地讲，应当是要求为W的线性函数，

即增长率为线性函数（中点处的增长率刚好为平均增长率）CopyrightbyLiXinliang21针对Euler方程的具体构造引入新变量：则：目的：使得F(w)是W二次齐函数（增长率为线性函数）f(U)不是U的二次齐函数二次齐函数！中点处的斜率即为平均斜率。Roe点Roe点为：增长率为线性函数！CopyrightbyLiXinliang22最终：其中如下计算：平均增长率（矩阵）含义：左、右两个状态点的某种平均（称为Roe平均，为密度加权平均）

该状态点对应的增长率（矩阵）为平均增长率（矩阵）

实际上是一种“等效平均”。效果优于简单的算数（或几何）平均。

三维情况下，还有其他量（如压力、温度、音速等）用这三个量计算（5）简单易记：CopyrightbyLiXinliang23Roe格式的计算步骤（半离散）已知n时刻所有网格点上的物理量，对于j点：1）利用差分格式计算UR,UL

2）采用Roe平均公式（5）计算Roe平均值3）将Jacobian矩阵进行特征分解:

计算4)计算5）计算6）计算空间导数7）时间推进，计算下一时间步的值。j-1jj+1与前文（第3，4讲）的形式相同，仅需把式中的密度、压力、速度等换成经过Roe平均的密度、压力、速度即可其中：CopyrightbyLiXinliang24可能出现导数不连续，可能引起数值振荡实际使用时可用如下函数代替——所谓“熵修正”实际上是在特征值0点周围增加了耗散Roe格式的优点：1）保持守恒性的同时，严格保证了特征方向2）便于推广到高精度格式——特征投影分裂中使用Roe平均即可（见本PPT第5页）。推广到高阶后，虽不再保证严格的特征方向，但仍优于采用算数平均方法。Roe格式的不足：

本身精度只有一阶；

推广到高阶后，特征方向无法严格保证；

推广到二维或三维后，特征方向无法严格保证，出现振荡。CopyrightbyLiXinliang25关于f(U)与f(W)深入讨论新变量虽然是“一次齐函数”但有变量在分母上干净的二次齐函数自变量W的线性函数！实质区别：是自变量的线性函数，而是自变量的非线性函数！自变量U的非线性函数使用U做自变量的优点：物理意义鲜明（质量密度、动量密度和能量密度），守恒性好使用W做自变量的优点：Jocabian矩阵为线性矩阵思考：如果在CFD计算中，使用W替换U做自变量会怎样？4阶3阶（TVD型）2阶CopyrightbyLiXinliang26

§

3时间推进方法

1.显格式推荐方法：Runge-Kutta法更高阶……CopyrightbyLiXinliang272.隐格式算法简介（以二维Euler方程为例）

1）原理介绍（1）时间离散后方案A：直接将（1）进行空间离散，得到Un+1的代数方程组困难：大型非线性方程组，求解困难方案B:设计一种迭代算法令两端同时添加显式项右端项，已知是已知项，可采用某种差分方法显式计算得到（对算法无限制）CopyrightbyLiXinliang28(2)已知，（2）为线性方程。（1）(2)是一个线性化过程

含义：先用显格式计算，再用隐格式计算修正量线性方程，离散求解离散后为大型带状方程组，求解计算量大LU-SGS原理：将矩阵分解为上、下三角阵，避免矩阵求逆运算隐式问题显式化未知量：显式部分隐式修正若隐式修正为0，则为显格式CopyrightbyLiXinliang29Step1:求解Step2:求解LU-SGSi,ji-1,ji+1,ji,j+1i,ji-1,ji+1,ji,j+1i,j+1i,j+12）LU-SGS方法原理简介：忽略右、上未知量忽略左、下未知量i,ji-1,ji+1,ji,j+1i,j+1原则上：（1）（2）两步反复迭代，直至收敛CopyrightbyLiXinliang30LU-SGS方法求解Euler方程a)将矩阵A,B分裂为了简化计算，通常采用L-F分裂矩阵分裂1阶迎风格式离散整理迭代收敛后，