program Main
implicit none
integer nx, ny
integer it
real*8 cfl,p0, t0, rho0, ma0,gamma,resOutput
real*8, allocatable :: x(:,:), y(:,:), dxdi(:,:), dydi(:,:), lengthI(:,:), dxdj(:,:), dydj(:,:), lengthJ(:,:), area(:,:), dt(:,:),rho(:,:), u(:,:), v(:,:), p(:,:), q(:,:,:), residual(:,:,:)
nx = 401
ny = 101
cfl = 0.1d0
ma0 = 0.6d0
p0 = 101325.0d0
t0 = 288.15d0
rho0 = p0 / t0 / 287.0d0
gamma = 1.4d00
! 分配内存
allocate( x(nx,ny), y(nx,ny),dxdi(nx,ny-1), dydi(nx,ny-1), lengthI(nx,ny-1), dxdj(nx-1,ny), dydj(nx-1,ny), lengthJ(nx-1,ny), area(nx-1,ny-1), dt(nx-1,ny-1),rho(-1:nx+1,-1:ny+1), u(-1:nx+1,-1:ny+1), v(-1:nx+1,-1:ny+1), p(-1:nx+1,-1:ny+1), q(4,nx-1,ny-1),residual(4,nx-1,ny-1) )
! 划分网格
call Grid
! 初始化流场
call Ini
! 打开残差文件
open(unit = 1, file = 'residual.plt')
write(1,*) 'variables=it,residual'
it = 0
resOutput = 1.0d0
! 循环计算
do while(resOutput > 1.0d-4 .and. it<200000)
! 边界条件
call CalBC
! 计算残差
call Advance
! 输出残差
resOutput = sum(abs(residual(1,:,:))) / (nx-1) / (ny-1)
it = it + 1
write(* ,"(I6,1X,E11.4)") it, resOutput
write(1 ,"(I6,1X,E11.4)") it, resOutput
end do
! 关闭残差文件
close(1)
! 输出数据
call Output
pause
contains
subroutine Grid
implicit none
integer i, j, i1, j1
real*8 temp
do j = 1, ny
do i = 1, nx
x(i,j) = (i-1.0d0) / (nx-1.0d0) * 4.0d0
temp = 0.0d0
if(x(i,j)>1.0d0 .and. x(i,j)<2.0d0) then
temp = dsqrt(3.0d0**2.0d0 - (x(i,j)-1.5d0)**2.0d0) - dsqrt(3.0d0**2.0d0 - 0.5**2.0d0)
endif
y(i,j) = temp + (j-1.0d0) / (ny-1.0d0) * (1.0d0-temp)
end do
end do
do j = 1, ny-1
do i = 1, nx
dxdi(i,j) = y(i,j+1) - y(i,j )
dydi(i,j) = x(i,j ) - x(i,j+1)
lengthI(i,j) = dsqrt( dxdi(i,j)*dxdi(i,j) + dydi(i,j)*dydi(i,j) )
dxdi(i,j) = dxdi(i,j) / lengthI(i,j)
dydi(i,j) = dydi(i,j) / lengthI(i,j)
end do
end do
do j = 1, ny
do i = 1, nx-1
dxdj(i,j) = y(i ,j) - y(i+1,j)
dydj(i,j) = x(i+1,j) - x(i ,j)
lengthJ(i,j) = dsqrt( dxdj(i,j)*dxdj(i,j) + dydj(i,j)*dydj(i,j) )
dxdj(i,j) = dxdj(i,j) / lengthJ(i,j)
dydj(i,j) = dydj(i,j) / lengthJ(i,j)
end do
end do
do j =1, ny-1
do i = 1, nx-1
i1 = i + 1
j1 = j + 1
area(i,j) = 0.5d0*( dabs( ( x(i ,j ) - x(i1,j ) )*y(i1,j1)&
+( x(i1,j ) - x(i1,j1) )*y(i ,j )&
+( x(i1,j1) - x(i ,j ) )*y(i1,j ) )&
+ dabs( ( x(i ,j ) - x(i1,j1) )*y(i ,j1)&
+( x(i1,j1) - x(i ,j1) )*y(i ,j )&
+( x(i ,j1) - x(i ,j ) )*y(i1,j1) ) )
end do
end do
end
subroutine Ini
implicit none
integer i, j
do j = -1, ny+1
do i = -1, nx+1
rho(i,j) = rho0
u(i,j) = ma0*dsqrt(gamma*p0/rho0)
v(i,j) = 0.0d0
p(i,j) = p0
end do
end do
do j = 1, ny-1
do i = 1, nx-1
q(1,i,j) = rho(i,j)
q(2,i,j) = rho(i,j)*u(i,j)
q(3,i,j) = rho(i,j)*v(i,j)
q(4,i,j) = p(i,j)/(gamma - 1.0d0) + 0.5d0*rho(i,j)*(u(i,j)*u(i,j)+v(i,j)*v(i,j))
end do
end do
end
subroutine CalBC
implicit none
integer i, j
real*8 vNorm, vTemp
do j = -1, ny+1
rho(0,j) = 2.0d0*rho0 - rho(1,j)
u(0,j) = 2.0d0*ma0 * dsqrt(gamma*p0/rho0) - u(1,j)
v(0,j) = 2.0d0*0.d0 - v(1,j)
p(0,j) = 2.0d0*p0 - p(1,j)
rho(-1,j) = 2.d0*rho(0,j) - rho(1,j)
u(-1,j) = 2.d0*u(0,j) - u(1,j)
v(-1,j) = 2.d0*v(0,j) - v(1,j)
p(-1,j) = 2.d0*p(0,j) - p(1,j)
end do
do j = -1, ny+1
if(ma0 > 1.0) then
rho(nx,j) = rho(nx-1,j)
u(nx,j) = u(nx-1,j)
v(nx,j) = v(nx-1,j)
p(nx,j) = p(nx-1,j)
else
rho(nx,j) = rho(nx-1,j)
u(nx,j) = u(nx-1,j)
v(nx,j) = v(nx-1,j)
p(nx,j) = p0
endif
rho(nx+1,j) = 2.d0*rho(nx,j) - rho(nx-1,j)
u(nx+1,j) = 2.d0*u(nx,j) - u(nx-1,j)
v(nx+1,j) = 2.d0*v(nx,j) - v(nx-1,j)
p(nx+1,j) = 2.d0*p(nx,j) - p(nx-1,j)
end do
do i = 1, nx-1
vNorm = u(i,ny-1)*dxdj(i,ny) + v(i,ny-1)*dydj(i,ny)
vTemp = 2.d0 * vNorm
rho(i,ny) = rho(i,ny-1)
u(i,ny) = u(i,ny-1) - vTemp*dxdj(i,ny)
v(i,ny) = v(i,ny-1) - vTemp*dydj(i,ny)
p(i,ny) = p(i,ny-1)
rho(i,ny+1) = 2.d0*rho(i,ny) - rho(i,ny-1)
u(i,ny+1) = 2.d0*u(i,ny) - u(i,ny-1)
v(i,ny+1) = 2.d0*v(i,ny) - v(i,ny-1)
p(i,ny+1) = 2.d0*p(i,ny) - p(i,ny-1)
end do
do i = 1, nx-1
vNorm = u(i,1)*dxdj(i,1) + v(i,1)*dydj(i,1)
vTemp = 2.d0 * vNorm
rho(i,0) = rho(i,1)
u(i,0) = u(i,1) - vTemp*dxdj(i,1)
v(i,0) = v(i,1) - vTemp*dydj(i,1)
p(i,0) = p(i,1)
rho(i,-1) = 2.d0*rho(i,0) - rho(i,1)
u(i,-1) = 2.d0*u(i,0) - u(i,1)
v(i,-1) = 2.d0*v(i,0) - v(i,1)
p(i,-1) = 2.d0*p(i,0) - p(i,1)
end do
end
subroutine Advance
implicit none
real*8::rhoL, rhoR, uL, uR, vL, vR, pL, pR
real*8::flux(4)
real*8,allocatable::fluxI(:,:,:), fluxJ(:,:,:)
integer::i,j
real*8 resMax, uBar, vBar, c
allocate(fluxI(4,nx,ny-1), fluxJ(4,nx-1,ny))
do j = 1, ny-1
do i = 1, nx
call MUSCL(rho(i-2,j), rho(i-1,j), rho(i,j), rho(i+1,j), rhoL, rhoR)
call MUSCL(u(i-2,j), u(i-1,j), u(i,j), u(i+1,j), uL, uR)
call MUSCL(v(i-2,j), v(i-1,j), v(i,j), v(i+1,j), vL, vR)
call MUSCL(p(i-2,j), p(i-1,j), p(i,j), p(i+1,j), pL, pR)
call StegerWarming1D(rhoL, uL*dxdi(i,j)+vL*dydi(i,j), -uL*dydi(i,j)+vL*dxdi(i,j), pL, rhoR, uR*dxdi(i,j)+vR*dydi(i,j), -uR*dydi(i,j)+vR*dxdi(i,j), pR, flux)
fluxI(1,i,j) = flux(1)
fluxI(2,i,j) = flux(2)*dxdi(i,j) - flux(3)*dydi(i,j)
fluxI(3,i,j) = flux(2)*dydi(i,j) + flux(3)*dxdi(i,j)
fluxI(4,i,j) = flux(4)
end do
end do
do j = 1,ny
do i = 1,nx-1
call MUSCL(rho(i,j-2), rho(i,j-1), rho(i,j), rho(i,j+1), rhoL, rhoR)
call MUSCL(u(i,j-2), u(i,j-1), u(i,j), u(i,j+1), uL, uR)
call MUSCL(v(i,j-2), v(i,j-1), v(i,j), v(i,j+1), vL, vR)
call MUSCL(p(i,j-2), p(i,j-1), p(i,j), p(i,j+1), pL,pR)
call StegerWarming1D(rhoL, uL*dxdj(i,j)+vL*dydj(i,j), -uL*dydj(i,j)+vL*dxdj(i,j), pL, rhoR, uR*dxdj(i,j)+vR*dydj(i,j), -uR*dydj(i,j)+vR*dxdj(i,j), pR, flux)
fluxJ(1,i,j) = flux(1)
fluxJ(2,i,j) = flux(2)*dxdj(i,j) - flux(3)*dydj(i,j)
fluxJ(3,i,j) = flux(2)*dydj(i,j) + flux(3)*dxdj(i,j)
fluxJ(4,i,j) = flux(4)
end do
end do
do j = 1, ny-1
do i = 1, nx-1
residual(1,i,j) = fluxI(1,i+1,j)*lengthI(i+1,j) - fluxI(1,i,j)*lengthI(i,j) + fluxJ(1,i,j+1)*lengthJ(i,j+1) - fluxJ(1,i,j)*lengthJ(i,j)
residual(2,i,j) = fluxI(2,i+1,j)*lengthI(i+1,j) - fluxI(2,i,j)*lengthI(i,j) + fluxJ(2,i,j+1)*lengthJ(i,j+1) - fluxJ(2,i,j)*lengthJ(i,j)
residual(3,i,j) = fluxI(3,i+1,j)*lengthI(i+1,j) - fluxI(3,i,j)*lengthI(i,j) + fluxJ(3,i,j+1)*lengthJ(i,j+1) - fluxJ(3,i,j)*lengthJ(i,j)
residual(4,i,j) = fluxI(4,i+1,j)*lengthI(i+1,j) - fluxI(4,i,j)*lengthI(i,j) + fluxJ(4,i,j+1)*lengthJ(i,j+1) - fluxJ(4,i,j)*lengthJ(i,j)
dt(i,j) = cfl*area(i,j)/((dsqrt(gamma*p(i,j)/rho(i,j))+dabs(u(i,j)*dxdi(i,j)+v(i,j)*dydi(i,j)))*lengthI(i,j)+(dsqrt(gamma*p(i,j)/rho(i,j))+dabs(u(i,j)*dxdj(i,j)+v(i,j)*dydj(i,j)))*lengthJ(i,j))
q(:,i,j) = q(:,i,j) - dt(i,j)*residual(:,i,j)/area(i,j)
rho(i,j) = q(1,i,j)
u(i,j) = q(2,i,j) / q(1,i,j)
v(i,j) = q(3,i,j) / q(1,i,j)
p(i,j) = (gamma - 1.0d0) * ( q(4,i,j) - 0.5d0*rho(i,j)*(u(i,j)*u(i,j)+v(i,j)*v(i,j)) )
end do
end do
end
subroutine MUSCL(u1,u2,u3,u4,uL,uR)
implicit none
real*8 u1, u2, u3, u4, uL, uR
uL = u2 + 0.5d0*minmod(u2 - u1,u3 - u2)
uR = u3 - 0.5d0*minmod(u3 - u2,u4 - u3)
end subroutine MUSCL
function minmod(a, b)
implicit none
real*8 a, b, minmod
if(a * b < 0.0d0) then
minmod = 0.0d0
elseif(abs(a) < abs(b)) then
minmod = a
else
minmod = b
end if
end function minmod
subroutine StegerWarming1D(rhoL, uL, vL, pL, rhoR, uR, vR, pR, flux)
implicit none
real*8, parameter:: eps = 0.1d0
real*8 rhoL, uL, vL, pL, rhoR, uR, vR, pR, flux(4), fluxP(4), fluxM(4), cL, cR
real*8 temp1, temp2, lamda1, lamda2, lamda3
temp1=2.d0*(gamma - 1.0d0)
temp2=(3.d0-gamma)/(2.d0*(gamma - 1.0d0))
cL = dsqrt(gamma*pL/rhoL)
cR = dsqrt(gamma*pR/rhoR)
lamda1 = (uL + dsqrt(uL**2 + eps**2)) * 0.5d0
lamda2 = (uL - cL + dsqrt((uL - cL)**2 + eps**2)) * 0.5d0
lamda3 = (uL + cL + dsqrt((uL + cL)**2 + eps**2)) * 0.5d0
fluxP(1) = temp1*lamda1 + lamda2 + lamda3
fluxP(2) = temp1*lamda1*uL + lamda2*(uL-cL) + lamda3*(uL+cL)
fluxP(3) = temp1*lamda1*vL + lamda2*vL + lamda3*vL
fluxP(4) = lamda1*(gamma - 1.0d0)*(uL*uL+vL*vL) + lamda2*((uL-cL)**2+vL*vL)*0.5d0 + lamda3*((uL+cL)**2+vL*vL)*0.5d0 + temp2*cL*cL*(lamda2+lamda3)
fluxP = fluxP * rhoL / (2.d0 * gamma)
lamda1 = (uR - dsqrt(uR**2 + eps**2)) * 0.5d0
lamda2 = (uR - cR - dsqrt((uR - cR)**2 + eps**2)) * 0.5d0
lamda3 = (uR + cR - dsqrt((uR + cR)**2 + eps**2)) * 0.5d0
fluxM(1) = temp1*lamda1 + lamda2 + lamda3
fluxM(2) = temp1*lamda1*uR + lamda2*(uR-cR) + lamda3*(uR+cR)
fluxM(3) = temp1*lamda1*vR + lamda2*vR + lamda3*vR
fluxM(4) = lamda1*(gamma - 1.0d0)*(uR*uR+vR*vR) + lamda2*((uR-cR)**2+vR*vR)*0.5d0 + lamda3*((uR+cR)**2+vR*vR)*0.5d0 + temp2*cR*cR*(lamda2+lamda3)
fluxM = fluxM * rhoR / (2.d0 * gamma)
flux = fluxP+fluxM
end subroutine StegerWarming1D
subroutine Output
implicit none
integer i, j
open(unit = 11, file = 'solution.plt')
write(11,*) 'variables="x","y","rho","u","v","p","Ma"'
write(11,*) 'zone i=', nx, ', j=', ny
do j = 1, ny
do i = 1, nx
WRITE(11,'(20(1X,e18.8))') x(i,j),y(i,j),0.25d0*(rho(i,j) + rho(i-1,j) + rho(i,j-1) + rho(i-1,j-1)),0.25d0*(u(i,j) + u(i-1,j) + u(i,j-1) + u(i-1,j-1)),0.25d0*(v(i,j) + v(i-1,j) + v(i,j-1) + v(i-1,j-1)),0.25d0*(p(i,j) + p(i-1,j) + p(i,j-1) + p(i-1,j-1)),0.25d0*(dsqrt(u(i,j)**2+v(i,j)**2)/dsqrt(gamma*p(i,j)/rho(i,j)) + dsqrt(u(i-1,j)**2+v(i-1,j)**2)/dsqrt(gamma*p(i-1,j)/rho(i-1,j)) + dsqrt(u(i,j-1)**2+v(i,j-1)**2)/dsqrt(gamma* p(i,j-1)/rho(i,j-1)) + dsqrt(u(i-1,j-1)**2+v(i-1,j-1)**2)/dsqrt(gamma*p(i-1,j-1)/rho(i-1,j-1)))
end do
end do
close(11)
end
end program Main