计算流体课程-02.ppt

《计算流体课程-02.ppt》由会员分享，可在线阅读，更多相关《计算流体课程-02.ppt（239页珍藏版）》请在三一文库上搜索。

1、计算流体力学讲义 第九讲 有限体积法（1）,知识点：,1,Copyright by Li Xinliang,有限体积法的基本概念 重构和反演 迎风型有限体积法Riemann求解器；Roe格式的新理解：近似Riemann解 多维迎风型有限体积法坐标旋转,Copyright by Li Xinliang,2,知识回顾,1. 差分方法的基本概念：,差分格式、修正方程、相容性、收敛性、稳定性、LAX等价定理,2. 精度分析、稳定性分析与分辨率分析（修正波数）,Taylor分析,Fourier分析,修正波数,激波捕捉格式 GVC, NND, Roe, Godnov, MUSCL, TVD, WENO,E

2、uler (N-S) 方程的通量分裂 逐点分裂、特征投影分裂 （建议使用Roe平均）,5. 隐格式求解的LU-SGS方法,要点： a. 引入差量，方程线性化 b. 单边差分，隐式代数方程显式（推进）化,以一维为例，多维可直接推广,方法1：直接隐式离散,直接求解,非线性方程组，计算量大,方法2,差量化,线性化,已知项,线化微分方程,Copyright by Li Xinliang,3,Copyright by Li Xinliang,4,求解思路：如果直接离散，得到线性代数方程组，仍需求解，计算量大（多维情况）,如果能单侧差分就好解了！,多对角方程组，不好解（多维情况）,中心（双侧）离散,如果单

3、侧离散,单侧离散，可推进求解，免受解方程组之苦。真简单,Copyright by Li Xinliang,5,可是，A有正有负，无法单侧差分化,还是个三对角的,奇思妙想：如果分成两个子步，各自用单侧值，就简单多了,强行单侧差分会不稳定的,近似LU分解,Step 1:,近似LU分解,Step 2:,均为递推求解 （两次扫描），免受解方程组之苦,j -1 - j,j+1 j,以上描述适用于求解定常问题，求解非定常问题该过程可用于内迭代。,迭代收敛后q趋于0，精度由右端项决定,Copyright by Li Xinliang,6, 9.1 有限体积法入门,有限体积法主要优势： 处理复杂网格,差分法处

4、理复杂外形 坐标变换,坐标变换函数必须足够光滑 否则损失精度,实际问题： 外形复杂， 光滑的结构网格生成困难,Copyright by Li Xinliang,7,9.1.1 有限体积法 的基本概念,实质： 把几何信息包含于离散过程中,虽然简单，但有助于建立基本概念,j-1 j j+1,j-1/2 j+1/2,1. 全离散型过程,含义： f在j+1/2点的值 （注意与差分法的区别）,在控制体上积分原方程,定义：,空间平均,时间平均,精确推导，不含误差,提示： 为区间内的空间及时间平均值，如果把它们理解为某点的值，会产生误差,Copyright by Li Xinliang,8,积分（精确）,重

5、构（Reconstruction),有限差分法的离散：数值微分过程 有限体积法的离散：数值积分过程,积分方程,离散化,反演（evolution),(1) 重构过程,A. 零阶重构，假设分片常数,j-1,B. 线性重构，假设分片线性函数,零阶重构与一阶重构示意图,j,j+1,or,or,或其他方法,C. 更高阶的重构例如: 分片二次函数 （PPM）, WENO等,重构是有限体积的空间离散化过程，有多种方法,Copyright by Li Xinliang,9,（2） 演化过程 （以线性方程为例）,需要得知时间演化信息，通常利用特征方程,若采用零阶重构：,则：,假设时间步长足够小,则方程为：,等价

6、于一阶迎风差分,Riemann解,Copyright by Li Xinliang,10,若采用线性重构,若,Warming-Beam,Lax-Wendroff,0阶重构 1阶精度 线性重构 2阶精度,一维均匀网格的有限体积法等价于有限差分法,Euler方程： 演化过程可通过Riemann解或近似Riemann解进行,Copyright by Li Xinliang,11,2. 半离散方法,全离散： 积分方程 代数方程 （守恒性好，但复杂） 半离散： 积分方程 常微分方程 （简便，便于使用R-K等成熟方法）,仅空间积分,f 在j+1/2点的值，仍需要使用周围点 进行插值,通常无法精确计算， 可

