利用Lax-wendroff两步差分格式求解一维激波管问题Fortran语言(沉鱼醉江湖制作).doc

5! Lax-wendroff1D.for !-!利用Lax-wendroff两步差分格式求解一维激波管问题! !- program Laxwendroff1D implicit double precision (a-h,o-z) parameter (M=1000) common /G_def/ GAMMA,PI,J,JJ,dL,TT,Sf dimension U(0:M+1,0:2),Uf(0:M+1,0:2) dimension Ef(0:M+1,0:2) GAMMA=1.4 !气体常数 PI=3.1415926 J=M !网格数 dL=2.0 !计算区域 TT=0.3 !总时间 Sf=0.8 !时间步长因子 call Init(U,dx) T=01 dt=CFL(U,dx) T=T+dt write(*,*)T=,T,dt=,dt call Laxwendroff_1D_Solver(U,Uf,Ef,dx,dt)call bound(U,dx) if(T.lt.TT)goto 1 call Output(U,dx) end!-!初始化!入口: 无；!出口: U, 已经给定的初始值，dx，网格宽度。!- subroutine Init(U,dx) implicit double precision (a-h,o-z) common /G_def/ GAMMA,PI,J,JJ,dL,TT,Sf dimension U(0:J+1,0:2)!-!左侧初始条件：1=1， u1=0， p1=1!右侧初始条件：2=0.125，u1=0， p1=0.1!- rou1=2.0 u1=0 p1=1.0 rou2=0.1 u2=0 p2=0.1 dx=dL/J !网格宽度!-! 薄膜放在激波管正中位置!- do 1 i=0,J/2 U(i,0)=rou1 U(i,1)=rou1*u1 U(i,2)=p1/(GAMMA-1)+0.5*rou1*u1*u11 continue do 2 i=J/2+1,J+1 U(i,0)=rou2 U(i,1)=rou2*u2 U(i,2)=p2/(GAMMA-1)+0.5*rou2*u2*u22 continue end!-!计算时间步长!入口: U， 当前物理量，dx, 网格宽度；!返回: 时间步长。!- double precision function CFL(U,dx) implicit double precision (a-h,o-z) common /G_def/ GAMMA,PI,J,JJ,dL,TT,Sf dimension U(0:J+1,0:2) dmaxvel=1e-10 do 10 i=1,J uu=U(i,1)/U(i,0) p=(GAMMA-1)*(U(i,2)-0.5*U(i,0)*uu*uu) vel=dsqrt(GAMMA*p/U(i,0)+dabs(uu) if(vel.gt.dmaxvel)dmaxvel=vel10 continue CFL=Sf*dx/dmaxvel end !-!一维 差分格式求解器!入口: U， 上一时刻U矢量，! Uf、Ef，临时变量，! dx，网格宽度，dt,，时间步长；!出口: U， 计算得到得当前时刻U矢量。!- subroutine Laxwendroff_1D_Solver(U,Uf,Ef,dx,dt) implicit double precision (a-h,o-z) common /G_def/ GAMMA,PI,J,JJ,dL,TT,Sf dimension U(0:J+1,0:2),Uf(0:J+1,0:2) dimension Ef(0:J+1,0:2) r=dt/dx dnu=0.25 do 20 i=1,J do 20 k=0,2!开关函数 q=dabs(dabs(U(i+1,0)-U(i,0)-dabs(U(i,0)-U(i-1,0)/(dabs(U(i+1,0)-U(i,0)+dabs(U(i,0)-U(i-1,0)+1e-10)!人工黏性项 Ef(i,k)=U(i,k)+0.5*dnu*q*(U(i+1,k)-2*U(i,k)+U(i-1,k) 20 continue do 21 k=0,2 do 21 i=1,J U(i,k)=Ef(i,k)21 continue call U2E(U,Ef,0,J+1) do 22 i=0,J do 22 k=0,2!U(n+1/2)(i+1/2) Uf(i,k)=0.5*(U(i+1,k)+U(i,k)-0.5*r*(Ef(i+1,k)-Ef(i,k) 22 continue !E(n+1/2)(i+1/2) call U2E(Uf,Ef,0,J) do 23 i=1,J do 23 k=0,2!U(n+1)(i) U(i,k)=U(i,k)-r*(Ef(i,k)-Ef(i-1,k) 23 continue end!-!边界条件!流动出口边界条件:NI,J=NI-1,J!入口: dx，网格宽度；!出口: U， 已经给定边界。!- subroutine bound(U,dx) implicit double precision (a-h,o-z) common /G_def/ GAMMA,PI,J,JJ,dL,TT,Sf dimension U(0:J+1,0:2)!左边界 do 30 k=0,2 U(0,k)=U(1,k)30 continue!右边界 do 31 k=0,2 U(J+1,k)=U(J,k)31 continue end!-!根据U计算E!入口: U，当前U矢量；!出口: E，计算得到的E矢量，! U、E定义见Euler方程组。!- subroutine U2E(U,E,is,in) implicit double precision (a-h,o-z) common /G_def/ GAMMA,PI,J,JJ,dL,TT,Sf dimension U(0:J+1,0:2),E(0:J+1,0:2) do 40 i=is,in uu=U(i,1)/U(i,0) p=(GAMMA-1)*(U(i,2)-0.5*U(i,1)*U(i,1)/U(i,0) E(i,0)=U(i,1) E(i,1)=U(i,0)*uu*uu+p E(i,2)=(U(i,2)+p)*uu40 continue end!-!输出结果, 用*.txt数据格式文件，再用Origin8.0画图!入口: U， 当前时刻U矢量，! dx，网格宽度；!出口: 无。!- subroutine Output(U,dx) implicit double precision (a-h,o-z) common /G_def/ GAMMA,PI,J,JJ,dL,TT,Sf dimension U(0:J+1,0:2) open(1,file=Lax-Wendroff_1D result.txt,status=unknown) do 50 i=0,J+1 rou=U(i,0