有人用FORTRAN语言计算过喷管的准一维流动？

骑士精魂

我想运行下网上下载的程序
发现没有结果啊！
希望高人指点啊
谢谢

sw_vehicle

俺也很想了解一下这方面的东西
不知道有没有高手指点哈

everestar

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!     拟一维喷管流动数值模拟 2010-06-29 李永刚          !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!定义各个变量，数组，文件

implicit real*8(a-h,q-z)
real(4) A(31)          !喷管离散面积数组，自然对数函数ALOG(x)自变量要求单精度实数
real(8) TL,x,dx,gama,dt,ddt,C,time,Rou(31),V(31),T(31),Rout(31),Vt(31),Tt(31)
real(8)        Roud(31),Vd(31),Td(31),Roudt(31),Vdt(31),Tdt(31),Ma(31),p(31)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!        TL:喷管总长；x:长度坐标；dx:长度方向增量，网格点间距                                                                        !
!        time:时间坐标；dt:时间增量；ddt:求dt用中间变量；nnend:总循环时间步数                                        !
!        gama:比热比；C：柯朗数；                                                                                                                                !
!        Rou(31),V(31),T(31),p(31),Ma(31):流体密度，速度，温度、压强、马赫数在t时刻的值                        !
!        Rout(31),Vt(31),Tt(31):流体密度，速度，温度的时间导数在t时刻的值，31个网格点                        !
!        Roud(31),Vd(31),Td(31):流体密度，速度，温度在t+dt时刻的值，31个网格点                                        !
!        Roudt(31),Vdt(31),Tdt(31):流体密度，速度，温度的时间导数在t+dt时刻的值，31个网格点                !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!                                                                                                                                                                                       
open(1,file='read.dat')  
read(1,*) TL,NodeN,gama,C,nnend

open(3,file='Rou.dat')         !存储密度
open(4,file='V.dat')         !存储速度
open(5,file='T.dat')         !存储温度
open(6,file='p.dat')         !存储压强
open(7,file='Ma.dat')         !存储马赫数
open(8,file='不同时间点的质量流量.dat')         !质量流量
open(9,file='定常解.dat')!存储定常值

!计算流场变量初始值，根据给定的函数
dx=TL/(NodeN-1)
x=0                                                                  
do i=1,NodeN
x=(i-1)*dx
A(i)=1+2.2*(x-1.5)**2
Rou(i)=1-0.3146*x
T(i)=1-0.2314*x
V(i)=(0.1+1.09*x)*sqrt(T(i))
end do

!无量纲化
do i=1,NodeN
A(i)=A(i)/A(16)
Rou(i)=Rou(i)/Rou(1)
T(i)=T(i)/T(1)
V(i)=V(i)/V(1)
x=0                                                                          
x=(i-1)*dx
end do

time=0
!输出初始值
write(3,200) time
write(4,200) time
write(5,200) time
write(6,200) time
write(7,200) time

do i=1,NodeN,1
if(i==NodeN)  then
        write(3,400) Rou(i)
        write(4,400) V(i)
        write(5,400) T(i)
        write(6,400) p(i)
        write(7,400) Ma(i)
else
        write(3,300) Rou(i)
        write(4,300) V(i)
        write(5,300) T(i)
        write(6,300) p(i)
        write(7,300) Ma(i)
end if
end do

!质量流量输出格式定义
write(8,430)
430 format(7X,'t',\)
do i=1,NodeN-1
write(8,440) i
440 format(I8,\)
end do
write(8,445) NodeN
445 format(I8)

!开始时间循环
do nn=1,nnend,1

!求t时刻中间网格点(2~30)流场变量的时间偏导数,用向前差分
do i=2,NodeN-1
Rout(i)=-V(i)*(Rou(i+1)-Rou(i))/dx&
                -Rou(i)*(V(i+1)-V(i))/dx&   
                -Rou(i)*V(i)*(ALOG(A(i+1))-ALOG(A(i)))/dx

Vt(i)=-V(i)*(V(i+1)-V(i))/dx&
                -(T(i+1)-T(i))/dx/gama&
                -T(i)*(Rou(i+1)-Rou(i))/Rou(i)/dx/gama

Tt(i)=-V(i)*(T(i+1)-T(i))/dx&
                -(gama-1)*T(i)*(V(i+1)-V(i))/dx&
                -(gama-1)*T(i)*V(i)*(ALOG(A(i+1))-ALOG(A(i)))/dx
end do

