|
发表于 2009-7-2 14:56:05
|
显示全部楼层
p2
c---
c left and right walls
c---
Do J=2,Ny
c O(1 ,J) = 2.0*(S(1,J) - S(2,J) )/Dxs
c O(Nx1,J) = 2.0*(S(Nx+1,J) - S(Nx,J))/Dxs
O(1 ,J) = (7.0*S(1 ,J)-8.0*S(2,J) +S(3, J))/(2.0*Dxs)
O(Nx1,J) = (7.0*S(Nx1,J)-8.0*S(Nx,J)+S(Nx-1,J))/(2.0*Dxs)
End Do
c--------------------------------
c Compute the velocity at the grid points
c by central differences
c--------------------------------
Do J=2,Ny
Do I=2,Nx
U(I,J) = (S(I,J+1)-S(I,J-1))/Dy2
V(I,J) = - (S(I+1,J)-S(I-1,J))/Dx2
End Do
End Do
c--------------------------------------------
c Iterate on Poisson's equation for the vorticity
c
c The source term is the rhs of equation (13.1.8)
c of Pozrikidis (1997)
c--------------------------------------------
Iter = 0
92 Continue
Iter = Iter + 1
c---
c Jacobi updating
c---
ErrO = 0.0D0 ! error register
Do J=2,Ny
Do I=2,Nx
source(I,J) = U(I,J)*(O(I+1,J)-O(I-1,J))/Dx2
+ + V(I,J)*(O(I,J+1)-O(I,J-1))/Dy2
source(I,J) = source(I,J)/rnu
res = O(I+1,J)+O(I-1,J)-2.0*O(I,J)
+ + beta*(O(I,J+1)+O(I,J-1)-2.0*O(I,J))
+ - Dxs * source(I,J)
Onew(I,J) = O(I,J) + rel_vort * res
If(abs(res).GT.ErrO) ErrO = abs(res)
End Do
End Do
c---
c replace
c---
Do J=2,Ny
Do I=2,Nx
O(I,J) = Onew(I,J)
End Do
End Do
If(Iter.le.Niter) Go to 92
c-----------------
c printing session
c-----------------
If(IseeV.eq.1) then
write (6,*)
write (6,*) " Vorticity"
write (6,*)
Do J=1,Ny1
write (6,101) (O(I,J),I=1,Nx1)
End Do
End If
write (6,103) Kter,ERRP,ERRO
c-----------------
c advance counters
c-----------------
Jter = Jter + 1
Kter = Kter + 1
If(Jter.le.Niterg) Go to 98
c------------------------
c End of batch iterations
c------------------------
Jter = 1
write (6,*)
write (6,*) " One more set of global iterations ?"
write (6,*)
write (6,*) " Enter 0 for No, 1 for YES "
write (6,*) " ---------------------------"
read (5,*) more
If(more.ne.0) Go to 98
c-------------------------
c-------------------------
c End of global iterations
c-------------------------
c-------------------------
99 Continue
c--------------------------------
c print the final stream function
c--------------------------------
open (4,file="cvt_sv.str")
write (4,*) Nx,Ny
write (4,*) Dx,Dy
Do j=1,Ny1
write (4,101) (S(i,j),i=1,Nx1)
End Do
write (4,*)
write (4,400) UUlid
write (4,401) Lx,Ly
write (4,402) Nx,Ny
write (4,403) visc
write (4,404) rho
write (4,405) rel_psi
write (4,406) rel_vort
write (4,407) Niter
write (4,408) Niterg
close (4)
c--------------------------------
c print the final vorticity field
c--------------------------------
open (4,file="cvt_sv.vort")
write (4,*) Nx,Ny
write (4,*) Dx,Dy
Do j=1,Ny1
write (4,101) (O(i,j),i=1,Nx1)
End Do
write (4,*)
write (4,400) UUlid
write (4,401) Lx,Ly
write (4,402) Nx,Ny
write (4,403) visc
write (4,404) rho
write (4,405) rel_psi
write (4,406) rel_vort
write (4,407) Niter
write (4,408) Niterg
close (4)
c--------------------------------------------------
c print the final upper wall vorticity distribution
c along with the local solution for Stokes flow
c--------------------------------------------------
open (4,file="cvt_sv.vortw")
write (6,*) " upper wall vorticity distribution"
write (6,*)
j=Ny1
write (4,*) Nx1
Do i=1,Nx1
xh = (i-1.0D0)/Nx
Ocorner = 0.0D0
If(xh.ne.0) then
Ocorner = - pi4/(pis-4.0D0) * UUlid/xh
End If
write (4,102) i,xh,O(i,j),Ocorner
write (6,102) i,xh,O(i,j),Ocorner
End Do
write (4,*) Null
close (4)
c-----------------------------
c Plot velocity vectors
c reduced by maximum magnitude
c for visualization
c-----------------------------
open (4,file="cvt_sv.vel")
c---
c print cavity walls
c---
write (4,104) Nfive
write (4,104) None,None,zero,zero
write (4,104) None,None,Lx ,zero
write (4,104) None,None,Lx ,Ly
write (4,104) None,None,zero,Ly
write (4,104) None,None,zero,zero
c-----
c Plot velocity vectors
c reduced by maximum magnitude
c for visualization
c----
Umax = 0.0D0
Do J=2,Ny
Do I=2,Nx
Umag = U(I,J)**2+V(I,J)**2
If(Umag.gt.Umax) Umax = Umag
End Do
End Do
Umax = Dsqrt(Umax)
write (6,*) " Maximum velocity: ", Umax
Do J=2,Ny
Do I=2,Nx
x1 = X(I,J)
y1 = Y(I,J)
Darx = 0.10*U(I,J)/Umax
Dary = 0.10*V(I,J)/Umax
call draw_arrow_2d
+
+ (x1,y1
+ ,Darx,Dary
+ ,vector
+ )
write (4,104) Nfive
Do l=1,5
write (4,104) J,I,vector(l,1),vector(l,2)
End Do
End Do
End Do
write (4,104) Null
write (4,*)
write (4,400) UUlid
write (4,401) Lx,Ly
write (4,402) Nx,Ny
write (4,403) visc
write (4,404) rho
write (4,405) rel_psi
write (4,406) rel_vort
write (4,407) Niter
close (4)
c-----
c Done
c-----
100 Format (2(1x,f15.10))
101 Format (200(1x,f10.5))
102 Format (1X,I4,3(1x,f15.8))
103 Format ("Iter: ",i4,", Max error in psi and omega ",
+ 2(1x,f15.10))
104 Format (2(1x,i2),20(1x,f10.5))
400 Format (1x," UUlid ",F10.5)
401 Format (1x," Lx,Ly ",F10.5,1x,F10.5)
402 Format (1x," Nx,Ny ",I3,1x,I3)
403 Format (1x," Viscosity ",f10.8)
404 Format (1x," Density ",f10.8)
405 Format (1x," rel_psi ",f10.8)
406 Format (1x," rel_vort ",f10.8)
407 Format (1x," Niter ",I5)
408 Format (1x," Niterg ",I5)
Stop
End |
|