****************************************** NON-LINEAR VISCOUS FLUD AROUND A PLATE/ TURBULENCES: ****************************************** program fluid character new character*25 d,fname real maxvalue,conv,init,R,totforce,omega,sum integer, parameter :: nx=300 integer, parameter :: ny=150 integer rnx,rny,nxp,nyp,display,ph,pd,dist,freq,choice integer, dimension(nx,ny)::xlab,ylab real, dimension(nx,ny)::psi,opsi,ozeta,mask,zeta real, dimension(nx,2*ny)::colour include 'pict_subs.h' 10 format('rescale_factor',2I3) 11 format('size_of_plate',2I3) 12 format('plate_distance',I3) 13 format('convergence',F8.4) 14 format('Reynolds_number',F8.4) 15 format('frequency_of_display',I3) 16 format('relaxation_parameter',F8.4) 17 format(5X,I5,21X,F12.4) CMF$ layout colour(:news, :news) CMF$ layout xlab(:news, :news) CMF$ layout ylab(:news, :news) CMF$ layout mask(:news, :news) CMF$ layout psi(:news, :news) CMF$ layout opsi(:news, :news) CMF$ layout ozeta(:news, :news) CMF$ layout zeta(:news, :news) c get parameter from file---------------------------------- c the parameters can be changed in the file 'parameter' before c implementing the program or by typing 'y' after execution 19 open(8,file='parameter',status='old') rewind 8 read(8,*) d,rnx,rny read(8,*) d,ph,pd read(8,*) d,dist read(8,*) d,conv read(8,*) d,R read(8,*) d,freq read(8,*) d,omega write(*,*) 'new parameter? (y,n) ' read(*,*) new if(new.eq.'y') then print*, '(a$)','Enter rescale factor for the image (x y): ' read*, rnx,rny print*, 'Enter size of plate (heigth depth) ' read*, ph,pd print*, 'Enter distance of plate from origin: ' read*, dist print*, 'Enter convergence criterium ' read*, conv print*, 'Enter Reynolds number ' read*, R print*, 'Enter frequency of updating display' read*, freq print*, 'Enter relaxation parameter ' read*, omega print*, 'Save new parameter? (y,n) ' read*, new if(new.eq.'y') then rewind 8 write(8,10) rnx,rny write(8,11) ph,pd write(8,12) dist write(8,13) conv write(8,14) R write(8,15) freq write(8,16) omega end if end if close(8) print*, 'write force values to file ? (Enter filename, or "no")' read*, fname print*, 'choose Display:' print*, '(1) stream function & vorticity' print*, '(2) stream function' print*, '(3) vorticity' read(*,*) choice c graphics-------------------------------------------------- nxp = rnx*nx nyp = rny*ny call FSR_select_display_menu(depth,nxp,nyp) call FSR_set_display_color_map("density") display = FSR_allocate_image_buffer(depth,nxp, nyp,.TRUE.) c making maskfile ------------------------------------------- c mask has an identifying number for each boundary c mask is zero (even) and 9 (odd) in the interior of the grid c and 10 at the plate. c xlab and ylab are two other maskfiles indicating the values of c x- and y-axis in terms of h. forall(i=1:nx,j=1:ny) xlab(i,j)=i forall(i=1:nx,j=1:ny) ylab(i,j)=ny+1-j where(mod(xlab+ylab,2).eq.1) mask=9 elsewhere mask=0 end where where(xlab.eq.1) mask=1 where(ylab.eq.ny) mask=2 where(xlab.eq.nx) mask=3 where(ylab.eq.1) mask=4 c placing the plate-------------------------------------------- where(xlab.gt.dist.and.xlab.le.dist+pd.and.ylab.le.ph) mask=10 where(xlab.eq.dist.and.ylab.le.ph) mask=7 where(xlab.gt.dist.and.xlab.le.dist+pd.and.ylab.eq.ph+1) mask=6 where(xlab.eq.dist+1+pd.and.ylab.le.ph) mask=5 c initial value for velocity---------------------------------- psi=ylab/2 c initial value for vorticity--------------------------------- zeta=0 c ensure that the boundary conditions are satisfied----------- call psi_boundary_cond(psi,mask,nx,ny) call zeta_boundary_cond(zeta,psi,mask,nx,ny) ozeta=zeta opsi=psi write(*,*) 'number of iteration Force on the plate' c gauss-seidel algorithem for interior of grid---------------- c and beginning of iteration-loop n=0 maxvalue=2*conv do 20 while(maxvalue.gt.conv) ! CONVERGENCECRITERIUM n=n+1 ! AND BEGIN OF ITERATION-LOOP where(mask.eq.0) psi = & ((cshift(psi,1,1)+cshift(psi,1,-1) & +cshift(psi,2,1)+cshift(psi,2,-1)-zeta)/4.) where(mask.eq.9) psi = & ((cshift(psi,1,1)+cshift(psi,1,-1) & +cshift(psi,2,1)+cshift(psi, 2,-1)-zeta)/4.) psi=omega*psi+(1.-omega)*opsi c ensure that the boundary condition--s are satisfied----------- call psi_boundary_cond(psi,mask,nx,ny) call zeta_boundary_cond(zeta,psi,mask,nx,ny) c obtain zeta------------------------------------------------- call eval_zeta(zeta,psi,mask,nx,ny,R) zeta=omega*zeta+(1.-omega)*ozeta call zeta_boundary_cond(zeta,psi,mask,nx,ny) c evaluate force on the plate---------------------------------- sum=0 do i=dist+2,dist+pd-1 sum=sum+zeta(i,ny-ph) end do sum=sum+(zeta(dist+1,ny-ph)+zeta(dist+pd,ny-ph))/2. totforce=1/(R*ph)*sum c write out the force function over iteration step------------- if (fname.ne.'no') then open(30,file=fname,status='unknown') write(30,*) n,totforce,sum end if c display stream function, vorticity and force----------------- c freq confirms how often the display is updated if (mod(n,freq).eq.0) then where(mask.eq.10) psi=-3 zeta=-.05 end where if (choice.eq.1) then forall(j=1:ny) colour(:,j)= x sqrt(abs(psi(:,j)+3)/(maxval(psi)+3)) forall(j=ny+1:2*ny) colour(:,j)= x sqrt(abs(zeta(:,2*ny+1-j)+.05)/(maxval(zeta)+.05)) else if (choice.eq.2) then forall(j=1:ny) colour(:,j)= x sqrt(abs(psi(:,j)+3)/(maxval(psi)+3)) forall(j=ny+1:2*ny) colour(:,j)= x sqrt(abs(psi(:,2*ny+1-j)+3)/(maxval(psi)+3)) else if(choice.eq.3) then forall(j=1:ny) colour(:,j)= x sqrt(abs(zeta(:,j)+.05)/(maxval(zeta)+.05)) forall(j=ny+1:2*ny) colour(:,j)= x sqrt(abs(zeta(:,2*ny+1-j)+.05)/(maxval(zeta)+.05)) end if call display_rescaled_image(nx,2*ny,nxp,nyp, x colour,display) write(*,17) n,totforce end if c monitor convergence ------------------------------------------ maxvalue=maxval(opsi-psi) opsi=psi ozeta=zeta 20 continue ! END OF ITERATION-LOOP c -------------------------------------------------------------- close(30) print*, 'convergencecriterium reached at step ',n print* 99 print*, 'once again, another display or quit (a,d,q)' read(*,*) new if (new.eq.'a') goto 19 if (new.eq.'d') then print*, 'choose Display:' print*, '(1) stream function & vorticity' print*, '(2) stream function' print*, '(3) vorticity' read(*,*) choice if (choice.eq.1) then forall(j=1:ny) colour(:,j)= x sqrt(abs(psi(:,j)+3)/(maxval(psi)+3)) forall(j=ny+1:2*ny) colour(:,j)= x sqrt(abs(zeta(:,2*ny+1-j)+.05)/(maxval(zeta)+.05)) else if (choice.eq.2) then forall(j=1:ny) colour(:,j)= x sqrt(abs(psi(:,j)+3)/(maxval(psi)+3)) forall(j=ny+1:2*ny) colour(:,j)= x sqrt(abs(psi(:,2*ny+1-j)+3)/(maxval(psi)+3)) else if(choice.eq.3) then forall(j=1:ny) colour(:,j)= x sqrt(abs(zeta(:,j)+.05)/(maxval(zeta)+.05)) forall(j=ny+1:2*ny) colour(:,j)= x sqrt(abs(zeta(:,2*ny+1-j)+.05)/(maxval(zeta)+.05)) end if call display_rescaled_image(nx,2*ny,nxp,nyp, x colour,display) end if goto 99 stop end c------SUBROUTINES-------------------------------------------- c compute boundary values for velocity------------------------ subroutine psi_boundary_cond(psi,mask,nx,ny) real, dimension(nx,ny)::psi,mask cmf$ layout psi(:news, :news) cmf$ layout mask(:news, :news) where(mask.eq.2) psi=cshift(psi,2,1)+1. ! G-BOUNDARY where(mask.eq.1) psi=cshift(psi,1,1) ! F-BOUNDARY where(mask.eq.3) psi=cshift(psi,1,-1) ! H-BOUNDARY c ABCDE-BOUNDARY where(mask.eq.4.or.mask.eq.5.or.mask.eq.6.or.mask.eq.7) psi=0 end subroutine c evaluate vorticity from given velocity--------------------- subroutine eval_zeta(zeta,psi,mask,nx,ny,R) real R real, dimension(nx,ny)::zeta,psi,mask cmf$ layout zeta(:news, :news) cmf$ layout psi(:news, :news) cmf$ layout mask(:news, :news) where(mask.eq.0) zeta= x (cshift(zeta,1, 1)+cshift(zeta,1,-1) x +cshift(zeta,2,-1)+cshift(zeta,2, 1))/4. x -R/16.*((cshift(psi, 2,-1)-cshift(psi, 2, 1)) x *(cshift(zeta,1, 1)-cshift(zeta,1,-1)) x -(cshift(psi, 1, 1)-cshift(psi, 1,-1)) x *(cshift(zeta,2,-1)-cshift(zeta,2, 1))) where(mask.eq.9) zeta = x (cshift(zeta,1, 1)+cshift(zeta,1,-1) x +cshift(zeta,2,-1)+cshift(zeta,2, 1))/4. x -R/16.*((cshift(psi, 2,-1)-cshift(psi, 2, 1)) x *(cshift(zeta,1, 1)-cshift(zeta,1,-1)) x -(cshift(psi, 1, 1)-cshift(psi, 1,-1)) x *(cshift(zeta,2,-1)-cshift(zeta,2, 1))) end subroutine c compute boundary values for vorticity---------------------- subroutine zeta_boundary_cond(zeta,psi,mask,nx,ny) real, dimension(nx,ny)::zeta,psi,mask cmf$ layout zeta(:news, :news) cmf$ layout psi(:news, :news) cmf$ layout mask(:news, :news) where(mask.eq.1.or.mask.eq.2.or.mask.eq.4) zeta=0 ! FGAE where(mask.eq.3) zeta=cshift(zeta,1,-1) ! H-BOUND where(mask.eq.5) zeta=(cshift(psi,1,1))*2. ! B-BOUND where(mask.eq.6) zeta=(cshift(psi,2,-1))*2. ! C-BOUND where(mask.eq.7) zeta=(cshift(psi,1,-1))*2. ! D-BOUND end subroutine