计算流体力学中两步Lax程序，可以运行，数据错误，求解 - Fortran Coder

[Fortran] 纯文本查看 复制代码
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