7、采用近似值 代替,等价于二阶中心差分,半离散,j-1 j j+1,j-1/2 j+1/2,重构,Copyright by Li Xinliang,12,9.1.2 一维Euler方程的迎风型有限体积法,j-1 j j+1,j-1/2 j+1/2,半离散,1. 重构,控制体积,j-1,j,j+1,左重构值,右重构值,选择不同的模板会得到不同的重构方案,向左偏的模板产生 向右偏的模板产生,差分法 同一点的导数可使用向前差分和向后差分，根据特征方向选择之,例如： 0阶重构 1阶单边重构,根据特征方向，选择左通量或右通量,途径1： FVS,途径2：FDS,Copyright by Li Xinlian

8、g,13,2. 分裂方法 （1）: FVS方法 （流通矢量分裂 逐点分裂）,具体方法： Steger-Warming 分裂 Lax-Friedrichs分裂 Van Leer分裂： Liou-Steffen分裂： （压力项与其他项分开， AUSM类格式的基础）,根据当地Mach数分裂,保证 的Jocabian阵特征值为正， 的为负,正通量： 向左偏斜重构； 负通量： 向右偏斜重构 偏重向上游,与迎风差分法类似： 网格基（或权重）偏重上游,差分、有限体积都可使用,一个参数，反映全部特征,Copyright by Li Xinliang,14,小知识： Liou-Steffen分裂,对流项,压力项

9、,思路： 决定特征的关键参数 当地Mach数,超音速，x-方向,超音速，x+方向,因此，对Mach数进行分裂更为简洁！,显然：,参考文献： Toro: Riemann Solvers and Numerical Methods for Fluid Dynamics, section 8.4.4 Liou: Ten Years in the making AUSM family, NASA TM-2001-210977,类似 Van Leer分裂，但压力单独处理,M,保证光滑过渡,M=1,Copyright by Li Xinliang,15,(3) FDS 方法 （通量差分分裂特征投影分裂）,

10、1. 利用精确Riemann解Godnov格式,目的：,j-1 j j+1,j-1/2 j+1/2,控制体积,j-1,j,j+1,左重构值,右重构值,1） 精确求解Riemann问题,2）,精度： 取决于重构的精度 （原则上可任意阶）,差分法：Godnov格式使用分片常数，精度1阶 有限体积法：先重构，再解Riemann问题，可高阶,精确Riemann解（见本讲座第2讲）需迭代求解，计算量大 - 近似Riemann解,整体思路： 先重构自变量（两种方案得到 ）， 再求解Riemann问题（或用FVS）得到通量的方法通常称为MUSCL方法。,Copyright by Li Xinliang,16

11、,差分法与有限体积法区别与联系（二阶迎风+FVS为例）,差分、有限体积,差分（通常做法）： 直接插值通量fi+1/2,有限体积：先插值自变量U，然后计算通量f：,先插值自变量，再计算通量的方法，称为MUSCL类方法。 是有限体积法的常用方法（差分法也可以用）,单侧重构，以避免跨过激波,还可使用FDS方法，重构后求解Riemann问题,当f=f(U) 连续时，对f插值与对U插值精度相同。,（称为数值流通量） 的含义,Copyright by Li Xinliang,17,重要概念澄清： 重构与插值,A. 有限差分法：,j+1/2,切线,j-1/2,j,j-1,注意： 与 f 在xj+1/2点的值

12、含义不同！,用周围几个点的值 计算 的过程称为“重构”，不能理解为用 来插值,记号 确实容易混淆，让人容易联想起 。记为 更好些,否则，最高只能达到2阶精度了！,是控制体内的平均值,（称为数值流通量） 的含义,Copyright by Li Xinliang,18,重要概念澄清： 重构与插值,B. 有限体积法：,j+1/2,j-1/2,确实为f在xj+1/2点的值 ！,通常做法： 1） 用 计算出 2）,u在xj+1/2点的值！,关键： 是用 计算 （称为重构） ，而不是用 计算 （是标准的插值）；否则最高也只能达到2阶精度。,19,概念：MUSCL与非MUSC类方法,j+1/2,切线,j-1

