anb1.f90

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

      program anb
      implicit none
      integer  lx,ly
      parameter(lx=120,ly=80)
      real*8   density,init_g
      real*8   omega_v,omega_g,mu,a,tao_v,tao_g !mu为粘度，a为导温系数，tao_v为质量分布函数的弛豫时间,
	                                  !tao_g为内能分布函数的弛豫时间
      real*8   accel_d
      integer  t_max
      real*8   r_rey
      integer  time
      logical  error
      logical  obst(lx,ly)
	  real*8  g_wall
      real*8  node(0:8,lx,ly),g(0:8,lx,ly)
      real*8  n_hlp(0:8,lx,ly),g_hlp(0:8,lx,ly),u_hlp_x(lx,ly),u_hlp_y(lx,ly),u_hlp_squ(lx,ly),u_hlp_c(0:8,lx,ly)
      real*8  vel
      write (6,*)
      write (6,*) 'anb 1.0 (2001-03-29)'
      write (6,*) 'Copyright (C) 1998-2001 Joerg Bernsdorf /'
      write (6,*) 'C&C Research Laboratories, NEC Europe Ltd.'
      write (6,*) 'anb comes with ABSOLUTELY NO WARRANTY;' 
      write (6,*) 'Revised by Bi Jincheng.'
      write (6,*) 'All rights reserved.'
      write (6,*)
      write (6,*) '****************************************************'
      write (6,*) '***              anb starting ...                ***'
      write (6,*) '****************************************************'
      write (6,*) '*** Precompiled for lattice size lx=',lx 
      write (6,*) '***                              ly=',ly 
      write (6,*) '****************************************************'
      write (6,*) '***'
      error=.false.
      call read_parametrs(error,t_max,density,init_g,accel_d,mu,a,r_rey,omega_v,omega_g,g_wall)
      if (error) goto 990
      call read_obstacles(error,obst,lx,ly)
      if (error) goto 990
      call init_density(lx,ly,density,init_g,node,g)
	  call init_u_hlp(lx,ly,u_hlp_x,u_hlp_y,u_hlp_c,u_hlp_squ)
      open(10,file='anb_qxm.out') 
      do 100 time=1, t_max
        if (time .ge. 10 .and. mod(time,t_max/10) .eq. 0) then
          call check_density(lx,ly,node,time)
        end if
        call redistribute(lx,ly,obst,node,g,accel_d,density,init_g)
        call propagate(lx,ly,node,n_hlp,g,g_hlp)
        call bounceback(lx,ly,obst,node,n_hlp)
		call thermal_bc(lx,ly,obst,g,g_wall,g_hlp)
        call 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)
        call write_velocity(lx,ly,time,obst,node,vel)  
  100 continue
      call write_results(lx,ly,obst,node,density,g)
      call comp_rey(lx,ly,obst,node,time,omega_v,density,r_rey,mu)
      goto 999
  990 write (6,*) '!!! error: program stopped during iteration =', time
      write (6,*) '!!!'
  999 continue
      close(10)
      write (6,*) '********************    end     ********************'
      stop
      end	

      subroutine read_parametrs(error,t_max,density,init_g,accel_d,mu,a,r_rey,omega_v,omega_g,g_wall)
      implicit none
      real*8  density,init_g,accel_d,mu,a,omega_v,omega_g,r_rey,g_wall
      integer  t_max
      logical  error
      open(5,file='anb.par')
      read(5,*,err=900)t_max
      read(5,*,err=900)density
	  read(5,*,err=900)init_g
      read(5,*,err=900)accel_d
      read(5,*,err=900)mu
      read(5,*,err=900)a
      read(5,*,err=900)g_wall
      read(5,*,err=900)r_rey
      close(5)      
	  omega_v=1.d0/(1.d0/2.d0+3.d0*mu)
	  omega_g=1.d0/(1.d0/2.d0+3.d0*a/2.d0)
      write (6,*) '*** Paramters read from file anb.par.'
      write (6,*) '***'
      goto 999
  900 write (6,*) '!!! Error reading file anb.par'
      write (6,*) '!!!'
      goto 990
  990 error=.true.
  999 continue
      return
      end

    subroutine read_obstacles(error,obst,lx,ly)
      implicit none
      integer  lx,ly
      logical  error,obst(lx,ly)
      integer  x,y
        do 10 y=1, ly
          do 10 x=1, lx
   10       obst(x,y)=.false.
     open(10,file='c:\bi\myfor\heattransfer\debug\anb.txt')
   20 continue
        read(10,*,end=50,err=900) x,y
       if (x .le. lx .and. y .le. ly) then
          obst(x,y)=.true.
       else
          write(6,*) '!!! Obstacle out of range, skipped'
          write(6,*) '!!! lx=', x, ' , ly=', y
          write(6,*) '!!!'
       end if
      goto 20
   50 continue
     close(10)
     write (6,*) '*** Geometry information read from file anb.txt.'
     write (6,*) '***'
     goto 999
  900 write (6,*) '!!! Error reading file anb1.obs'
     write (6,*) ' '
      goto 990
  990 error=.true.
  999 continue
     error=.false.
     return
      end

      subroutine init_density(lx,ly,density,init_g,node,g)
      implicit none
      integer  lx,ly
      real*8  density,init_g,node(0:8,lx,ly),g(0:8,lx,ly)
      integer  x,y
      real*8  t_0_d,t_1_d,t_2_d,t_0_g,t_1_g,t_2_g
      t_0_d=density* 4.d0/9.d0
      t_1_d=density/ 9.d0
      t_2_d=density/36.d0
	  t_0_g=0.d0
	  t_1_g=init_g/6.d0
	  t_2_g=init_g/12.d0
      do 10 x=1, lx
        do 10 y=1, ly
          node(0,x,y)=t_0_d
          node(1,x,y)=t_1_d
          node(2,x,y)=t_1_d
          node(3,x,y)=t_1_d
          node(4,x,y)=t_1_d
          node(5,x,y)=t_2_d
          node(6,x,y)=t_2_d
          node(7,x,y)=t_2_d
          node(8,x,y)=t_2_d
 
          g(0,x,y)=t_0_g
          g(1,x,y)=t_1_g
          g(2,x,y)=t_1_g
          g(3,x,y)=t_1_g
          g(4,x,y)=t_1_g
          g(5,x,y)=t_2_g
          g(6,x,y)=t_2_g
          g(7,x,y)=t_2_g
          g(8,x,y)=t_2_g         
   10 continue
      return
      end

	  subroutine init_u_hlp(lx,ly,u_hlp_x,u_hlp_y,u_hlp_c,u_hlp_squ)
	  implicit none
	  integer  lx,ly
	  integer  x,y,i
	  real*8   u_hlp_x(lx,ly),u_hlp_y(lx,ly),u_hlp_c(0:8,lx,ly),u_hlp_squ(lx,ly)
	  do 102 x=1,lx
	    do 102 y=1,ly
		  u_hlp_x(x,y)=0.d0
		  u_hlp_y(x,y)=0.d0
		  u_hlp_squ(x,y)=0.d0
		  do 102 i=0,8
		  u_hlp_c(i,x,y)=0.d0
	102 continue
	  return
	  end

      subroutine check_density(lx,ly,node,time)  
      implicit none
      integer  lx,ly,time
      real*8  node(0:8,lx,ly)
      integer  x,y,n
      real*8 n_sum
      n_sum=0.d0
        do 10 y=1, ly
          do 10 x=1, lx
            do 10 n=0, 8
   10         n_sum=n_sum+node(n,x,y)
      write(6,*) '*** Iteration number=', time
      write(6,*) '*** Integral density=', n_sum
      write(6,*) '***'
      return
      end

      subroutine redistribute(lx,ly,obst,node,g,accel_d,density,init_g)
      implicit none
      integer  lx,ly
      logical  obst(lx,ly)
      real*8   node(0:8,lx,ly),g(0:8,lx,ly),accel_d,density,init_g
      integer  y,i,j
      real*8  t_1,t_2,g_1,g_2,ggg
      t_1=density*accel_d/ 9.d0
      t_2=density*accel_d/36.d0
      ggg=5.d0
	  g_1=ggg/6.d0
	  g_2=ggg/12.d0
   do 10 y=1, ly
    if(.not.obst(1,y).and.node(3,1,y)-t_1>=0..and.node(6,1,y)-t_2 >= 0..and.node(7,1,y)-t_2 >= 0.) then
          node(1,1,y)=node(1,1,y)+t_1
          node(3,1,y)=node(3,1,y)-t_1
          node(5,1,y)=node(5,1,y)+t_2
          node(6,1,y)=node(6,1,y)-t_2
          node(7,1,y)=node(7,1,y)-t_2
          node(8,1,y)=node(8,1,y)+t_2
     endif
	! if(.not.obst(1,y))then
  !        do 101 i=1,4
