二阶迎风型差分格式求解一维可压缩黏性流动问题Fortran语言

Fortran77语言源程序
! UpwindTVD_1D.f
!----------------------------------------------------------------------- ! 二阶迎风型TVD差分格式求解一维可压缩黏性流动问题
!（Fortran77语言版本）
!----------------------------------------------------------------------- program UPWIND_TVD_1D
implicit real*8 (a-h,o-z)
parameter(mx=201,Tt=5.0 )
dimension Q(3,mx),Qold(3,mx) ! Q: [rou, rou*u, E]
dimension rou(mx),u(mx),p(mx),T(mx),E(mx)
real*8 Ma
common/para_def/Sf,Re,Ma,pr,gama,cp,cv,dx,dt,im
im=mx
call Initialize(Q)
time=0.0
n=0
1 n=n+1
time=time+dt
do i=1,im
do l=1,3
Qold(l,i)=Q(l,i)
end do
end do
call UpwindTVD_1D_Solver(Q)
if(n/10000*10000.eq.n)then
write(*,20)n,time,error(Q,Qold)
call OutPut(Q)
end if
20 format(1x,i10,' time=',e16.9,' error=',e16.9,1x)
if((time.lt.Tt) .and. (error(Q,Qold).gt.1.0e-8))then
goto 1
end if
call OutPut(Q)
write(*,*) 'Program end !'
end
!--------------------------------------------------------------------
subroutine Initialize(Q)
implicit real*8 (a-h,o-z)
real*8 Ma
common/para_def/Sf,Re,Ma,pr,gama,cp,cv,dx,dt,im
dimension rou(im),u(im),p(im),T(im),E(im),Q(3,im)
Re=50.0 ! 雷诺数
Ma=2.0 ! 马赫数
pr=0.72 ! 普朗特数
gama=1.4 ! 气体常数
cv=1.0/(gama*(gama-1.0)*Ma**2) ! 比定容热容
cp=gama*cv ! 比定压热容
dt=1.0e-6
dx=1.0/(1.0*(im-1)) ! 空间步长
xl=0.0
xr=1.0
rou0=1.0
u0=1.0
T0=1.0
temp=(gama+1.)/(gama-1.)
u1=(2./(gama-1.)+Ma**2)/(temp*Ma**2)
T1=(2.*gama/(gama+1)*Ma**2-1./temp)*(1./temp+2./(gama+1.)/Ma**2) do i=1,im
rou(i)=rou0
xi=(i-1)*dx
u(i)=(xi-xl)/(xr-xl)*(u1-u0)+u0
T(i)=(xi-xl)/(xr-xl)*(T1-T0)+T0
P(i)=rou(i)*T(i)/(gama*Ma**2)
E(i)=rou(i)*(cv*T(i)+0.5*u(i)**2)
end do
do i=1,im
Q(1,i)=rou(i)
Q(2,i)=rou(i)*u(i)
Q(3,i)=E(i)
end do
end
!-------------------------------------------------------------------- subroutine Q2U(Q,rou,u,p,T,E,a)
implicit real*8 (a-h,o-z)
real*8 Ma
common/para_def/Sf,Re,Ma,pr,gama,cp,cv,dx,dt,im
dimension Q(3,im),rou(im),u(im)
dimension p(im),T(im),E(im),a(im)
do i=1,im
rou(i)=Q(1,i)
u(i)=Q(2,i)/Q(1,i)
E(i)=Q(3,i)
T(i)=(E(i)/rou(i)-0.5*u(i)**2)/cv
p(i)=rou(i)*T(i)/(gama*Ma**2)
a(i)=dsqrt(T(i))/Ma
end do
end
!-------------------------------------------------------------------- subroutine UpwindTVD_1D_Solver(Q)
implicit real*8 (a-h,o-z)
real*8 Ma,lamda,minmod
common/para_def/Sf,Re,Ma,pr,gama,cp,cv,dx,dt,im
dimension Q(3,im),Qold(3,im),dq(3,im),rou(im),u(im),p(im),T(im) dimension E(im),a(im),f(3,im),fwave(3,im),fr(3),fl(3)
dimension Gv(3,im),g(3,im),lamda(3,im),sigma(3,im)
dimension alfa(3,im),Rn(3,3,im),R(3,3,im),beta(3,im)
call Q2U(Q,rou,u,p,T,E,a)
do i=1,im
f(1,i)=rou(i)*u(i)
f(2,i)=rou(i)*u(i)**2+p(i)
f(3,i)=u(i)*(E(i)+p(i))
end do
do i=1,im-1
lamda(1,i)=(u(i)+u(i+1))/2.0
lamda(2,i)=(u(i)-a(i)+u(i+1)-a(i+1))/2.0
lamda(3,i)=(u(i)+a(i)+u(i+1)+a(i+1))/2.0
end do
do i=1,im-1
do l=1,3
dq(l,i)=(Q(l,i+1)-Q(l,i))
sigma(l,i)=qz(lamda(l,i))/2.0 ! 定常情况
end do
end do
do i=1,im-1
ut=(u(i)+u(i+1))/2.0 !ut 临时变量
at=(a(i)+a(i+1))/2.0 !at 临时变量
R(1,1,i)=1.0
R(1,2,i)=1.0
R(1,3,i)=1.0
R(2,1,i)=ut
R(2,2,i)=ut-at
R(2,3,i)=ut+at
R(3,1,i)=ut*ut/2.0
R(3,2,i)=at*at/(gama-1.0)+ut*ut/2.0-ut*at
R(3,3,i)=at*at/(gama-1.0)+ut*ut/2.0+ut*at
Rn(1,1,i)=1.0-(gama-1.0)*ut*ut/2.0/at/at
Rn(1,2,i)=(gama-1.0)*ut/at/at
Rn(1,3,i)=-(gama-1.0)/at/at
Rn(2,1,i)=ut/2.0/at+(gama-1.0)*ut*ut/4.0/at/at
Rn(2,2,i)=-(gama-1.0)*ut/2./at/at-1.0/2.0/at
Rn(2,3,i)=(gama-1.0)/2.0/at/at
Rn(3,1,i)=-ut/2.0/at+(gama-1.0)*ut*ut/4.0/at/at
Rn(3,2,i)=(1.0-gama)*ut/2./at/at+1.0/2.0/at
Rn(3,3,i)=(gama-1.0)/2./at/at
end do
do i=1,im
do l=1,3
alfa(l,i)=0.0
fwave(l,i)=0.0
g(l,i)=0.0
end do
end do
do i=1,im-1
do l=1,3
do k=1,3
alfa(l,i)=alfa(l,i)+Rn(l,k,i)*dq(k,i)
end do
end do
end do
do i=2,im-1
do l=1,3
g(l,i)=minmod(sigma(l,i)*alfa(l,i),sigma(l,i-1)*alfa(l,i-1)) end do
end do
do i=1,im-1
do l=1,3
if(alfa(l,i).ne.0.0)then
beta(l,i)=(g(l,i+1)-g(l,i))/alfa(l,i)
else
beta(l,i)=0.0
end if
end do
end do
do i=1,im-1
do l=1,3
do k=1,3
fwave(l,i)=fwave(l,i)+0.5*R(l,k,i)*(g(k,i)+g(k,i+1)
* -qz(lamda(k,i)+beta(k,i))*alfa(k,i))
end do
end do
do i=2,im-1
do l=1,3
fl(l)= 0.5*(f(l,i-1)+f(l,i)) +fwave(l,i-1)
fr(l)= 0.5*(f(l,i+1)+f(l,i)) +fwave(l,i)
Q(l,i)=Q(l,i)-dt/dx*(fr(l)-fl(l))
end do
end do
! 处理黏性项
do i=2,im-1
Gv(1,i)=0.0
Gv(2,i)=(T(i)*(u(i+1)-2.*u(i)+u(i-1))+(T(i+1)-T(i-1))*
& (u(i+1)-u(i-1))/4.)*4./3./re/dx**2
Gv(3,i)=((T(i)*(u(i+1)-2.*u(i)+u(i-1))*u(i)+(T(i+1)-T(i-1))*
& (u(i+1)-u(i-1))*u(i)/4.+T(i)*(u(i+1)-u(i-1))**2/4.)*4./3. & +(T(i)*(T(i+1)-2.*T(i)+T(i-1))+(T(i+1)-T(i-1))**2/4.)
& *cp/pr)/re/dx**2
end do
do i=2,im-1
do l=1,3
Q(l,i)=Q(l,i)+Gv(l,i)*dt
end do
end do
! 边界条件
Q(1,im)=Q(1,im-1)
Q(2,im)=Q(1,im-1)*u(im)
Q(3,im)=Q(1,im-1)*(cv*T(im)+0.5*u(im)**2)
end
!-------------------------------------------------------------------- function error(Q,Qold)
implicit real*8 (a-h,o-z)
real*8 Ma,err,error
common/para_def/Sf,Re,Ma,pr,gama,cp,cv,dx,dt,im
dimension Q(3,im),Qold(3,im)
err=1.0e-10
do i=1,im
do l=1,3
if(Qold(l,i).ne.0.0)then
ddq=abs((Q(l,i)-Qold(l,i))/Qold(l,i))
else
ddq=abs(Q(l,i)-Qold(l,i))
if(err.lt.ddq) err=ddq
end do
end do
error=err
end
!-------------------------------------------------------------------- function qz(x)
real*8 x,qz
if (abs(x).gt.0.1) then
qz=abs(x)
else
qz=(x**2/0.1+0.1)/2.
endif
end function qz
!---------------------------------------------------------------------------- function minmod(x,y)
real*8 x,y,minmod
if(x*y.le.0.0) then
minmod=0.0
else if(abs(x).gt.abs(y)) then
minmod=y
else
minmod=x
endif
end function minmod
!--------------------------------------------------------------------
subroutine OutPut(Q)
implicit real*8 (a-h,o-z)
real*8 Ma
common/para_def/Sf,Re,Ma,pr,gama,cp,cv,dx,dt,im
dimension Q(3,im)
open(1,file='result.dat')
write(1,*)'variables=x,rou,u,p,e'
do i=1,im
x=(i-1)*dx
rou=Q(1,i)
u=Q(2,i)/rou
E=Q(3,i)
T=(E/rou-0.5*u*u)/cv
p=rou*T/(gama*Ma**2)
write(1,10)x,rou,u,p,cv*T
end do
close(1)
10 format(1x,f5.3,1x,3(e15.8,1x))
end
------------------------------------------------------------------------ ----------------------------------------------------------------------- -----------------------------------------------------------------------。