13、/2,j-1,差分,有限体积,方法1 （非MUSCL类）： 直接利用周围几个点的函数值 或 ）直接计算 （或 ）,如何计算 或 ？,方法2 （MUSCL类）： 利用周围几个点的自变量值 （或 ）计算出 （或 ）； 然后再计算 （或 ）,当f=f(u)是连函数时，二者精度相同,f的误差与u的误差同阶,Copyright by Li Xinliang,20,2. 近似Riemann解 例： Roe格式,与差分法的Roe格式形式相同,理解： 近似Riemann解（Euler方程 常系数线性化解）,u,f(u),uL,uR,uRoe,利用Roe平均， 刚好是左右两点间的平均增长率，实现了常系数线性化。

14、,常系数双曲方程组，易解！,思路： 用平均增长率矩阵 取代瞬时增长率矩阵A,不但实现了线性化，而且实现了常系数化。 利用二次齐函数的性质，可找到了Roe点（即Roe平均点），该点处的增长率刚好等于平均增长率。,Roe平均,常系数化,线性化,常系数线性单波方程的Riemann问题太简单了,21,常系数方程组的Riemann问题,解耦了的单波方程，有精确解,初值,Copyright by Li Xinliang,22,解为,线性化条件 （并利用齐函数性质）,与差分法的Roe格式相同,还有各种其他类型的近似Riemann解（今后介绍）,Copyright by Li Xinliang,23,9.1.

15、3 多维问题的有限体积法,二维问题,一维Riemann问题，坐标选取不当，变为 “二维”Riemann问题,x,y,差分法： 独立计算 只考虑各自的特征方向,由于非线性， 实际（二维）特征方向并非x,y方向特征量的线性组合。 特征方向计算不严格，带来误差,差分方法： 多维情况，特征理论复杂，通常x,y方向独立计算,转化为 x方向与y方向的两个一维问题,逐点分裂,特征投影分裂,完全按照一维情况独立处理,局部坐标旋转？ 差分算法设计造成局部旋转困难,差分法 的多维 处理方法,1. 小知识： 差分方法如何处理高维问题的 ？优缺点？,优点： 简单 缺点：特征方向计算不准,Copyright by Li

16、 Xinliang,24,2. 二维有限体积方法的离散过程,在以某节点为中心的控制体上积分,i,j,k,非结构网格的控制体,i+1,j,i-1,j,i,j+1,i,j-1,k3,k1,k2,k4,k5,结构网格的控制体,x,y,n,体积平均,控制体边界垂直于节点连线（也可选其他方式）,垂直平分线,n,1) 建立控制体,2) 在控制体上积分,离散方程,需要重构： 由节点上平均值 给出函数分布，最终给出通量,表示第m个界面上的值,1） 重构 两种不同的重构方案，向左偏及向右偏。 给出两种结果： 及,Copyright by Li Xinliang,25,i,j,i+1,j,i-1,j,i,j+1,

17、i,j-1,n,左重构,右重构,2） 由左右重构得到的自变量： 和 给出通量 方案A: FVS 方案B: 解Riemann问题 （常用）,3. 二维迎风型有限体积法,例如： 0阶重构：,线性重构:,用i, i-1点的值 插i+1/2点的值 （网格剧烈变化时，应当用实际坐标插值）,用i, i+1点的值 插i+1/2点的值,x,y,看似二维Riemann问题，其实是一维的，坐标旋转一下就行了,Copyright by Li Xinliang,26,x,y,x,y,（通常）进行坐标旋转,旋转q角后的坐标系 (x,y),性质： Euler方程的旋转不变性,形式上完全不变 （仅需把u,v,x,y换成u,

