计算流体力学SOD激波管(共14页).docx
精选优质文档-倾情为你奉上一、 题目如图所示初始时刻在x=0左右两侧的气体密度和压力存在间断，计算t=0.1时刻管道内的工质物理参数分布。二、 格式介绍本次作业采用迎风格式计算，使用Steger-Warming和Roe格式分别计算了流场分布。其中Steger-Warming采用了一阶和二阶迎风格式，Roe格式采用一阶迎风格式。此外对于Steger-Warming格式使用了Van-leer限制器和Min-mod限制器来提高分辨率。对于Roe格式尝试使用MUSCL重构，但是没能调试完成。 1. Steger-Warming格式Euler方程Ut+f(U)x=0U=uE fU=uu2+pu(E+p)Ut+AUx=0对流通量为AU，A是U的函数，无法确定其正负因此无法直接对对流通量使用迎风格式进行差分。将A相似对角化后，A=S-1S。矢通量分裂格式是通过将对角阵的特征值分裂成正特征值和负特征值之和，进而对每个点得到正通量和负通量，再分别对正负通量进行迎风差分。令特征值k±=k±k2+22其中是一个常数，存在的目的是为了使k在接近0的时候变化连续。将正负特征值代入对角阵得到±，进而得到A±以及f±=A±U，在每个网格节点上对该点计算一个正通量计算一个负通量，总通量等于正负通量之和：fj=fj+fj-差分时分别对f+x和f-x使用迎风格式进行差分。本次作业使用一阶时间向前推进算法进行时间推进求解，流通矢量的差分分别使用了一阶和二阶迎风差分，得到的计算公式为Ujn+1=Ujn-tx(fj+12n+-fj-12n+fj+12n-fj-12n-)其中，对于一阶迎风格式fj+12+=fj+ fj+12-=fj+1-对于二阶迎风格式fj+12+=fj+fj+-fj-1+2 fj+12-=fi+1-+fj+1-fj+2-2由于二阶格式精度高但是会存在一定的数值振荡，因此可以在格式中加入限制器来提高格式的分辨率。2. 限制器本次作业选用了Van-Leer和Min-mod两种限制器对Steger-Warming格式进行了优化。限制器本质上是用二阶格式对一阶格式进行优化，即fj+12=fj+修正函数使用二阶迎风和二阶中心同时修正，当两个二阶函数的修正都为正或者都为负时，选用其中一个对一阶迎风格式进行修正，而当二者修正方向相反时，则不修正。这样得到的结果就是在流场变化不剧烈的时候，二阶格式发挥作用降低耗散，在流场变化剧烈的时候，一阶格式发挥作用防止数值振荡。具体的格式为fj+12+=fj+r1/2+fj+1+-fj+2 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，即不修正。3. Roe格式Roe格式是一种通量差分分裂格式，具体的做法是将Euler方程变为Ut+AUx=0A=fU如果能得到fU的平均值，以其作为A的值，就能把原方程化为一个常系数的方程，进而实现通量的分裂。Roe的方法是将U进行整理，使其变化为W：W=1uHA关于W是一个二次齐次函数，根据二次齐次函数的性质，在WR,WL上中点处导数的值就是导数的平均值，因此可以用这种方法得到A的平均值A。然后在UR,UL上使用Roe平均参数计算A，在UR,UL上即可得知f1+12=fUR+fUL-12|A(UR,UL)|(UR,UL)AUR,UL=S-1|S|即为所有特征值取绝对值的对角阵。可以看到上式在特征值为正时为正通量，特征值为负时，为负通量，实现了矢量的分裂。三、 结果分析1. Steger-Warming格式 一阶和二阶Steger-Warming格式的结果如上，可以看到一阶格式的耗散相对比较严重，但是二阶格式虽然耗散减小，数值振荡却非常明显。从上图可以看到，加了限制器以后二阶格式的数值振荡明显消失了，此外Van-Leer限制器对于数值振荡的减小相比于Min-mod限制器更加优秀。2. Roe格式 一阶Roe格式也没有数值振荡，但是耗散也还是比较大，使用MUSCL重构来获得UR和UL的值可以提高分辨率，但是本次作业的调试尚未完成。四、 源代码#include <stdio.h>#include <stdlib.h>#include <math.h>const double gama=1.4,x0=-1.0,x1=1.0;const int dot_x=501;double fi(double);double mmd(double,double);int main()double delt_x,delt_t;double predot_x,vlodot_x,rhodot_x,sondot_x,u1dot_x,u2dot_x,u3dot_x;int i=0,j=0,k=0;delt_x=(x1-x0)/(dot_x-1);delt_t=0.0001;/*初始化物理量*/for(i=0;i<dot_x;i+)vloi=0.0;if(i*delt_x-1<0|i*delt_x-1=0)prei=1.0;rhoi=1.0;elseprei=0.1;rhoi=0.125;soni=sqrt(gama*prei/rhoi);u1i=rhoi;u2i=rhoi*vloi;u3i=prei/(gama-1)+0.5*rhoi*vloi*vloi;/*Steger-Warming求解*/ double eps_sw=0;double F_sw1_pdot_x,F_sw2_pdot_x,F_sw3_pdot_x,F_sw1_ndot_x,F_sw2_ndot_x,F_sw3_ndot_x,lam3,lam_p3,lam_n3;double r_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;j<dot_x;j+)lam0=vloj;lam1=vloj-sonj;lam2=vloj+sonj;for(k=0;k<3;k+)lam_pk=(lamk+sqrt(lamk*lamk+eps_sw*eps_sw)/2;lam_nk=(lamk-sqrt(lamk*lamk+eps_sw*eps_sw)/2;F_sw1_pj=(rhoj/(2*gama)*(2*(gama-1)*lam_p0+lam_p1+lam_p2);F_sw1_nj=(rhoj/(2*gama)*(2*(gama-1)*lam_n0+lam_n1+lam_n2);F_sw2_pj=(rhoj/(2*gama)*(2*(gama-1)*lam_p0*vloj+lam_p1*(vloj-sonj)+lam_p2*(vloj+sonj);F_sw2_nj=(rhoj/(2*gama)*(2*(gama-1)*lam_n0*vloj+lam_n1*(vloj-sonj)+lam_n2*(vloj+sonj);F_sw3_pj=(rhoj/(2*gama)*(gama-1)*lam_p0*vloj*vloj+0.5*lam_p1*(vloj-sonj)*(vloj-sonj)+0.5*lam_p2*(vloj+sonj)*(vloj+sonj)+(3-gama)*(lam_p1+lam_p2)*sonj*sonj/(2*(gama-1);F_sw3_nj=(rhoj/(2*gama)*(gama-1)*lam_n0*vloj*vloj+0.5*lam_n1*(vloj-sonj)*(vloj-sonj)+0.5*lam_n2*(vloj+sonj)*(vloj+sonj)+(3-gama)*(lam_n1+lam_n2)*sonj*sonj/(2*(gama-1);/一阶迎风格式/*for(j=1;j<dot_x-1;j+)u1j=u1j-(delt_t/delt_x)*(F_sw1_pj-F_sw1_pj-1+F_sw1_nj+1-F_sw1_nj);u2j=u2j-(delt_t/delt_x)*(F_sw2_pj-F_sw2_pj-1+F_sw2_nj+1-F_sw2_nj);u3j=u3j-(delt_t/delt_x)*(F_sw3_pj-F_sw3_pj-1+F_sw3_nj+1-F_sw3_nj);/二阶迎风格式for(j=2;j<dot_x-2;j+)u1j=u1j-(0.5*delt_t/delt_x)*(3*F_sw1_pj-4*F_sw1_pj-1+F_sw1_pj-2-F_sw1_nj+2+4*F_sw1_nj+1-3*F_sw1_nj);u2j=u2j-(0.5*delt_t/delt_x)*(3*F_sw2_pj-4*F_sw2_pj-1+F_sw2_pj-2-F_sw2_nj+2+4*F_sw2_nj+1-3*F_sw2_nj);u3j=u3j-(0.5*delt_t/delt_x)*(3*F_sw3_pj-4*F_sw3_pj-1+F_sw3_pj-2-F_sw3_nj+2+4*F_sw3_nj+1-3*F_sw3_nj);*/限制器 for(j=2;j<dot_x-2;j+)r_p1R=(F_sw1_pj-F_sw1_pj-1)/(F_sw1_pj+1-F_sw1_pj+eps);r_p1L=(F_sw1_pj-1-F_sw1_pj-2)/(F_sw1_pj-F_sw1_pj-1+eps);r_n1R=(F_sw1_nj+1-F_sw1_nj+2)/(F_sw1_nj-F_sw1_nj+1+eps);r_n1L=(F_sw1_nj-F_sw1_nj+1)/(F_sw1_nj-1-F_sw1_nj+eps);u1j=u1j-(delt_t/delt_x)*(F_sw1_pj+fi(r_p1R)*0.5*(F_sw1_pj+1-F_sw1_pj)-F_sw1_pj-1-fi(r_p1L)*0.5*(F_sw1_pj-F_sw1_pj-1)+F_sw1_nj+1+fi(r_n1R)*0.5*(F_sw1_nj-F_sw1_nj+1)-F_sw1_nj-fi(r_n1L)*0.5*(F_sw1_nj-1-F_sw1_nj);r_p2R=(F_sw2_pj-F_sw2_pj-1)/(F_sw2_pj+1-F_sw2_pj+eps);r_p2L=(F_sw2_pj-1-F_sw2_pj-2)/(F_sw2_pj-F_sw2_pj-1+eps);r_n2R=(F_sw2_nj+1-F_sw2_nj+2)/(F_sw2_nj-F_sw2_nj+1+eps);r_n2L=(F_sw2_nj-F_sw2_nj+1)/(F_sw2_nj-1-F_sw2_nj+eps);u2j=u2j-(delt_t/delt_x)*(F_sw2_pj+fi(r_p2R)*0.5*(F_sw2_pj+1-F_sw2_pj)-F_sw2_pj-1-fi(r_p2L)*0.5*(F_sw2_pj-F_sw2_pj-1)+F_sw2_nj+1+fi(r_n2R)*0.5*(F_sw2_nj-F_sw2_nj+1)-F_sw2_nj-fi(r_n2L)*0.5*(F_sw2_nj-1-F_sw2_nj);r_p3R=(F_sw3_pj-F_sw3_pj-1)/(F_sw3_pj+1-F_sw3_pj+eps);r_p3L=(F_sw3_pj-1-F_sw3_pj-2)/(F_sw3_pj-F_sw3_pj-1+eps);r_n3R=(F_sw3_nj+1-F_sw3_nj+2)/(F_sw3_nj-F_sw3_nj+1+eps);r_n3L=(F_sw3_nj-F_sw3_nj+1)/(F_sw3_nj-1-F_sw3_nj+eps);u3j=u3j-(delt_t/delt_x)*(F_sw3_pj+fi(r_p3R)*0.5*(F_sw3_pj+1-F_sw3_pj)-F_sw3_pj-1-fi(r_p3L)*0.5*(F_sw3_pj-F_sw3_pj-1)+F_sw3_nj+1+fi(r_n3R)*0.5*(F_sw3_nj-F_sw3_nj+1)-F_sw3_nj-fi(r_n3L)*0.5*(F_sw3_nj-1-F_sw3_nj); /mmd/*for(j=1;j<dot_x-1;j+)u1j=u1j-(delt_t/delt_x)*(F_sw1_pj+0.5*mmd(F_sw1_pj-F_sw1_pj-1,F_sw1_pj+1-F_sw1_pj)-F_sw1_pj-1-0.5*mmd(F_sw1_pj-1-F_sw1_pj-2,F_sw1_pj-F_sw1_pj-1)+F_sw1_nj+1-0.5*mmd(F_sw1_nj-F_sw1_nj+1,F_sw1_nj+1-F_sw1_nj+2)-F_sw1_nj+0.5*mmd(F_sw1_nj-1-F_sw1_nj,F_sw1_nj-F_sw1_nj+1);u2j=u2j-(delt_t/delt_x)*(F_sw2_pj+0.5*mmd(F_sw2_pj-F_sw2_pj-1,F_sw2_pj+1-F_sw2_pj)-F_sw2_pj-1-0.5*mmd(F_sw2_pj-1-F_sw2_pj-2,F_sw2_pj-F_sw2_pj-1)+F_sw2_nj+1-0.5*mmd(F_sw2_nj-F_sw2_nj+1,F_sw2_nj+1-F_sw2_nj+2)-F_sw2_nj+0.5*mmd(F_sw2_nj-1-F_sw2_nj,F_sw2_nj-F_sw2_nj+1);u3j=u3j-(delt_t/delt_x)*(F_sw3_pj+0.5*mmd(F_sw3_pj-F_sw3_pj-1,F_sw3_pj+1-F_sw3_pj)-F_sw3_pj-1-0.5*mmd(F_sw3_pj-1-F_sw3_pj-2,F_sw3_pj-F_sw3_pj-1)+F_sw3_nj+1-0.5*mmd(F_sw3_nj-F_sw3_nj+1,F_sw3_nj+1-F_sw3_nj+2)-F_sw3_nj+0.5*mmd(F_sw3_nj-1-F_sw3_nj,F_sw3_nj-F_sw3_nj+1); */for(j=0;j<dot_x;j+)rhoj=u1j;vloj=u2j/u1j;prej=(gama-1)*(u3j-0.5*u2j*u2j/u1j);sonj=sqrt(gama*prej/rhoj);/*Roe格式*/*double UR1dot_x,UL1dot_x,UR2dot_x,UL2dot_x,UR3dot_x,UL3dot_x;/UR和UL double f1_rhdot_x,f1_lhdot_x,f2_rhdot_x,f2_lhdot_x,f3_rhdot_x,f3_lhdot_x;/f(i+1/2)和f(i-1/2) double rhoRdot_x,rhoLdot_x,vloRdot_x,vloLdot_x,entRdot_x,entLdot_x,preRdot_x,preLdot_x;double u1_avdot_x,u2_avdot_x,u3_avdot_x;/Roe平均参数 double rho_avdot_x,vlo_avdot_x,ent_avdot_x,egy_avdot_x,pre_avdot_x,son_avdot_x;double s133,s233,a33;double temp1,temp2,temp3,lam1,lam2,lam3,s,eps=1e-6,delt_p,delt_n;int l=0,m=0,n=0;for(i=1;i*delt_t<=0.1;i+)/计算f(j+1/2)/一阶迎风 for(j=0;j<dot_x-1;j+)UR1j=u1j+1;UR2j=u2j+1;UR3j=u3j+1;UL1j=u1j;UL2j=u2j;UL3j=u3j;/MUSCL重构for(j=2;j<dot_x-2;j+)delt_n=u1j+1-u1j;delt_p=u1j+2-u1j+1;s=(2*delt_n*delt_p+eps)/(delt_n*delt_n+delt_p*delt_p+eps);UR1j=u1j+1+0.25*s*(1-s/3)*delt_n+(1+s/3)*delt_p);delt_n=u2j+1-u2j;delt_p=u2j+2-u2j+1;s=(2*delt_n*delt_p+eps)/(delt_n*delt_n+delt_p*delt_p+eps);UR2j=u2j+1+0.25*s*(1-s/3)*delt_n+(1+s/3)*delt_p);delt_n=u3j+1-u3j;delt_p=u3j+2-u3j+1;s=(2*delt_n*delt_p+eps)/(delt_n*delt_n+delt_p*delt_p+eps);UR3j=u3j+1+0.25*s*(1-s/3)*delt_n+(1+s/3)*delt_p);delt_n=u1j-u1j-1;delt_p=u1j+1-u1j;s=(2*delt_n*delt_p+eps)/(delt_n*delt_n+delt_p*delt_p+eps);UL1j=u1j+0.25*s*(1-s/3)*delt_n+(1+s/3)*delt_p);delt_n=u2j-u2j-1;delt_p=u2j+1-u2j;s=(2*delt_n*delt_p+eps)/(delt_n*delt_n+delt_p*delt_p+eps);UL1j=u2j+0.25*s*(1-s/3)*delt_n+(1+s/3)*delt_p);delt_n=u3j-u3j-1;delt_p=u3j+1-u3j;s=(2*delt_n*delt_p+eps)/(delt_n*delt_n+delt_p*delt_p+eps);UL3j=u3j+0.25*s*(1-s/3)*delt_n+(1+s/3)*delt_p); for(j=2;j<dot_x-2;j+)rhoRj=UR1j;rhoLj=UL1j;vloRj=UR2j/UR1j;vloLj=UL2j/UL1j;entRj=(UR3j+(gama-1)*(UR3j-0.5*UR2j*UR2j/UR1j)/UR1j;entLj=(UL3j+(gama-1)*(UL3j-0.5*UL2j*UL2j/UL1j)/UL1j;preRj=(gama-1)*(UR3j-0.5*UR2j*UR2j/UR1j);preLj=(gama-1)*(UL3j-0.5*UL2j*UL2j/UL1j);for(j=2;j<dot_x-2;j+)rho_avj=(sqrt(rhoRj)+sqrt(rhoLj)/2)*(sqrt(rhoRj)+sqrt(rhoLj)/2);vlo_avj=(1/sqrt(rho_avj)*(sqrt(rhoLj)*vloLj+sqrt(rhoRj)*vloRj)*0.5;ent_avj=(1/sqrt(rho_avj)*(sqrt(rhoLj)*entLj+sqrt(rhoRj)*entRj)*0.5;u1_avj=rho_avj;u2_avj=rho_avj*vlo_avj;u3_avj=(rho_avj*ent_avj+0.5*(gama-1)*(rho_avj*vlo_avj*vlo_avj)/gama;pre_avj=(gama-1)*(u3_avj-0.5*u2_avj*u2_avj/u1_avj);son_avj=sqrt(gama*pre_avj/rho_avj);for(j=2;j<dot_x-2;j+)lam1=vlo_avj;lam2=vlo_avj-son_avj;lam3=vlo_avj+son_avj;s100=(0.5*vlo_avj*vlo_avj-son_avj*son_avj/(gama-1)*fabs(lam1);s101=(-vlo_avj)*fabs(lam1);s102=fabs(lam1);s110=(-vlo_avj-(gama-1)*vlo_avj*vlo_avj/(son_avj*2)*fabs(lam2);s111=(1+(gama-1)*vlo_avj/son_avj)*fabs(lam2);s112=(-(gama-1)/son_avj)*fabs(lam2);s120=(-vlo_avj+(gama-1)*vlo_avj*vlo_avj/(son_avj*2)*fabs(lam3);s121=(1-(gama-1)*vlo_avj/son_avj)*fabs(lam3);s122=(gama-1)/son_avj)*fabs(lam3);s200=-(gama-1)/(son_avj*son_avj);s201=-1/(2*son_avj);s202=1/(2*son_avj);s210=-vlo_avj*(gama-1)/(son_avj*son_avj);s211=-(vlo_avj-son_avj)/(2*son_avj);s212=(vlo_avj+son_avj)/(2*son_avj);s220=-(gama-1)*vlo_avj*vlo_avj/(2*son_avj*son_avj);s221=-(1/(2*son_avj)*(vlo_avj*vlo_avj*0.5+son_avj*son_avj/(gama-1)-vlo_avj*son_avj); s222=(1/(2*son_avj)*(vlo_avj*vlo_avj*0.5+son_avj*son_avj/(gama-1)+vlo_avj*son_avj);for(l=0;l<3;l+)for(m=0;m<3;m+)alm=0; for(l=0;l<3;l+)for(m=0;m<3;m+)for(n=0;n<3;n+)alm=alm+s2ln*s1nm;temp1=UR1j-UL1j;temp2=UR2j-UL2j;temp3=UR3j-UL3j; f1_rhj=0.5*(rhoRj*vloRj+rhoLj*vloLj)-0.5*(a00*temp1+a01*temp2+a02*temp3);f2_rhj=0.5*(rhoRj*vloRj*vloRj+preRj+rhoLj*vloLj*vloLj+preLj)-0.5*(a10*temp1+a11*temp2+a12*temp3);f3_rhj=0.5*(vloRj*(UR3j+preRj)+vloLj*(UL3j+preLj)-0.5*(a20*temp1+a21*temp2+a22*temp3);/计算f(i-1/2) /一阶迎风 for(j=1;j<dot_x;j+)UR1j=u1j;UR2j=u2j;UR3j=u3j;UL1j=u1j-1;UL2j=u2j-1;UL3j=u3j-1;/MUSCL重构for(j=2;j<dot_x-2;j+)delt_n=u1j-u1j-1;delt_p=u1j+1-u1j;s=(2*delt_n*delt_p+eps)/(delt_n*delt_n+delt_p*delt_p+eps);UR1j=u1j+0.25*s*(1-s/3)*delt_n+(1+s/3)*delt_p);delt_n=u2j-u2j-1;delt_p=u2j+1-u2j;s=(2*delt_n*delt_p+eps)/(delt_n*delt_n+delt_p*delt_p+eps);UR2j=u2j+0.25*s*(1-s/3)*delt_n+(1+s/3)*delt_p);delt_n=u3j-u3j-1;delt_p=u3j+1-u3j;s=(2*delt_n*delt_p+eps)/(delt_n*delt_n+delt_p*delt_p+eps);UR3j=u3j+0.25*s*(1-s/3)*delt_n+(1+s/3)*delt_p);delt_n=u1j-1-u1j-2;delt_p=u1j-u1j-1;s=(2*delt_n*delt_p+eps)/(delt_n*delt_n+delt_p*delt_p+eps);UL1j=u1j-1+0.25*s*(1-s/3)*delt_n+(1+s/3)*delt_p);delt_n=u2j-1-u2j-2;delt_p=u2j-u2j-1;s=(2*delt_n*delt_p+eps)/(delt_n*delt_n+delt_p*delt_p+eps);UL2j=u2j-1+0.25*s*(1-s/3)*delt_n+(1+s/3)*delt_p);delt_n=u3j-1-u3j-2;delt_p=u3j-u3j-1;s=(2*delt_n*delt_p+eps)/(delt_n*delt_n+delt_p*delt_p+eps);UL3j=u3j-1+0.25*s*(1-s/3)*delt_n+(1+s/3)*delt_p); for(j=2;j<dot_x-2;j+)rhoRj=UR1j;rhoLj=UL1j;vloRj=UR2j/UR1j;vloLj=UL2j/UL1j;entRj=(UR3j+(gama-1)*(UR3j-0.5*UR2j*UR2j/UR1j)/UR1j;entLj=(UL3j+(gama-1)*(UL3j-0.5*UL2j*UL2j/UL1j)/UL1j;preRj=(gama-1)*(UR3j-0.5*UR2j*UR2j/UR1j);preLj=(gama-1)*(UL3j-0.5*UL2j*UL2j/UL1j);for(j=2;j<dot_x-2;j+)rho_avj=(sqrt(rhoRj)+sqrt(rhoLj)/2)*(sqrt(rhoRj)+sqrt(rhoLj)/2);vlo_avj=(1/sqrt(rho_avj)*(sqrt(rhoLj)*vloLj+sqrt(rhoRj)*vloRj)*0.5;ent_avj=(1/sqrt(rho_avj)*(sqrt(rhoLj)*entLj+sqrt(rhoRj)*entRj)*0.5;u1_avj=rho_avj;u2_avj=rho_avj*vlo_avj;u3_avj=(rho_avj*ent_avj+0.5*(gama-1)*(rho_avj*vlo_avj*vlo_avj)/gama;pre_avj=(gama-1)*(u3_avj-0.5*u2_avj*u2_avj/u1_avj);son_avj=sqrt(gama*pre_avj/rho_avj);for(j=2;j<dot_x-2;j+)lam1=vlo_avj;lam2=vlo_avj-son_avj;lam3=vlo_avj+son_avj;s100=(0.5*vlo_avj*vlo_avj-son_avj*son_avj/(gama-1)*fabs(lam1);s101=(-vlo_avj)*fabs(lam1);s102=fabs(lam1);s110=(-vlo_avj-(gama-1)*vlo_avj*vlo_avj/(son_avj*2)*fabs(lam2);s111=(1+(gama-1)*vlo_avj/son_avj)*fabs(lam2);s112=(-(gama-1)/son_avj)*fabs(lam2);s120=(-vlo_avj+(gama-1)*vlo_avj*vlo_avj/(son_avj*2)*fabs(lam3);s121=(1-(gama-1)*vlo_avj/son_avj)*fabs(lam3);s122=(gama-1)/son_avj)*fabs(lam3);s200=-(gama-1)
