基于非结构网格二维欧拉方程的求解方法

基于非结构网格
二维Euler方程的Jameson
求解方法
姓名：
学号：
摘要
本文介绍了基于CFD理论的求解二维可压缩流Euler方程的Jameson中心格式方法。

在空间离散上采用的是有限体积法，时间上采用的是四步显式Runge －Kutta迭代求解。

人工耗散项为守恒变量的二阶和四阶差分项。

边界条件采用的是无反射边界条件，并采用当地时间步长进行加速收敛。

最后对NACA0012翼型划分了三角形，并应用本文程序进行数值模拟，结果较为理想。

关键字：CFD，Jameson中心格式，Euler方程，有限体积法
Abstract
A method for the numerical solution of the two-dimensional Euler equations has been developed. The cell-centred symmetric finite-volume spatial discretisation is applied in a general formulation. The integration in time, to a steady-state solution, is performed using an explicit, four-stage Runge－Kutta procedure. The artificial dissipation is constructed as a blending of second and fourth differences of the conserved variables. And in the boundary, there is none of the outgoing waves are reflected back into the computational domain. An acceleration technique called local time stepping is used. At last, standard test cases for both subsonic and supersonic flows have been used to validate the method.
Key words：CFD, Jameson method,Euler equations, finite-volume
第一章引言
在工程应用的推动下，计算流体力学随着计算机技术的发展和计算格式的不断更新而迅猛发展。

在航空航天领域，CFD已经与地面试验和飞行试验共同构成了飞行器设计，飞行器性能分析和飞行器空气动力学设计的三大工具。

由此可见，CFD数值模拟对于实际问题的求解是十分重要的。

在目前，实际流动问题的解多是通过求解Navier-Stokes方程来获得的，因为Navier-Stokes方程能够反映流体的质量守恒、动量守恒和能量守恒的规律。

而由于本文针对的是二维无粘可压缩流，所以可以将Navier-Stokes方程简化为二维Euler方程，因此本文研究的是Euler方程的求解。

计算流体力学常用的求解方法有有限体积法、有限单元法和差分方法。

本文所采用的求解方法是Jameson求解方法，该求解方法是中心格式的有限体积法。

这种方法有两个特点：1、在空间的离散上，内场边的通量是通过将左右单元中心处的通量值取平均获得的；2、为了更好的捕捉间断点和提高格式的稳定性，在计算通量时加入了人工耗散项，这人工耗散项是由守恒变量的二阶和四阶差分项组成的。

Jameson格式可以方便地求解以非结构存储方式存储的网格。

实际求解过程中，对初始化的全流场进行直接的时间积分直至收敛以获得稳定的数值解。

本文的计算是基于非结构网格进行的，在时间上则通过四步Runge-Kutta迭代来获得定常解，边界条件则是采用无反射边界条件，并且采用了加速收敛技术。

最后，分别对NACA0012翼型的跨音速和超音速绕流进行了数值模拟。

本文第二章是对Jameson求解方法的理论推导，第三章是提供算例和结果进行验证，第四章是总结与展望。

第二章 方法论述
本章具体介绍了二维Euler 方程的Jameson 求解方法，Jameson 格式是一种有限体积的求解方法。

这种方法就是对每个网格单元进行积分，然后通过计算得到流场结果的方法。

为此，本章首先从二维Euler 方程的控制方程着手，对其进行积分，得到积分方程。

然后分别介绍了时间和空间上的离散，人工耗散项的计算，以及边界的处理。

2.1 控制方程
由第一章的介绍可知控制方程是Euler 方程，对于二维的控制体的积分域为Ω，边界为S 。