18、v,x,y即可）,其中：旋转矩阵,旋转 q 角,矩阵表示,Copyright by Li Xinliang,27,i,j,i+1,j,i-1,j,i,j+1,i,j-1,左重构,右重构,局部坐标系,x,y,x,y,旋转q角后的坐标系 (x,y),习题： 设 n 为平行x轴的向量，试证明：,证明：,坐标旋转，标量不变 向量的模不变,Copyright by Li Xinliang,28,i,j,i+1,j,i-1,j,i,j+1,i,j-1,左重构,右重构,于是：,x,y,x,y,旋转q角后的坐标系 (x,y),其中下标m表示控制体第m个面（线）， 表示该面的面积 （长度）,于是，问题转化为求控

19、制面上的,这个量有两个重构方案,方法1： FVS: 方法2： 需要求解Riemann问题，旋转后，转化为“扩展的”1维Riemann问题,Copyright by Li Xinliang,29,解释： “扩展的”一维Riemann问题,x,y,x,y,旋转q角后的坐标系 (x,y),问题本身是一维的 ： 所有变量都只沿着x方向分布，沿y方向均匀 允许有y方向的速度v (比纯一维问题多一个变量）,v的存在对流动的一维性质无任何影响,举例： Sod激波管问题（一维）。 如果在沿y方向匀速运动的坐标系中观察，则方程为“扩展的一维问题”，但不影响其一维性质,坐标系沿y方向匀速运动,x,y,可用精确Ri

20、emann解，也可用Roe等近似解,Copyright by Li Xinliang,30,二维迎风型有限体积法求解步骤 1） 对n时刻的平均量 进行重构，给出控制面上的左、右重构值 ， 2）将以上值旋转到（每个）控制面法向的局部坐标系下： 3） 求解上述“扩展的”一维Riemann问题，给出后续时刻控制面上的值 4）利用积分型方程： 计算下一时刻的平均量,i,j,i+1,j,i-1,j,i,j+1,i,j-1,左重构,右重构,0阶重构：,1阶重构 （线性重构）：,更复杂的重构（WENO等）,下标m指的是第m个控制面上的值,Copyright by Li Xinliang,31,知识回顾： R

21、iemann问题精确解,Riemann问题,问题描述： 初始时刻，物理量分布存在单个间断；间断两侧物理量为常数。,求解思路： 采用积分方程 单个间断，且间断两侧物理量为常数情况下： 积分方程转化为代数方程,代数方程： 质量、动量、能量守恒,计算出 ， 将 与这三个值进行比较，判断会产生的情况。具体见下图：,Copyright by Li Xinliang,32,Riemann问题的具体计算步骤 （全流场）,1. 判断可能会出现的情况（五种情形之一）,a. 定义函数,b. 进行判断,情况5,情况4,情况3,情况1,情况5,情况4,情况2,情况1,单调增函数，性质很好,Copyright by L

22、i Xinliang,33,2. 求解中心区的压力和速度,单未知数的代数方程，迭代求解（例如Newton法，F(p)性质好，求解不困难）,3. 确定中心区接触间断两侧的密度 以及左、右波传播的速度 a. 左波为激波的情况（情况1,3）,b. 左波为稀疏波的情况 （情况2，4，5）,中心区 接触间断左侧的物理量,膨胀波的波头及波尾速度,激波的传播速度,对于情况（5），波尾速度为：,中心区为真空，音速 无定义，改由该式计算,Copyright by Li Xinliang,34,c. 右波为激波的情况（情况1，2）,中心区 接触间断右侧的物理量,b. 右波为稀疏波的情况 （情况2，4，5）,4.

23、计算稀疏波区域的值（如果有稀疏波的话）,a. 左稀疏波 b. 右稀疏波,情况2，4,情况5：,注意： 教科书32页c的公式有误！,Copyright by Li Xinliang,35,有限体积法 “扩展的”Riemann问题的计算方法 （中心线x=0处）,迎风型有限体积法，需要求解“扩展的”一维Riemann问题,x,y,物理问题分析： 所有物理量均沿x方向一维分布，沿y方向均匀分布。,仅需计算t时刻x=0处各物理量的值,v和w跟随流体运动，相当于“被动标量”,穿过激波及稀疏波，切向速度不变,Copyright by Li Xinliang,36,求解t时刻x=0处物理量的具体步骤,Step

24、 1: 求解 得到中心区压力 Step 2: 计算中心区的速度 Step 3: 根据 及 判断会出现哪种情况（五种情况之一） Step 4: 根据具体情况（左、右波是激波还是膨胀波）计算出中心区接触间断两侧的密度 及 Step 5: 如果中心区x方向速度0, 则中心线（x=0）处的密度及切向速度为接触间断右侧的值，否则为接触间断左侧的值。,具体公式见本PPT33-34页， 本页中上标“L”和“R”分别对应原先的下标“1”和“2”,x,t,v和w没有给计算过程带来任何麻烦,先无视v和w的存在，求解标准1维Riemann问题。再根据u的符号确定中心线的v和w,Copyright by Li Xin

25、liang,37,作业9.1 求解 “坐标旋转的Sod激波管问题”,x,y,物理问题描述：如右图示，有一条与x轴夹角120的直线，左右两侧充满同种介质的无粘完全气体。初始时刻，左右两侧的气体状态为： 左侧： ，右侧 试计算t=0.14时刻的流动分布。,计算要求： 1） 计算域 网格 2）初值设置如图所示； 3）空间离散采用有限体积法计算。采用线性重构（见上页）。时间推进采用3阶Runge-Kutta方法 4) 分别采用FVS方法及Roe-FDS （又称Roe- Riemann近似解，见本PPT 22页）两种不同的方法计算通量。 5） 绘制出t=0.14时刻密度、速度u,v及压力的二维分布。 6

