I wrote a D2Q9 fortran code for horizontal poiseuille flow that not include any force term. Now I want to develop it for vertical flow and I have to implement force term in my code. Does anybody know the differences between these situations? What should I change in my code? It would be greatly appreciated if anybody could help me.

Presumably you have implemented boundary conditions to apply Poiseuille (pressure-driven) flow by specifying densities at the inlet and outlet? If so, you should not need to include this type of boundary if you are applying a constant body force instead: in this situation, you can just apply periodic boundary conditions at the inlet and outlet.

In any case, there are a few methods by which you can apply forces to your lattice Boltzmann fluid.

One method is to simply modify the velocity used when calculating local equilibrium distribution functions for the collision process: Martys and Chen (Phys Rev E, 53, 743-750 (1996)) suggest adding tau*F/rho to the fluid velocity (where tau is the relaxation time, rho is the fluid density and F is the force). Apart from adding this term to the velocity, no other modification would be needed.

Another method is to devise an additional forcing term for each lattice link that can be added to the collision operator. One example of this is Kupershtokh and Medvedev’s Equal Difference Method (EDM) scheme (J Electrostatics, 64, 581-585 (2006)), which defines the forcing term as a difference in local equilibrium distribution functions between the fluid at a velocity increased by F/rho and the fluid at its measured velocity.

There are also forcing methods that make use of both a velocity correction and an additional forcing term. One well-established method is given by Guo and co-workers (Phys Rev E, 65, 046308 (2002)), which corrects the velocity used in local equilibrium distribution functions and applies a forcing term that depends on the relaxation time (tau).

These forcing methods vary according to computational efficiency and accuracy: the Martys/Chen approach is comparatively quick but may be inaccurate for larger forces, while Guo’s method is very accurate but also more computationally intensive. How you code them up will depend on your pre-existing code, e.g. EDM could be coded easily by calling a subroutine to calculate the local equilibrium distribution function twice - once each for the adjusted and original velocities - and finding the difference between the resulting values. They are all comparatively straightforward, however, so it shouldn’t be difficult to try each method and see which work the best for your system.

Thanks for your helpful reply. I think my problem is in boundary conditions as you mentioned. Now, may I ask you what should I choose for initial values of rho in gravitational Poiseuille flow? Should I assign a constant value, like rho0, for all the lattice nodes?

I am trying to code the exact difference scheme in Matlab…
Do you have any code avaiable for multiphase system which uses pesudo Forcing term…?
I would like to have an idea of how to code Multiphase system in LBM