超松弛迭代法求解泊松方程

计算流作业

1. 题目

对15cm?12cm的矩形平板，假设板表面绝热，仅在四条边有热流通过，且导热系数为常数。设一条15cm边温度保持为100?，其余三条边温度为20?。计算板内稳定的温度场并绘出等温线。采用迭代法计算时，采用不同的松弛因子，并比较收敛速度，观察在不同的网格尺度下，最佳松弛因子的变化。

2. 理论基础

超松弛SOR迭代法迭代公式为：

?1?im?im,j?(1??),j??2(1??2)m?12m?12(?im?im?1,j??i?1,j??(,j?1??i,j?1)?x?i,j)

3. 结果分析

计算时取x,y方向网格数分别为(4,5),(8,10),为1.0,1.0475,1.095,1.95。控制固定余量为10?5 3.1 温度场

网格尺度为dx?dy?0.5cm时，最佳松弛因子下的平板温度场分布结果基本上是对称的。

(32,40)，松弛因子取值范围

3.2 最佳松弛因子-网格尺度

显然，随着网格尺度增加，最佳松弛因子逐渐减小，基本上呈现指数规律。

4. 附件

!初始温度分布

subroutine initial(tt_in,nx_in,ny_in) implicit none integer::i,j

integer,intent(in)::nx_in,ny_in !横纵网格数

real(kind=8),dimension(0:nx_in,0:ny_in)::tt_in !温度矩阵

tt_in=0.

tt_in(0,0:ny_in)=20.

tt_in(nx_in,0:ny_in)=100. tt_in(1:(nx_in-1),0)=20.

tt_in(1:(nx_in-1),ny_in)=20. do i=1,nx_in-1 do j=1,ny_in-1 tt_in(i,j)=60. end do end do

end subroutine initial !超松弛迭代法

subroutine sor(tt_sor,nx_sor,ny_sor,ome,be)

implicit none integer::i,j

integer,intent(in)::nx_sor,ny_sor !横纵网格数

real(kind=8),dimension(0:nx_sor,0:ny_sor)::tt_sor !温度 real(kind=8),intent(in)::ome,be !松弛因子,网格参数h real(kind=8)::tem1

do i=1,nx_sor-1 do j=1,ny_sor-1

tem1=tt_sor(i+1,j)+tt_sor(i-1,j)+be**2*(tt_sor(i,j+1)+tt_sor(i,j-1)) tt_sor(i,j)=(1-ome)*tt_sor(i,j)+ome*tem1/2./(1+be**2) end do end do end subroutine sor

program poisson implicit none integer::i,j,k,l

integer,parameter::numk=8 !计算时取的网格大小种类数 real(kind=8),allocatable,dimension(:,:)::tem !温度 integer::numx,numy !横纵网格数

real(kind=8),parameter::lx=12.,ly=15. !计算域长宽 integer,parameter::nw=20 !计算时omega取值的数目

integer,dimension(nw)::ncell !存储同一omega下迭代需要的步数 real(kind=8)::om,bt !松弛因子,网格参数h

real(kind=8)::dx,dy !横纵方向上的网格参数 integer::ct !中间变量

real(kind=8)::w_opt !最佳松弛因子

integer,parameter::status=-1 !状态变量.1-计算不同网格尺度下的最佳松弛因子

!-1-计算特定尺度，取最佳松弛因子时得到的最终温度分布

open(unit=15,file='poisson.dat') open(unit=16,file='poisson2.txt') write(16,'(\

write(16,'(\

write(16,'(\ do l=1,8

numx=4*l numy=5*l dx=lx/numx dy=ly/numy bt=dx/dy

allocate(tem(0:numx,0:numy)) ct=10**10 do k=1,nw

om=1.+0.95*dble(k)/dble(nw)

call caclulate_w(tem,numx,numy,om,ncell(k),bt)

write(16,'(1x,d15.5,3x,d15.5,3x,i10)')l,om,ncell(k) if(ncell(k)

if(status==1) then

write(15,'(1x,d15.5,3x,d15.5)')dx,w_opt else if(status==-1 .and. l==6) then write(15,'(\

write(15,'(\

write(15,'(\

call caclulate_w(tem,numx,numy,w_opt,ct,bt) do j=0,numy do i=0,numx

write(15,*)i*dx,j*dy,tem(i,j) end do end do end if

deallocate(tem) end do end program poisson

!计算特定网格和松弛因子下的温度

subroutine caclulate_w(tt,nx,ny,omega,nncell,beta) implicit none

integer::i,j

integer,intent(in)::nx,ny !横纵网格数 real(kind=8)::error=1.d-5 !误差限

integer,intent(out)::nncell !存储迭代次数

real(kind=8)::omega,etem,beta !etem-两次计算间的平均误差

real(kind=8),dimension(0:nx,0:ny)::tem0 !存储上一次计算得到的温度矩阵

real(kind=8),dimension(0:nx,0:ny)::tt !温度矩阵

etem=1.d10

nncell=0

call initial(tt,nx,ny)

do while(etem>error) nncell=nncell+1 tem0=tt

call sor(tt,nx,ny,omega,beta) etem=0. do i=1,nx-1 do j=1,ny-1

etem=etem+(tt(i,j)-tem0(i,j))**2 end do end do

etem=sqrt(etem) end do

end subroutine caclulate_w