26、） 绘制出t=0.14时刻x轴（垂直于初始间断面，见右图）上的密度、沿x方向的速度及压力的一维分布，并同精确解进行比较。,x,y,计算流体力学讲义 第十讲 有限体积法（2）,知识点：,38,Copyright by Li Xinliang,近似Riemann求解器 HLL / HLLC 中心型有限体积法人工粘性 具体问题的计算 翼型绕流 (边界条件、具体解法、粘性项处理),Copyright by Li Xinliang,39,知识回顾,在以某节点为中心的控制体上积分,i,j,k,非结构网格的控制体,i+1,j,i-1,j,i,j+1,i,j-1,k3,k1,k2,k4,k5,结构网格的控制体

27、,x,y,n,体积平均,控制体边界垂直于节点连线（也可选其他方式）,垂直平分线,n,1) 建立控制体,2) 在控制体上积分,离散方程,重构： 由节点上平均值 给出函数分布，最终给出通量,表示第m个界面上的值,1. 有限体积法的离散过程重构,1） 重构 两种不同的重构方案，向左偏及向右偏。 给出两种结果： 及,Copyright by Li Xinliang,40,i,j,i+1,j,i-1,j,i,j+1,i,j-1,n,左重构,右重构,2） 由左右重构得到的自变量： 和 给出通量 方案A: FVS 方案B: 解Riemann问题 （常用）,例如： 0阶重构：,线性重构:,用i, i-1点的值

28、 插i+1/2点的值 （网格剧烈变化时，应当用实际坐标插值）,用i, i+1点的值 插i+1/2点的值,x,y,看似二维Riemann问题，其实是一维的，坐标旋转一下就行了,2. 迎风型有限体积法,（称为数值流通量） 的含义,Copyright by Li Xinliang,41,重要概念澄清： 重构与插值,A. 有限差分法：,j+1/2,切线,j-1/2,j,j-1,注意： 与 f 在xj+1/2点的值含义不同！,用周围几个点的值 计算 的过程称为“重构”，不能理解为用 来插值,记号 确实容易混淆，让人容易联想起 。记为 更好些,否则，最高只能达到2阶精度了！,是控制体内的平均值,（称为数值

29、流通量） 的含义,Copyright by Li Xinliang,42,重要概念澄清： 重构与插值,B. 有限体积法：,j+1/2,j-1/2,确实为f在xj+1/2点的值 ！,通常做法： 1） 用 计算出 2）,u在xj+1/2点的值！,关键： 是用 计算 （称为重构） ，而不是用 计算 （是标准的插值）；否则最高也只能达到2阶精度。,43,概念：MUSCL与非MUSC类方法,j+1/2,切线,j-1/2,j-1,差分,有限体积,方法1 （非MUSCL类）： 直接利用周围几个点的函数值 或 ）直接计算 （或 ）,如何计算 或 ？,方法2 （MUSCL类）： 利用周围几个点的自变量值 （或

