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

jujuju987 · 发表于 2016-5-24 20:30:09

[Fortran] 纯文本查看 复制代码

001Program TwoStep_Lax
002Implicit None
003 
004Real(8),Parameter:: lendom=10.0,Tini=300.0
005Real(8),Parameter:: PiniL=10.0,PiniR=1.0
006Real(8),Parameter:: uini=0.0,gamma=1.4,R=287.0
007Real(8),Parameter:: CFL=0.98,timeOut=0.006
008Integer,Parameter:: nnod=101
009Integer:: i,j,n
010Real(8):: dx,DiniL,DiniR,Hini,aini,dt,Eini,time,comp1,comp2
011 
012Real(8):: Den(nnod),Tem(nnod),Pre(nnod),Vel(nnod)
013Real(8):: M(nnod),Ene(nnod),Ent(nnod),a(nnod)
014Real(8):: U(3,nnod),F(3,nnod),Un(3,nnod)
015 
016 
017Open(Unit=10,File="TwostepLaxUn1.TXT")
018Open(Unit=11,File="TwostepLaxF11.TXT")
019Open(Unit=12,File="TwostepLaxF12.TXT")
020Open(Unit=13,File="TwostepLaxF13.TXT")
021Open(Unit=14,File="TwostepLaxF21.TXT")
022Open(Unit=15,File="TwostepLaxF22.TXT")
023Open(Unit=16,File="TwostepLaxF23.TXT")
024Open(Unit=17,File="TwostepLaxUn2.TXT")
025Open(Unit=19,File="TwostepLaxDen.TXT")
026Open(Unit=20,File="TwostepLaxVel.TXT")
027Open(Unit=21,File="TwostepLaxPre.TXT")
028Open(Unit=22,File="TwostepLaxM.TXT")
029 
030 
031dx=lendom/(nnod-1)
032DiniL=100000.0*PiniL/(R*Tini)
033DiniR=100000.0*PiniR/(R*Tini)
034Eini=R*Tini/(gamma-1.0)
035Hini=Eini+R*Tini
036aini=(gamma*R*Tini)**0.5
037dt=CFL*dx/aini
038 
039 
040Den(1:nnod)=0.0
041Tem(1:nnod)=0.0
042Pre(1:nnod)=0.0
043Vel(1:nnod)=0.0
044M(1:nnod)=0.0
045Ene(1:nnod)=0.0
046Ent(1:nnod)=0.0
047a(1:nnod)=0.0
048 
049 
050U(1,1:50)=DiniL
051U(1,51:nnod)=DiniR
052U(2,1:nnod)=0.0
053U(3,1:50)=DiniL*Eini
054U(3,51:nnod)=DiniR*Eini
055 
056 
057F(1,1:nnod)=0.0
058F(2,1:50)=PiniL
059F(2,51:nnod)=PiniR
060F(3,1:nnod)=0.0
061 
062Un(1:3,2:100)=0
063Un(1:3,1)=U(1:3,1)
064Un(1:3,nnod)=U(1:3,nnod)
065 
066 
067time=0.0
068Do While (time<timeOut)
069  time=time+dt
070 
071  Do i=1,3
072    Do j=2,100
073      Un(i,j)=(U(i,j+1)+U(i,j))*0.5-dt/dx*(F(i,j+1)-F(i,j))*0.5
074      Write(10,"('Un(',I1,',',I3,')',E17.10)") i,j,Un(i,j)
075    End Do
076  End Do
077 
078  Do j=1,nnod
079    Den(j)=Un(1,j)
080    Vel(j)=Un(2,j)/Den(j)
081    Ene(j)=Un(3,j)/Den(j)
082    Tem(j)=(Ene(j)-0.5*Vel(j)**2)*(gamma-1.0)/R
083    Pre(j)=R*Den(j)*Tem(j)
084    Ent(j)=Ene(j)+R*Tem(j)
085    a(j)=(gamma*R*Tem(j))**0.5
086    M(j)=Vel(j)/a(j)
087 
088    F(1,j)=Den(j)*Vel(j)
089    F(2,j)=Den(j)*Vel(j)**2+Pre(j)
090    F(3,j)=Den(j)*Vel(j)*Ent(j)
091    Write(11,*) F(1,j)
092    Write(12,*) F(2,j)
093    Write(13,*) F(3,j)
094  End Do
095 
096  Do i=1,3
097    Do j=2,100
098      Un(i,j)=U(i,j)-dt/dx*(F(i,j)-F(i,j-1))
099      Write(17,"('Un(',I1,',',I3,')',E17.10)") i,j,Un(i,j)
100    End Do
101  End Do
102 
103  Do j=1,nnod
104    Den(j)=Un(1,j)
105    Vel(j)=Un(2,j)/Den(j)
106    Ene(j)=Un(3,j)/Den(j)
107    Tem(j)=(Ene(j)-0.5*Vel(j)**2)*(gamma-1.0)/R
108    Pre(j)=R*Den(j)*Tem(j)
109    Ent(j)=Ene(j)+R*Tem(j)
110    a(j)=(gamma*R*Tem(j))**0.5
111    M(j)=Vel(j)/a(j)
112 
113    F(1,j)=Den(j)*Vel(j)
114    F(2,j)=Den(j)*Vel(j)**2+Pre(j)
115    F(3,j)=Den(j)*Vel(j)*Ent(j)
116    Write(14,*) F(1,j)
117    Write(15,*) F(2,j)
118    Write(16,*) F(3,j)
119  End Do
120 
121  comp1=0.0
122  Do j=1,nnod
123    comp2=a(j)+Vel(j)
124    If (comp2>comp1) comp1=comp2
125  End Do
126  dt=CFL*dx/comp1
127 
128  U(1:3,1:nnod)=Un(1:3,1:nnod)
129 
130  n=n+1
131 
132End Do
133 
134Do j=1,nnod
135  Write(19,*)Den(j)
136  Write(20,*)Vel(j)
137  Write(21,*)Pre(j)
138  Write(22,*)M(j)
139End Do