用有限体积方法求解欧拉方程

有限体积法求解二维可压缩Euler方程

——计算流体力学课程大作业

老师：夏健、刘学强

学生：徐锡虎

学号：SQ***********

日期：2010年2月5日

目录

一、内容摘要 (2)

二、流动控制方程 (2)

三、有限体积法的空间离散 (2)

四、人工耗散 (3)

五、时间离散 (4)

六、边界条件 (5)

七、计算结果 (8)

八、结论与展望 (11)

参考文献 (11)

一、内容摘要

本文通过运用JAMESON 有限体积法求解了二维定常和非定常可压缩Euler 方程。程序实现语言为C++。其中，使用的网格是三角形非结构网格。在时间推进上使用的是四步龙—库塔推进格式。推进的时间步长取的是当地的时间步长。为了消除迭代误差、round-off 等误差，本文采用了添加人工耗散项的办法。另外，本文计算了NACA0012翼型在跨音速下不同迎角的情况，并与fluent 软件的计算结果进行了比较，来验证程序的准确性。

二、流动控制方程

守恒形式的Euler 方程：

0=-+Ω∂∂

⎰⎰Ω

Gdx Fdy wd t S

（1） 其中x 和y 代表笛卡儿坐标系。W 是守恒变量。

⎥⎥⎥⎥

⎦

⎤⎢⎢⎢⎢⎣⎡=E V U W ρρρρ (2)

F,G 表示通量

⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡+=UH UV P U U F ρρρρ2, ⎥⎥⎥⎥

⎦

⎤⎢⎢⎢⎢⎣⎡+=VH P V UV V G ρρρρ2

(3) ρ,P , H 和E 表示密度，压强，单元总焓和单元总能量。U,V 表示笛卡儿坐标系下

的速度矢量。这些量由理想气体的单位体积的总能量和总焓相互联系。

2/122）（）（V U P E ++-=ργρ (4)

P E H +=ρρ (5)

三、有限体积法的空间离散

计算域被划分为互不重叠的单元。在每个单元运用守恒形式的Euler 方程。由于每个单

元相对于时间都是不变的，所以等式（1）可以写成：

⎰⎰Ω

Ω

--=∂∂d Gdx Fdy t

W

S )

( （8）

其中Ω和S 是单元的体积和边界。W 是单元的平均值。

在对上述方程进行时间离散前，先对空间进行离散，则方程（6）可以写为：

k k k Q

dt

dW Ω-= （9）

其中k Ω表示第k 个单元的体积，k W 是第k 个单元的守恒变量。k Q 表示第k 个单元的通量。方程（7）的右边项可以写成：

∑=∆-∆=

kedges

i i

k x G y F Q 1

)

( （10）

其中 a b i a b i y y y x x x -=∆-=∆, （11） （8）式中的求和是对第k 个单元的所有边进行的。守恒参数的量是单元中心值，在求通量时，第I 条边的守恒参数值是用左右单元的平均来表示的：

2/)(W p k i W W += （12）

引入变量：

i i i i i x V y U Z ∆-∆= （13）

则第k 单元的Euler 方程可以写为：

∑=⎥

⎥

⎥

⎥

⎦⎤⎢⎢⎢⎢⎣⎡∆-∆+Ω-=kedges i i

k k H Z x P V Z y P U Z Z dt dW 11ρρρρ （14） 在本文中，采用的是JAMENSON 有限体积法，为了减少存储的相关信息的量，其存储的

方式选择的是按边存储的方法。在存储的每条边的信息中，包含了这条边的边号，左右单元号和边的端点。在计算通量时采用按边循环的方式：

do I=1,nedge

k=connmatrix(I,1) a=connmatrix(I,2) b=connmatrix(I,3) p=connmatrix(I,4)

flux=function(k,a,b,p) sum(k)=sum(k)+flux sum(p)=sum(p)-flux end do

这里给出的是FOTRAN 语言的形式，我编写采用的是C ，具体表现在上交的程序中。 在计算时间步长、人工耗散项等也可用象这样按边循环，从此处我们可以看出求解时与单元的形状无关。

四、人工耗散

人工粘性模型对方法的成功应用起着关键作用，人工粘性抑制解在激波附近的振荡，又阻尼解在光滑区域内的高阶误差，对解的线性稳定和收敛于定态是很重要的。本文在方程（14）的右端加入了人工耗散项，如对于单元k ，其表达式可以表示为：

()

k k k k D Q dt

dW Ω--= （15）

在有限体积法中，耗散项的公式可以表示为：

∑∑==+

=

kedges

i i

kedges

i i

k d

d

D 1

)4(1

)2( （16）

其中：

()()i

k p

i

i i

i k p i i i W W

d

W W d 2

2

)4()4()2()2(∇-∇=-=εαεα （17）

其中的I 表示单元k 和p 的公共边，2

∇定义为：

()∑=-=

∇kedges i k j

k W W

W 1

2

（18）

上面的j 表示与k 相邻的单元。

∑∑==+-=

kedges i i

k p

kedges

i i

k p

k P P

P P

1

1)()(ν （19）

)

,0max(),max()2()

4()4()2()2(i

i

i k p i k

k ε

ε

ννε-== （20）

其中的量)4()

2(,k k

的范围是：0.121,12561)2()4(<<<

在计算时发现上面方法得到的人工耗散项并不太适合。其在光滑区域耗散项太大，而在大剃度区域又显得太小，为了弥补上面的不足，作下面的修改：

k

p k p i P P P P +-=

ν （21）

自适应系数为：

)

,0max()2()

4()4()2()2(i

i

i i k

k ε

ε

νε-== （22）

尺度系数为：

()22y x c x V y U i ∆+∆+∆-∆=α （23）

其中的U,V 表示边上的值，C 表示当地声速。

五、时间离散

方程最后的稳定解是通过时间上的迭代得到的，可以写为：

k k

R dt

dW = （24） 右边项的表达式为：