program duct include 'ductin' c c duct computes the fully developed viscous laminar flow in a c rectangular channel of aspect ratio, b/a (= bar), by fem (linear c elements) and fdm (3pt centered differences). c sqr is used to obtain the iterative solution c c input parameters: call pread c write some parameters: call pwrite c determine grid and some parameters needed for finite vol method: call grid c determine exact solution: call vel c start solution: call initcon c iterate using sor it = 0 do while( (it .lt. niter) .and. (rms .gt. eps) ) call sor it = it+1 end do if (rms .gt. eps) write(16,15) niter,rms if (rms .le. eps) write(16,16) it,rms 15 format(//,' Convergence NOT achieved in',i5,' steps',5x, + ' rms=',e12.5,//) 16 format(//,' Convergence achieved in',i5,' steps',5x, + ' rms=',e12.5,//) c compare solution with exact call check write(14) me,it,rms,om close(14) c stop end c ------------ : ------------ subroutine pread include 'ductin' c input data: c parameters explained in input file duct.dat open(11,file='duct.dat') open(16,file='duct.out') read(11,1) me,nterm,ipr,niter read(11,2) bar,eps,om 1 format(i8) 2 format(e10.3) close(11) return end c ------------ : ------------ subroutine pwrite include 'ductin' open(14,file='duct.bin',form='unformatted') write(14) nx,ny open(6,file='duct.out') if (me .eq. 1) write(16,3) if (me .eq. 2) write(16,4) write(16,5) nx,ny,nterm,bar,niter,eps,om 3 format(' duct flow by fem: linear elements') 4 format(' duct flow by fdm: 3pt centerd differences') 5 format(' nx=', i3,' ny=', i3,' nterm=', i3,' b/a=',f5.2, + /,' max iter=',i4, + ' conv tol=',e10.3,' om=',f6.3,//) c return end c ------------ : ------------ subroutine grid include 'ductin' integer j,k c uniform grid dx = 2./(nx-1) dxs = dx*dx dy = 2./(ny-1) dys = dy*dy bas = bar*bar/dxs par1 = bas + 1./dys par2 = (2./dys - bas)/3. par3 = (2.*bas - 1./dys)/3. c c set exact solution and initial solutions c do 10 j = 1,nx 10 x(j) = -1. + (j-1)*dx do 20 k = 1,ny 20 y(k) = -1. + (k-1)*dy write(16,7) (x(k),k=1,nx) write(16,8) (y(k),k=1,ny) write(16,9) write(14) x,y 7 format(' x=', 20f8.4) 8 format(' y=', 20f8.4) 9 format(/) return end c ------------ : ------------ subroutine vel include 'ductin' real xan(nx),yan(ny),con,cx,cy,sig,den,dum integer i,j,k,l c c for given x, y computes axial velocity in channel c and computes the flow rate c con = (8./pi/pi)**2 do 10 l = 1,nx do 10 k = 1,ny xan(l) = 0.5*pi*x(l) yan(k) = 0.5*pi*y(k) fexact(l,k) = 0. 10 continue do 110 l = 1,nx do 110 k = 1,ny do 100 i = 1,nterm,2 cx = cos(i*xan(l)) do 100 j = 1,nterm,2 cy = cos(j*yan(k)) ij = (i+j)/2 -1 if (ij .eq. 0) sig = 1. if (ij .gt. 0) sig = (-1.)**ij den = i*j*( (i*bar)**2 + j*j ) dum = sig*cx*cy/den fexact(l,k) = fexact(l,k) + dum 100 continue fexact(l,k) = con*fexact(l,k) 110 continue flowex = 0. do 120 i = 1,nterm,2 do 120 j = 1,nterm,2 ij = (i+j)/2 -1 if (ij .eq. 0) sig = 1. if (ij .gt. 0) sig = (-1.)**ij den = i*j*( (i*bar)**2 + j*j ) flowex = flowex + 1./den/j/i 120 continue flowex = (16./pi/pi)*con*flowex if (ipr .ne. 0) write(16,6) do 200 l = 1,nx 200 if (ipr .ne. 0) write(16,7) (fexact(l,k),k=1,ny) if (ipr .eq. 0) write(16,8) 6 format(' fexact=',/) 7 format(41f8.4) 8 format(/) return end c ------------ : ------------ subroutine initcon include 'ductin' integer l,k do 10 l = 1,nx do 10 k = 1,ny c f(l,k) = fexact(l,k) f(l,k) = 0 10 continue rms = 2.*eps return end c ------------ : ------------ subroutine sor include 'ductin' integer j,jm,jp,k,km,kp real wdd,dum,sum c c iterate using sor c if ( me .eq. 1 ) then sum = 0. do 13 j = 2,nx-1 jm = j-1 jp = j+1 do 13 k = 2,ny-1 km = k-1 kp = k+1 dum = 1. + par1/6.*( f(jm,kp)+ f(jp,kp)+f(jm,km)+f(jp,km) ) dum = dum + par2*( f(j,kp)+ f(j,km) ) + + par3*( f(jm,k)+f(jp,k) ) wdd = 0.75*dum/par1 dif = wdd - f(j,k) sum = sum + dif*dif f(j,k) = f(j,k) + om*dif 13 continue end if if ( me .ne. 1 ) then sum = 0. do 14 j = 2,nx-1 jm = j-1 jp = j+1 do 14 k = 2,ny-1 km = k-1 kp = k+1 dum = 1. + bas*(f(jm,k)+f(jp,k)) + (f(j,kp)+ f(j,km))/dys wdd = 0.5*dum/par1 dif = wdd - f(j,k) sum = sum + dif*dif f(j,k) = f(j,k) + om*dif 14 continue end if rms = sqrt((sum/nx)/ny) return end c ------------ : ------------ subroutine check include 'ductin' integer j,k,nhx,nhy real dif,sum c compare solution with exact sum = 0. flow = 0. write(16,12) do 19 j = 2,nx-1 do 17 k = 2,ny-1 dif = f(j,k) - fexact(j,k) sum = sum + dif*dif flow = flow + f(j,k) 17 continue if (ipr .ne. 0) write(16,18) (f(j,k),k=2,ny-1) 19 continue flow = flow*dx*dy rms = sqrt((sum/nx)/ny) write(16,20) rms nhx = nx/2+1 nhy = ny/2+1 write(16,21) f(nhx,nhy),fexact(nhx,nhy) write(16,22) flow,flowex write(14) f,fexact write(14) flow,flowex 12 format(' f=') 18 format(41f8.4) 20 format(/,' difference between computed and exact solution,', + ' sol rms=',e12.5,//) 21 format(' flow at center: f=',f8.4,' fexact(c/l) =',f8.4,//) 22 format(' comp. flow rate=',e12.5,' exact flow rate=',e12.5,//) return end