喷管内流场计算程序
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
喷管内流场计算程序
!本程序采用三种格式对Buckley-Leverett方程进行求解
!计算过程中采用人工粘性进行处理
!name,name1是用于进行变文件名输出数据的字符串参数
!n,m分别表示空间网格节点和选择哪种计算方法
!uN,SN分别表示前一时刻的速度、人工粘性值
!u,FN分别表示这一时刻的速度,前一时刻对流项的函数值
!dx,dt,time分别表示空间尺度、时间尺度和总计算时间
!Cx分别表示人工粘性系数
program main
implicit none
character(len=15) :: name,name1
integer :: i,n=201,m
real(kind=8) :: uN(201),SN(201),FN(201)
real(kind=8) :: u(201),AN(201)
real(kind=8) :: dx,dt,time,t,Cx
!给定输入参数，对于不同的边界条件需要修改
dx=2.0/(n-1)
t=0.0
time=0.4
dt=0.0001
Cx=0.006
m=2
!给定初始时刻给定的速度值，不同边界条件时需要修改
do i=1,n
if(-1.0+(i-1)*dx<=0.0.and.-1.0+(i-1)*dx>=-0.5)then uN(i)=1.0
else
uN(i)=0.0
end if
end do
!选择方法进行计算
if(m==1) then
name1="Lax_Friedrichs"
do while(t<=time)
t=t+dt
do i=1,n
FN(i)=4.0*uN(i)**2/(4*uN(i)**2+(1-uN(i))**2) end do
do i=2,n-1
SN(i)=Cx*(uN(i+1)-2*uN(i)+uN(i-1))
u(i)=(uN(i+1)+uN(i-1))/2.0-dt*(FN(i+1)-FN(i-1))/(2.0*dx)+SN(i) end do
u(1)=uN(1)
u(n)=uN(n)
do i=1,n
uN(i)=u(i)
end do
end do
else if(m==2) then
name1="Lax_Wendroff"
do while(t<=time)
t=t+dt
do i=1,n
FN(i)=4.0*uN(i)**2/(4*uN(i)**2+(1-uN(i))**2)
end do
do i=1,n-1
if(uN(i)==uN(i+1)) then
AN(i)=1.0
else
AN(i)=(FN(i+1)-FN(i))/(uN(i+1)-uN(i))
end if
end do
do i=2,n-1
SN(i)=Cx*(uN(i+1)-2*uN(i)+uN(i-1))
u(i)=uN(i)-dt*(FN(i+1)-FN(i-1))/(dx*2.0)+(AN(i)*dt/dx)**2*(uN(i+1)-& uN(i))/2.0-(AN(i-1)*dt/dx)**2*(uN(i)-uN(i-1))/2.0+SN(i)
end do
u(1)=uN(1)
u(n)=uN(n)
do i=1,n
uN(i)=u(i)
end do
end do
else
name1="upwind"
do while(t<=time)
t=t+dt
do i=1,n
FN(i)=4.0*uN(i)**2/(4*uN(i)**2+(1-uN(i))**2)
end do
do i=1,n-1
if(uN(i)==uN(i+1)) then
AN(i)=1.0
else
AN(i)=(FN(i+1)-FN(i))/(uN(i+1)-uN(i))
end if
end do
do i=2,n-1
SN(i)=Cx*(uN(i+1)-2*uN(i)+uN(i-1))
u(i)=uN(i)-dt*(FN(i+1)-FN(i-1))/(dx*2.0)+abs(AN(i)*dt/dx)*(uN(i+1)-& uN(i))/2.0-abs(AN(i-1)*dt/dx)*(uN(i)-uN(i-1))/2.0+SN(i)
end do
u(1)=uN(1)
u(n)=uN(n)
do i=1,n
uN(i)=u(i)
end do
end do
end if
write(name,"(A15)") name1
open(10,file=''//trim(adjustl(name))//'.dat')
do i=1,n
!输出速度值，对于不同边界条件需要修改
write(10,"(f9.5,2x,f9.5)") -1.0+(i-1)*dx,u(i)
end do
stop
end