Fortran语言源程序_flow_两步差分格式_二维斜激波在平面刚壁上反射问题的数值解_compressible_

program Lax_Wendroff2D implicit double precision (a-h,o-z) parameter (M1=400,M2=100) common /G_def/ GAMA,PI,Jx,Jy,dLx,dLy,TT,Sf dimension U(0:M1+1,0:M2+1,0:3),Uf(0:M1+1,0:M2+1,0:3) dimension EFf(0:M1+1,0:M2+1,0:3) !气体常数 GAMA=1.4 PI=3.1415926 !网格数 Jx=M1 Jy=M2 !计算区域 dLx=4.0 dLy=1.0 !总时间 TT=3.0 !时间步长因子 Sf=0.8 call Init(U,dx,dy) T=0 1 dt=CFL(U,dx,dy) T=T+dt write(*,*)'T=',T,'dt=',dt call Lax_Wendroff_2DSolver(U,Uf,EFf,dx,dy,dt) if(T.lt.TT)goto 1 call Output(U,dx,dy) end !------------------------------------------------------- !计算时间步长 !入口: U， 当前物理量，dx、dy, 网格宽度； !返回: 时间步长。 !------------------------------------------------------- double precision function CFL(U,dx,dy) implicit double precision (a-h,o-z) common /G_def/ GAMA,PI,Jx,Jy,dLx,dLy,TT,Sf dimension U(0:Jx+1,0:Jy+1,0:3) dmaxvel=1e-10 do 10 i=1,Jx do 10 j=1,Jy uu=U(i,j,1)/U(i,j,0) vv=U(i,j,2)/U(i,j,0) p=(GAMA-1)*(U(i,j,3)-0.5*U(i,j,0)*(uu*uu+vv*vv)) vel=dsqrt(GAMA*p/U(i,j,0))+dsqrt(uu*uu+vv*vv) if(vel.gt.dmaxvel)dmaxvel=vel 10 continue CFL=Sf*dmin1(dx,dy)/dmaxvel end !------------------------------------------------------- !初始化 !入口: 无； !出口: U, 已经给定的初始值，dx、dy, 网格宽度 !------------------------------------------------------- subroutine Init(U,dx,dy) implicit double precision (a-h,o-z) common /G_def/ GAMA,PI,Jx,Jy,dLx,dLy,TT,Sf dimension U(0:Jx+1,0:Jy+1,0:3) !初始条件 rou1=1.0 u1=2.9 v1=0 p1=0.71429 rou2=1.69997 u2=2.61934 v2=-0.50632 p2=1.52819 dx=dLx/Jx dy=dLy/Jy !入射激波角度 theta=29*PI/180 do 20 j=0,Jy+1 do 20 i=0,Jx+1 tx=(1.0-j*dy)/tan(theta) if(i*dx.le.tx)then U(i,j,0)=rou1 U(i,j,1)=rou1*u1 U(i,j,2)=rou1*v1 U(i,j,3)=p1/(GAMA-1)+rou1*(u1*u1+v1*v1)/2 else U(i,j,0)=rou2 U(i,j,1)=rou2*u2 U(i,j,2)=rou2*v2 U(i,j,3)=p2/(GAMA-1)+rou2*(u2*u2+v2*v2)/2 endif 20 continue end !------------------------------------------------------- !边界条件 !入口: dx、dy， 网格宽度； !出口: U，已经给定边界。 !------------------------------------------------------- subroutine bound(U,dx,dy) implicit double precision (a-h,o-z) common /G_def/ GAMA,PI,Jx,Jy,dLx,dLy,TT,Sf dimension U(0:Jx+1,0:Jy+1,0:3) !初始条件 rou1=1.0 u1=2.9 v1=0 p1=0.71429 rou2=1.69997 u2=2.61934 v2=-0.50632 p2=1.52819 !左边界 do 30 j=0,Jy+1 U(0,j,0)=rou1 U(0,j,1)=rou1*u1 U(0,j,2)=rou1*v1 U(0,j,3)=p1/(GAMA-1)+rou1*(u1*u1+v1*v1)/2 30 continue !右边界 do 31 j=0,Jy+1 do 31 k=0,3 U(Jx+1,j,k)=U(Jx,j,k) 31 continue !上边界 do 32 i=0,Jx+1 U(i,Jy+1,0)=rou2 U(i,Jy+1,1)=rou2*u2 U(i,Jy+1,2)=rou2*v2 U(i,Jy+1,3)=p2/(GAMA-1)+rou2*(u2*u2+v2*v2)/2 32 continue !下边界 do 33 i=0,Jx+1 U(i,0,0)=U(i,1,0) U(i,0,1)=U(i,1,1) U(i,0,2)=0 U(i,0,3)=U(i,1,3) 33 continue end !------------------------------------------------------- !根据U计算E !入口: U， 当前U矢量； !出口: E， 计算得到的E矢量， ! U、E定义见Euler方程组。 !------------------------------------------------------- subroutine U2E(U,E,is,in,js,jn) implicit double precision (a-h,o-z) common /G_def/ GAMA,PI,Jx,Jy,dLx,dLy,TT,Sf dimension U(0:Jx+1,0:Jy+1,0:3),E(0:Jx+1,0:Jy+1,0:3) do 40 i=is,in do 40 j=js,jn uu=U(i,j,1)/U(i,j,0) vv=U(i,j,2)/U(i,j,0) p=(GAMA-1)*(U(i,j,3) $ -0.5*(U(i,j,1)*U(i,j,1)+U(i,j,2)*U(i,j,2))/U(i,j,0)) E(i,j,0)=U(i,j,1) E(i,j,1)=U(i,j,0)*uu*uu+p E(i,j,2)=U(i,j,0)*uu*vv E(i,j,3)=(U(i,j,3)+p)*uu 40 continue end !------------------------------------------------------- !根据U计算F !入口: U， 当前U矢量； !出口: E，计算得到的F矢量， ! U、F定义见Euler方程组。 !------------------------------------------------------- subroutine U2F(U,F,is,in,js,jn) implicit double precision (a-h,o-z) common /G_def/ GAMA,PI,Jx,Jy,dLx,dLy,TT,Sf dimension U(0:Jx+1,0:Jy+1,0:3),F(0:Jx+1,0:Jy+1,0:3) do 50 i=is,in do 50 j=js,jn uu=U(i,j,1)/U(i,j,0) vv=U(i,j,2)/U(i,j,0) p=(GAMA-1)*(U(i,j,3) $ -0.5*(U(i,j,1)*U(i,j,1)+U(i,j,2)*U(i,j,2))/U(i,j,0)) F(i,j,0)=U(i,j,2) F(i,j,1)=U(i,j,0)*uu*vv F(i,j,2)=U(i,j,0)*vv*vv+p F(i,j,3)=(U(i,j,3)+p)*vv 50 continue end !------------------------------------------------------- ! 二维 差分格式x方向的差分算子 !入口: U，上一时刻U矢量，Uf、Ef，临时变量 ! dx，x方向的网格宽度，dt， 时间步长； !出口: U， 计算得到的当前时刻U矢量。 !------------------------------------------------------- subroutine LLx(U,Uf,Ef,dx,dt) implicit double precision (a-h,o-z) common /G_def/ GAMA,PI,Jx,Jy,dLx,dLy,TT,Sf dimension U(0:Jx+1,0:Jy+1,0:3),Uf(0:Jx+1,0:Jy+1,0:3) dimension Ef(0:Jx+1,0:Jy+1,0:3) r=dt/dx dnu=3.0*r*(1-3.0*r) do 60 i=1,Jx do 60 j=0,Jy+1 do 60 k=0,3 !开关函数 q=dabs(dabs(U(i+1,j,0)-U(i,j,0))-dabs(U(i,j,0)-U(i-1,j,0))) $ /(dabs(U(i+1,j,0)-U(i,j,0))+dabs(U(i,j,0)-U(i-1,j,0))+1e-10) !人工黏性项 Ef(i,j,k)=U(i,j,k)+0.5*dnu*q*(U(i+1,j,k)-2*U(i,j,k)+U(i-1,j,k)) 60 continue do 61 j=0,Jy+1 do 61 k=0,3 do 61 i=1,jx U(i,j,k)=Ef(i,j,k) 61 continue call U2E(U,Ef,0,Jx+1,0,Jy+1) do 63 i=0,Jx do 63 j=0,Jy+1 do 63 k=0,3 !U(n+1/2)(i+1/2)(j) Uf(i,j,k)=0.5*(U(i+1,j,k)+U(i,j,k)) $ -0.5*r*(Ef(i+1,j,k)-Ef(i,j,k)) 63 continue !E(n+1/2)(i+1/2)(j) call U2E(Uf,Ef,0,Jx,0,Jy+1) do 64 i=1,Jx do 64 j=0,Jy+1 do 64 k=0,3 !U(n+1)(i)(j) U(i,j,k)=U(i,j,k)-r*(Ef(i,j,k)-Ef(i-1,j,k)) 64 continue end !------------------------------------------------------- ! 二维 差分格式y方向的差分算子 !入口: U，上一时刻U矢量，Uf、Ff，临时变量， ! dy， y方向的网格宽度，dt， 时间步长； !出口: U，计算得到的当前时刻U矢量。 !------------------------------------------------------- subroutine LLy(U,Uf,Ff,dy,dt) implicit double precision (a-h,o-z) common /G_def/ GAMA,PI,Jx,Jy,dLx,dLy,TT,Sf dimension U(0:Jx+1,0:Jy+1,0:3),Uf(0:Jx+1,0:Jy+1,0:3) dimension Ff(0:Jx+1,0:Jy+1,0:3) r=dt/dy dnu=3.0*r*(1-3.0*r) do 70 i=0,Jx+1 do 70 j=1,Jy do 70 k=0,3 q=dabs(dabs(U(i,j+1,0)-U(i,j,0))-dabs(U(i,j,0)-U(i,j-1,0))) $ /(dabs(U(i,j+1,0)-U(i,j,0))+dabs(U(i,j,0)-U(i,j-1,0))+1e-10) !人工黏性项 Ff(i,j,k)=U(i,j,k)+0.5*dnu*q*(U(i,j+1,k)-2*U(i,j,k)+U(i,j-1,k)) 70 continue do 71 i=0,Jx+1 do 71 k=0,3 do 71 j=1,Jy U(i,j,k)=Ff(i,j,k) 71 continue call U2F(U,Ff,0,Jx+1,0,Jy+1) do 73 i=0,Jx+1 do 73 j=0,Jy do 73 k=0,3 !U(n+1/2)(i)(j+1/2) Uf(i,j,k)=0.5*(U(i,j+1,k)+U(i,j,k)) $ -0.5*r*(Ff(i,j+1,k)-Ff(i,j,k)) 73 continue !F(n+1/2)(i)(j+1/2) call U2F(Uf,Ff,0,Jx+1,0,Jy) do 74 i=0,Jx+1 do 74 j=1,Jy do 74 k=0,3 !U(n+1)(i)(j) U(i,j,k)=U(i,j,k)-r*(Ff(i,j,k)-Ff(i,j-1,k)) 74 continue end !-------------------------------------------