30、）计算出 （或 ）； 然后再计算 （或 ）,当f=f(u)是连函数时，二者精度相同,f的误差与u的误差同阶,Copyright by Li Xinliang,44,10.1 近似Riemann解,迎风型有限体积法通常需要求解Riemann问题； 精确Riemann解计算量大,10.1.1 HLL Riemann近似求解器 （Harten, Lax & van Leer）,思路： 在控制体上积分，分析积分方程,在控制体 上积分Euler方程,Ref.： E. F. Toro: Riemann Solvers and Numerical Methods for Fluid Dynamics, Sp

31、ringer, 2009 (Third Edition),控制体内质量（动量、能量）的减少等于流出控制面的通量,Copyright by Li Xinliang,45,若控制体空间足够大（或时间跨度足够小）， 扰动波未达到控制体的边界（如图）,未扰动的值,边界上，未受扰动，保持常数,把积分域分成三段； 中间为受扰动段,为扰动波传播的速度（激波或稀疏波的波头传播速度）,左、右波的速度假设已知,Copyright by Li Xinliang,46,未扰动段，不影响积分，可不考虑,未扰动区，保持常值UL UR,消项,T时刻，扰动区内的平均值,T时刻的流动,以平均值 代替分布值,Copyright

32、by Li Xinliang,47,思路： 以扰动区内的平均值，代替瞬时值（ 扰动区很小，误差也不大）,含义： T时刻扰动区内的平均值,于是，HLL的近似Riemann解为：,HLL近似解是 左右波都是激波，中心区域均匀分布的一个近似模型,双激波近似,假设： 1）左右波都是激波 2）中心区无接触间断 （近似）,显然， 假设中心区无接触间断是不合理的（造成方程数多于未知数个数），但可作为一种近似 （类似最小二乘解）,Copyright by Li Xinliang,48,利用HLL近似Riemann解计算中心线上的通量,在控制体 （图中红色部分）积分Euler方程，得：,控制体积 （图中绿色部分

33、）积分,二者计算结果相同,带入 表达式,Copyright by Li Xinliang,49,最终界面的HLL通量如下：,原理： 1） 双激波近似， 假设初始间断演化成两道激波，之间物理量为常数； 2）利用积分关系式，确定了两道激波（左行激波与右行激波）之间的物理量； 3） 利用积分关系式，计算了穿过x=0界面的流通量,优点： 计算简单 缺点： 没有考虑到接触间断。 该现象（初始间断产生两道激波、中间没有接触间断）不可能发生。 Riemann问题是三波问题，该近似为两波问题，假设过强。,Copyright by Li Xinliang,50,10.1.2 HLLC Riemann近似求解器

34、（Toro）,发展了HLL近似解，用三波模型来近似 （如图）,三波近似， 左、中、右波的波速为 ， ， 中间波为接触间断,HLLC 近似解,T 时刻的流动状态,Copyright by Li Xinliang,51,根据积分关系，可知,红色区域积分可得,蓝色区域积分可得,假设中间波为接触间断：,利用积分关系计算接触间断的速度及其左右的物理量,如果把ZL,ZR也当成未知数，可精确求解 （6方程6未知数，Riemann精确解），但计算复杂。,R-H关系式； 弱解定义式,含义： 控制体内质量的增加等于,控制体1 控制体2,Copyright by Li Xinliang,52,问题最终描述：求解超定

35、方程 4个未知数，6个方程,求解方程组：,其中：,假设波速ZR,ZL已知，未知数 4个,与精确Riemann解的区别： 本近似解假设左、右波速度已知，因此，未知数只有4个 ； 精确Riemann解认为ZL,ZR也是未知数，因此6个方程，6个未知数,6个方程，4个未知数，怎么解？,6个方程,求解思路 只使用其中的4个方程 方案（1） ： 利用的连续性方程及动量方程（共4个），解出全部4个未知数,求解方案不能太复杂，复杂度不能超过求解原6个未知数的方程（Riemann精确解）； 原则务必给出解析解,Copyright by Li Xinliang,53,近似解： 假设左、右波速已知，简化了计算 （

