二维非定常热传导问题的有限体积法数值模拟求解(fortan)

program main
!有限体积参数
integer i,j
integer,parameter:: imax=52,jmax=52 !内部结点个数为imax-2,jmax-2,分别加上边界条件，共有imax和jmax个
real,parameter:: lx=1.0,ly=1.0 !x和y方向总长度
real,parameter:: dx=lx/(imax-2),dy=ly/(imax-2) !x和y方向的步长
real,parameter:: pc=1.0,k=1.0 !k：扩散系数
real t_total,dt !计算总时间和时间步长



!边界条件
real,parameter:: tem_w=20,tem_e=20,tem_s=20,tem_n=20 !边界温度

!网格划分
real aw(imax,jmax),ae(imax,jmax),an(imax,jmax),as(imax,jmax),tem(imax,jmax),tem_old(imax,jmax),&
ap_old(imax,jmax),ap(imax,jmax),sp(imax,jmax),su(imax,jmax)

real x(imax),y(jmax)
real delta,t
delta=0.001
t=0
t_total=0.1
dt=0.0001
!节点坐标
x(1)=0.0
x(2)=0.5*dx
x(imax)=lx
do i=3,imax-1
x(i)=x(i-1)+dx
enddo
y(1)=0.0
y(2)=0.5*dy
y(imax)=ly
do j=3,jmax-1
y(j)=y(j-1)+dy
enddo

!定义温度的初始条件和边界条件
do i=1,imax
do j=1,jmax
tem(i,j)=0.0
tem_old(i,j)=100.0
enddo
enddo
do j=1,jmax
tem(1,j)=tem_w
tem_old(1,j)=tem_w
tem(imax,j)=tem_e
tem_old(imax,j)=tem_e
enddo
do i=1,imax
tem(i,1)=tem_s
tem_old(i,1)=tem_s
tem(i,jmax)=tem_n
tem_old(i,jmax)=tem_n
enddo

!计算离散方程的各项系数aw,ae,as,an,ap_old,sp,su,ap
do i=2,imax-1
do j=2,jmax-1
aw(i,j)=k*dy/dx
ae(i,j)=k*dy/dx
as(i,j)=k*dx/dy
an(i,j)=k*dx/dy
ap_old(i,j)=pc*dy*dx/dt
sp(i,j)=0.0
su(i,j)=0.0
enddo
enddo
! left side wall
do j=2,jmax-1
aw(2,j)=0
ae(2,j)=k*dy/dx
as(2,j)=k*dx/dy
an(2,j)=k*dx/dy
ap_old(2,j)=pc*dy*dx/dt
sp(2,j)=-2*k*dy/dx
su(2,j)=2*k*tem_w
enddo
! right side wall
do j=2,jmax-1
aw(imax-1,j)=k*dy/dx
ae(imax-1,j)=0
as(imax-1,j)=k*dx/dy
an(imax-1,j)=k*dx/dy
ap_old(imax-1,j)=pc*dy*dx/dt
sp(imax-1,j)=-2*k*dy/dx
su(imax-1,j)=2*k*tem_e
enddo
!bottom side wall
do i=2,imax-1
aw(i,2)=k*dy/dx
ae(i,2)=k*dy/dx
as(i,2)=0
an(i,2)=k*dx/dy
ap_old(i,2)=pc*dy*dx/dt
sp(i,2)=-2*k*dy/dx
su(i,2)=2*k*tem_s
enddo
!top side wall
do i=2,imax-1
aw(i,jmax-1)=k*dy/dx
ae(i,jmax-1)=k*dy/dx
as(i,jmax-1)=k*dx/dy
an(i,jmax-1)=0
ap_old(i,jmax-1)=pc*dy*dx/dt
sp(i,jmax-1)=-2*k*dy/dx
su(i,jmax-1)=2*k*tem_n
enddo
! special nodes
! (2,2)点
aw(2,2)=0
ae(2,2)=k*dy/dx
as(2,2)=0
an(2,2)=k*dx/dy
ap_old(2,2)=pc*dy*dx/dt
sp(2,2)=-4*k
su(2,2)=2*k*(tem_w+tem_s)
! (2,jamx-1)点
aw(2,jmax-1)=0
ae(2,jmax-1)=k*dy/dx
as(2,jmax-1)=k*dx/dy
an(2,jmax-1)=0
ap_old(2,jmax-1)=pc*dy*dx/dt
sp(2,jmax-1)=-4*k
su(2,jmax-1)=2*k*(tem_w+tem_n)
! (imax-1,2)点
aw(imax-1,2)=k*dy/dx
ae(imax-1,2)=0
as(imax-1,2)=0
an(imax-1,2)=k*dx/dy
ap_old(imax-1,2)=pc*dy*dx/dt
sp(2,2)=-4*k

su(2,2)=2*k*(tem_e+tem_s)
! (imax-1,jmax-1)点
aw(imax-1,2)=k*dy/dx
ae(imax-1,2)=0
as(imax-1,2)=k*dx/dy
an(imax-1,2)=0
ap_old(imax-1,2)=pc*dy*dx/dt
sp(2,2)=-4*k
su(2,2)=2*k*(tem_e+tem_n)
do i=2,imax-1
do j=2,jmax-1
ap(i,j)=aw(i,j)+ae(i,j)+as(i,j)+an(i,j)+ap_old(i,j)-sp(i,j)
enddo
enddo

!sor 迭代算法计算温度值
10 do i=2,imax-1
do j=2,jmax-1
tem(i,j)=(aw(i-1,j)*tem(i-1,j)+ae(i+1,j)*tem(i+1,j)+as(i,j-1)*tem(i,j-1)+an(i,j+1)*tem(i,j+1)+ap_old(i,j)*tem_old(i,j)+su(i,j))/ap(i,j)
enddo
enddo
do i=2,imax-1
do j=2,jmax-1
if(abs(tem(i,j)-tem_old(i,j))/tem(i,j)goto 20
else
tem_old(i,j)=tem(i,j) !不满足精度要求,继续迭代
goto 10
end if
enddo
enddo
20 do i=2,imax-1
do j=2,jmax-1
t=t+dt
if(ttem_old(i,j)=tem(i,j)
goto 10
else
goto 30 !到达预定时间要求后，输出温度值。
endif
enddo
enddo
30 do i=1,imax
do j=1,jmax
write(*,*) x(i),y(j),tem(i,j),ap(i,j)
enddo
enddo


open (1,file='temprature.dat',status='replace')
do i=1,imax
do j=1,jmax
write(1,*) x(i),y(j),tem(i,j)
enddo
enddo
close(1)
end