Bounceback Conditions

Lattice Boltzmann Method LBM: Algorithms
#1

Hello,

I’ve been given an LBM code and I cannot understand what kind of bounceback condition is this:


         if (obst(x,y,z)) then
c
c...........rotate all ensities and write back to node

            node(1,x,y,z) = n_hlp(3,x,y,z)

            node(2,x,y,z) = n_hlp(4,x,y,z)

            node(3,x,y,z) = n_hlp(1,x,y,z)

            node(4,x,y,z) = n_hlp(2,x,y,z)

            node(5,x,y,z) = n_hlp(6,x,y,z)

            node(6,x,y,z) = n_hlp(5,x,y,z)

            node(7,x,y,z) = n_hlp(9,x,y,z)

            node(8,x,y,z) = n_hlp(10,x,y,z)

		  node(9,x,y,z) = n_hlp(7,x,y,z)

		  node(10,x,y,z) = n_hlp(8,x,y,z)

            node(11,x,y,z) = n_hlp(13,x,y,z)

            node(12,x,y,z) = n_hlp(14,x,y,z)

            node(13,x,y,z) = n_hlp(11,x,y,z)

            node(14,x,y,z) = n_hlp(12,x,y,z)

            node(15,x,y,z) = n_hlp(17,x,y,z)

            node(16,x,y,z) = n_hlp(18,x,y,z)

            node(17,x,y,z) = n_hlp(15,x,y,z)

		  node(18,x,y,z) = n_hlp(16,x,y,z)

I was reading boundary review and I couldn’t find it. I think it is not really well implemented, have you got any idea?

Thank you

Antonio

#2

This is bounceback. But I cannot tell you exactly which one. I assume it is fullway bounceback where you temporary store the populations on the wall nodes.

#3

Dear Timm,

thank you for the answer