36、求解线性方程） 精确解： 左、右波速都是压力的函数，方程复杂（非线性），需迭代求解,线性方程 动量积分关系式,求出,解出,利用压力相等，求出解出速度,利用质量及动量方程,Copyright by Li Xinliang,54,最终HLLC （方案1）的表达式为：,HLLC 近似解 (方案1),Copyright by Li Xinliang,55,HLLC 近似解( 方案2),根据方案(1)算出u*,p*,近似解方案（2）与方案（1） 的区别 方案（1）只使用了质量及动量积分关系（未使用能量积分关系）； 方案（2）在最后一步也使用了能量积分关系,HLLC （方案2）给出的最终形式：,Copyr

37、ight by Li Xinliang,56,HLLC 近似解 （方案3）,假设压力时左、右两个解的平均,方案（1）令这两个解相等，以此计算Z* 方案（3）不要求他们相等，而以平均值代替,最终形式 （近似解方案3）：,6个方程，4个未知数，近似解方法很多，希望读者思考，开发出更好的方案,Copyright by Li Xinliang,57,10.1.3 波速的估算,HLL 及HLLC 均假设ZL, ZR已知，实际上它们仍需要估算,准确计算ZL,ZR，实际是计算Riemann精确解，计算量大,方法1： 直接估算,1a: 假设以声速传播 （Davis）,竟然假设激波以声速传播，太OUT了,小常识

38、： 激波的传播速度 激波相对于波前介质以超声速传播，相对于波后介质以亚声速传播 ； 弱激波（Ma趋近于1）以声速传播。,1b: 左、右两种状态声速的平均 (Davis, Einfeldt),要平均吗？ 用Roe平均吧，效果可好了。,激波速度介于波后（相对）声速与波前（相对）声速之间，平均是个好思路,Copyright by Li Xinliang,58,Roe 平均：,1c. Roe平均的修正 （Einfeldt）,方法2： 基于压力的波速估算法 （Toro）,已知中心区压力，容易计算波传播速度,中心区的估算：,Copyright by Li Xinliang,59,10.1.4 “扩展一维R

39、iemann问题”的HLLC近似解,特点： 切向速度v,w 的处理方法与被动标量相同； 被动标量场穿过激波值不变,Copyright by Li Xinliang,60,10.1.5 HLLC 方法计算通量的步骤,Step1： 估算左、中、右波的速度 （选用压力估算，还可采用直接估算）,Step 2: 计算HLLC型通量 （选用近似方案1， 还可选用其他两种方法）,左、右激波的速度,中间波（接触间断）的速度,中心区压力的估算,Copyright by Li Xinliang,61,在有限体积法中的具体使用方法：,控制体上积分，得到积分方程,j-1 j j+1 j+2,j-1/2 j+1/2,S

40、tep 1: 重构， 用 计算,相当于差分方法 的差分格式； 方法非常丰富,例如：,Step 2: 求通量,把 代人上页的公式，得到：,相当于差分法的通量分裂（属于FDS类），方法很多,Copyright by Li Xinliang,62,10.2 中心型有限体积法,中心型及迎风型方法的对比,为了使计算稳定，中心型格式通常添加人工粘性,Copyright by Li Xinliang,63,j-1 j j+1 j+2,j-1/2 j+1/2,控制体上积分,采用两侧的值重构j+1/2点的通量， 例如：,非MUSCL型,MUSCL型,注意： 不要把重构理解为插值，否则精度不会超过2阶,添加人工粘

41、性,二阶人工粘性,四阶人工粘性,Copyright by Li Xinliang,64,Jameson 人工粘性,j-1 j j+1 j+2,j-1/2 j+1/2,二阶和四阶人工粘性,光滑区： 间断区,经验常数，可调整人工粘性大小,Copyright by Li Xinliang,65,10. 3 具体算例讲解RAE2822翼型绕流,算例： RAE2822超临界翼型绕流的数值模拟,RAE2822翼型及压力测量装置,问题描述： 二维跨声速流动吹过RAE2822 翼型，试计算流场及翼型表面的压力分布。,流动参数： Mach 0.729, 攻角（AoA）= 2.31， Re= 6.5106 （基于

