二维Euler方程的Jameson求解方法

基于非结构网格二维Euler方程的Jameson求解方法姓名：陈皓学号：0501211日期：2006年5月21目录中英文摘要．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．3符号说明．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．4引言．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．5方法论述．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．62.1控制方程．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．62.2空间离散．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．72.3人工耗散项．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．82.4时间离散．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．112.5边界条件．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．11算例分析．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．123.1马赫数M=0.84,攻角α=0°．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．123.2马赫数M=0.8,攻角α=1.25°．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．123.3马赫数M=2.0,攻角α=10°．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．12附图．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．13结论．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．21参考文献．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．．21摘要本文介绍了基于非结构网格求解二维Euler方程的Jameson中心格式方法。在空间离散上采用的是有限体积法，时间上采用的是四步显式Runge－Kutta迭代求得最后的定常解。人工耗散项为守恒变量的二阶和四阶差分项。边界条件采用的是无反射边界条件，并采用当地时间步长进行加速收敛。最后对翼型绕流进行了数值模拟，结果比较令人满意。关键词：Jameson中心格式，Euler方程，非结构网格，有限体积法AbstractAmethodforthenumericalsolutionofthetwo-dimensionalEulerequationsonunstructuredgridshasbeendeveloped.Thecell-centredsymmetricfinite-volumespatialdiscretisationisappliedinageneralformulation.Theintegrationintime,toasteady-statesolution,isperformedusinganexplicit,four-stageRunge－Kuttaprocedure.Theartificialdissipationisconstructedasablendingofsecondandfourthdifferencesoftheconservedvariables.Andintheboundary,thereisnoneoftheoutgoingwavesarereflectedbackintothecomputationaldomain.Anaccelerationtechniquecalledlocaltimesteppingisused.Atlast,standardtestcasesforbothsubsonicandsupersonicflowshavebeenusedtovalidatethemethod.Keyword:Jamesonmethod,Eulerequations,unstructuredgrid，finite-volume符号说明CL,CD升力、阻力系数C音速Dk离散的耗散通量积分d(2),d(4)耗散通量E单位体积总能量F，G流场通量H单位体积总焓k(2),k(4)离散通量系数P静压Qk离散通量积分Rk残值SΩ域的边界t时间U,V笛卡尔坐标系下的速度分量W守恒变量矢量x,y笛卡尔坐标系α比例因子γ比热比△t时间步长隐式残值光顺系数(2),(4)离散项系数v激波传感器ρ密度Ω积分域第一章引言在目前，实际流动问题的解多是通过求解Navier-Stokes方程来获得的，因为Navier-Stokes方程能够反映流体的质量守恒、动量守恒和能量守恒的规律。而由于本文针对的是二维无粘可压缩流，所以可以将Navier-Stokes方程简化为二维Euler方程，因此本文研究的是Euler方程的求解。所采用的Euler方程的求解方法是Jameson求解方法，这种方法有两个特点：1、在空间的离散上，内场边的通量是通过将左右单元中心处的通量值取平均获得的；2、为了更好的捕捉间断点和提高格式的稳定性，在计算通量时加入了人工耗散项，这人工耗散项是由守恒变量的二阶和四阶差分项组成的。本文的计算是基于非结构网格进行的，在时间上则通过四步Runge-Kutta迭代来获得定常解，边界条件则是采用无反射边界条件，并且采用了加速收敛技术。最后，分别对NACA0012翼型的跨音速和超音速绕流进行了数值模拟。第二章方法论述本章具体介绍了二维Euler方程的Jameson求解方法，Jameson格式是一种有限体积的求解方法。这种方法就是对每个网格单元进行积分，然后通过计算得到流场结果的方法。为此，本章首先从二维Euler方程的控制方程着手，对其进行积分，得到积分方程。然后分别介绍了时间和空间上的离散，人工耗散项的计算，以及边界的处理，最后介绍了加速收敛技术。2.1控制方程由第一章的介绍可知控制方程是Euler方程，对于二维的控制体的积分域为Ω，边界为S。则Euler方程可以改写为：（2-1）其中：X,Y是笛卡尔坐标系，W为守恒变量矢量：（2-2）F,G为流量矢量：，（2-3），P，H和E分别为密度，压强，单位体积的总焓和单位体积的总能。U,V分别为速度在X,Y方向的分量。且对于理想气体，总焓和总能可以表示为：（2-4）（2-5）2.2空间离散将计算域划分为有限个互不重叠的单元，并将积分守恒方程应用于每个单元，由于各个单元的面积不随时间变化，所以可将（2-1）式改写为：（2-6）其中：Ω和S均是指每个单元，W是单元的平均值。该方程对空间的离散是应用有限体积方法，然后对得到的半离散方程（时间的离散并未完成）按时间步长推进从而得到精确解。在进行时间离散之前，对方程（2-6）空间离散得到一个关于时间的常微分方程：（2-7）其中：Ωk是第k个单元的面积，Wk是守恒变量矢量，是通量积分的离散近似值，可写成：（2-8）其中：,(2-9)这是第k个单元所有边界的累计总和。采用中心格式的有限体积法，把守恒通量都控制在单元中心上，流体流经第i条边的值由两个相邻单元中心点（k和p）的平均值决定：（2-10）下标的意义如下图：图1单元中心，边和网格顶点设∶（2-11）则Euler方程对k单元的离散解，方程（2-7）可写为：（2-12）上述的一般形式既能用于结构网格，也可以用于非结构网格，唯一的不同点是：前一种情形中方程（2-12）的形式更紧凑，因为前者在采用曲线坐标时，每个单元由一个二阶矩阵来定义。非结构网格中，相邻的网格点和网格单元间并没有对应关系，所以非结构网格需要一个连接矩阵，这个连接矩阵中储存了所有的必须信息（例如网格顶点的坐标，边数，单元数等等）。2.3人工耗散项如上所述的中心格式是不包含耗散项的，所以任何误差（离散误差，循环误差等等）都不会衰减，最后的定常解可能会出现振荡。为了减少这些振荡，在方程（2-7）的右边加了人工耗散项，对于k单元，方程变为：（2-13）在目前的工作中，采用了Jameson的方法，耗散项Dk取守恒变量Wk的二阶和四阶差分项的组合。含四阶差分项的部分加在流场域中的平滑部分，但是在激波区被关掉。此时，打开含二阶差分项的部分，从而减少激波区的振荡，这个值可以非常大。这一开关函数是由基于当地压强的二阶差分的激波感应器控制的。在非结构网格中，Dk可以表示为：(2-14)其中二阶耗散项和四阶耗散项都是单元边界上的耗散通量之和：，(2-15)表示第i条边上的二阶耗散通量，表示四阶耗散通量：(2-16)(2-17)其中指数I表示单元k和p的分界边，被定义为：(2-18)ε(2)与ε(4)为二阶和四阶耗散项的自适应系数。下式表示了当地压强的二阶差分项：(2-19)由此确定的自适应系数为：（2-20）其中k(2)和k(4)是两个经验常数，一般取1/256<k(4)<1/32和1/2<k(2)<1.0。由上述公式可以看出，开关函数是受压力梯度控制的。虽然格式的精确性被保持下来，但是我们对上述的人工粘性项并不是完全满意。在上述非结构网格中，人工粘性项的大小对于平滑区域太高了，而对于大梯度的区域又太低了。在目前的工作中，我们采用了如下的简单方法：构造每一条边的激波感应器和比例因子，仅仅使用其两个相邻单元k和p的流动变量：(2-21)由此确定的自适应系数变为：(2-22)比例因子αi取沿单元边界的Jacobian矩阵F/W和G/W的最大特征值：(2-23)其中U，V和c是边界上的平均值，c取当地音速。在人工粘性项中，二阶耗散项的作用是抑制解在激波附近的振荡,在流场中压力梯度不大的区域内此项作用很小。四阶耗散项的作用是抑制高频振荡，并使解趋于定态。总之，在加入人工耗散项后，可以更好的捕捉激波间断点和使格式稳定。2.4时间离散定常解是由对常微分方程进行时间积分获得的，此方程可以写成：（2-24）方程的右侧的值表示每个单元k中心的残值：（2-25）对于方程（2-24）的积分可以采用显式的四步格式来完成，由于对于定常解时间的精确并不重要，所以选择此格式只是因为它的稳定和衰减的特性。目前采用的是以下的格式：W(0)=WnW(m)=W(0)+αm△tR(m-1)form=1to4Wn+1=W(4)(2-26)式中n是当时的时间步数，n+1是新的时间步数：（2-27）其中的系数是：,,,（2-28）为了减少计算时间，只在第一步时计算耗散项D，然后在接下来的各步中D是一个常值。对于一个标量模型方程，这种做法改变了格式的稳定区，但是精度和收敛特性被保留了下来。以上格式在应用于非定常Euler方程时，能保持稳定的CFL数最大可以取到2。这个显式格式最大的缺点是最大的时间步长受到限制，因为稳定域受到限制。而且，对于多维的方程组，最大时间步长只可以用近似的方法获得。对于固定形状的网格，采用以下的表达式：（2-29）2.5边界条件在流场中，需要考虑的边界条件包括远场边界条件和物面边界条件。在这一小节中将分别介绍这两种边界条件。2.5.1物面边界条件对于无粘流体，在物面边界上，应该满足流动方向与物面相切，即物面的法向速度分量为零，可知。这就是物面无穿透边界条件。因此，X,Y动量方向的通量不是零，而与压强有关。2.5.2远场边界条件用有限体积法进行空间离散时，只能取一个有限远的边界作为远场边界，而在绕物体的实际流动中并不存在这个边界，物体所产生的扰动可以认为被传到无穷远而没有反射，即无反射边界条件。这里，采用一维Riemann不变量来处理远场边界条件。令，分别为边界外法向速度和切向速度，根据特征线理论，Riemann不变量可写为：（2-30）Riemann不变量的取值可分为以下四种情况：1)亚音速入流（－a<<0）R1,R2,R3取来流值，R4取内场值。2)亚音速出流（0<<a）R1,R2,R4取内场值，R3取来流值。3)超音速入流（<-a）所有物理量都取来流值。4)超音速出流（>a）所有物理量都取内场值。其他所有物理量都由以上四个不变量求出。第三章算例分析本文所选择的进行算例分析的是NACA0012翼型，并采用以三角形为基本单元的非结构网格作为计算网格，图2为计算所用的非结构网格示意图。网格总节点数为5306，总单元数为10382,总边数为15688。马赫数M=0.84,攻角α=0°图3～6给出了这个算例的计算结果（其中包括：上下表面压力分布图，以及压力、密度和马赫数的分布云图）。图3给出了翼型表面的压力分布曲线。从图上可以看出，在翼型65%弦长处，上下表面出现了一对上下对称的激波。计算所得的:-0.0001，0.0379。3.2马赫数M=0.8,攻角α=1.25°图7～10给出了此算例的计算结果（同上）。图7给出了翼型表面的压力分布曲线。从图上可以明显看出，在翼型的60%弦长附近的上表面处出现一道激波，在约35%弦长附近的下表面上又出现一道弱激波。计算所得的：0.2764，0.0159。3.3马赫数M=2.0,攻角α=10°图11～14给出了此算例的计算结果（同上）。图11给出了翼型表面的压力分布曲线。从图上可以明显看出，在翼型前缘上下表面出现了一对上下不对称的脱体激波，并且在翼型的头部前端出现了跨音速区。计算所得的：0.3854，0.0913。附图：图2计算所用的非结构网格示意图图3NACA0012翼型表面压力分布,Ma=0.84,alpha=0图4NACA0012绕流密度云图,Ma=0.84,Alpha=0图5NACA0012绕流压力云图,Ma=0.84,Alpha=0图6NACA0012绕流马赫数图,Ma=0.84,Alpha=0图7NACA0012绕流残值收敛图,Ma=0.84,Alp