Hi,
In LBM algorithm, the calculation of macrovariables is placed before the collision/streaming step.
I have tried replacing the macrovariables at the last. But i end up with an error…
May anyone can help why cant macrovariables comes at the last?

It is mandatory to use this order of routines if you really want to get reliable values.

the normal order is:
initialization (vel,dens,feq)
boundary treatment (bounce-backs and most of open boundaries)
compute macroscopic properties (velocity and density)
compute f_eq (with uptaded velocity and density)
collision
streaming, and back to boundary treatment

After streaming step, the macroscopic values need to be recalculated because all the nodes get different f’s values from the neighbour nodes.

So if you want to have the correct velocity and density, you need to compute before them in order to get a correct fEq, in order to have a correct collision operator, in order to propagate again the f’s, and so on…

Boundary methods such as Zou and He are “on node” methods, meaning that the LB stream and collide algorithm is applied everywhere. So if you have a domain that runs from i=1 to i=nx in the horizontal direction and j=1 to j=ny in the vertical direction then you apply the collision operator to all node {1,nx)x{1,ny}. The method must find the unknown distribution functions at a boundary. That is, the distribution functions at the wall moving into the flow (eg f2, f5 and f6 at a south wall, using the standard labelling of the lattice directions). Your algorithm is probably doing the collision step first and then streaming. After the streaming you don’t know the functions that are pointing into the domain, so you will apply the boundary conditions here (eg, at j=1 for a south wall after streaming). At the next time step you perform collisions again, and you must do this everywhere. If you don’t do collisions at j=1 then your LBE won’t be able to reproduce the Navier-Stokes equations near the wall. One way to see this is to consider the Chapman-Enkog expansion: no collisions means no viscous diffusion and you won’t have the velocity gradients at the wall correct (which are present in the non equilibrium parts of the f_i). It will make no difference if you do further collisions to the particles that have streamed to j=0, for example, because these are outside of the domain.

The reason you are not seeing any discernible difference in Poiseuille flow is either because you are applying boundaries at j=1 and doing a further collision at j=0 or, if you are applying the boundary conditions at j=0 and then deciding not to collide here, then the flow is just too simple (or you’re not looking in the right place). Mass and momentum are conserved by collisions so it shouldn’t matter if you calculate velocity with pre or post collisional distributions (and hence why your velocities are looking so similar). The difference will be in the stress field and the higher order moments. The stress tangential to the wall should be zero (because the velocity is constant). You can check this, if you wish (it can be found from the non equilibrium part of the momentum flux, ie Pi_xx^(non-eq) at a horizontal no slip boundary), although I imagine it will be ok because of the no slip condition. It might be more enlightening to look at the shear stress (non equilibrium part of Pi_xy) at the wall so see if this is correct (note that you have an analytic solution for this). There is also the stress normal to the wall. This is zero is Poiseuile flow because the velocity is one dimensional but I imagine it will not be the same for wall-collisions and not-wall collisions in other, more complicated flows with a 2D velocity field that depends on both x and y.

Another thing to note is that, for Poiseuille flow, the LBE should give you the exact solution on a very coarse mesh. Since there is no x-dependence you can write your code in 1D (for force driven periodic flow) and use only three grid points in the y-direction (that is, j=1,2,3, with no slip conditions applied at j=1 and j=3). You can check that both wall collsions and not wall collisions cases do this. Also, you must make sure that your variables are define correctly. If you have a flow in which you have a characteristic length L=1 (eg, the height of your channel is 1 in the non dimensional system) and you have N points with boundaries at j=1 and j=N then your grid spacing is dx=1/(N-1). If you study a richer flow, such as lid driven cavity for example, I imagine you’ll notice a decrease in both accuracy and rate of convergence.