42、弦长及来流值）,背景： NASA NPARC Alliance Verification and Validation Archive 的标准算例（Bechmark）之一， 有可靠的实验数据、网格及多个数值结果。 成为目前CFD程序验证的标准模型之一。,网站：http:/www.grc.nasa.gov/WWW/wind/valid/validation.html （用Google 搜索 RAE2822即可）,Copyright by Li Xinliang,66,有限差分法的计算结果 5阶及7阶迎风差分+Steger-Warming FVS+ 3阶Runge-Kutta BL 湍流模型,网格

43、：NPARC 网站提供， 369*65,流动为湍流，二维计算需使用湍流模型； 湍流模型将在后续课程介绍，目前先按层流处理。并假设流动定常,Copyright by Li Xinliang,67,1. 控制方程、网格及边界条件 控制方程 无量纲N-S方程,无粘和粘性项,计算网格,边界条件： 外边界： 无穷远来流， 出流条件 固壁： 无滑移绝热边界,无穷远来流条件,出流条件,网格： 从NASA网站下载,Copyright by Li Xinliang,68,边界条件的处理方法：,方案1 ：严格按照特征方向 （精确，但复杂）,外边界1,外边界2,边界,对于每个边界单元上讨论特征方向， 在法方向局部一

44、维无粘化,边界处的物理量,在边界法向的坐标系中，控制方程为近似为“扩展的1维Euler方程”,求解 ， 计算边界值，单边差分（重构）,Copyright by Li Xinliang,69,计算模型：4个特征波，分别以速度 在该一维直线上传播，根据传播方向独立讨论。 若传播方向指向边界内，则给定流动参数 ；否则，用内部点计算,边界,方案2： 近似处理,外边界1 ： 流场远离翼型，扰动很弱，认为 （无论亚、超均按超音速处理） 外边界2： 按超音速出流处理,外边界1,外边界2,近似 模型：,对于内流，有问题,固壁边界条件,物理条件：,补充计算条件（常用）：,Copyright by Li Xinl

45、iang,70,2. 空间离散,无粘通量,粘性通量,控制体上积分：,注： 选取控制体的方法很多,本计算采用的控制体节点中心型,控制体的多种选择方式,Copyright by Li Xinliang,71,2. 空间离散,无粘通量,粘性通量,（1） 无粘通量的计算,旋转到xy坐标系,（迎风）重构控制体边界处的左、右值,线性重构：平均值 =中心点值,求解“扩展”Riemann问题，得到控制体边界m处的通量,建议：使用HLL 或HLLC型通量 （公式见本PPT 12、17页）,q 是x轴与x轴的夹角,Copyright by Li Xinliang,72,具体步骤: 无粘通量,讨论以i,j 为中心的

46、控制体,在 所在的控制体表面（记为m）上,1） 重构自变量的左、右值,2）进行坐标旋转，得到坐标系（x,y）下的自变量,q 是x轴与x轴的夹角,3） 求解Riemann问题，获得界面m上的通量,m,使用HLL 或HLLC型通量， 见12，17页,4） 在其他3个控制面（图中虚线所示）上，同样求出通量,5）4个面的通量求和，得到该控制体的总无粘通量,Copyright by Li Xinliang,73,（2） 粘性通量的计算 采用中心格式,可用周围角点的平均,关键问题： 如何计算 ？,方案1： 借鉴有限差分法，引入坐标变换,具体表达式见: 任玉新等 计算流体力学基础 130-131页及127页,转化到计算坐标系,通常，令,计算坐标 实际上就是下标（i,j）,（1）,Copyright by Li Xinliang,74,于是：,坐标变换系数的计算：,控制体的体积,半点上的Jocabian系数（体积加权）插值得到,完全的有限差分法， 坐标变换Jocabian系数的计算是全局的； 本方法中坐标变换Jocabian系数的计算是局部的，只用到周围几个点的坐标。Jocab