则Euler 方程可以改写为：
⎰⎰Ω=-+Ω∂∂
s
Gdx Fdy Wd t 0)( （2-1） 其中：X,Y 是笛卡尔坐标系，W 为守恒变量矢量：
⎥⎥⎥⎥
⎦
⎤⎢⎢⎢⎢⎣⎡=E v u W ρρρρ （2-2） F,G 为流量矢量：
⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡+=UH UV P U U F ρρρρ2 ⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎣⎡+=VH P V UV V G ρρρρ2 （2-3） ρ，P ，H 和 E 分别为密度，压强，单位体积的总焓和单位体积的总能。

U,V 分别为速度在X,Y 方向的分量。

且对于理想气体，总焓和总能可以表示为：
2/)()1/(22V U P E ++-=ργρ （2-4）
P E H +=ρρ （2-5）
将计算域划分为有限个互不重叠的单元，并将积分守恒方程应用于每个单元，由于各个单元的面积不随时间变化，所以可将（2-1）式改写为：
⎰⎰Ω
Ω--=∂∂d Gdx Fdy t W S )( （2-6） 该方程对空间的离散是应用有限体积方法，然后对得到的半离散方程（时间的离散并未完成）按时间步长推进从而得到精确解。

2.2 空间离散
在进行时间离散之前，对方程（2-6）空间离散得到一个关于时间的常微分方程：
k k k Q dt
dW Ω-=/ （2-7） 其中：Ωk 是第k 个单元的面积，W k 是守恒变量矢量，k Q 是通量积分的离散近似值，可写成：
∑=∆-∆=
kedges
i i k x G y F Q 1)( （2-8）
其中： a b i x x x -=∆a b i y y y -=∆ （2-9）
这是第k 个单元所有边界的累计总和。

采用中心格式的有限体积法，把守恒通量都控制在单元中心上，流体流经第i 条边的值由两个相邻单元中心点（k 和p ）的平均值决定：
2/)(p k i W W W += （2-10）
下标的意义如下图：
图1 单元中心，边和网格顶点
设∶ i i i i i x V y U Z ∆-∆= （2-11） 则Euler 方程对k 单元的离散解，方程（2-7）可写为：
∑=⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎣⎡∆-∆+Ω-=kedges i k k H Z x P V Z y P U Z Z dt dW 11ρρρρ （2-12） 上述的一般形式既能用于结构网格，也可以用于非结构网格，唯一的不同点是：前一种情形中方程（2-12）的形式更紧凑，因为前者在采用曲线坐标时，每个单元由一个二阶矩阵来定义。

非结构网格中，相邻的网格点和网格单元间并没有对应关系，所以非结构网格需要一个连接矩阵，这个连接矩阵中储存了所有的必须信息（例如网格顶点的坐标，边数，单元数等等）。

2.3 人工耗散项
如上所述的中心格式是不包含耗散项的，所以任何误差（离散误差，循环误差等等）都不会衰减，最后的定常解可能会出现振荡。

为了减少这些振荡，在方程（2-7）的右边加了人工耗散项，对于k 单元，方程变为：
k k k k D Q dt
dW Ω--=/)( （2-13） 在目前的工作中，采用了Jameson 的方法，耗散项D k 取守恒变量W k 的二阶和四阶差分项的组合。

含四阶差分项的部分加在流场域中的平滑部分，但是在激波区被关掉。

此时，打开含二阶差分项的部分，从而减少激波区的振荡，这个
值可以非常大。

这一开关函数是由基于当地压强的二阶差分的激波感应器控制的。

在非结构网格中，D k 可以表示为：
)4()2(i i k D D D += (2-14)
其中二阶耗散项和四阶耗散项都是单元边界上的耗散通量之和：
∑==kedges K i i d D
1)2()2(，∑==kedges K i i d D 1)4()4( (2-15) )2(i d 表示第i 条边上的二阶耗散通量，)4(i d 表示四阶耗散通量：
i k p i i i W W d )()2()2(-=εα (2-16)
i
k p i i i W W d )(22)4()4(∇-∇-=εα (2-17) 其中指数I 表示单元k 和p 的分界边，2∇被定义为：
∑=-=∇kedges
j k j k W W
W 12)( (2-18)
ε(2) 与ε(4)为二阶和四阶耗散项的自适应系数。

构造每一条边的激波感应器和比例因子，仅仅使用其两个相邻单元k 和p 的流动变量：
k P k
P i P P P P v +-= (2-19)
由此确定的自适应系数变为：
i i v k )2()2(=ε
i i i k ),0max()2()4()4(εε-= (2-20)
比例因子i α取沿单元边界的Jacobian 矩阵∂F/∂W 和∂G/∂W 的最大特征值：
)(22y x c x V y U i ∆+∆+∆-∆=α (2-21)
其中U ，V 和c 是边界上的平均值，c 取当地音速。

在人工粘性项中，二阶耗散项的作用是抑制解在激波附近的振荡,在流场中压力梯度不大的区域内此项作用很小。

四阶耗散项的作用是抑制高频振荡，并使解趋于定态。

2.4 时间离散
定常解是由对常微分方程进行时间积分获得的，此方程可以写成：
k k R dt
dW = （2-22） 方程的右侧的值表示每个单元k 中心的残值：
k k k k D Q R Ω--=/)( （2-23）
对于方程（2-24）的积分可以采用显式的四步格式来完成，由于对于定常解时间的精确并不重要，所以选择此格式只是因为它的稳定和衰减的特性。

目前采用的是以下的格式：
W (0)=W n
W (m)=W (0)+αm △tR (m-1) for m=1 to 4
W n+1=W (4) (2-24)
式中n 是当时的时间步数，n+1是新的时间步数：
Ω--=/)()0()()(D Q R m m （2-25）
其中的系数是：
411=α,312=α,213=α,14=α （2-26）
为了减少计算时间，只在第一步时计算耗散项D ，然后在接下来的各步中D 是一个常值。

对于一个标量模型方程，这种做法改变了格式的稳定区，但是精度和收敛特性被保留了下来。

以上格式在应用于非定常Euler 方程时，能保持稳定的CFL 数最大可以取到22。

这个显式格式最大的缺点是最大的时间步长受到限制，因为稳定域受到限制。

而且，对于多维的方程组，最大时间步长只可以用近似的方法获得。

对于固定形状的网格，采用以下的表达式：
1k
k kedges i i i i i i t CFL U y V x c =Ω∆=⋅∆-∆+∑
（2-27）
2.5 边界条件
在流场中，需要考虑的边界条件包括远场边界条件和物面边界条件。

在这一小节中将分别介绍这两种边界条件。

2.5.1 物面边界条件
对于无粘流体，在物面边界上，应该满足流动方向与物面相切，即物面的法向速度分量为零，可知 0=⋅n v 。

这就是物面无穿透边界条件。

此时屋面边界流量为：
(00)T i flux P y P x =∆-∆ （2-28）
2.5.2 远场边界条件
用有限体积法进行空间离散时，只能取一个有限远的边界作为远场边界，而在绕物体的实际流动中并不存在这个边界，物体所产生的扰动可以认为被传到无穷远而没有反射，即无反射边界条件。

这里，采用一维Riemann 不变量来处理远场边界条件。

令n q ，t q 分别为边界外法向速度和切向速度，根据特征线理论，Riemann 不变量可写为：
t q R =1
s R =2
1
23--=K a q R n 124-+=K a q R n （2-29）
Riemann 不变量的取值可分为以下四种情况：
1) 亚音速入流（－a<n q <0）
R 1,R 2,R 3取来流值，R 4取内场值。

2) 亚音速出流（0<n q <a ）
R 1,R 2,R 4取内场值，R 3取来流值。

3) 超音速入流（n q <-a ）
所有物理量都取来流值。

4) 超音速出流（
q>a）
n
所有物理量都取内场值。

其他所有物理量都由以上四个不变量求出。