Is there any realtion why the simulation after a bit that the mass decreases suddenly indefinitely increases (at least for 100,000 iterations?

Thank you

Antonio

#4

Antonio, could you please rewrite the sentence in a way I can understand it? :wink:
And could you give more information on the behavior of the mass?

#5

Dear Timm,

sorry for the strange answer…

Looking at picture: http://yfrog.com/1niterationsp you can see that the desntiy, defined as the sum of the node all over the domain, at the beginning goes down and then increases. I made simulations up to 100,000 iterations and the mass keeps rising with constant gradient.

Is there any realtion between this behaviour and the fullway bounceback?

Thank you

Antonio

#6

Antonio, bounceback is mass-conserving (if there is no bug in the code). Halfway bounceback is exactly mass-conserving, and fullway bounceback on average (because you temporarily store the populations on the walls, and there may be small fluctuations). The behavior I see in the picture seems to be not related to bounceback. What kind of flow do you have? What are the boundary conditions at inlet/outlet? How do you initialize your velocity and density?
Timm

#7

Dear Timm,

the code has not been made by me, that’s why I did not investigate too much.
However I’m simulationg a micro-reactor http://yfrog.com/5omicroreactorp which:

dimensions D_PHY=300micrometers; D_LB = 70 lu
viscosity nu_LB=0.0175 lu^2/ts
dx=4.286E-6[m/lu], dt=0.3367E-6[s/ts]
densit RHO_PHY=1002[kg/m^3], RHO_LB=1[mass/lu^3]
Reynolds 100

Inlet conditions are defined by imposing the velocity at the inlet, since it’s known the velocity, for each node of the inlet is given u_x=0.05 lu/ts


      t_0 = 1.d0 / 3.d0
      t_1 = 1.d0 / 18.d0
      t_2 = 1.d0 / 36.d0
      c_squ = 1.d0 / 3.d0	
        lxt = 1 
        if(time .eq. 1) lxt = lx         
         do 2 z = 1, lz 
	    do 2 y = 1, ly
		 do 2 x = 1, lxt
           d_loc = 0.d0
            do 20 i = 0, 18 
              d_loc = d_loc + node(i,x,y,z)
   20       continue
           u_x=0.05
		   u_y	 =  0.0 
		   u_z	 =  0.0
		  
		    u_n(1) =   u_x
            u_n(2) =         u_y
            u_n(3) = - u_x
            u_n(4) =       - u_y
            u_n(5) =   u_z 
            u_n(6) =       - u_z
            u_n(7) =   u_x + u_y
            u_n(8) = - u_x + u_y
		    u_n(9) = - u_x - u_y
            u_n(10) =   u_x - u_y
		    u_n(11) =   u_x - u_z
            u_n(12) = - u_x - u_z
		    u_n(13) = - u_x + u_z
            u_n(14) =   u_x + u_z
            u_n(15) =   u_z + u_y
            u_n(16) = - u_z + u_y
            u_n(17) = - u_z - u_y
            u_n(18) =   u_z - u_y

		  c_squ = 1.d0 / 3.d0 
		  u_squ = u_x * u_x + u_y * u_y + u_z * u_z

c...........equilibrium densities
c...........zero velocity density

            ro=d_loc

            node(0,x,y,z) = t_0  * ro *(1.d0 - u_squ / (2.d0 * c_squ))

c...........axis speeds (factor: t_1)

            do 3 i = 1,6 
            node(i,x,y,z) = t_1 * ro *(1.d0 + u_n(i) / c_squ 
     &               + u_n(i) ** 2.d0 / (2.d0 * c_squ ** 2.d0) 
     &               - u_squ / (2.d0 * c_squ))
c
    3       continue 

c...........diagonal speeds (factor: t_2)           

            do 4 i= 7, 18
            node(i,x,y,z) = t_2  * ro *(1.d0 + u_n(i) / c_squ 
     &               + u_n(i) ** 2.d0 / (2.d0 * c_squ ** 2.d0) 
     &               - u_squ / (2.d0 * c_squ))

    4       continue 
    2  continue

here above is explained the initialisation (all the domain at 0.05) and the inlet conditions. (I changed this a bit, since the mass flow in the inner tube is less, I initialise the inner tube to a slower velocity and all the rest at 0.05. So for a part of the geometry the velocity is not 0.05 but 0.04)

For the oulet is used a convective boundary condition:


      t_0 = 1.d0 / 3.d0
      t_1 = 1.d0 / 18.d0
      t_2 = 1.d0 / 36.d0
      c_squ = 1.d0 / 3.d0	
	x = lx

! u_curr: at boundary
! u_bdy1: at boundary - 1 (in x)
! u_prev: at boundary - 1 (in time)


c.....first compute the mean outflow velocity, Uc

	Uc = 0.0d0
	num = 0
	do 1 y = 1, ly
	 do 1 z = 1, lz
	  if (.not.obst(x,y,z)) then
		Uc = Uc + u_curr(y,z)
		num = num + 1
	  endif
    1	continue 	
	Uc = Uc/num !if (num.gt.0) check not needed


c.....compute the new velocities (based on convective BC)
	do 2 y = 1, ly 
	 do 2 z = 1, lz
	  if (.not.obst(x,y,z)) then 		
		u_curr(y,z) = (u_prev(y,z) + Uc*u_bdy1(y,z))/(1.d0+Uc)
		u_prev(y,z) = u_curr(y,z)		
		v_curr(y,z) = (v_prev(y,z) + Uc*v_bdy1(y,z))/(1.d0+Uc)
		v_prev(y,z) = v_curr(y,z)			
          w_curr(y,z) = (w_prev(y,z) + Uc*w_bdy1(y,z))/(1.d0+Uc)
		w_prev(y,z) = w_curr(y,z)
	  endif	  
    2	continue 
      
c
c.....now re-assign the densities at the outflow boundary
	do 3 y = 1, ly
	 do 3 z = 1, lz
	  if (.not.obst(x,y,z)) then
            d_loc = 0.d0
            do i = 0, 18 
              d_loc = d_loc + n_hlp(i,x,y,z)
		  enddo
c
             u_x = u_curr(y,z)
             u_y = v_curr(y,z)
             u_z = w_curr(y,z)
		
c
c...........square velocity
c
            u_squ = u_x * u_x + u_y * u_y + u_z * u_z
c
c...........n- velocity compnents (n = lattice node connection vectors)
c...........this is only necessary for clearence, and only 3 speeds would
c...........be necessary
c
		  u_n(1) =   u_x
            u_n(2) =         u_y
            u_n(3) = - u_x
            u_n(4) =       - u_y
            u_n(5) =   u_z 
            u_n(6) =       - u_z
            u_n(7) =   u_x + u_y
            u_n(8) = - u_x + u_y
		  u_n(9) = - u_x - u_y
            u_n(10) =   u_x - u_y
		  u_n(11) =   u_x - u_z
            u_n(12) = - u_x - u_z
		  u_n(13) = - u_x + u_z
            u_n(14) =   u_x + u_z
            u_n(15) =   u_z + u_y
            u_n(16) = - u_z + u_y
            u_n(17) = - u_z - u_y
            u_n(18) =   u_z - u_y
c

c...........equilibrium densities
c...........this can be rewritten to improve computational performance
c...........considerabely !
c
c...........zero velocity density
c
            n_equ(0) = t_0  * d_loc*(1.d0 - u_squ / (2.d0 * c_squ))

c
c...........axis speeds (factor: t_1)
c
            do 11 i = 1, 6
            n_equ(i) = t_1 * d_loc*(1.d0 + u_n(i) / c_squ 
     &               + u_n(i) ** 2.d0 / (2.d0 * c_squ ** 2.d0) 
     &               - u_squ / (2.d0 * c_squ))
 
   11       continue
c
c...........diagonal speeds (factor: t_2)           
c
            do 12 i = 7, 18 
            n_equ(i) = t_2  * d_loc*(1.d0 + u_n(i) / c_squ 
     &               + u_n(i) ** 2.d0 / (2.d0 * c_squ ** 2.d0) 
     &               - u_squ / (2.d0 * c_squ))
c
   12       continue        
c
           do i = 0, 18
              node(i,x,y,z) =n_equ(i) 
              !nod_conv(i,y,z)= n_equ(i)
            enddo
	   endif
    3	continue
c
	return
	end

That’s the inlet and outlet condition, for the rest in order go:

redistribute (initialisation and inlet), streaming, bounceback, relaxation, convective_BC.

Thanks for the time and the patience.

Do not hesitate to make more questions if something is not clear.

I thought was normal this behaviour of the LBM

Antonio

#8

Dear Antonio & Tomm

I read this text and I see:

“Looking at picture… the desntiy, defined as the sum of the node all over the domain, at the beginning goes down and then increases. … and the mass keeps rising with constant gradient.”

I have this problem in some of my works too.
did you find reason of that?
is it possible to explain for me?

Thank you

seyedmehdi