附录B一维可压缩黏性流动问题的数值解法与计算程序

附录B 一维可压缩黏性流动问题的数值解法与计算程序一维可压缩黏性流动是气体动力学中最经典的黏性流动问题，对它采用迎风型TVD 差分算法进行数值求解。

同时，为了初学者入门和练习方便,这里给出了由C 语言和Fortran77语言编写的、计算一维可压缩黏性流动问题的计算程序，供大家学习参考。

B-1 利用二阶迎风型TVD 差分格式求解一维可压缩黏性流动问题1.一维可压缩黏性流动问题在两端开口管道中充满了可压缩黏性流体。

当黏性流体以超声速从左向右运动时，一定会在管道中形成一道正激波，如图B.1所示。

111u p ρ、、和222u p ρ、、分别为激波波前和波后的参数。

该问题可简化为一维可压缩黏性流动问题。

当数值解达到稳定时，在管道中可求解得到一道稳定的激波。

2.基本方程组、初始条件和边界条件设流体是黏性流体。

一维可压缩黏性流动问题，在数学上可以用一维可压缩黏性流动N-S 方程组来描述。

量纲为一的一维N-S 方程组为：t x x∂∂∂+=∂∂∂u f s（B.1）其中()122304,,343p s u u u u p s Re x E E p u s C K u u T Re x Pr Rr x ρρμρρμ⎛⎫ ⎪ ⎪⎛⎫⎛⎫⎡⎤⎪∂ ⎪ ⎪⎢⎥==+== ⎪ ⎪ ⎪⎢⎥∂ ⎪ ⎪ ⎪⎢⎥+⎣⎦⎝⎭⎝⎭⎪∂∂ ⎪+ ⎪∂∂⎝⎭u f s（B.1a ） 图B.1一维可压缩黏性流动问题示意图()2221,,21,,,1p v p v v c Tp E C T u Pr M k c U L Re c a M c M μρργργμγγ∞∞∞∞∞∞∞∞⎛⎫==+= ⎪⎝⎭====- （B.1b ）其中u p ρ、、和E 分别是量纲为一的密度、速度、压力和单位体积总能，s 为流体的黏性项。

Pr 为普朗特数（此处公式中p c k μ∞∞、、是有量纲量），Re 为雷诺数，v c 为比定容热容，p c 为比定压热容，,p v c c 是量纲为一的量，称为气体绝热指数，a 为当地声速。

求解区域为10x ≤≤。

取2,50,0.72, 1.4,M Re Pr k T γμ∞======。

初始条件：在0t =时刻，(),01x ρ=，其他物理量采用线性插值得到。

边界条件：左边界()0x =处：()()()0,0,0,1t u t T t ρ=== （B.2）右边界()1x =处：()()()()12222,0211,1121121,1111x x t xM u t M T t M M ργγγγγγγγγγ=∞∞∞∞∂=∂+-=+-⎛⎫⎛⎫--=-+ ⎪ ⎪ ⎪++++⎝⎭⎝⎭（B.3）3.二阶精度迎风型TVD 差分格式一维N-S 方程组中的对流项采用Harten 的二阶精度迎风型TVD 差分格式：11122n nnn jjj j r t x ++-⎛⎫∂=--+∆ ⎪∂⎝⎭suu h h （B.4） 111122212j j j j j ++++⎛⎫=++ ⎪⎝⎭h f f R Φ （B.4a ）其中t r x∆=∆。

