计算流体力学SOD激波管

1一、题目如图所示初始时刻在x=0左右两侧的气体密度和压力存在间断，计算t=0.1时刻管道内的工质物理参数分布。二、格式介绍本次作业采用迎风格式计算，使用Steger-Warming和Roe格式分别计算了流场分布。其中Steger-Warming采用了一阶和二阶迎风格式，Roe格式采用一阶迎风格式。此外对于Steger-Warming格式使用了Van-leer限制器和Min-mod限制器来提高分辨率。对于Roe格式尝试使用MUSCL重构，但是没能调试完成。1.Steger-Warming格式Euler方程𝜕𝑈𝜕𝑡+𝜕𝑓(𝑈)𝜕𝑥=0𝑈=[𝜌𝜌𝑢𝐸]𝑓(𝑈)=[𝜌𝑢𝜌𝑢2+𝑝𝑢(𝐸+𝑝)]𝜕𝑈𝜕𝑡+𝜕𝐴𝑈𝜕𝑥=0对流通量为AU，A是U的函数，无法确定其正负因此无法直接对对流通量使用迎风格式进行差分。将A相似对角化后，A=𝑆−1Λ𝑆。矢通量分裂格式是通过将对角阵Λ的特征值分裂成正特征值和负特征值之和，进而对每个点得到正通量和负通量，再分别对正负通量进行迎风差分。令特征值𝜆𝑘±=𝜆𝑘±√𝜆𝑘2+𝜀22其中𝜀是一个常数，存在的目的是为了使𝜆𝑘在接近0的时候变化连续。将正负特征值代入对角阵得到𝛬±，进而得到𝐴±以及𝑓±=𝐴±𝑈，在每个网格节点上对该点计算一个正通量计算一个负通量，总通量等于正负通量之和：𝑓𝑗=𝑓𝑗++𝑓𝑗−差分时分别对𝜕𝑓+𝜕𝑥和𝜕𝑓−𝜕𝑥使用迎风格式进行差分。本次作业使用一阶时间向前推进算法进行时间推进求解，流通矢量的差分分别使用了一阶和二阶迎风差分，得到的计算公式为𝑈𝑗𝑛+1=𝑈𝑗𝑛−∆𝑡∆𝑥(𝑓𝑗+12𝑛+−𝑓𝑗−12𝑛++𝑓𝑗+12𝑛−−𝑓𝑗−12𝑛−)其中，对于一阶迎风格式𝑓𝑗+12+=𝑓𝑗+𝑓𝑗+12−=𝑓𝑗+1−2对于二阶迎风格式𝑓𝑗+12+=𝑓𝑗++𝑓𝑗+−𝑓𝑗−1+2𝑓𝑗+12−=𝑓𝑖+1−+𝑓𝑗+1−−𝑓𝑗+2−2由于二阶格式精度高但是会存在一定的数值振荡，因此可以在格式中加入限制器来提高格式的分辨率。2.限制器本次作业选用了Van-Leer和Min-mod两种限制器对Steger-Warming格式进行了优化。限制器本质上是用二阶格式对一阶格式进行优化，即𝑓𝑗+12=𝑓𝑗+𝜑修正函数使用二阶迎风和二阶中心同时修正，当两个二阶函数的修正都为正或者都为负时，选用其中一个对一阶迎风格式进行修正，而当二者修正方向相反时，则不修正。这样得到的结果就是在流场变化不剧烈的时候，二阶格式发挥作用降低耗散，在流场变化剧烈的时候，一阶格式发挥作用防止数值振荡。具体的格式为𝑓𝑗+12+=𝑓𝑗++𝜑(𝑟1/2+)𝑓𝑗+1+−𝑓𝑗+2𝑓𝑗+12−=𝑓𝑗+1−+𝜑(𝑟1/2−)𝑓𝑗−−𝑓𝑗+1−2𝑟1/2+=𝑓𝑗+−𝑓𝑗−1+𝑓𝑗+1+−𝑓𝑗++𝜀𝑟1/2−=𝑓𝑗+1−−𝑓𝑗+2−𝑓𝑗−−𝑓𝑗+1−+𝜀其中𝜀是一个低于要求计算精度的小量，防止出现分母为0的情况，本次计算取10−6。对于Van-Leer限制器𝜑(𝑟)=|𝑟|+𝑟1+𝑟对于Min-mid限制器𝜑(𝑟)=𝑀𝑖𝑛𝑚𝑜𝑑(𝑟,1)𝑀𝑖𝑛𝑚𝑜𝑑函数具体的形式是当迎风格式与中心格式都为正或都为负时，取其中绝对值较小的一个作为修正，当二者正负不同时，𝜑(𝑟)=0，即不修正。3.Roe格式Roe格式是一种通量差分分裂格式，具体的做法是将Euler方程变为𝜕𝑈𝜕𝑡+𝐴𝜕𝑈𝜕𝑥=0𝐴=𝜕𝑓𝜕𝑈如果能得到𝜕𝑓𝜕𝑈的平均值，以其作为A的值，就能把原方程化为一个常系数的方程，进而实现通量的分裂。Roe的方法是将U进行整理，使其变化为W：𝑊=√𝜌[1𝑢𝐻]A关于W是一个二次齐次函数，根据二次齐次函数的性质，在[𝑊𝑅,𝑊𝐿]上中点处导数的值就是导数的平均值，因此可以用这种方法得到A的平均值𝐴̃。然后在[𝑈𝑅,𝑈𝐿]上使用Roe平均参数计算𝐴̃，在[𝑈𝑅,𝑈𝐿]上即可得知𝑓1+12=[𝑓(𝑈𝑅)+𝑓(𝑈𝐿)]−12|𝐴̃(𝑈𝑅,𝑈𝐿)|(𝑈𝑅,𝑈𝐿)3|𝐴̃(𝑈𝑅,𝑈𝐿)|=𝑆−1|Λ|𝑆|Λ|即为所有特征值取绝对值的对角阵。可以看到上式在特征值为正时为正通量，特征值为负时，为负通量，实现了矢量的分裂。三、结果分析1.Steger-Warming格式4一阶和二阶Steger-Warming格式的结果如上，可以看到一阶格式的耗散相对比较严重，但是二阶格式虽然耗散减小，数值振荡却非常明显。5从上图可以看到，加了限制器以后二阶格式的数值振荡明显消失了，此外Van-Leer限制器对于数值振荡的减小相比于Min-mod限制器更加优秀。2.Roe格式6一阶Roe格式也没有数值振荡，但是耗散也还是比较大，使用MUSCL重构来获得𝑈𝑅和𝑈𝐿的值可以提高分辨率，但是本次作业的调试尚未完成。四、源代码#includestdio.h#includestdlib.h#includemath.hconstdoublegama=1.4,x0=-1.0,x1=1.0;constintdot_x=501;doublefi(double);doublemmd(double,double);intmain(){doubledelt_x,delt_t;doublepre[dot_x],vlo[dot_x],rho[dot_x],son[dot_x],u1[dot_x],u2[dot_x],u3[dot_x];inti=0,j=0,k=0;delt_x=(x1-x0)/(dot_x-1);delt_t=0.0001;/*初始化物理量*/for(i=0;idot_x;i++){vlo[i]=0.0;if(i*delt_x-10||i*delt_x-1==0){pre[i]=1.0;rho[i]=1.0;}else{pre[i]=0.1;rho[i]=0.125;}son[i]=sqrt(gama*pre[i]/rho[i]);u1[i]=rho[i];u2[i]=rho[i]*vlo[i];7u3[i]=pre[i]/(gama-1)+0.5*rho[i]*vlo[i]*vlo[i];}/*Steger-Warming求解*/doubleeps_sw=0;doubleF_sw1_p[dot_x],F_sw2_p[dot_x],F_sw3_p[dot_x],F_sw1_n[dot_x],F_sw2_n[dot_x],F_sw3_n[dot_x],lam[3],lam_p[3],lam_n[3];doubler_p1R,r_p1L,r_n1R,r_n1L,r_p2R,r_p2L,r_n2R,r_n2L,r_p3R,r_p3L,r_n3R,r_n3L,eps=1e-6;for(i=1;i*delt_t=0.1;i++){for(j=0;jdot_x;j++){lam[0]=vlo[j];lam[1]=vlo[j]-son[j];lam[2]=vlo[j]+son[j];for(k=0;k3;k++){lam_p[k]=(lam[k]+sqrt(lam[k]*lam[k]+eps_sw*eps_sw))/2;lam_n[k]=(lam[k]-sqrt(lam[k]*lam[k]+eps_sw*eps_sw))/2;}F_sw1_p[j]=(rho[j]/(2*gama))*(2*(gama-1)*lam_p[0]+lam_p[1]+lam_p[2]);F_sw1_n[j]=(rho[j]/(2*gama))*(2*(gama-1)*lam_n[0]+lam_n[1]+lam_n[2]);F_sw2_p[j]=(rho[j]/(2*gama))*(2*(gama-1)*lam_p[0]*vlo[j]+lam_p[1]*(vlo[j]-son[j])+lam_p[2]*(vlo[j]+son[j]));F_sw2_n[j]=(rho[j]/(2*gama))*(2*(gama-1)*lam_n[0]*vlo[j]+lam_n[1]*(vlo[j]-son[j])+lam_n[2]*(vlo[j]+son[j]));F_sw3_p[j]=(rho[j]/(2*gama))*((gama-1)*lam_p[0]*vlo[j]*vlo[j]+0.5*lam_p[1]*(vlo[j]-son[j])*(vlo[j]-son[j])+0.5*lam_p[2]*(vlo[j]+son[j])*(vlo[j]+son[j])+((3-gama)*(lam_p[1]+lam_p[2])*son[j]*son[j]/(2*(gama-1))));F_sw3_n[j]=(rho[j]/(2*gama))*((gama-1)*lam_n[0]*vlo[j]*vlo[j]+0.5*lam_n[1]*(vlo[j]-son[j])*(vlo[j]-son[j])+0.5*lam_n[2]*(vlo[j]+son[j])*(vlo[j]+son[j])+((3-gama)*(lam_n[1]+lam_n[2])*son[j]*son[j]/(2*(gama-1))));}//一阶迎风格式/*for(j=1;jdot_x-1;j++){u1[j]=u1[j]-(delt_t/delt_x)*(F_sw1_p[j]-F_sw1_p[j-1]+F_sw1_n[j+1]-F_sw1_n[j]);u2[j]=u2[j]-(delt_t/delt_x)*(F_sw2_p[j]-F_sw2_p[j-1]+F_sw2_n[j+1]-F_sw2_n[j]);u3[j]=u3[j]-(delt_t/delt_x)*(F_sw3_p[j]-F_sw3_p[j-1]+F_sw3_n[j+1]-F_sw3_n[j]);}//二阶迎风格式for(j=2;jdot_x-2;j++){8u1[j]=u1[j]-(0.5*delt_t/delt_x)*(3*F_sw1_p[j]-4*F_sw1_p[j-1]+F_sw1_p[j-2]-F_sw1_n[j+2]+4*F_sw1_n[j+1]-3*F_sw1_n[j]);u2[j]=u2[j]-(0.5*delt_t/delt_x)*(3*F_sw2_p[j]-4*F_sw2_p[j-1]+F_sw2_p[j-2]-F_sw2_n[j+2]+4*F_sw2_n[j+1]-3*F_sw2_n[j]);u3[j]=u3[j]-(0.5*delt_t/delt_x)*(3*F_sw3_p[j]-4*F_sw3_p[j-1]+F_sw3_p[j-2]-F_sw3_n[j+2]+4*F_sw3_n[j+1]-3*F_sw3_n[j]);}*///限制器for(j=2;jdot_x-2;j++){r_p1R=(F_sw1_p[j]-F_sw1_p[j-1])/(F_sw1_p[j+1]-F_sw1_p[j]+eps);r_p1L=(F_sw1_p[j-1]-F_sw1_p[j-2])/(F_sw1_p[j]-F_sw1_p[j-1]+eps);r_n1R=(F_sw1_n[j+1]-F_sw1_n[j+2])/(F_sw1_n[j]-F_sw1_n[j+1]+eps);r_n1L=(F_sw1_n[j]-F_sw1_n[j+1])/(F_sw1_n[j-1]-F_sw1_n[j]+eps);u1[j]=u1[j]-(delt_t/delt_x)*(F_sw1_p[j]+fi(r_p1R)*0.5*(F_sw1_p[j+1]-F_sw1_p[j])-F_sw1_