anb1.f90

anb 版的LBM程序 This code was written to show beginners in a simple and c short way the relevant proce

          if (obst(x,y))then
            g(1,x,y)=g_1 !g_hlp(3,x,y)+g_1
            g(2,x,y)=g_1 !g_hlp(4,x,y)+g_1
            g(3,x,y)=g_1 !g_hlp(1,x,y)+g_1
            g(4,x,y)=g_1 !g_hlp(2,x,y)+g_1
            g(5,x,y)=g_2 !g_hlp(7,x,y)+g_2
            g(6,x,y)=g_2 !g_hlp(8,x,y)+g_2
            g(7,x,y)=g_2 !g_hlp(5,x,y)+g_2
            g(8,x,y)=g_2 !g_hlp(6,x,y)+g_2
          end if
   101 continue
      return
      end

      subroutine propagate(lx,ly,node,n_hlp,g,g_hlp)
      implicit none
      integer  lx,ly
      real*8   node(0:8,lx,ly),n_hlp(0:8,lx,ly),g(0:8,lx,ly),g_hlp(0:8,lx,ly)
      integer  x,y,x_e,x_w,y_n,y_s
      do 10 x=1, lx
        do 10 y=1, ly
          y_n=mod(y,ly)+1
          x_e=mod(x,lx)+1
          y_s=ly-mod(ly+1-y, ly)
          x_w=lx-mod(lx+1-x, lx)
          n_hlp(0,x  ,y  )=node(0,x,y)
          n_hlp(1,x_e,y  )=node(1,x,y)
          n_hlp(2,x  ,y_n)=node(2,x,y)
          n_hlp(3,x_w,y  )=node(3,x,y)
          n_hlp(4,x  ,y_s)=node(4,x,y)
          n_hlp(5,x_e,y_n)=node(5,x,y)
          n_hlp(6,x_w,y_n)=node(6,x,y)
          n_hlp(7,x_w,y_s)=node(7,x,y)
          n_hlp(8,x_e,y_s)=node(8,x,y)
    10 continue
	  do 101 x=1, lx
        do 101 y=1, ly
          y_n=mod(y,ly)+1
          x_e=mod(x,lx)+1
          y_s=ly-mod(ly+1-y, ly)
          x_w=lx-mod(lx+1-x, lx)  
       if(x==1)then
	      g_hlp(0,x  ,y  )=g(0,x,y)
          g_hlp(1,x_e,y  )=g(1,x,y)
          g_hlp(2,x  ,y_n)=g(2,x,y)
          g_hlp(3,x_w,y  )=0
          g_hlp(4,x  ,y_s)=g(4,x,y)
          g_hlp(5,x_e,y_n)=g(5,x,y)
          g_hlp(6,x_w,y_n)=0
          g_hlp(7,x_w,y_s)=0
          g_hlp(8,x_e,y_s)=g(8,x,y)
	      else if(x==lx)then
          g_hlp(0,x  ,y  )=g(0,x,y)
          g_hlp(1,x_e,y  )=0
          g_hlp(2,x  ,y_n)=g(2,x,y)
          g_hlp(3,x_w,y  )=g(3,x,y)
          g_hlp(4,x  ,y_s)=g(4,x,y)
          g_hlp(5,x_e,y_n)=0
          g_hlp(6,x_w,y_n)=g(6,x,y)
          g_hlp(7,x_w,y_s)=g(7,x,y)
          g_hlp(8,x_e,y_s)=0
		  else
		  g_hlp(0,x  ,y  )=g(0,x,y)
          g_hlp(1,x_e,y  )=g(1,x,y)
          g_hlp(2,x  ,y_n)=g(2,x,y)
          g_hlp(3,x_w,y  )=g(3,x,y)
          g_hlp(4,x  ,y_s)=g(4,x,y)
          g_hlp(5,x_e,y_n)=g(5,x,y)
          g_hlp(6,x_w,y_n)=g(6,x,y)
          g_hlp(7,x_w,y_s)=g(7,x,y)
          g_hlp(8,x_e,y_s)=g(8,x,y)
           endif
   101 continue
  
      return
      end

      subroutine relaxation(density,omega_v,lx,ly,node,n_hlp,g,g_hlp,u_hlp_x,u_hlp_y,u_hlp_c,u_hlp_squ,obst)
      implicit none
      integer  lx,ly
      logical  obst(lx,ly)
      real*8   density,omega_g,omega_v
	  real*8   node(0:8,lx,ly),n_hlp(0:8,lx,ly),g(0:8,lx,ly),g_hlp(0:8,lx,ly),q(0:8,lx,ly)
	  real*8   u_hlp_x(lx,ly),u_hlp_y(lx,ly),u_hlp_squ(lx,ly),u_hlp_c(0:8,lx,ly)
      integer  x,y,i
      real*8  c_squ,t_0,t_1,t_2,u_x,u_y,u_n(8),u_squa(lx,ly)
	  real*8  n_equ(0:8),g_equ(0:8),u_squ,d_loc,g_loc,u(lx,ly)
      t_0=4.d0/9.d0
      t_1=1.d0/9.d0
      t_2=1.d0/36.d0
      c_squ=1.d0/3.d0
      do 10 x=1, lx
        do 10 y=1, ly
          if (.not. obst(x,y)) then
            d_loc=0.d0 
         	g_loc=0.d0
            do 20 i=0, 8 
             d_loc=d_loc+n_hlp(i,x,y)
   20       continue
            u_x=(n_hlp(1,x,y)+n_hlp(5,x,y)+n_hlp(8,x,y)-(n_hlp(3,x,y)+n_hlp(6,x,y)+n_hlp(7,x,y)))/d_loc
            u_y=(n_hlp(2,x,y)+n_hlp(5,x,y)+n_hlp(6,x,y)-(n_hlp(4,x,y)+n_hlp(7,x,y)+n_hlp(8,x,y)))/d_loc
            u_squ=u_x*u_x+u_y*u_y
            u_n(1)=  u_x
            u_n(2)=      u_y
            u_n(3) =-u_x
            u_n(4)=     -u_y
            u_n(5)=  u_x+u_y
            u_n(6) =-u_x+u_y
            u_n(7) =-u_x-u_y
            u_n(8)=  u_x-u_y
            n_equ(0)=t_0*d_loc*(1.d0-u_squ/(2.d0*c_squ))
            n_equ(1)=t_1*d_loc*(1.d0+u_n(1)/c_squ+u_n(1)**2.d0/(2.d0*c_squ**2.d0)-u_squ/(2.d0*c_squ))
            n_equ(2)=t_1*d_loc*(1.d0+u_n(2)/c_squ+u_n(2)**2.d0/(2.d0*c_squ**2.d0)-u_squ/(2.d0*c_squ))
            n_equ(3)=t_1*d_loc*(1.d0+u_n(3)/c_squ+u_n(3)**2.d0/(2.d0*c_squ**2.d0)-u_squ/(2.d0*c_squ))
            n_equ(4)=t_1*d_loc*(1.d0+u_n(4)/c_squ+u_n(4)**2.d0/(2.d0*c_squ**2.d0)-u_squ/(2.d0*c_squ))
            n_equ(5)=t_2*d_loc*(1.d0+u_n(5)/c_squ+u_n(5)**2.d0/(2.d0*c_squ**2.d0)-u_squ/(2.d0*c_squ))
            n_equ(6)=t_2*d_loc*(1.d0+u_n(6)/c_squ+u_n(6)**2.d0/(2.d0*c_squ**2.d0)-u_squ/(2.d0*c_squ))
            n_equ(7)=t_2*d_loc*(1.d0+u_n(7)/c_squ+u_n(7)**2.d0/(2.d0*c_squ**2.d0)-u_squ/(2.d0*c_squ))
            n_equ(8)=t_2*d_loc*(1.d0+u_n(8)/c_squ+u_n(8)**2.d0/(2.d0*c_squ**2.d0)-u_squ/(2.d0*c_squ))
         
            u_squa(x,y)=u_squ
            do 201 i=0, 8 
             g_loc=g_loc+g_hlp(i,x,y)
   201       continue

			g_equ(0)=-2.d0*g_loc*u_squ/3.d0
			g_equ(1)=g_loc*(3.d0/2.d0+3.d0*u_n(1)/2.d0+9.d0*u_n(1)**2.d0/2.d0-3.d0*u_squ/2.d0)/9.d0
			g_equ(2)=g_loc*(3.d0/2.d0+3.d0*u_n(2)/2.d0+9.d0*u_n(2)**2.d0/2.d0-3.d0*u_squ/2.d0)/9.d0
			g_equ(3)=g_loc*(3.d0/2.d0+3.d0*u_n(3)/2.d0+9.d0*u_n(3)**2.d0/2.d0-3.d0*u_squ/2.d0)/9.d0
			g_equ(4)=g_loc*(3.d0/2.d0+3.d0*u_n(4)/2.d0+9.d0*u_n(4)**2.d0/2.d0-3.d0*u_squ/2.d0)/9.d0
			g_equ(5)=g_loc*(3.d0+6.d0*u_n(5)+9.d0*u_n(5)**2.d0/2.d0-3.d0*u_squ/2.d0)/36.d0
			g_equ(6)=g_loc*(3.d0+6.d0*u_n(6)+9.d0*u_n(6)**2.d0/2.d0-3.d0*u_squ/2.d0)/36.d0
			g_equ(7)=g_loc*(3.d0+6.d0*u_n(7)+9.d0*u_n(7)**2.d0/2.d0-3.d0*u_squ/2.d0)/36.d0
			g_equ(8)=g_loc*(3.d0+6.d0*u_n(8)+9.d0*u_n(8)**2.d0/2.d0-3.d0*u_squ/2.d0)/36.d0
            
			do 30 i=0, 8
			  if(i.eq.0)then
			  q(i,x,y)=-(u_hlp_x(x,y)*u_x+u_hlp_y(x,y)*u_y)+u_hlp_squ(x,y)
			  else
			  q(i,x,y)=u_n(i)-u_hlp_c(i,x,y)-(u_hlp_x(x,y)*u_x+u_hlp_y(x,y)*u_y)+u_hlp_squ(x,y)
              endif
			  node(i,x,y)=n_hlp(i,x,y)+omega_v*(n_equ(i)-n_hlp(i,x,y))
			  g(i,x,y)=g_hlp(i,x,y)+omega_g*(g_equ(i)-g_hlp(i,x,y))-node(i,x,y)*q(i,x,y)
			  u_hlp_x(x,y)=u_x
			  u_hlp_y(x,y)=u_y
			  if(i.ne.0) u_hlp_c(i,x,y)=u_n(i)
   30       continue
              u_hlp_squ(x,y)=u_squa(x,y)
          end if
   10 continue
      return
      end

      subroutine write_velocity(lx,ly,time,obst,node,vel)   
      implicit none
      integer  lx,ly,time
      logical  obst(lx,ly)
      real*8   node(0:8,lx,ly),vel
      integer  x,y,i,n_free
      real*8  u_x,d_loc
      x=int(float(lx)/2.d0)
      n_free=0
      u_x=0.d0
      do 10 y=1, ly
        if(.not. obst(x,y)) then
          d_loc=0.d0
          do 20 i=0, 8 
            d_loc=d_loc+node(i,x,y)
   20     continue
          u_x=u_x+(node(1,x,y)+node(5,x,y)+node(8,x,y)-(node(3,x,y)+node(6,x,y)+node(7,x,y)))/d_loc
          n_free=n_free+1
        end if
   10 continue
      vel=u_x/float(n_free)
      write(10,*) time, vel
      return
      end

     subroutine  write_results(lx,ly,obst,node,density,g)
      implicit none
      integer  lx,ly
      real*8  node(0:8,lx,ly),density,g(0:8,lx,ly)
      logical  obst(lx,ly)
      integer  x,y,i,obsval
      real*8  u_x,u_y,d_loc,press,c_squ,pres(lx,ly),flux(lx,ly),xvel(lx,ly),in_e,gco(lx,ly)
      c_squ=1.d0/3.d0
      open(11,file='anb.dat')
      write(11,*) 'VARIABLES=X, Y, VX, VY, PRESS, OBST' 
      write(11,*) 'ZONE I=', lx, ', J=', ly, ', F=POINT'
      do 10 x=1, lx
        do 10 y=1, ly
          if (obst(x,y)) then 
            obsval=1
            u_x=0.d0
            u_y=0.d0
            press=density*c_squ
          else
            d_loc=0.d0
            do 20 i=0, 8 
              d_loc=d_loc+node(i,x,y)
   20       continue
            u_x=(node(1,x,y)+node(5,x,y)+node(8,x,y)-(node(3,x,y)+node(6,x,y)+node(7,x,y)))/d_loc
            u_y=(node(2,x,y)+node(5,x,y)+node(6,x,y)-(node(4,x,y)+node(7,x,y)+node(8,x,y)))/d_loc
            press=d_loc*c_squ
            obsval=0
          end if
		  in_e=0.d0
		  do 201 i=0,8
		    in_e=g(i,x,y)+in_e
		201 continue
		  gco(x,y)=in_e
          pres(x,y)=press	
		  flux(x,y)=u_x*d_loc				
		  xvel(x,y)=u_x						 
          write(11,'(x,2i6,4x,5(e15.6,2x))')x,y,u_x,u_y,press,in_e,obsval
   10 continue
      close(11)
	  open(12,file='pressure.dat')
	  write(12,'(120(x,f10.5))')pres
	  close(12)
      open(13,file='flux.dat')
	  write(13,'(120(x,f10.5))')flux
	  close(13)
      open(14,file='xvel.dat')
	  write(14,'(120(x,f10.5))')xvel
	  close(14)
      open(15,file='energe.dat')
	  write(15,'(120(x,f10.5))')gco
	  close(15)
      return
      end

      subroutine comp_rey(lx,ly,obst,node,time,omega_v,density,r_rey,mu)
      implicit none
      integer  lx,ly,time
      real*8  density,node(0:8,lx,ly),omega_v,r_rey
      logical  obst(lx,ly)
      real*8  vel,mu,rey
      call write_velocity(lx,ly,time,obst,node,vel)    
      !visc=1.d0/6.d0*(2.d0/omega-1.d0)
      rey=vel*r_rey/mu
	  open(13,file='rev.txt')
      write(13,*) '*** viscosity=', mu
      write (13,*) '*** average velocity=', vel
      write (13,*) '*** Reynolds number=', rey
	  close(13)
      write (6,*) '*** Calculations finished, results:'
      write (6,*) '***'
      write (6,*) '*** viscosity=', mu
      write (6,*) '*** average velocity=', vel
      write (6,*) '*** Reynolds number=', rey
      write (6,*) '***'
      write (6,*) '*** In the file anb.dat, you can find local'
      write (6,*) '*** information about the simulated flow.'
      write (6,*) '***'
      write (6,*) '*** In the file anb_qx.out, you can find the average'
      write (6,*) '*** flow velocity plotted as a function of time.'
      write (6,*) '***'
      return
      end