[推荐]经典PISO算法程序

hyacinth · 发表于 2005-8-2 10:53

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？我要加入

x

C#####################################################################C
C ME58:269 COMPUTATIONAL FLUID DYNAMICS AND HEAT TRANSFER C
C C
C C.-L. Lin C
C#####################################################################C
c234567890123456789012345678901234567890123456789012345678901234567890
program final_project
include 'param_piso.f'
open(unit=2,file='piso_out.plt')
open(unit=4,file='piso_par.dat')
call init
if (lunst) open(unit=5,file='piso_con.plt')
if (lmoni) open(unit=3,file='piso_mon.plt')
c------------------------------ Beginning of time loop
1 continue
time=time+dt
itim=itim+1
c===================================================================
c Outer iteration
c===================================================================
do iter=1,itermax
if (lcalcv) call veloc
if (lcalcp) then
call press
call press2
end if
if (lcalct) call temps
c------------------------------ check whether u(m)=u(m-1)=u(n+1);
if (lmoni) write(3,999) iter, (resu(i),i=1,4)
999 format(i5,4(1pe12.5,1x))
source=max(resu(1),resu(2),resu(3),resu(4))
if (source.lt.srmin) go to 20
if (source.gt.srmax) then
print *,'iter=', iter,' diveragence with source=',source
stop
end if
end do
20 continue
c===================================================================
c check whether u(n+1)=u(n), steady solutions
c===================================================================
if (lunst) then
rest(1)=0.
rest(2)=0.
rest(3)=0.
do i=2,nxm
do ij=li(i)+2,li(i)+nym
rest(1)=rest(1)+abs(u(ij)-uo(ij))
rest(2)=rest(2)+abs(v(ij)-vo(ij))
rest(3)=rest(3)+abs(t(ij)-to(ij))
end do
end do
rest(1)=rest(1)*fnxy
rest(2)=rest(2)*fnxy
rest(3)=rest(3)*fnxy
write(5,999) itim, time, (rest(i),i=1,3)
end if
if (itim.gt.itimmax) go to 30
c===================================================================
c assign new solutions <u, v, t> to <uo, vo, to>
c and start a new iteration if necessary.
c===================================================================
do ij=1,nxy
uo(ij)=u(ij)
vo(ij)=v(ij)
to(ij)=t(ij)
end do
c-- begins a new time step
go to 1
c===================================================================
30 continue
c------------------------------ output
write(2,*) 'VARIABLES="x" "y" "u" "v" "p" "t"'
write(2,*) 'ZONE F=POINT,I=',nxm2,',J=',nym2
do i=2,nxm
do j=2,nym
ij=li(i)+j
write(2,9999) (i-2)*dx,(j-2)*dy,u(ij),v(ij),p(ij),t(ij)
9999 format(2(f8.4,1x),4(E12.5,1x))
end do
end do
stop
end
C#####################################################################C
subroutine init
include 'param_piso.f'
C#####################################################################C
c234567890123456789012345678901234567890123456789012345678901234567890
time=0.0
itim=0
c
c parameters: rho, density; vis, viscosity; prm, Prandtl number
c grav, gravity; beta, coef of thermal expansion.
c
rho=1.0
vis=1.e-3
prm=0.1
alpha=vis/prm
grav=0.
c beta=1/rho*(d\rho/dt). No minus sign. Be careful about the source term
c in y-momentum equation.
beta=0.
th=0.
tc=0.
tref=0.
ulid=1.
xl=1.0
yl=1.0
dt=1.E20
ipiso=1
if (ipiso.eq.0) then
c -----SIMPLE
urf(1)=0.8
urf(2)=0.8
urf(3)=0.2
urf(4)=0.9
else
c -----PISO
urf(1)=0.8
urf(2)=0.8
urf(3)=1.0
urf(4)=1.0
end if
dx = xl/float(nxm2)
dy = yl/float(nym2)
vol = dx*dy
vdydx=vis*dy/dx
vdxdy=vis*dx/dy
adydx=vdydx/prm
adxdy=vdxdy/prm
lcalcv=.true.
lcalcp=.true.
lcalct=.false.
lmoni =.false.
lunst =.false.
apt=0.
if (lunst) apt=vol/dt
srmax=1.e+3
srmin=1.e-4
fnxy=1./float(nxm2)/float(nym2)
do i=1,4
resu(i)=0.
end do
do i=1,nx
li(i)=(i-1)*ny
end do
c===================================================================
c dimensionless parameter
c===================================================================
re=ulid*yl/vis
rei=1./re
gr=grav*beta*(th-tc)*yl**3/vis**2
ra=gr*prm
grre=gr/re/re
prre=1./re/pr
c234567890123456789012345678901234567890123456789012345678901234567890
write(4,777)re,gr,ra,prm
777 format('Re =',1pe12.5,' Gr=',1pe12.5,' Ra=',1pe12.5,
& ' Pr=',1pe12.5)
c===================================================================
c initial fields
c===================================================================
uin=0.
vin=0.
pin=0.
tin=0.
do ij=1,nxy
u(ij) =uin
v(ij) =vin
p(ij) =pin
t(ij) =tin
uo(ij)=uin
vo(ij)=vin
to(ij)=tin
end do
do ij=1,nxy
ae1(ij)=0.
aw1(ij)=0.
an1(ij)=0.
as1(ij)=0.
ae2(ij)=0.
aw2(ij)=0.
an2(ij)=0.
as2(ij)=0.
ae3(ij)=0.
aw3(ij)=0.
an3(ij)=0.
as3(ij)=0.
uprime(ij)=0.
vprime(ij)=0.
utilde(ij)=0.
vtilde(ij)=0.
flxe(ij)=0.
flxw(ij)=0.
flxn(ij)=0.
flxs(ij)=0.
flxe2(ij)=0.
flxw2(ij)=0.
flxn2(ij)=0.
flxs2(ij)=0.
su(ij)=0.
sv(ij)=0.
apu(ij)=0.
apv(ij)=0.
end do
do i=2,nxm
do j=2,nym
ij=li(i)+j
flxe(ij)=+0.5*(u(ij)+u(ij+ny))*dy
flxw(ij)=-0.5*(u(ij)+u(ij-ny))*dy
flxn(ij)=+0.5*(v(ij)+v(ij+1) )*dx
flxs(ij)=-0.5*(v(ij)+v(ij-1) )*dx
end do
end do
c====================================================================
c assign values to fictitious nodes according to boundary conditions.
c====================================================================
do i=2,nxm
ij=li(i)+ny
u(ij)=2.*ulid-u(ij-1)
c uo(ij)=u(ij)
end do
do j=2,nym
c isothermal wall
ij=li(1)+j
t(ij)=2.*tc-t(ij+ny)
c to(ij)=t(ij)
ij=li(nx)+j
t(ij)=2.*th-t(ij-ny)
c to(ij)=t(ij)
end do
c
return
end
C#####################################################################C
subroutine veloc
include 'param_piso.f'
C#####################################################################C
c234567890123456789012345678901234567890123456789012345678901234567890
c
do i=2,nxm
do j=2,nym
ij=li(i)+j
dpx(ij)=(0.5*(p(ij+ny)+p(ij))-0.5*(p(ij)+p(ij-ny)))/dx
dpy(ij)=(0.5*(p(ij+1) +p(ij))-0.5*(p(ij)+p(ij-1) ))/dy
end do
end do
do i=2,nxm
do ij=li(i)+2,li(i)+nym
cemin=min(flxe(ij),0)
cwmin=min(flxw(ij),0)
cnmin=min(flxn(ij),0)
csmin=min(flxs(ij),0)
cemax=max(flxe(ij),0)
cwmax=max(flxw(ij),0)
cnmax=max(flxn(ij),0)
csmax=max(flxs(ij),0)
ae1(ij)=cemin-vdydx
aw1(ij)=cwmin-vdydx
an1(ij)=cnmin-vdxdy
as1(ij)=csmin-vdxdy
ae2(ij)= ae1(ij)
aw2(ij)= aw1(ij)
an2(ij)= an1(ij)
as2(ij)= as1(ij)
c234567890123456789012345678901234567890123456789012345678901234567890
su(ij)= apt*uo(ij) - dpx(ij)*vol
& +cemin*u(ij+ny)+cemax*u(ij)-flxe(ij)*0.5*(u(ij+ny)+u(ij))
& +cwmin*u(ij-ny)+cwmax*u(ij)-flxw(ij)*0.5*(u(ij-ny)+u(ij))
& +cnmin*u(ij+1) +cnmax*u(ij)-flxn(ij)*0.5*(u(ij+1) +u(ij))
& +csmin*u(ij-1) +csmax*u(ij)-flxs(ij)*0.5*(u(ij-1) +u(ij))
sv(ij)= apt*vo(ij) - dpy(ij)*vol
& +cemin*v(ij+ny)+cemax*v(ij)-flxe(ij)*0.5*(v(ij+ny)+v(ij))
& +cwmin*v(ij-ny)+cwmax*v(ij)-flxw(ij)*0.5*(v(ij-ny)+v(ij))
& +cnmin*v(ij+1) +cnmax*v(ij)-flxn(ij)*0.5*(v(ij+1) +v(ij))
& +csmin*v(ij-1) +csmax*v(ij)-flxs(ij)*0.5*(v(ij-1) +v(ij))
end do
end do
c
c------ Buoyancy force
c
if (lcalct) then
do i=2,nxm
do ij=li(i)+2,li(i)+nym
sv(ij)=sv(ij) - beta*grav*(t(ij)-tref)*vol
end do
end do
end if
call bounduv
c===================================================================
c ap for x-momentum
c===================================================================
do i=2,nxm
do ij=li(i)+2,li(i)+nym
ap(ij)=apt-ae1(ij)-aw1(ij)-an1(ij)-as1(ij)
end do
end do
C
C.....PRINT INITIAL COEFFICIENT IF DESIRED
C
c ------ adjust coefficients of boundary nodes for x-momentum equation
c as a result of shear stress
do i=2,nxm
c ------ South
ij=li(i)+2
ap(ij)=ap(ij)+2.*vdxdy
c ------ North
ij=li(i)+nym
ap(ij)=ap(ij)+2.*vdxdy
end do
c
c===================================================================
c.....under-relaxation, solving equation system for u-velocity
c ------ store apu=1/ap from u-momentum for pressure correction
c===================================================================
c
do i=2,nxm
do ij=li(i)+2,li(i)+nym
ap(ij)=ap(ij)/urf(1)
su(ij)=su(ij)+(1.-urf(1))*ap(ij)*u(ij)
apu(ij)=1./ap(ij)
end do
end do
call sipsolm(nx,ny,ae1,aw1,an1,as1,ap,u,su,resu(1))
c===================================================================
c ap for y-momentum
c===================================================================
do i=2,nxm
do ij=li(i)+2,li(i)+nym
ap(ij)=apt-ae2(ij)-aw2(ij)-an2(ij)-as2(ij)
end do
end do
c ------ adjust coefficients of boundary nodes for y-momentum equation
c as a result of shear stress
do j=2,nym
c ------ East
ij=li(nxm)+j
ap(ij)=ap(ij)+2.*vdydx
c ------ West
ij=li(2)+j
ap(ij)=ap(ij)+2.*vdydx
end do
c
c===================================================================
c.....under-relaxation, solving equation system for v-velocity
c ------ store apv=1/ap from v-momentum for pressure correction
c===================================================================
c
do i=2,nxm
do ij=li(i)+2,li(i)+nym
ap(ij)=ap(ij)/urf(2)
sv(ij)=sv(ij)+(1.-urf(2))*ap(ij)*v(ij)
apv(ij)=1./ap(ij)
end do
end do
call sipsolm(nx,ny,ae2,aw2,an2,as2,ap,v,sv,resu(2))
return
end
C#####################################################################C
subroutine press
include 'param_piso.f'
C#####################################################################C
c234567890123456789012345678901234567890123456789012345678901234567890
do i=2,nxm
do ij=li(i)+2,li(i)+nym
flxe(ij)=+0.5*(u(ij)+u(ij+ny))*dy
flxw(ij)=-0.5*(u(ij)+u(ij-ny))*dy
flxn(ij)=+0.5*(v(ij)+v(ij+1) )*dx
flxs(ij)=-0.5*(v(ij)+v(ij-1) )*dx
end do
end do
c===================================================================
c assign fluxes adjacent to fictitious CVs.
c===================================================================
do i=2,nxm
ij=li(i)+2
flxs(ij)=0.
ij=li(i)+nym
flxn(ij)=0.
end do
do j=2,nym
ij=li(2) +j
flxw(ij)=0.
ij=li(nxm)+j
flxe(ij)=0.
end do
sumo=0.
do i=2,nxm
do ij=li(i)+2,li(i)+nym
sumo=sumo+(flxe(ij)+flxw(ij)+flxn(ij)+flxs(ij))**2
end do
end do
sumo=sqrt(sumo*fnxy)
vole=vol
do i=2,nxm2
do ij=li(i)+2,li(i)+nym
apue=0.5*(apu(ij)+apu(ij+ny))
finc=dy*vole*apue
& *((p(ij+ny)-p(ij))/dx-0.5*(dpx(ij)+dpx(ij+ny)))
flxe(ij)=flxe(ij)-finc
flxw(ij+ny)=-flxe(ij)
ae3(ij)=-dy*dy*apue
aw3(ij+ny)=ae3(ij)
end do
end do
voln=vol
do i=2,nxm
do ij=li(i)+2,li(i)+nym2
apvn=0.5*(apv(ij)+apv(ij+1))
finc=dx*voln*apvn
& *((p(ij+1)-p(ij))/dy-0.5*(dpy(ij)+dpy(ij+1)))
flxn(ij)=flxn(ij)-finc
flxs(ij+1)=-flxn(ij)
an3(ij)=-dx*dx*apvn
as3(ij+1)=an3(ij)
end do
end do
C sum=0.
do i=2,nxm
do ij=li(i)+2,li(i)+nym
ap(ij)=-(ae3(ij)+aw3(ij)+an3(ij)+as3(ij))
su(ij)=-(flxe(ij)+flxw(ij)+flxn(ij)+flxs(ij))
pp(ij)=0.
C sum=sum+su(ij)
end do
end do
C print *,'sum=',sum
C if (iter.eq.2) then
C CALL PRINT(flxe,' Fe ')
C CALL PRINT(flxn,' Fn ')
C CALL PRINT(flxw,' Fw ')
C CALL PRINT(flxs,' Fs ')
C CALL PRINT( dpx,' DPX ')
C CALL PRINT( dpy,' DPY ')
C CALL PRINT( apu,' APU ')
C CALL PRINT( apv,' APV ')
C CALL PRINT( su,' SU ')
C STOP
C endif
call sipsolm(nx,ny,ae3,aw3,an3,as3,ap,pp,su,resu(3))
C
c===================================================================
c enforce dpp/dn=0 at boundary nodes.
c no correction is performed on these points.
c use pp(fictitious)=-pp(boundary) ==> linear interpolation at
c boundary faces (pp_f+pp_b)/2=0 ==> no correction on u' is required.
c===================================================================
do i=2,nxm
c ... south
ij=li(i)+1
pp(ij)=pp(ij+1)
c pp(ij)=2.*pp(ij+1)-pp(ij+2)
c ... north
ij=li(i)+ny
pp(ij)=pp(ij-1)
c pp(ij)=2.*pp(ij-1)-pp(ij-2)
end do
do j=2,nym
c ... west
ij=li(1)+j
pp(ij)=pp(ij+ny)
c pp(ij)=2.*pp(ij+ny)-pp(ij+2*ny)
c ... east
ij=li(nx)+j
pp(ij)=pp(ij-ny)
c pp(ij)=2.*pp(ij-ny)-pp(ij-2*ny)
end do
c===================================================================
c if (iter.eq.2) then
c CALL PRINT(AE,' AE ')
c CALL PRINT(AW,' AW ')
c CALL PRINT(AN,' AN ')
c CALL PRINT(AS,' AS ')
c CALL PRINT(AP,' AP ')
c CALL PRINT(PP,' PP ')
c CALL PRINT(SU,' SU ')
c STOP
c endif
c===================================================================
c correct fluxes at interior cell face
c===================================================================
do i=2,nxm2
do j=2,nym
ij=li(i)+j
flxe(ij)=flxe(ij)+ae3(ij)*(pp(ij+ny)-pp(ij))
flxw(ij+ny)=-flxe(ij)
end do
end do
do i=2,nxm
do j=2,nym2
ij=li(i)+j
flxn(ij)=flxn(ij)+an3(ij)*(pp(ij+1)-pp(ij))
flxs(ij+1)=-flxn(ij)
end do
end do
c===================================================================
c Check continuity
c===================================================================
sumn=0.
do i=2,nxm
do ij=li(i)+2,li(i)+nym
sumn=sumn+(flxe(ij)+flxw(ij)+flxn(ij)+flxs(ij))**2
end do
end do
sumn=sqrt(sumn*fnxy)
print *,' sumn1=', sumn
c===================================================================
c correct pressure and velocities at cell center.
c pp(ij) is at cell center, calculate pp at cell face (ppe, ppw...)
c then take d(pp)/dx and d(pp)/dy to correct u, v, p at cell center.
c===================================================================
C
C.....VALUE OF P' AT REFERENCE LOCATION TO BE SUBTRACTED FROM ALL P'
C
IJPREF=LI(IPR)+JPR
PPO=PP(IJPREF)
c print *,'ipr=', ipr,'jpr=', jpr,'ppo=', ppo
c
do i=2,nxm
do ij=li(i)+2,li(i)+nym
ppe=0.5*(pp(ij+ny)+pp(ij))
ppw=0.5*(pp(ij-ny)+pp(ij))
ppn=0.5*(pp(ij+1) +pp(ij))
pps=0.5*(pp(ij-1) +pp(ij))
uprime(ij)=-dy*(ppe-ppw)*apu(ij)
vprime(ij)=-dx*(ppn-pps)*apv(ij)
u(ij)=u(ij)+uprime(ij)
v(ij)=v(ij)+vprime(ij)
p(ij)=p(ij)+urf(3)*(pp(ij)-ppo)
end do
end do
c===================================================================
c assign pressure at fictitious cell.
c===================================================================
do i=2,nxm
c ... south
ij=li(i)+1
p(ij)=p(ij+1)
c p(ij)=2.*p(ij+1)-p(ij+2)
c ... north
ij=li(i)+ny
p(ij)=p(ij-1)
c p(ij)=2.*p(ij-1)-p(ij-2)
end do
do j=2,nym
c ... west
ij=li(1)+j
p(ij)=p(ij+ny)
c p(ij)=2.*p(ij+ny)-p(ij+2*ny)
c ... east
ij=li(nx)+j
p(ij)=p(ij-ny)
c p(ij)=2.*p(ij-ny)-p(ij-2*ny)
end do
return
end
C#####################################################################C
subroutine press2
include 'param_piso.f'
C#####################################################################C
c234567890123456789012345678901234567890123456789012345678901234567890
do i=2,nxm
do ij=li(i)+2,li(i)+nym
utilde(ij)=-apu(ij)*(ae1(ij)*uprime(ij+nj)+aw1(ij)*uprime(ij-nj)
+ +an1(ij)*uprime(ij+ 1)+as1(ij)*uprime(ij- 1))
vtilde(ij)=-apv(ij)*(ae2(ij)*vprime(ij+nj)+aw2(ij)*vprime(ij-nj)
+ +an2(ij)*vprime(ij+ 1)+as2(ij)*vprime(ij- 1))
end do
end do
do i=2,nxm
do ij=li(i)+2,li(i)+nym
flxe2(ij)=+0.5*(utilde(ij)+utilde(ij+ny))*dy
flxw2(ij)=-0.5*(utilde(ij)+utilde(ij-ny))*dy
flxn2(ij)=+0.5*(vtilde(ij)+vtilde(ij+1) )*dx
flxs2(ij)=-0.5*(vtilde(ij)+vtilde(ij-1) )*dx
end do
end do
c===================================================================
c assign fluxes adjacent to fictitious CVs.
c===================================================================
do i=2,nxm
ij=li(i)+2
flxs2(ij)=0.
ij=li(i)+nym
flxn2(ij)=0.
end do
do j=2,nym
ij=li(2) +j
flxw2(ij)=0.
ij=li(nxm)+j
flxe2(ij)=0.
end do
C sum=0.
do i=2,nxm
do ij=li(i)+2,li(i)+nym
ap(ij)=-(ae3(ij)+aw3(ij)+an3(ij)+as3(ij))
su(ij)=-(flxe2(ij)+flxw2(ij)+flxn2(ij)+flxs2(ij))
pp(ij)=0.
C sum=sum+su(ij)
end do
end do
C print *,'sum=',sum
C if (iter.eq.2) then
C CALL PRINT(flxe,' Fe ')
C CALL PRINT(flxn,' Fn ')
C CALL PRINT(flxw,' Fw ')
C CALL PRINT(flxs,' Fs ')
C CALL PRINT( dpx,' DPX ')
C CALL PRINT( dpy,' DPY ')
C CALL PRINT( apu,' APU ')
C CALL PRINT( apv,' APV ')
C CALL PRINT( su,' SU ')
C STOP
C endif
call sipsolm(nx,ny,ae3,aw3,an3,as3,ap,pp,su,resu(3))
C
c===================================================================
c enforce dpp/dn=0 at boundary nodes.
c no correction is performed on these points.
c use pp(fictitious)=-pp(boundary) ==> linear interpolation at
c boundary faces (pp_f+pp_b)/2=0 ==> no correction on u' is required.
c===================================================================
do i=2,nxm
c ... south
ij=li(i)+1
pp(ij)=pp(ij+1)
c pp(ij)=2.*pp(ij+1)-pp(ij+2)
c ... north
ij=li(i)+ny
pp(ij)=pp(ij-1)
c pp(ij)=2.*pp(ij-1)-pp(ij-2)
end do
do j=2,nym
c ... west
ij=li(1)+j
pp(ij)=pp(ij+ny)
c pp(ij)=2.*pp(ij+ny)-pp(ij+2*ny)
c ... east
ij=li(nx)+j
pp(ij)=pp(ij-ny)
c pp(ij)=2.*pp(ij-ny)-pp(ij-2*ny)
end do
c===================================================================
c if (iter.eq.2) then
c CALL PRINT(AE,' AE ')
c CALL PRINT(AW,' AW ')
c CALL PRINT(AN,' AN ')
c CALL PRINT(AS,' AS ')
c CALL PRINT(AP,' AP ')
c CALL PRINT(PP,' PP ')
c CALL PRINT(SU,' SU ')
c STOP
c endif
c===================================================================
c correct fluxes at interior cell face
c===================================================================
do i=2,nxm2
do j=2,nym
ij=li(i)+j
flxe(ij)=flxe(ij)+flxe2(ij)+ae3(ij)*(pp(ij+ny)-pp(ij))
flxw(ij+ny)=-flxe(ij)
end do
end do
do i=2,nxm
do j=2,nym2
ij=li(i)+j
flxn(ij)=flxn(ij)+flxn2(ij)+an3(ij)*(pp(ij+1)-pp(ij))
flxs(ij+1)=-flxn(ij)
end do
end do
c===================================================================
c Check continuity
c===================================================================
sumn=0.
do i=2,nxm
do ij=li(i)+2,li(i)+nym
sumn=sumn+(flxe(ij)+flxw(ij)+flxn(ij)+flxs(ij))**2
end do
end do
sumn=sqrt(sumn*fnxy)
print *,' sumn2=', sumn
c===================================================================
c correct pressure and velocities at cell center.
c pp(ij) is at cell center, calculate pp at cell face (ppe, ppw...)
c then take d(pp)/dx and d(pp)/dy to correct u, v, p at cell center.
c===================================================================
C
C.....VALUE OF P' AT REFERENCE LOCATION TO BE SUBTRACTED FROM ALL P'
C
IJPREF=LI(IPR)+JPR
PPO=PP(IJPREF)
c print *,'ipr=', ipr,'jpr=', jpr,'ppo=', ppo
c
do i=2,nxm
do ij=li(i)+2,li(i)+nym
ppe=0.5*(pp(ij+ny)+pp(ij))
ppw=0.5*(pp(ij-ny)+pp(ij))
ppn=0.5*(pp(ij+1) +pp(ij))
pps=0.5*(pp(ij-1) +pp(ij))
u(ij)=u(ij)+utilde(ij)-dy*(ppe-ppw)*apu(ij)
v(ij)=v(ij)+vtilde(ij)-dx*(ppn-pps)*apv(ij)
p(ij)=p(ij)+urf(3)*(pp(ij)-ppo)
end do
end do
c===================================================================
c assign pressure at fictitious cell.
c===================================================================
do i=2,nxm
c ... south
ij=li(i)+1
p(ij)=p(ij+1)
c p(ij)=2.*p(ij+1)-p(ij+2)
c ... north
ij=li(i)+ny
p(ij)=p(ij-1)
c p(ij)=2.*p(ij-1)-p(ij-2)
end do
do j=2,nym
c ... west
ij=li(1)+j
p(ij)=p(ij+ny)
c p(ij)=2.*p(ij+ny)-p(ij+2*ny)
c ... east
ij=li(nx)+j
p(ij)=p(ij-ny)
c p(ij)=2.*p(ij-ny)-p(ij-2*ny)
end do
return
end
C#####################################################################C
subroutine temps
include 'param_piso.f'
C#####################################################################C
c234567890123456789012345678901234567890123456789012345678901234567890
c
do i=2,nxm
do ij=li(i)+2,li(i)+nym
cemin=min(flxe(ij),0)
cwmin=min(flxw(ij),0)
cnmin=min(flxn(ij),0)
csmin=min(flxs(ij),0)
cemax=max(flxe(ij),0)
cwmax=max(flxw(ij),0)
cnmax=max(flxn(ij),0)
csmax=max(flxs(ij),0)
ae1(ij)=cemin-adydx
aw1(ij)=cwmin-adydx
an1(ij)=cnmin-adxdy
as1(ij)=csmin-adxdy
ap(ij)=apt-ae1(ij)-aw1(ij)-an1(ij)-as1(ij)
c234567890123456789012345678901234567890123456789012345678901234567890
su(ij)= apt*to(ij)
& +cemin*t(ij+ny)+cemax*t(ij)-flxe(ij)*0.5*(t(ij+ny)+t(ij))
& +cwmin*t(ij-ny)+cwmax*t(ij)-flxw(ij)*0.5*(t(ij-ny)+t(ij))
& +cnmin*t(ij+1) +cnmax*t(ij)-flxn(ij)*0.5*(t(ij+1) +t(ij))
& +csmin*t(ij-1) +csmax*t(ij)-flxs(ij)*0.5*(t(ij-1) +t(ij))
end do
end do
call boundt
c
c===================================================================
c.....under-relaxation, solving equation system for temperature
c===================================================================
c
do i=2,nxm
do ij=li(i)+2,li(i)+nym
ap(ij)=ap(ij)/urf(4)
su(ij)=su(ij)+(1.-urf(4))*ap(ij)*t(ij)
end do
end do
call sipsolm(nx,ny,ae1,aw1,an1,as1,ap,t,su,resu(4))
c===================================================================
c assign values at fictitious nodes to satisfy boundary conditions.
c===================================================================
do i=2,nxm
c south: adiabatic wall
ij=li(i)+1
t(ij)=t(ij+1)
c north: adiabatic wall
ij=li(i)+ny
t(ij)=t(ij-1)
end do
do j=2,nym
c west: T_c isothermal wall
ij=li(1) +j
t(ij)=2*tc-t(ij+ny)
c east: T_h isothermal wall
ij=li(nx)+j
t(ij)=2*th-t(ij-ny)
end do
return
end
C#####################################################################C
subroutine bounduv
include 'param_piso.f'
C#####################################################################C
c234567890123456789012345678901234567890123456789012345678901234567890
c normal stress du/dx=0, dv/dy and shear stress du/dy, dv/dx
c East: (du/dx)_e=0
do j=2,nym
ij=li(nxm)+j
ae1(ij)=0.
ae2(ij)=0.
end do
c West: (du/dx)_w=0
do j=2,nym
ij=li(2)+j
aw1(ij)=0.
aw2(ij)=0.
end do
c South: absorb as into ap, u_s=2*ubottom-u_p, ubottom=0.
do i=2,nxm
ij=li(i)+2
as1(ij)=0.
as2(ij)=0.
end do
c North: absorb an into ap, u_n=2*ulid-u_p
do i=2,nxm
ij=li(i)+nym
an1(ij)=0.
an2(ij)=0.
su(ij)=su(ij)+2*vdxdy*ulid
end do
return
end
C#####################################################################C
subroutine boundt
include 'param_piso.f'
C#####################################################################C
c234567890123456789012345678901234567890123456789012345678901234567890
c East: isothermal wall; absorb ae into ap, t_e=2*th-t_p
do j=2,nym
ij=li(nxm)+j
ae1(ij)=0.
ap(ij)=ap(ij)+adydx
su(ij)=su(ij)+2*adydx*th
end do
c West: isothermal wall; absorb aw into ap, t_w=2*tc-t_p
do j=2,nym
ij=li(2)+j
aw1(ij)=0.
ap(ij)=ap(ij)+adydx
su(ij)=su(ij)+2*adydx*tc
end do
c South: adiabatic wall
do i=2,nxm
ij=li(i)+2
as1(ij)=0.
ap(ij)=ap(ij)-adxdy
end do
c North: adiabatic wall
do i=2,nxm
ij=li(i)+nym
an1(ij)=0.
ap(ij)=ap(ij)-adxdy
end do
return
end
C#####################################################################C
SUBROUTINE SIPSOLM(MI,MJ,GE,GW,GN,GS,GP,FI,Q,RES0)
C#####################################################################C
C REVISED FOR ME58:269 FINAL PROJECT
C
C MI: # of actual + fictitious CVs in the x-dir.
C MJ: # of actual + fictitious CVs in the y-dir.
C
C GW*FI(K-MJ)+GS*FI(K-1)+GP*FI(K)+GN*FI(K+1)+GE*FI(K+MJ)=Q(K)
C
C FI=INITIAL GUESS. FOR EXAMPLE,
c USE OLD U FOR U-COMPONENT VELOCITY.
C
C-------------------------------------------------------------
C This is the ILU solver after Stone; see Sect. 5.3.4 for
C a description of the algorithm.
C
C M. Peric, Institut fuer Schiffbau, Hamburg, 1995
C=============================================================
PARAMETER (MAXITSIP=20, MX=102, MY=102, MXY=MX*MY)
REAL LW,LS,LPR,RES0
DIMENSION GE(1),GW(1),GN(1),GS(1),GP(1),Q(1)
DIMENSION FI(MXY),UN(MXY),UE(MXY),RES(MXY)
DIMENSION LW(MXY),LS(MXY),LPR(MXY)
DIMENSION LI(MX)
IF (MI.GT.MX) THEN
PRINT *,'MI IS GREATER THAN MX IN SR. SIPSOLM'
STOP
ENDIF
ALFA=0.92
RESMAX=0.2
MIM=MI-1
MJM=MJ-1
MIJ=MI*MJ
DO I=1,MI
LI(I)=(I-1)*MJ
END DO
C
C.....CALCULATE ELEMENTS OF [L] AND [U] MATRICES
C
DO I=2,MIM
DO IJ=LI(I)+2,LI(I)+MJM
LW(IJ)=GW(IJ)/(1.+ALFA*UN(IJ-MJ))
LS(IJ)=GS(IJ)/(1.+ALFA*UE(IJ-1))
P1=ALFA*LW(IJ)*UN(IJ-MJ)
P2=ALFA*LS(IJ)*UE(IJ-1)
LPR(IJ)=1./(GP(IJ)+P1+P2-LW(IJ)*UE(IJ-MJ)-LS(IJ)*UN(IJ-1))
UN(IJ)=(GN(IJ)-P1)*LPR(IJ)
UE(IJ)=(GE(IJ)-P2)*LPR(IJ)
END DO
END DO
C
C.....CALCULATE RESIDUAL AND AUXILLIARY VECTORS; INNER ITERATION LOOP
C
DO N=1,MAXITSIP
C
RESN=0.
DO I=2,MIM
DO IJ=LI(I)+2,LI(I)+MJM
RES(IJ)=Q(IJ)-GP(IJ)*FI(IJ)-GN(IJ)*FI(IJ+1)-
* GS(IJ)*FI(IJ-1)-GE(IJ)*FI(IJ+MJ)-GW(IJ)*FI(IJ-MJ)
RESN=RESN+ABS(RES(IJ))
RES(IJ)=(RES(IJ)-LS(IJ)*RES(IJ-1)-LW(IJ)*RES(IJ-MJ))*LPR(IJ)
END DO
END DO
IF(N.EQ.1) RES0=RESN
C
C.....CALCULATE INCREMENT AND CORRECT VARIABLE
C
DO I=MIM,2,-1
DO IJ=LI(I)+MJM,LI(I)+2,-1
RES(IJ)=RES(IJ)-UN(IJ)*RES(IJ+1)-UE(IJ)*RES(IJ+MJ)
FI(IJ)=FI(IJ)+RES(IJ)
END DO
END DO
C
C.....CONVERGENCE CHECK
C
RSM=RESN/(RES0+1.E-20)
IF(RSM.LT.RESMAX) RETURN
END DO
C
PRINT *,'EXCEEDS MAXITSIP IN SR. SIPSOL, RSM = ',RSM
PRINT *,'RESMAX=', RESMAX, 'MAXITSIP=', MAXITSIP
C
RETURN
END
C
C#####################################################################C
C#############################################################
SUBROUTINE PRINT(PHI,HEDPHI)
C#############################################################
C This routine prints 2D array in an easy to read format.
C--------------------------------------------------------------
include 'param_piso.f'
PARAMETER (NPHI=4)
DIMENSION PHI(NXY)
CHARACTER*6 HEDPHI
C--------------------------------------------------------
C
WRITE(3,20) HEDPHI
NL=(NX-1)/12+1
C
DO L=1,NL
IS=(L-1)*12+1
IE=MIN(NX,L*12)
WRITE(3,21) (I,I=IS,IE)
WRITE(3,22)
C
DO J=NY,1,-1
WRITE(3,23) J,(PHI(LI(I)+J),I=IS,IE)
END DO
END DO
C
20 FORMAT(2X,26('*-'),5X,A6,5X,26('-*'))
21 FORMAT(3X,'I = ',I3,11I10)
22 FORMAT(2X,'J')
23 FORMAT(1X,I3,1P12E10.2)
C
RETURN
END

复制代码

[ 本帖最后由 yejet 于 2006-11-19 19:39 编辑 ]

dgl0611 · 发表于 2009-1-2 16:07

太伟大了，楼主我也是学传热的，有机会交流一下，173251077

账号		自动登录	找回密码
密码			我要加入

[Fortran] [推荐]经典PISO算法程序

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。