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

! 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=0
1 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*u1 1 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*u2 2 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=vel
10 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)*uu
40 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)
uu=U(i,1)/rou
p=(GAMMA-1)*(U(i,2)-0.5*U(i,0)*uu*uu)
write(1,51)i*dx,rou,uu,p,U(i,2)
50 continue
close(1)
51 format(D20.10,D20.10,D20.10,D20.10,D20.10)
end。