向量Φ在第l 个特征方向上分量lΦ为： 111112222l l l llj j j j j j g g Q a Φγα+++++⎛⎫=+-+ ⎪⎝⎭ （B.4b ）11112222,l l l l l j j j j j g min mod σασα++--⎛⎫= ⎪⎝⎭（B.4c ）111212212,00,0l l j jl l j lj j l j g g ααγα+++++⎧-≠⎪⎪=⎨⎪=⎪⎩ （B.4d ） ()1111222l ll j j j R u α-+++=⨯∆（B.4e ） ()()2,,0.050.51,2z z Q z z z εεεεε⎧≥⎪=≤≤⎛⎫⎨+≤ ⎪⎪⎝⎭⎩ （B.4f ）2112212121212l l i i li l i Q r Q λλσλ++++⎧⎡⎤⎛⎫-⎪⎢⎥ ⎪⎢⎥⎪⎝⎭⎣⎦=⎨⎛⎫⎪ ⎪⎪⎝⎭⎩非定常定常,,（B.4g ）流通量矢量f 的非线性Jiacobian 系数矩阵-1=A R ΛR 为：()()222320103312232112u uu aa u u u γγγγγγγγ⎛⎫⎪ ⎪∂- ⎪==--- ⎪∂ ⎪-- ⎪-+ ⎪--⎝⎭f A u u （B.5）非线性Jiacobian 系数矩阵A 的特征值为：123123000000,,u u a u aλλλλλλ⎛⎫⎪= ⎪ ⎪⎝⎭==-=+Λ （B.6）非线性Jiacobian 系数矩阵A 的右特征矢量-1、R R 为：211112uu a u a u h ua h ua ⎛⎫⎪⎪=-+ ⎪ ⎪-+ ⎪⎝⎭R（B.7）()()()()()()222221222222211112111124222111124222u ua a a u u u aa a a a u u u a a a aa γγγγγγγγγ-⎛⎫-----⎪⎪ ⎪---=+-- ⎪ ⎪ ⎪----+-+ ⎪ ⎪⎝⎭R （B.8）一维N-S 方程组中的黏性项t x∂∆∂s采用二阶精度中心差分格式。

4.计算结果分析采用用C 语言和Fortran77语言对一维可压缩黏性流动问题编制了计算程序，并对雷诺数50Re =的流动进行了计算，计算结果如图B.2和图 B.3所示。

图B.2和图 B.3是计算达到稳定后激波间断位置和密度()ρ、速度()u 、压力()p 和单位质量内能()e 的分布。

由上述计算结果中可以看出，采用二阶精度迎风型TVD 差分格式计算一维可压缩黏性流动问题得到的数值解和经典文献中的结果是完全一致的。

计算结果表明，迎风型TVD 差分格式能够精确地捕捉激波间断，计算效果较好。

由于本问题中黏性较大(50)Re =，所以计算得到的激波比较光滑，有一定的宽度。

一维可压缩黏性流动问题的解是连续、光滑的。

