C This program solves the 2-D Steady Incompressible Navier-Stokes C equations with using Successive Over Relaxation (SOR) method C Visit http://www.cavityflow.com C C s(i,j) ==> streamfunction variable C v(i,j) ==> vorticity variable C x(i) ==> x-coordinate C y(j) ==> y-coordinate C dh ==> grid spacing C Re ==> Reynolds Number program main implicit double precision (a-h,o-z) parameter(N=128) common / flow variables / >s(0:N,0:N),v(0:N,0:N), >s_old(0:N,0:N),v_old(0:N,0:N) common / geometry / >x(0:N),y(0:N) write(*,*) "Input Re" read(*,*) Re Write(*,*) "Re=",Re beta=0.5 dh=1.0d0/dble(N) do 1 k=0,N x(k)=dble(k)/dble(N) y(k)=dble(k)/dble(N) 1 continue do 999 iteration=1,10000 do 99 iter=1,5000 do 3 i=0,N do 3 j=0,N s_old(i,j)=0.0 v_old(i,j)=0.0 3 continue C Iterate on Streamfunction do 4 i=1,N-1 do 4 j=1,N-1 s(i,j)=beta*(0.25d0*( > s(i-1,j)+s(i+1,j)+s(i,j-1)+s(i,j+1) >+dh**2.*v(i,j) > ) > )+(1.d0-beta)*s(i,j) 4 continue C Calculate Vorticity at the Boundaries using Jensen's formula do 5 i=1,N-1 v(i,0)=0.5d0*(-8.d0*s(i,1)+s(i,2))/dh**2. v(i,N)=0.5d0*(-8.d0*s(i,N-1)+s(i,N-2))/dh**2.-3.d0/dh 5 continue do 6 j=1,N-1 v(0,j)=0.5d0*(-8.d0*s(1,j)+s(2,j))/dh**2. v(N,j)=0.5d0*(-8.d0*s(N-1,j)+s(N-2,j))/dh**2. 6 continue C Calculate Vorticity at the corners C NOTE:We do not need these points for interior solution C We only need these points for plotting purposes C For these boundary conditions please refer to: C T. Stortkuhl, C. Zenger, S. Zimmer, "An Asymptotic Solution for C the Singularity at the Angular Point of the Lid Driven Cavity", C International Journal of Numerical Methods for Heat & Fluid Flow C 1994, Vol 4, pp 47--59 v(0,0)=(-(s(1,1))/(3.d0*dh**2.)-(0.5d0*v(1,0)+0.5d0*v(0,1) >+0.25d0*v(1,1))/(9.d0))*9.d0 v(N,0)=(-(s(N-1,1))/(3.d0*dh**2.)-(0.5d0*v(N-1,0)+0.5d0*v(N,1) >+0.25d0*v(N-1,1))/(9.d0))*9.d0 v(N,N)=(-0.5d0/dh-(s(N-1,N-1))/(3.d0*dh**2.)-(0.5d0*v(N-1,N) >+0.5d0*v(N,N-1)+0.25d0*v(N-1,N-1))/(9.d0))*9.d0 v(0,N)=(-0.5d0/dh-(s(1,N-1))/(3.d0*dh**2.)-(0.5d0*v(1,N) >+0.5d0*v(0,N-1)+0.25d0*v(1,N-1))/(9.d0))*9.d0 C Iterate on Vorticity do 7 i=1,N-1 do 7 j=1,N-1 v(i,j)=beta*(0.25d0*( > v(i-1,j)+v(i+1,j)+v(i,j-1)+v(i,j+1) >-0.25d0*Re*(s(i,j+1)-s(i,j-1))*(v(i+1,j)-v(i-1,j)) >+0.25d0*Re*(s(i+1,j)-s(i-1,j))*(v(i,j+1)-v(i,j-1)) > ) > )+(1.d0-beta)*v(i,j) 7 continue 99 continue res1f_A=0.d0 res1g_A=0.d0 res2f_A=0.d0 res2g_A=0.d0 res3f_A=0.d0 res3g_A=0.d0 do 8 i=1,N-1 do 8 j=1,N-1 res1f_B=abs( > (s(i-1,j)-2.d0*s(i,j)+s(i+1,j))/dh**2. >+(s(i,j-1)-2.d0*s(i,j)+s(i,j+1))/dh**2. >+v(i,j) ) res1g_B=abs( > (v(i-1,j)-2.d0*v(i,j)+v(i+1,j))/(Re*dh**2.) >+(v(i,j-1)-2.d0*v(i,j)+v(i,j+1))/(Re*dh**2.) >-((s(i,j+1)-s(i,j-1))/(2.0d0*dh))*((v(i+1,j)-v(i-1,j))/(2.0d0*dh)) >+((s(i+1,j)-s(i-1,j))/(2.0d0*dh))*((v(i,j+1)-v(i,j-1))/(2.0d0*dh)) > ) res1f_A=max(res1f_A,res1f_B) res1g_A=max(res1g_A,res1g_B) res2f_B=abs(s(i,j)-s_old(i,j)) res2g_B=abs(v(i,j)-v_old(i,j)) res2f_A=max(res2f_A,res2f_B) res2g_A=max(res2g_A,res2g_B) res3f_B=abs((s(i,j)-s_old(i,j))/s_old(i,j)) res3g_B=abs((v(i,j)-v_old(i,j))/v_old(i,j)) res3f_A=max(res3f_A,res3f_B) res3g_A=max(res3g_A,res3g_B) 8 continue write(*,*) iteration,iter,res1f_A,res1g_A C Output the residuals c open(2,file='residuals.txt',ACCESS='APPEND') c write(2,*) 'iteration number=',iteration*50000 c write(2,*) 'res1f=',res1f_A,' res1g=',res1g_A c write(2,*) 'res2f=',res2f_A,' res2g=',res2g_A c write(2,*) 'res3f=',res3f_A,' res3g=',res3g_A c close(2) C Output the variables c open(3,file='outputfile.dat') c rewind(3) c do 9 i=0,N c do 9 j=0,N c write(3,"(4e14.6)") i*dh,j*dh,s(i,j),v(i,j) c 9 continue c close(3) C Condition to stop the code if((res1f_A.lt.1.d-9).AND.(res1g_A.lt.1.d-9)) then goto 10 endif 999 continue 10 continue open(4,file='final.dat') rewind(3) do i=0,N do j=0,N write(4,"(4e14.6)") i*dh,j*dh,-s(i,j),v(i,j) end do write(4,*) end do close(4) stop end