可压缩二维无粘流动_二维_欧拉方程_有限差分_MacCormack_Bump

可压缩二维无粘流动

摘要本题利用欧拉方程求解可压缩二维无粘流动，并将其与Numeca Fine/Turbo 的计算结果对比。流道由上平板固壁和带有凸起的下固壁组成，进口给定总温、总压和速度方向，出口给定压力。自编代码求解时，基于有限差分方法，利用MacCormack 格式对控制方程进行离散，根据黎曼不变量和边界条件由内层网格数据外推获得边界数据。文中给出了计算收敛残差历史、密度、速度、压力、马赫数和熵分布，并将其和Numeca 计算结果对比，分析自编代码计算结果的合理性和误差来源。

关键词二维；欧拉方程；有限差分；MacCormack ；Bump

1 问题提出

该问题是经典的Bump 计算问题[1]，如图1所示，上壁为平板，下壁带有凸起，均为滑移边界。进口为轴向进气，且给定总参数为0280T K =和50 1.110p Pa =⨯，出口为5110out p Pa =⨯。

图1准一维管道示意图

本题的分析思路：首先，建立计算域中的主控方程，然后根据MacCormack 格式对方程进行离散，最后通过边界条件和黎曼不变量确定边界数据。收敛条件为相邻时间步的压力差的最大值小于610-Pa 。 2模型建立

物理域中的主控方程为二维欧拉方程，如式(1)所示。将物理域x-y 变换到计算域ξ-η，控制方程变为式(2)，J 为坐标变换的雅克比行列式，y η、x η、y ξ和x ξ均为物理域坐标对计算域的偏导数。

22

0,,,u v u u p uv where v uv v p t x y

E Hu Hv ρρρρρρρρρρρρ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥

+∂∂∂⎢⎥⎢⎥⎢⎥++= ===⎢⎥⎢⎥⎢⎥+∂∂∂⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦

Q F G Q F G (1)

0,,,where J y x y x t ηηξξξη

∂∂∂++= ==- =-+∂∂∂Q'F'G'

Q'Q F'F G G'F G (2)

未知物理量为,,,,,p u v H E ρ共6个，因此为了方程组的封闭西需要补充两个方程。补充方程如下：

()22111, 1.422p E u v γργ⎛

⎫=---= ⎪⎝

⎭

(3)

p

H E ρ

=+

(4)

至此，已实现方程组的封闭，可以进行求解。 3网格划分

计算网格已经给定，分析网格数据可知，第一个数字表示网格块的数目，本题中是1；第二到四个数字表示xyz 三个方向的节点数，以L1为例，分别是65、17和1，说明网格为长方向有65个节点、宽方向有17个节点的单层网格；剩余的数字均为网格坐标，前65×17个数是所有点的x 坐标，中间65×17个数是所有点的y 坐标，最后的是z 坐标（本题中全是0）。计算时将物理域转换成矩形的计算域ξ-η，取Δξ=Δη=1。 5差分离散和边界条件

MacCormack 格式广泛应用于求解流动方程，本题中其预估和校正格式分别如式(5)和式(6)，式中i 和j 分别对应η和ξ节点顺序。

1.,,1,1,,,1,,11,,1,22()()11

[(2)(2)]

n n n n n n i j i j i j i j i j i j n n n n n n i j i j i j i j i j i j t t

t ξη

αξη++++-+-∆∆=-

---+∆∆∆-++-+∆∆Q'Q'F'F'G'G'Q'Q'Q'Q'Q'Q'

(5)

111111

,,,,,1,1,1111111,,1,,1

,,1

2

2

1[()()

222(

)]

n n n n n n n i j i j i j i j i j i j i j n n n n n n i j

i j

i j

i j i j

i j t t t ξηαηξ++++++--+++++++-+-∆∆=+----∆∆-+-++∆+

∆∆Q'Q'Q'F'F'G'G'Q'

Q'

Q'

Q'

Q'

Q'

(6)

以上离散格式只能用于计算非边界网格点，边界网格点数据需由相邻非边界点外推得到。

判断本问题为二维亚音入流，因此在进口处需要保持相邻两层的黎曼不变量相同，再结合总温、总压和速度方向，计算出,,,,,p u v H E ；出口给法类似于进口，不过出口

只需一个已知量（如压力）即可；在固壁处，让其压力和密度和相邻内层的相等，将内层速度的法向分量去掉，赋于边界，作为边界处速度。

6计算结果

以L1为例展示计算结果，计算时人工粘性为0，如图2所示。经过约5万步计算达到收敛，在第4万步时计算曾被暂停，优化代码后继续进行。

图2L1计算结果

由于是亚音无粘流动，因此流动过程等熵。考虑其物理过程，因为几何完全对称，进口给定总问总压和速度方向，出口给定静压，故以中心（x = 0.5）为对称轴的标量对

x 10

4

收敛历史

步数

l o g (残差)

密度分布(kg/m 3)

1.23

1.241.251.261.271.281.291.31.31速度分布

(m/s)

90

100

110120

130

140

150

压力分布(Pa)

9.4

9.59.69.79.89.91010.110.210.3

10.44

马赫数分布

0.28

0.30.320.34

0.36

0.38

0.40.420.440.46熵分布(J/K)

8012

8012.5

8013

8013.58014

8014.5