计算流体力学中科院力学所第9讲-有限体积法.ppt
计算流体力学讲义 第九讲 有限体积法（1）李新亮；力学所主楼219；82543801,知识点：,1,讲义、课件上传至(流体中文网）-“流体论坛”-“CFD基础理论”讲课录像及讲义上传至网盘,Copyright by Li Xinliang,有限体积法的基本概念 重构和反演迎风型有限体积法Riemann求解器；Roe格式的新理解：近似Riemann解多维迎风型有限体积法坐标旋转,Copyright by Li Xinliang,2,知识回顾,1.差分方法的基本概念：,差分格式、修正方程、相容性、收敛性、稳定性、LAX等价定理,2.精度分析、稳定性分析与分辨率分析（修正波数）,Taylor分析,Fourier分析,修正波数,激波捕捉格式 GVC,NND,Roe,Godnov,MUSCL,TVD,WENO,Euler(N-S)方程的通量分裂 逐点分裂、特征投影分裂（建议使用Roe平均）,5.隐格式求解的LU-SGS方法,要点：a.引入差量，方程线性化 b.单边差分，隐式代数方程显式（推进）化,以一维为例，多维可直接推广,方法1：直接隐式离散,直接求解,非线性方程组，计算量大,方法2,差量化,线性化,已知项,线化微分方程,Copyright by Li Xinliang,3,Copyright by Li Xinliang,4,求解思路：如果直接离散，得到线性代数方程组，仍需求解，计算量大（多维情况）,如果能单侧差分就好解了！,多对角方程组，不好解（多维情况）,中心（双侧）离散,如果单侧离散,单侧离散，可推进求解，免受解方程组之苦。真简单,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 有限体积法入门,有限体积法主要优势：处理复杂网格,差分法处理复杂外形 坐标变换,坐标变换函数必须足够光滑 否则损失精度,实际问题：外形复杂，光滑的结构网格生成困难,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,积分（精确）,重构（Reconstruction),有限差分法的离散：数值微分过程有限体积法的离散：数值积分过程,积分方程,离散化,反演（evolution),(1)重构过程,A.零阶重构，假设分片常数,j-1,B.线性重构，假设分片线性函数,零阶重构与一阶重构示意图,j,j+1,or,or,或其他方法,C.更高阶的重构例如:分片二次函数（PPM）,WENO等,重构是有限体积的空间离散化过程，有多种方法,Copyright by Li Xinliang,9,（2）演化过程（以线性方程为例）,需要得知时间演化信息，通常利用特征方程,若采用零阶重构：,则：,假设时间步长足够小,则方程为：,等价于一阶迎风差分,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点的值，仍需要使用周围点 进行插值,通常无法精确计算，可采用近似值 代替,等价于二阶中心差分,半离散,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 Xinliang,13,2.分裂方法（1）:FVS方法（流通矢量分裂 逐点分裂）,具体方法：Steger-Warming 分裂 Lax-Friedrichs分裂 Van Leer分裂：Liou-Steffen分裂：（压力项与其他项分开，AUSM类格式的基础）,根据当地Mach数分裂,保证 的Jocabian阵特征值为正，的为负,正通量：向左偏斜重构；负通量：向右偏斜重构 偏重向上游,与迎风差分法类似：网格基（或权重）偏重上游,差分、有限体积都可使用,一个参数，反映全部特征,Copyright by Li Xinliang,14,小知识：Liou-Steffen分裂,对流项,压力项,思路：决定特征的关键参数 当地Mach数,超音速，x-方向,超音速，x+方向,因此，对Mach数进行分裂更为简洁！,显然：,参考文献：Liou:Ten Years in the making AUSM family,NASA TM-2001-210977,类似 Van Leer分裂，但压力单独处理,M,保证光滑过渡,M=1,Copyright by Li Xinliang,15,(3)FDS 方法（通量差分分裂特征投影分裂）,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,差分法与有限体积法区别与联系（二阶迎风+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点的值含义不同！,用周围几个点的值 计算 的过程称为“重构”，不能理解为用 来插值,记号 确实容易混淆，让人容易联想起。记为 更好些,否则，最高只能达到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/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平均，刚好是左右两点间的平均增长率，实现了常系数线性化。,常系数双曲方程组，易解！,思路：用平均增长率矩阵 取代瞬时增长率矩阵A,不但实现了线性化，而且实现了常系数化。利用二次齐函数的性质，可找到了Roe点（即Roe平均点），该点处的增长率刚好等于平均增长率。,Roe平均,常系数化,线性化,常系数线性单波方程的Riemann问题太简单了,21,常系数方程组的Riemann问题,解耦了的单波方程，有精确解,初值,Copyright by Li Xinliang,22,解为,线性化条件（并利用齐函数性质）,与差分法的Roe格式相同,还有各种其他类型的近似Riemann解（今后介绍）,Copyright by Li Xinliang,23,9.1.3 多维问题的有限体积法,二维问题,一维Riemann问题，坐标选取不当，变为“二维”Riemann问题,x,y,差分法：独立计算只考虑各自的特征方向,由于非线性，实际（二维）特征方向并非x,y方向特征量的线性组合。特征方向计算不严格，带来误差,差分方法：多维情况，特征理论复杂，通常x,y方向独立计算,转化为 x方向与y方向的两个一维问题,逐点分裂,特征投影分裂,完全按照一维情况独立处理,局部坐标旋转？差分算法设计造成局部旋转困难,差分法的多维处理方法,1.小知识：差分方法如何处理高维问题的？优缺点？,优点：简单缺点：特征方向计算不准,Copyright by Li 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,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,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个面（线），表示该面的面积（长度）,于是，问题转化为求控制面上的,这个量有两个重构方案,方法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,可用精确Riemann解，也可用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,知识回顾：Riemann问题精确解,Riemann问题,问题描述：初始时刻，物理量分布存在单个间断；间断两侧物理量为常数。,求解思路：采用积分方程 单个间断，且间断两侧物理量为常数情况下：积分方程转化为代数方程,代数方程：质量、动量、能量守恒,计算出，将 与这三个值进行比较，判断会产生的情况。具体见下图：,Copyright by Li Xinliang,32,Riemann问题的具体计算步骤（全流场）,1.判断可能会出现的情况（五种情形之一）,a.定义函数,b.进行判断,情况5,情况4,情况3,情况1,情况5,情况4,情况2,情况1,单调增函数，性质很好,Copyright by Li 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.计算稀疏波区域的值（如果有稀疏波的话）,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 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 Xinliang,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）绘制出t=0.14时刻x轴（垂直于初始间断面，见右图）上的密度、沿x方向的速度及压力的一维分布，并同精确解进行比较。,x,y,
