计算流体力学SOD激波管.docx

《计算流体力学SOD激波管.docx》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。

1、一、题目如图所示初始时刻在x=0左右两侧的气体密度和压力存在间断，计算t=0.1时刻管道内的工质物理参数分布。二、格式介绍本次作业采用迎风格式计算，使用Steger-Warming和Roe格式分别计算了流场分布。其中Steger-Warming采用了一阶和二阶迎风格式，Roe格式采用一阶迎风格式。此外对于Steger-Warming格式使用了Van-leer限制器和Min-mod限制器来提高分辨率。对于Roe格式尝试使用MUSCL重构，但是没能调试完成。1.Steger-Warming格式Euler方程∂U∂t+∂f(U)∂x=0U

2、=ρρuEfU=ρuρu2+pu(E+p)∂U∂t+∂AU∂x=0对流通量为AU，A是U的函数，无法确定其正负因此无法直接对对流通量使用迎风格式进行差分。将A相似对角化后，A=S-1ΛS。矢通量分裂格式是通过将对角阵Λ的特征值分裂成正特征值和负特征值之和，进而对每个点得到正通量和负通量，再分别对正负通量进行迎风差分。令特征值λk±=λk±λk2+ε22其中ε是一个常数，存在的目的是为了使λk在接近0的时候变化连续。将正负特征值代入对角阵得到Λ±，进而得到A±以及f±=A±U，在每个网格节点上对该点计算一个正通量计算一个负通量，总通量

3、等于正负通量之和：fj=fj++fj-差分时分别对∂f+∂x和∂f-∂x使用迎风格式进行差分。本次作业使用一阶时间向前推进算法进行时间推进求解，流通矢量的差分分别使用了一阶和二阶迎风差分，得到的计算公式为Ujn+1=Ujn-∆t∆x(fj+12n+-fj-12n++fj+12n--fj-12n-)其中，对于一阶迎风格式fj+12+=fj+fj+12-=fj+1-对于二阶迎风格式fj+12+=fj++fj+-fj-1+2fj+12-=fi+1-+fj+1--fj+2-2由于二阶格式精度高但是会存在一定的数值振荡，因此可以在格式中加入限

4、制器来提高格式的分辨率。1.限制器本次作业选用了Van-Leer和Min-mod两种限制器对Steger-Warming格式进行了优化。限制器本质上是用二阶格式对一阶格式进行优化，即fj+12=fj+φ修正函数使用二阶迎风和二阶中心同时修正，当两个二阶函数的修正都为正或者都为负时，选用其中一个对一阶迎风格式进行修正，而当二者修正方向相反时，则不修正。这样得到的结果就是在流场变化不剧烈的时候，二阶格式发挥作用降低耗散，在流场变化剧烈的时候，一阶格式发挥作用防止数值振荡。具体的格式为fj+12+=fj++φr1/2+fj+1+-fj+2

5、fj+12-=fj+1-+φr1/2-fj--fj+1-2r1/2+=fj+-fj-1+fj+1+-fj++εr1/2-=fj+1--fj+2-fj--fj+1-+ε其中ε是一个低于要求计算精度的小量，防止出现分母为0的情况，本次计算取10-6。对于Van-Leer限制器φr=r+r1+r对于Min-mid限制器φr=Minmod(r,1)Minmod函数具体的形式是当迎风格式与中心格式都为正或都为负时，取其中绝对值较小的一个作为修正，当二者正负不同时，φr=0，即不修正。2.Roe格式Roe格式是一种通量差分分裂格式，具体的做法是

6、将Euler方程变为∂U∂t+A∂U∂x=0A=∂f∂U如果能得到∂f∂U的平均值，以其作为A的值，就能把原方程化为一个常系数的方程，进而实现通量的分裂。Roe的方法是将U进行整理，使其变化为W：W=ρ1uHA关于W是一个二次齐次函数，根据二次齐次函数的性质，在[WR,WL]上中点处导数的值就是导数的平均值，因此可以用这种方法得到A的平均值A。然后在[UR,UL]上使用Roe平均参数计算A，在[UR,UL]上即可得知f1+12=fUR+fUL-12

7、A(UR,UL)

8、(UR,UL)AUR,UL=S-1

9、Λ

10、S

11、Λ

12、即为所有特征值取绝

13、对值的对角阵。可以看到上式在特征值为正时为正通量，特征值为负时，为负通量，实现了矢量的分裂。一、结果分析1.Steger-Warming格式一阶和二阶Steger-Warming格式的结果如上，可以看到一阶格式的耗散相对比较严重，但是二阶格式虽然耗散减小，数值振荡却非常明显。从上图可以看到，加了限制器以后二阶格式的数值振荡明显消失了，此外Van-Leer限制器对于数值振荡的减小相比于Min-mod限制器更加优秀。1.Roe格式一阶Roe格式也没有数值振荡，但是耗散也还是比较大，使用MUSCL重构来获得UR和UL的值可以提高分辨率，但

14、是本次作业的调试尚未完成。一、源代码#include#include#includeconstdoublegama=1.4,x0=-1.0,x1=1.0;constintdot_x