!		     g(i,1,y)=g_1
!		101  continue
!		  do 102 j=5,8
!		     g(j,1,y)=g_2
!		102  continue
     if(.not.obst(1,y))then      !.and.g(3,1,y)-g_1>=0..and.g(6,1,y)-g_2 >= 0..and.g(7,1,y)-g_2 >= 0.)
          g(1,1,y)=g(1,1,y)+g_1
          g(5,1,y)=g(5,1,y)+g_2
          g(8,1,y)=g(8,1,y)+g_2
     end if
   10 continue	   
      return
      end

      subroutine bounceback(lx,ly,obst,node,n_hlp)
      implicit none
      integer  lx,ly
      logical  obst(lx,ly)
      real*8   node(0:8,lx,ly),n_hlp(0:8,lx,ly)
      integer  x,y
      do 101 x=1, lx
        do 101 y=1, ly
          if (obst(x,y))then
            node(1,x,y)=n_hlp(3,x,y)
            node(2,x,y)=n_hlp(4,x,y)
            node(3,x,y)=n_hlp(1,x,y)
            node(4,x,y)=n_hlp(2,x,y)
            node(5,x,y)=n_hlp(7,x,y)
            node(6,x,y)=n_hlp(8,x,y)
            node(7,x,y)=n_hlp(5,x,y)
            node(8,x,y)=n_hlp(6,x,y)
          end if
   101 continue
      return
      end

	  subroutine thermal_bc(lx,ly,obst,g,g_wall,g_hlp)
      implicit none
      integer  lx,ly
      logical  obst(lx,ly)
      real*8   g(0:8,lx,ly),g_wall,g_1,g_2,g_hlp(0:8,lx,ly)
      integer  x,y
	  g_1=g_wall/6.d0
      g_2=g_wall/12.d0
      do 101 x=1, lx
        do 101 y=1, ly