(7.6.1)--6.8Jameson中心有限体积法(1)-讲义.pdf
计算流体力学基础讲义 1 第六章 简单的计算流体力学技术入门知识 6 6.8 8 JamesonJameson 中心有限体积法 首先我们给出本节内容的原始文献和参考教材。Jameson，Schmidt 和 Turkel 在 1981 年发表了著名的论文：Numerical Solution of the Euler Equation by Finite Volume Methods with Runge-Kutta Time Stepping Schemes.AIAA Paper 81-1259,1981.该论文开创性地提出了有限体积法框架下的一种高效鲁棒的数值模拟方法。本节内容主要参考德国宇航院 Kroll 教授采用 Jameson 方法的公开发行报告 DFVLR-FB 87-41:N.Kroll,Solution of two-dimensional Euler equations experience with a finite volume code.部分评述参考自北京航空航天大学阎超教授的教材:计算流体力学方法及应用（第 6 章），阎超编著，北京航空航天大学出版社。6.8.1 引言 CFD 一直沿着两个大不相同的方向发展，这个大不相同主要是指处理间断方法的不同：一个是认为流动处处“间断”，也就是把离散后的每一个网格同别的网格的交界面视为“间断点”，即构成“Riemann”问题，并采用不同的方法求解这些 Riemann 问题，形成了不同的计算格式，如Godunov 格式、Van Leer 格式、Roe 格式等。我们将在第七章讨论这一方向的具体处理方法。另一个则相反，是处处采用连续性处理方法，把间断处理成梯度很大的连续性特殊情况，利用格式的数值耗散效应把间断进行“光滑”处理，形成数值求解方法，此处理方法的典型代表就是“人工粘性”方法。6.7 节我们介绍了属于此类方法的 MacCormack 格式的应用。本节我们讨论该类方法的最典型代表：Jameson 中心格式有限体积法。1981 年，Jameson，Schmidt 和 Turkel 提出的显式中心格式有限体积法被称为 CFD 发展史上的一个突破，为纪念这三位科学家的贡献，中心有限体积格式也被称为 JST 格式。该方法经过不断改进和完善，成为目前 CFD 的主要计算方法之一，已经在航空领域亚、跨、超声速流动计算和叶轮机械等领域取得了巨大成功。计算流体力学基础讲义 2 有限体积方法求解的是积分形式的控制方程，即针对选定的有限控制体模型应用三大定律得到的控制方程。按照不同的离散方式，即待求解流动参数的存储位置，有限体积法分为格心格式和格点格式两大类。下图是结构化网格体系下，格心格式有限控制体的示意图。有限控制体单元如ABCD 所示，待求解未知数的位置在 ABCD 的形心。作为对比，下图给出结构化网格体系下格点格式有限控制体的示意图。有限控制体单元 ABCD由四个小网格单元组成，未知数的位置在网格点上。ABCDii+1j+1ji+1,ji,j+1i,ji-1,j格格心格式有限控制体示意图心格式有限控制体示意图Ui,ji,j-1ABCDii+1j+1j格点格式有限控制体示意图格点格式有限控制体示意图i-1j-1Ui,j计算流体力学基础讲义 3 针对守恒积分形式的控制方程组直接进行空间离散时，会很自然地选择有限体积法。由于有限体积法是将控制方程直接在物理空间进行离散，因此避免了从物理空间向计算空间的转换，并使控制方程离散时的物理概念简单明了，更能体现数值模拟的特点。此外有限体积法还具有良好的守恒性，以及适合具有复杂几何边界的流场问题求解等优点。这里我们再次强调一下，有限体积法是在物理空间直接离散控制方程，无需像有限差分法那样必须转换到计算平面的均匀网格，因此，可以采用非结构网格，也就是说，有限控制体单元可以是任意形状的。如下图给出的非结构三角形网格，我们可以直接在物理平面针对一个三角形单元进行数值离散，将离散未知数取在三角形单元的形心，如右边这个放大的三角形 ABC 所示，即格心格式离散。也可以将离散未知数取在网格顶点，有限控制体单元由围绕此顶点的 6 个三角形单元组成，即格点格式的有限体积单元如右边放大的多边形 ABCDEF 所示。下面我们给出本节的路线图，以翼型定常绕流的 Euler 方程数值求解为例，第 2 小节我们给出积分形式的控制方程，第 3 小节学习中心格式的有限体积空间离散。第 4 小节学习如何处理边界条件，第 5 小节介绍人工耗散项，第 6 小节学习时间推进格式，第 7 小节学习加速收敛措施。ABC格心格式CEABDF格点格式非结构三角形网格示意图（格心格式和格点格式）计算流体力学基础讲义 4 图 6.8.1 有限体积方法的路线图 6.8.2 积分形式控制方程 在笛卡儿坐标系下，二维、非定常、可压缩 Euler 方程的守恒微分形式控制方程为（6.8.1）式：UFG0 txy+=(6.8.1)其中解矢量 U，通量项 F 和 G 如（6.8.2）所示:UuvE=；222vF G+()()22uupuvvuv pVVeupuevpv+=+；(6.8.2)大写 E 代表单位质量的内能加动能：222()=2(1)2VpuvEe+=+(6.8.3)取任意控制体 V，对（6.8.1）进行体积分，可得(6.8.4)式：计算流体力学基础讲义 5 UFG0 VdVtxy+=（）(6.8.4)如果控制体 V 固定于空间，则U/t 的体积分中的微分符号/t 可以移到体积分符号之外，如（6.8.5）所示：UUVVdVdVtt=(6.8.5)对通量项应用高斯定理，（FGxy+）的体积分可以转化为沿控制体封闭表面的面积分：FGF VVdVndSxy+=（）（6.8.6）其中F为 x 方向通量 F 和 y 方向通量 G 的矢量相加：F=FGij+我们称为F通量张量。将 F 和 G 的具体表达式代入，可得通量张量表达式（6.8.7）：2222+()()22uvupuvFijvuvpVVeupuevpv+=+（6.8.7）即：2222+()+F=()()()22uivjup iuvjvuivp jVVeu ipuievpv j+（6.8.8）引入速度矢量 qu iv j=+，pHE=+，可得通量张量F：+F=+quq pivq pjHq （6.8.9）计算流体力学基础讲义 6 这样，我们就推导出了积分形式的二维、非定常、可压缩 Euler 方程（6.8.10）式：UF0VVdVndSt+（6.8.10）为使方程封闭，还需要（6.8.12）状态方程和内能计算公式（6.8.13）：pRT=（6.8.12）vec T=（6.8.13）6.8.3 有限体积空间离散（格心格式-cell centered）下面我们进入第 3 小节，介绍积分形式控制方程的空间离散。以二维翼型绕流问题为例，我们可以生成如图 6.8.2 所示的 O 型结构化网格，如果采用最常用的格心格式有限体积法，则每一个离散的有限控制体单元是相邻等和等线组成的四边形，如图中的四边形 ABCD，未知数取在四边形的形心。图 6.8.2 有限体积空间离散的绕翼型的网格示意图 对于格心格式有限体积法，待求解未知数是取在有限控制体的形心上的，即解矢量 U 以及位置，以及相应的有限控制体单元可用一对（i，j）标号指定。另外，每个有限控制体单元的顶点可以用一对标号（i，j）来指定。采用格心格式有限体积法需要注意的是，未知数点的坐标位置和网格点坐标位置不同。xyABCD计算流体力学基础讲义 7 图 6.8.2A 格心格式有限控制体单元未知数取值位置示意图 取(i,j)单元（ABCD）为有限控制体，顶点 A 的坐标为（1,1,ijijxy+，）；顶点 B 的坐标为（1,11,1+ijijxy+，）；顶点 C 的坐标为（,1,1i ji jxy+，）；顶点 D 的坐标为（,i ji jxy，）；图 6.8.2B 格心格式有限控制体的外法向矢量示意图 这样，每个表面的外法向及其面积大小都可以确定出来，各表面外法向方向如红色箭头所示。控制体的体积大小，对于二维情况就是四边形 ABCD 的面积，由公式（6.8.14）计算求得：,1,1,11,1,1,11,0.5()()()()i jiji ji jijiji ji jijVxxyyyyxx+=（6.8.14）ABCDii+1j+1ji+1,ji,j+1i,ji-1,j格心格式有限控制格心格式有限控制体体单元单元i,j-1MNRSUi，jUi，j+1Ui-1，jUi+1，jUi，j-1(i,j)(x,y)i+1,j(x,y)i+1,j+1(x,y)i,j+1(x,y)i,j1+,2()ijn dSABCD1-,2()ijn dS1,-2()i jn dS1,+2()i jn dS计算流体力学基础讲义 8 令S矢量的大小表示控制体表面的面积，方向为控制体表面沿 i 和 j 增长方向的法向方向，则由下面左、右图的对比，可以看出，CD 面的外法向方向的ndS的面积分为负的1-2ijS，。DA 面的外法向方向的ndS面积分为负的1-2ijS，。图 6.8.2C 有限控制体的面矢量S的定义 容易求得，S在有限控制体单元各表面的计算公式，如公式（6.8.15a、b、c、d）。1,+1,+1-,2()()i jiji ji jijSyyixxj=+，（6.8.15a）1+1,+11+1,+1,+1+,2()()ijijijijijSyyixxj+=+，（6.8.15b）1,1,1,2()()i jijiji ji jSyyixxj+=+（6.8.15c）1,11,11,1,1,2()()i jijiji ji jSyyixxj+=+（6.8.15d）控制体单元的体积和各表面面积值及矢量方向值确定了，就可以进行控制方程的离散了。首先我们看一下体积分项的离散表达，即UVdVt的数值离散表达。,Ui j位于控制体单元的中心，当有限体积单元趋于很小时，可近似认为,Ui j就是该单元的体积平均值，即（6.8.16）式：,1UUi,ji ji jVdVV=（6.8.16）进一步可写成（6.8.17）式：,U=Ui,ji ji jVdV V （6.8.17）(i,j)(x,y)i+1,j(x,y)i+1,j+1(x,y)i,j+1(x,y)i,j1,2i jS+1+,2ijS1-,2ijS1,-2i jSABCD(i,j)(x,y)i+1,j(x,y)i+1,j+1(x,y)i,j+1(x,y)i,j1+,2()ijn dSABCD1-,2()ijn dS1,-2()i jn dS1,+2()i jn dS计算流体力学基础讲义 9 这样，U 的体积分的时间偏导数就等于,Ui ji jV对时间求偏导：,U(U)i ji ji jVdVVtt=(6.8.18)因为有限控制体为固定于空间不随时间变化的，所以有：,U(U)=(U)i ji ji ji ji jVddVVVttdt=(6.8.19)由（6.8.19）可见，关于时间的偏微分项通过数值离散变为常微分项，说明积分形式的控制方程的时间离散和空间离散完全分开，这种方法被称为半离散方法。对有限控制体单元进行离散(半离散），得(6.8.20)式:(),0Ui ji ji jdVQdt+=(6.8.20)由（6.8.20）式可以看出，该方程代表了关于解矢量 U 的常微分方程，只有空间方向离散了，这是我们称其为半离散的原因。(6.8.20)式中矢量,i jQ是控制体表面相关项积分的数值表达：1111,2222(F)-(F)(F)-(F)ijijiji ji jQSSSS+=+(6.8.21)（6.8.21）式各通量项的计算是类似的，例如：通量1,2(F)i jS的四个分量如公式（6.8.22）所示：()111,22211111,222221,211111,22222111,222()(F)=()i ji ji ji ji ji ji ji ji ji ji ji ji ji ji ji ji jqSuqSpSiSvqSpSjHqS+（6.8.22）图 6.8.2D 给出了离散控制方程中涉及到的变量位置。解矢量 U 的离散点取在实心圆点处，通量的面积分计算在点处。计算流体力学基础讲义 10 图 6.8.2D 离散解矢量 U 的位置和通量计算时变量的位置图示 因此，计算控制体表面F时，其对应的解矢量 U 的值（分数下标的物理量1,2Uij+1-2Uij，等）需要取相邻控制体单元的平均值，Jameson 直接用如公式（6.8.23）和（6.8.24）所示的算术平均值：1,2,1,1FF(UU)2iji jij+=+（6.8.23）1,2,11FF(UU)2i ji ji j=+（6.8.24）本次课就到这里，谢谢大家。,i j（）1,ij+（）-1,ij（）1-,2ij（）,-1i j（）,1i j+（）1,2+ij（）12i j+（,）1-2i j（,）解解矢量矢量U U的位置的位置通量通量的位置的位置F n
