Velocity boundary conditions & mass conservation


I’ve just finished the implementation of the velocity boundary conditions (regularized method) for a 3D Poiseuille flow (D3Q15). The results are very nice. However, the total mass in my system increases a bit. At the beginning I have exactly 64000, but after 10000 iterations the mass is ca. 64100, i. e. I observe an increase of 0.16%.
Now one could argue that due to the velocity field the fluid is slightly compressed (my maximum Mach number is 0.01), which would explain the gain of mass.
After 4000 iterations the system is in a steady state, which means that the mass should remain constant from thereon. This is however not the case. The gain rate seems to become constant after those 4000 iterations.
Now I have two questions:

  1. Is this rate of mass increase a problem? Or is the finite accuracy of the computer responsible? The relaxation parameter is around 0.85.
  2. If it is a problem, how can I find out what the source is?

Thanks for your explanations!

Hi Timm,
the mass non conservation in your flow is a common “feature” of the usual boundary conditions and is not usually a problem (it becomes problematic only if the densities become huge).

Since we are solving the incompressible Navier-Stokes equations you don’t care about the absolute value of the density, but rather about the pressure gradient (which is proportional to the density gradient) and therefore what you really want is that the density gradients remain correct. To check if the everything is allright you should run accuracy tests and see if you recover the “second order” accuracy of the LBM. And of course check if the pressure gradient is correct (and remains correct) when you have reached the steady state.

Hello Orestis,

after some thinking it has become clear why there is a mass increase. Let us assume that we have a straight rectangular channel with constant velocity at the inlet and at the outlet. Both boundary velocities shall be the same, u_B (where the velocity is perpendicular to the inlet / outlet plane). Then the density at the inlet is proportional to 1 / (1 - u_B) and the density at the outlet to 1 / (1 + u_B). Since the velocities at both boundaries are the same, but the densities are not, there is a rough estimate for the mass increase:
dM = [1 / (1 - u_B) - 1 / (1 + u_B)] * u_B * LY * LZ * #t
LY and LZ are the lattice dimensions of the channel cross section and #t is the number of the iterations. I have found that this formula reproduces the actual mass increase by a constant factor. My above estimate is larger than the simulated increase by a factor of 4.6. Since the total mass of the system scales with LX * LY * LZ, the mass increase scales with 1 / LX. A finer mesh therefore reduces the rate of mass increase. On the other hand a reduction of the Mach number reduces this rate, since the density difference at inlet and outlet is smaller than. Obviously, the mass increase should linearly scale with the simulation time.
Isn’t it now possible to modify the inlet / outlet velocities in such a way that the mass flow at the inlet is equal to the mass flow at the outlet?
u_in / (1 - u_in) = u_out / (1 + u_out)
Or is this a bad approach?

Thank you,

I think that you are perfectly right. And that your approach for correcting the mass loss could be fine. But just for curiosity, why is it so important to you to have exact mass conservation?

Hi Orestis,

actually I am not really interested in exact mass conservation. If I keep the Mach number at around 0.01 the mass increase is less than 1 % after the complete simulation. The velocity error is in the same range. This way my simulation is much better than that of the corresponding experiment which I want to simulate.
I was just interested in a method to keep the mass under control. But I think that its exact conservation is really of no relevance in my case.


Hello Tim,

Could you give me your program source code. I would be easier to study the existing code instead of from zero. My e-mail is
Thank you

Hi Noraz,

my code is a specialized one using the message passing interface and written in c++. If you are looking for a simple code for starting with LB simulations, my code is not well suited for this purpose. Maybe it would be better to have a look at the openLB code first.


Hello Timm,
this will not really answer your question but you could read the following paper with a “regularized mass conserving BC”.

Enhanced, mass-conserving closure scheme for lattice Boltzmann equation hydrodynamics. A, Hollis, I, Halliday, C. M. Care, J. Phys. A: Math. Gen. 29, (2006), 10589-10601.

You should notice that to our opinion (we discussed it with Jonas) this BC suits very well for no slip walls but not for in/outlets, because the mass conservation in this paper is enforced by saying that all the mass leaving the domain on a BC node must reenter on the same BC node. Which to our opinion is not physical. But this BC should be really fine for no-slip walls.


Hi Orestis,

thanks for the hint. I have downloaded the paper and will read it today. By the way, you have a PM.



There is another issue connected to the mass increase using velocity boundary conditions.
Orestis said that the mass increase is not problematic as long as the pressure gradient is correct. The absolute value of the density does not influence the velocity, since

u_i = \sum_k f_k c_ki / \sum_k f_k

and the amplitudes of the f_k have no effect here.
I have run a simulation with a large Mach number (Ma = 0.23). As I have written in an earlier post in this thread, the mass increase grows very fast if the Mach number is large. I then have compared the simulated velocity profile with the theoretical profile, and the mean relative deviation is 1 %. This means that the heavy mass increase of more than 36 % does not really influence the velocity.
But on the other hand, I have calculated the shear stress profile which obeys

\sigma_{ij} \propto \sum_k f_k^neq c_ki c_kj

Here, I do not divide by \sum_k f_k which means that the mass increase directly influences the values of sigma. In fact, my result is that the mass increase of 36 % leads to an error of around 35 % in sigma. This cannot be an coincidence.
Of course this effect can be compensated. But I wonder what is the correct way.

  1. Divide all values of sigma by the mean density of my entire simulation
  2. Divide all values of sigma by the local density
    Both approaches will result in a corrected expression for sigma, but which one is more accurate? I tend to take the second approach.

Thank again,

To put the question into another form:
Is it good to rescale my distribution functions f_i at every time step like this?

f_i \to f_i / \bar \rho

where \bar \rho is the mean density of my entire lattice. I wonder if this would influence my second order accuracy.
Does anybody know an answer?