图 B.2 Fortran77语言程序得到的结果 图 B.3 C 语言程序得到的结果B-2 一维可压缩黏性流动问题的数值计算源程序1.C语言源程序// UpwindTVD_1D.c//------------------------------------------------------------------------ // 二阶迎风型TVD差分格式求解一维可压缩黏性流动问题（C语言版本）//------------------------------------------------------------------------- #include "stdio.h"#include "math.h"#define gama 1.4#define Tt 5.0#define im 201 //网格数//全局变量：double Q[3][im],Qold[3][im] ; //Q: [rou, rou*u, E]double rou[im],u[im],p[im],T[im],E[im],a[im] ;double Pr,Re,cv,cp,Ma,dx,dt ;//-------------------------------------------------------------------void initial(){double xl,xr,x;double ul,Tl,ur,Tr;//进出口的u，T值int i;dx=1.0/(im-1) ;dt=1.0e-6 ;Pr=0.72 ;Re=50.0 ;Ma=2.0 ;cv=1.0/(gama*(gama-1.0)*Ma*Ma) ;cp=gama*cv ;xl=0.0 ;xr=1.0 ;ul=1.0 ;Tl=1.0 ;ur=(2.0/(gama-1.0)+Ma*Ma)/((gama+1.0)/(gama-1.0)*Ma*Ma) ;Tr=2.0*gama/(gama+1.0)*Ma*Ma-(gama-1.0)/(gama+1.0) ;Tr=Tr*( (gama-1.0)/(gama+1.0) + 2.0/(gama+1.0)/Ma/Ma ) ;for(i=0; i<=im-1; i++){x=i*dx ;rou[i]=1.0 ;u[i]=(x-xl)/(xr-xl)*(ur-ul)+ul ;T[i]=(x-xl)/(xr-xl)*(Tr-Tl)+Tl ;p[i]=rou[i]*T[i]/(gama*Ma*Ma) ;E[i]=rou[i]*(cv*T[i]+0.5*u[i]*u[i]) ;}for(i=0; i<im; i++){Q[0][i]=rou[i] ;Q[1][i]=rou[i]*u[i] ;Q[2][i]=E[i] ;}}//------------------------------------------------------------------- double qz(double x){double result ;if (fabs(x)>0.1){result=fabs(x) ;}else{result=(x*x/0.1+0.1)/2.0 ;}return result ;}//------------------------------------------------------------------- double minmod(double x,double y){double result ;if(x*y<=0.0){result=0.0 ;}else if(fabs(x)>fabs(y)){result=y ;}else if(fabs(x)<fabs(y)){result=x ;}return result ;}void Q2U(){int i ;for(i=0; i<im; i++){rou[i]=Q[0][i] ;u[i]=Q[1][i]/rou[i] ;E[i]=Q[2][i] ;T[i]=(E[i]/rou[i]-0.5*u[i]*u[i])/cv ;p[i]=rou[i]*T[i]/(gama*Ma*Ma) ;a[i]=sqrt(T[i])/Ma ;}}//------------------------------------------------------------------- void UpwindTVD_1D_Solver(){int i,l,k ;double dq[3][im],lamda[3][im],sigma[3][im],f[3][im] ;double fl[3],fwave[3][im],fr[3],g[3][im],Gv[3][im] ;double alfa[3][im],beta[3][im],Rn[3][3][im],R[3][3][im] ;double ut,at ; //临时变量Q2U() ;for(i=0; i<im; i++){f[0][i]=rou[i]*u[i] ;f[1][i]=rou[i]*u[i]*u[i]+p[i] ;f[2][i]=u[i]*(E[i]+p[i]) ;}for(i=0; i<=im-2; i++){lamda[0][i]=(u[i]+u[i+1])/2.0 ;lamda[1][i]=(u[i]-a[i]+u[i+1]-a[i+1])/2.0 ;lamda[2][i]=(u[i]+a[i]+u[i+1]+a[i+1])/2.0 ;}for(i=0; i<=im-2; i++)for(l=0; l<=2; l++){dq[l][i]=Q[l][i+1]-Q[l][i] ;sigma[l][i]=qz(lamda[l][i])/2.0 ;}//计算左右特征向量for(i=0; i<=im-2; i++){ut=(u[i]+u[i+1])/2.0 ;at=(a[i]+a[i+1])/2.0 ;R[0][0][i]=1.0 ;R[0][1][i]=1.0 ;R[0][2][i]=1.0 ;R[1][0][i]=ut ;R[1][1][i]=ut-at ;R[1][2][i]=ut+at ;R[2][0][i]=ut*ut/2.0 ;R[2][1][i]=at*at/(gama-1.0)+ut*ut/2.0-ut*at ;R[2][2][i]=at*at/(gama-1.0)+ut*ut/2.0+ut*at ;Rn[0][0][i]=1.0-(gama-1.0)*ut*ut/2.0/at/at ;Rn[0][1][i]=(gama-1.0)*ut/at/at ;Rn[0][2][i]=-(gama-1.0)/at/at ;Rn[1][0][i]=ut/2.0/at+(gama-1.0)*ut*ut/4.0/at/at ;Rn[1][1][i]=-(gama-1.0)*ut/2./at/at-1.0/2.0/at ;Rn[1][2][i]=(gama-1.0)/2.0/at/at ;Rn[2][0][i]=-ut/2.0/at+(gama-1.0)*ut*ut/4.0/at/at ;Rn[2][1][i]=(1.0-gama)*ut/2./at/at+1.0/2.0/at ;Rn[2][2][i]=(gama-1.0)/2./at/at ;}for(i=0; i<=im-1; i++)for(l=0; l<=2; l++){alfa[l][i]=0.0 ;fwave[l][i]=0.0 ;g[l][i]=0.0 ;}//计算alfafor(i=0; i<=im-2; i++)for(l=0; l<=2; l++)for(k=0; k<=2; k++){alfa[l][i]=alfa[l][i]+Rn[l][k][i]*dq[k][i] ;}// 计算gfor(i=1; i<=im-2; i++)for(l=0; l<=2; l++){g[l][i]=minmod(sigma[l][i]*alfa[l][i],sigma[l][i-1]*alfa[l][i-1]) ;}for(i=0; i<=im-2; i++)for(l=0; l<=2; l++){if(alfa[l][i]!=0.0){beta[l][i]=(g[l][i+1]-g[l][i])/alfa[l][i] ;}else{beta[l][i]=0.0 ;}}for(i=0; i<=im-2; i++)for(l=0; l<=2; l++)for(k=0; k<=2; k++){fwave[l][i]=fwave[l][i]+0.5*R[l][k][i]*(g[k][i]+g[k][i+1]-qz(lamda[k][i]+beta[k][i])*alfa[k][i] ) ;}for(i=1; i<=im-2; i++)for(l=0; l<=2; l++){fl[l]=0.5*(f[l][i-1]+f[l][i])+fwave[l][i-1] ;fr[l]=0.5*(f[l][i+1]+f[l][i])+fwave[l][i] ;Q[l][i]=Q[l][i]-dt/dx*(fr[l]-fl[l]) ;}//黏性项的处理for(i=1; i<=im-2; i++){Gv[0][i]=0.0 ;Gv[1][i]=(T[i]*(u[i+1]-2.*u[i]+u[i-1])+(T[i+1]-T[i-1])*(u[i+1]-u[i-1])/4.)*4./3./Re/(dx*dx) ;Gv[2][i]=((T[i]*(u[i+1]-2.*u[i]+u[i-1])*u[i]+(T[i+1]-T[i-1])* (u[i+1]-u[i-1])*u[i]/4.+T[i]*(u[i+1]-u[i-1])*(u[i+1]-u[i-1])/4.)*4./3.+(T[i]*(T[i+1]-2.*T[i]+T[i-1])+(T[i+1]-T[i-1])*(T[i+1]-T[i-1])/4.)*cp/Pr)/Re/(dx*dx) ;}for(i=1; i<=im-2; i++)for(l=0; l<=2; l++){Q[l][i]=Q[l][i]+Gv[l][i]*dt ;}//边界条件Q[0][im-1]=Q[0][im-2] ;Q[1][im-1]=Q[0][im-1]*u[im-1] ;Q[2][im-1]=Q[0][im-1]*(cv*T[im-1]+0.5*u[im-1]*u[im-1]) ;}//------------------------------------------------------------------- double error(double Q1[3][im],double Q2[3][im]){int i,l ;double err,ddu ;err=1.0e-10 ;for(i=0; i<im; i++)for(l=0; l<3 ; l++){if(Q2[l][i]!=0.0){ddu=fabs((Q2[l][i]-Q1[l][i])/Q2[l][i]) ;}else{ddu=fabs(Q2[l][i]-Q1[l][i]) ;}if(err<ddu) err=ddu ;}return err ;}//------------------------------------------------------------------- void output(){double x,rou1,u1,E1,p1,T1 ; //临时变量int i ;FILE *fp ;fp=fopen("result.dat","w") ;fprintf(fp,"variables=x,rou,u,p,e\n");for(i=0; i<im; i++){x=i*dx ;rou1=Q[0][i] ;u1=Q[1][i]/rou1 ;E1=Q[2][i] ;T1=(E1/rou1-0.5*u1*u1)/cv ;p1=rou1*T1/(gama*Ma*Ma) ;fprintf(fp,"%15.8e %15.8e %15.8e %15.8e %15.8e\n",x,rou1,u1,p1,cv*T1 ) ;}fclose(fp) ;}//-------------------------------------------------------------------main(){double t,err ;int n,i,l ;initial() ;t=0.0 ;n=0 ;err=1.0e-7 ;while(t<Tt && err>1.0e-8){t=t+dt ;n=n+1 ;for(i=0; i<im; i++)for(l=0; l<3 ; l++){Qold[l][i]=Q[l][i] ;}UpwindTVD_1D_Solver() ;err=error(Q,Qold) ;if(n/10000*10000==n){printf("%10d time: %16.9e error: %16.9e\n",n,t,err) ;output() ;}}output() ;printf("Program end !\n") ;}-----------------------------------------2. Fortran77语言源程序! UpwindTVD_1D.f!----------------------------------------------------------------------- ! 二阶迎风型TVD差分格式求解一维可压缩黏性流动问题!（Fortran77语言版本）!----------------------------------------------------------------------- program UPWIND_TVD_1Dimplicit real*8 (a-h,o-z)parameter(mx=201,Tt=5.0 )dimension Q(3,mx),Qold(3,mx) ! Q: [rou, rou*u, E]dimension rou(mx),u(mx),p(mx),T(mx),E(mx)real*8 Macommon/para_def/Sf,Re,Ma,pr,gama,cp,cv,dx,dt,imim=mxcall Initialize(Q)time=0.0n=01 n=n+1time=time+dtdo i=1,imdo l=1,3Qold(l,i)=Q(l,i)end doend docall UpwindTVD_1D_Solver(Q)if(n/10000*10000.eq.n)thenwrite(*,20)n,time,error(Q,Qold)call OutPut(Q)end if20 format(1x,i10,' time=',e16.9,' error=',e16.9,1x)if((time.lt.Tt) .and. (error(Q,Qold).gt.1.0e-8))thengoto 1end ifcall OutPut(Q)write(*,*) 'Program end !'end!-------------------------------------------------------------------- subroutine Initialize(Q)implicit real*8 (a-h,o-z)real*8 Macommon/para_def/Sf,Re,Ma,pr,gama,cp,cv,dx,dt,imdimension rou(im),u(im),p(im),T(im),E(im),Q(3,im)Re=50.0 ! 雷诺数Ma=2.0 ! 马赫数pr=0.72 ! 普朗特数gama=1.4 ! 气体常数cv=1.0/(gama*(gama-1.0)*Ma**2) ! 比定容热容cp=gama*cv ! 比定压热容dt=1.0e-6dx=1.0/(1.0*(im-1)) ! 空间步长xl=0.0xr=1.0rou0=1.0u0=1.0T0=1.0temp=(gama+1.)/(gama-1.)u1=(2./(gama-1.)+Ma**2)/(temp*Ma**2)T1=(2.*gama/(gama+1)*Ma**2-1./temp)*(1./temp+2./(gama+1.)/Ma**2) do i=1,imrou(i)=rou0xi=(i-1)*dxu(i)=(xi-xl)/(xr-xl)*(u1-u0)+u0T(i)=(xi-xl)/(xr-xl)*(T1-T0)+T0P(i)=rou(i)*T(i)/(gama*Ma**2)E(i)=rou(i)*(cv*T(i)+0.5*u(i)**2)end dodo i=1,imQ(1,i)=rou(i)Q(2,i)=rou(i)*u(i)Q(3,i)=E(i)end doend!-------------------------------------------------------------------- subroutine Q2U(Q,rou,u,p,T,E,a)implicit real*8 (a-h,o-z)real*8 Macommon/para_def/Sf,Re,Ma,pr,gama,cp,cv,dx,dt,imdimension Q(3,im),rou(im),u(im)dimension p(im),T(im),E(im),a(im)do i=1,imrou(i)=Q(1,i)u(i)=Q(2,i)/Q(1,i)E(i)=Q(3,i)T(i)=(E(i)/rou(i)-0.5*u(i)**2)/cvp(i)=rou(i)*T(i)/(gama*Ma**2)a(i)=dsqrt(T(i))/Maend doend!-------------------------------------------------------------------- subroutine UpwindTVD_1D_Solver(Q)implicit real*8 (a-h,o-z)real*8 Ma,lamda,minmodcommon/para_def/Sf,Re,Ma,pr,gama,cp,cv,dx,dt,imdimension Q(3,im),Qold(3,im),dq(3,im),rou(im),u(im),p(im),T(im)dimension E(im),a(im),f(3,im),fwave(3,im),fr(3),fl(3) dimension Gv(3,im),g(3,im),lamda(3,im),sigma(3,im)dimension alfa(3,im),Rn(3,3,im),R(3,3,im),beta(3,im)call Q2U(Q,rou,u,p,T,E,a)do i=1,imf(1,i)=rou(i)*u(i)f(2,i)=rou(i)*u(i)**2+p(i)f(3,i)=u(i)*(E(i)+p(i))end dodo i=1,im-1lamda(1,i)=(u(i)+u(i+1))/2.0lamda(2,i)=(u(i)-a(i)+u(i+1)-a(i+1))/2.0lamda(3,i)=(u(i)+a(i)+u(i+1)+a(i+1))/2.0end dodo i=1,im-1do l=1,3dq(l,i)=(Q(l,i+1)-Q(l,i))sigma(l,i)=qz(lamda(l,i))/2.0 ! 定常情况end doend dodo i=1,im-1ut=(u(i)+u(i+1))/2.0 !ut 临时变量at=(a(i)+a(i+1))/2.0 !at 临时变量R(1,1,i)=1.0R(1,2,i)=1.0R(1,3,i)=1.0R(2,1,i)=utR(2,2,i)=ut-atR(2,3,i)=ut+atR(3,1,i)=ut*ut/2.0R(3,2,i)=at*at/(gama-1.0)+ut*ut/2.0-ut*atR(3,3,i)=at*at/(gama-1.0)+ut*ut/2.0+ut*atRn(1,1,i)=1.0-(gama-1.0)*ut*ut/2.0/at/atRn(1,2,i)=(gama-1.0)*ut/at/atRn(1,3,i)=-(gama-1.0)/at/atRn(2,1,i)=ut/2.0/at+(gama-1.0)*ut*ut/4.0/at/atRn(2,2,i)=-(gama-1.0)*ut/2./at/at-1.0/2.0/atRn(2,3,i)=(gama-1.0)/2.0/at/atRn(3,1,i)=-ut/2.0/at+(gama-1.0)*ut*ut/4.0/at/atRn(3,2,i)=(1.0-gama)*ut/2./at/at+1.0/2.0/atRn(3,3,i)=(gama-1.0)/2./at/atend dodo l=1,3alfa(l,i)=0.0fwave(l,i)=0.0g(l,i)=0.0end doend dodo i=1,im-1do l=1,3do k=1,3alfa(l,i)=alfa(l,i)+Rn(l,k,i)*dq(k,i)end doend doend dodo i=2,im-1do l=1,3g(l,i)=minmod(sigma(l,i)*alfa(l,i),sigma(l,i-1)*alfa(l,i-1)) end doend dodo i=1,im-1do l=1,3if(alfa(l,i).ne.0.0)thenbeta(l,i)=(g(l,i+1)-g(l,i))/alfa(l,i)elsebeta(l,i)=0.0end ifend doend dodo i=1,im-1do l=1,3do k=1,3fwave(l,i)=fwave(l,i)+0.5*R(l,k,i)*(g(k,i)+g(k,i+1)* -qz(lamda(k,i)+beta(k,i))*alfa(k,i))end doend doend dodo i=2,im-1do l=1,3fl(l)= 0.5*(f(l,i-1)+f(l,i)) +fwave(l,i-1)fr(l)= 0.5*(f(l,i+1)+f(l,i)) +fwave(l,i)Q(l,i)=Q(l,i)-dt/dx*(fr(l)-fl(l))end doend do! 处理黏性项Gv(1,i)=0.0Gv(2,i)=(T(i)*(u(i+1)-2.*u(i)+u(i-1))+(T(i+1)-T(i-1))*& (u(i+1)-u(i-1))/4.)*4./3./re/dx**2Gv(3,i)=((T(i)*(u(i+1)-2.*u(i)+u(i-1))*u(i)+(T(i+1)-T(i-1))*& (u(i+1)-u(i-1))*u(i)/4.+T(i)*(u(i+1)-u(i-1))**2/4.)*4./3. & +(T(i)*(T(i+1)-2.*T(i)+T(i-1))+(T(i+1)-T(i-1))**2/4.)& *cp/pr)/re/dx**2end dodo i=2,im-1do l=1,3Q(l,i)=Q(l,i)+Gv(l,i)*dtend doend do! 边界条件Q(1,im)=Q(1,im-1)Q(2,im)=Q(1,im-1)*u(im)Q(3,im)=Q(1,im-1)*(cv*T(im)+0.5*u(im)**2)end!-------------------------------------------------------------------- function error(Q,Qold)implicit real*8 (a-h,o-z)real*8 Ma,err,errorcommon/para_def/Sf,Re,Ma,pr,gama,cp,cv,dx,dt,imdimension Q(3,im),Qold(3,im)err=1.0e-10do i=1,imdo l=1,3if(Qold(l,i).ne.0.0)thenddq=abs((Q(l,i)-Qold(l,i))/Qold(l,i))elseddq=abs(Q(l,i)-Qold(l,i))end ifif(err.lt.ddq) err=ddqend doend doerror=errend!-------------------------------------------------------------------- function qz(x)real*8 x,qzif (abs(x).gt.0.1) thenqz=abs(x)elseqz=(x**2/0.1+0.1)/2.endifend function qz!---------------------------------------------------------------------------- function minmod(x,y)real*8 x,y,minmodif(x*y.le.0.0) thenminmod=0.0else if(abs(x).gt.abs(y)) thenminmod=yelseminmod=xendifend function minmod!--------------------------------------------------------------------subroutine OutPut(Q)implicit real*8 (a-h,o-z)real*8 Macommon/para_def/Sf,Re,Ma,pr,gama,cp,cv,dx,dt,imdimension Q(3,im)open(1,file='result.dat')write(1,*)'variables=x,rou,u,p,e'do i=1,imx=(i-1)*dxrou=Q(1,i)u=Q(2,i)/rouE=Q(3,i)T=(E/rou-0.5*u*u)/cvp=rou*T/(gama*Ma**2)write(1,10)x,rou,u,p,cv*Tend doclose(1)10 format(1x,f5.3,1x,3(e15.8,1x))end------------------------------------------------------------------------ ----------------------------------------------------------------------- -----------------------------------------------------------------------。