!求时间步长dt,根据中间网格点
ddt=0
dt=C*(dx/(V(2)+sqrt(T(2))))
do i=3,NodeN-1
ddt=C*(dx/(V(i)+sqrt(T(i))))
if (ddt<dt)        then
dt=ddt
endif
end do

!计算中间网格点(包括1点，因为要用1点处的预估值计算t+dt时刻偏导数)流场变量t+dt时刻预估值
do i=1,NodeN-1
Roud(i)=Rou(i)+Rout(i)*dt
Vd(i)=V(i)+Vt(i)*dt
Td(i)=T(i)+Tt(i)*dt
end do
Roud(1)=Rou(1)
!Vd(1)=V(1)
Td(1)=T(1)
!网格点1的预估值如何来取是个小问题，这里认为Rou(1)、T(1)为不变量，永远都是初始值，
!V(1)为变量,需计算其预估值

!计算t+dt时刻中间网格点流场变量的时间导数，用预估值向后差分表示空间导数
do i=2,NodeN-1
Roudt(i)=-Vd(i)*(Roud(i)-Roud(i-1))/dx &
                -Roud(i)*(Vd(i)-Vd(i-1))/dx   &
                -Roud(i)*Vd(i)*(ALOG(A(i))-ALOG(A(i-1)))/dx

Vdt(i)=-Vd(i)*(Vd(i)-Vd(i-1))/dx &
                -(Td(i)-Td(i-1))/dx/gama   &
                -Td(i)*(Roud(i)-Roud(i-1))/Roud(i)/dx/gama

Tdt(i)=-Vd(i)*(Td(i)-Td(i-1))/dx &
                -(gama-1)*Td(i)*(Vd(i)-Vd(i-1))/dx   &
                -(gama-1)*Td(i)*Vd(i)*(ALOG(A(i))-ALOG(A(i-1)))/dx
end do

!完成校正步计算
do i=2,NodeN-1

!计算平均时间导数
Rout(i)=0.5*(Rout(i)+Roudt(i))
Vt(i)=0.5*(Vt(i)+Vdt(i))
Tt(i)=0.5*(Tt(i)+Tdt(i))

!计算t+dt时刻流场变量校正值
Rou(i)=Rou(i)+Rout(i)*dt
V(i)=V(i)+Vt(i)*dt
T(i)=T(i)+Tt(i)*dt
end do

!计算边界网格点在t+dt时刻的流场变量，线性外插
V(1)=2*V(2)-V(3)
V(NodeN)=2*V(NodeN-1)-V(NodeN-2)
Rou(NodeN)=2*Rou(NodeN-1)-Rou(NodeN-2)
T(NodeN)=2*T(NodeN-1)-T(NodeN-2)

!计算压强和马赫数
do i=1,NodeN
p(i)=Rou(i)*T(i)
Ma(i)=V(i)/sqrt(T(i))
end do

time=time+dt
!输出结果
write(3,200) time
write(4,200) time
write(5,200) time
write(6,200) time
write(7,200) time
200 format(F8.4,\)

do i=1,NodeN,1
if(i==NodeN)  then
        write(3,400) Rou(i)
        write(4,400) V(i)
        write(5,400) T(i)
        write(6,400) p(i)
        write(7,400) Ma(i)
else
        write(3,300) Rou(i)
        write(4,300) V(i)
        write(5,300) T(i)
        write(6,300) p(i)
        write(7,300) Ma(i)
end if
300 format(F8.3,\)
400 format(F8.3)
end do

!在屏幕显示时间进程
print*,nn,time

!输出特定时间点的质量流量
if ((nn.eq.50).or.(nn.eq.100).or.(nn.eq.150).or.(nn.eq.200).or.(nn==1000)) then
write(8,449) nn
449 format(I8,\)
do i=1,NodeN-1
write(8,450) Rou(i)*V(i)*A(i)
450 format(F8.3,\)
end do
write(8,460) Rou(NodeN)*V(NodeN)*A(NodeN)
460 format(F8.3)
end if

end do        !完成一时间步循环

!输出稳态值
do i=1,NodeN
write(9,500) i,Rou(i),V(i),T(i),p(i),Ma(i),Rou(i)*V(i)*A(i)
500 format(I4,6F8.3)
end do

end

fifor

anderson的计算流体力学里的程序吗？

Heidiadalheid

fifor 发表于 2010-10-10 10:21

                               
登录/注册后可看大图

anderson的计算流体力学里的程序吗？

从代码上看好像是板凳自己写的程序

wlzs10

不懂