[其他行业算法] 计算流体力学中两步Lax程序，可以运行，数据错误，求解

Program TwoStep_Lax
Implicit None

Real(8),Parameter:: lendom=10.0,Tini=300.0
Real(8),Parameter:: PiniL=10.0,PiniR=1.0
Real(8),Parameter:: uini=0.0,gamma=1.4,R=287.0
Real(8),Parameter:: CFL=0.98,timeOut=0.006
Integer,Parameter:: nnod=101
Integer:: i,j,n
Real(8):: dx,DiniL,DiniR,Hini,aini,dt,Eini,time,comp1,comp2

Real(8):: Den(nnod),Tem(nnod),Pre(nnod),Vel(nnod)
Real(8):: M(nnod),Ene(nnod),Ent(nnod),a(nnod)
Real(8):: U(3,nnod),F(3,nnod),Un(3,nnod)


Open(Unit=10,File="TwostepLaxUn1.TXT")
Open(Unit=11,File="TwostepLaxF11.TXT")
Open(Unit=12,File="TwostepLaxF12.TXT")
Open(Unit=13,File="TwostepLaxF13.TXT")
Open(Unit=14,File="TwostepLaxF21.TXT")
Open(Unit=15,File="TwostepLaxF22.TXT")
Open(Unit=16,File="TwostepLaxF23.TXT")
Open(Unit=17,File="TwostepLaxUn2.TXT")
Open(Unit=19,File="TwostepLaxDen.TXT")
Open(Unit=20,File="TwostepLaxVel.TXT")
Open(Unit=21,File="TwostepLaxPre.TXT")
Open(Unit=22,File="TwostepLaxM.TXT")


dx=lendom/(nnod-1)
DiniL=100000.0*PiniL/(R*Tini)
DiniR=100000.0*PiniR/(R*Tini)
Eini=R*Tini/(gamma-1.0)
Hini=Eini+R*Tini
aini=(gamma*R*Tini)**0.5
dt=CFL*dx/aini


Den(1:nnod)=0.0
Tem(1:nnod)=0.0
Pre(1:nnod)=0.0
Vel(1:nnod)=0.0
M(1:nnod)=0.0
Ene(1:nnod)=0.0
Ent(1:nnod)=0.0
a(1:nnod)=0.0


U(1,1:50)=DiniL
U(1,51:nnod)=DiniR
U(2,1:nnod)=0.0
U(3,1:50)=DiniL*Eini
U(3,51:nnod)=DiniR*Eini


F(1,1:nnod)=0.0
F(2,1:50)=PiniL
F(2,51:nnod)=PiniR
F(3,1:nnod)=0.0

Un(1:3,2:100)=0
Un(1:3,1)=U(1:3,1)
Un(1:3,nnod)=U(1:3,nnod)


time=0.0
Do While (time<timeOut)
  time=time+dt

  Do i=1,3
    Do j=2,100
      Un(i,j)=(U(i,j+1)+U(i,j))*0.5-dt/dx*(F(i,j+1)-F(i,j))*0.5
      Write(10,"('Un(',I1,',',I3,')',E17.10)") i,j,Un(i,j)
    End Do
  End Do

  Do j=1,nnod
    Den(j)=Un(1,j)
    Vel(j)=Un(2,j)/Den(j)
    Ene(j)=Un(3,j)/Den(j)
    Tem(j)=(Ene(j)-0.5*Vel(j)**2)*(gamma-1.0)/R
    Pre(j)=R*Den(j)*Tem(j)
    Ent(j)=Ene(j)+R*Tem(j)
    a(j)=(gamma*R*Tem(j))**0.5
    M(j)=Vel(j)/a(j) 

    F(1,j)=Den(j)*Vel(j)
    F(2,j)=Den(j)*Vel(j)**2+Pre(j)
    F(3,j)=Den(j)*Vel(j)*Ent(j)
    Write(11,*) F(1,j)
    Write(12,*) F(2,j)
    Write(13,*) F(3,j)
  End Do

  Do i=1,3
    Do j=2,100
      Un(i,j)=U(i,j)-dt/dx*(F(i,j)-F(i,j-1))
      Write(17,"('Un(',I1,',',I3,')',E17.10)") i,j,Un(i,j)
    End Do
  End Do

  Do j=1,nnod
    Den(j)=Un(1,j)
    Vel(j)=Un(2,j)/Den(j)
    Ene(j)=Un(3,j)/Den(j)
    Tem(j)=(Ene(j)-0.5*Vel(j)**2)*(gamma-1.0)/R
    Pre(j)=R*Den(j)*Tem(j)
    Ent(j)=Ene(j)+R*Tem(j)
    a(j)=(gamma*R*Tem(j))**0.5
    M(j)=Vel(j)/a(j)

    F(1,j)=Den(j)*Vel(j)
    F(2,j)=Den(j)*Vel(j)**2+Pre(j)
    F(3,j)=Den(j)*Vel(j)*Ent(j)
    Write(14,*) F(1,j)
    Write(15,*) F(2,j)
    Write(16,*) F(3,j)
  End Do

  comp1=0.0
  Do j=1,nnod
    comp2=a(j)+Vel(j)
    If (comp2>comp1) comp1=comp2
  End Do
  dt=CFL*dx/comp1

  U(1:3,1:nnod)=Un(1:3,1:nnod)

  n=n+1 

End Do

Do j=1,nnod
  Write(19,*)Den(j)
  Write(20,*)Vel(j)
  Write(21,*)Pre(j)
  Write(22,*)M(j)
End Do