|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
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 编辑 ] |
|