Density loss when adding force term

Lattice Boltzmann Method LBM: Theory
#1

Hi to all,

I am a beginner to LBM, trying to set up a simulation of a Poiseuille flow using a D2Q9 lattice and periodic boundary conditions at the in/outlet. The problem is now the following: Whenever I add any kind of force term, be it an external body force or a pressure gradient, my local density starts to decay to 0 exponentially with time. Has anyone yet experienced such (mis-)behaviour and maybe knows a way around it?

thanks in advance for any help,
Philippe

#2

Dear Philippe,

obviously, there is something strange happening. A few months ago I added a new forcing term, and the mass was increasing. However, there was a flaw in my class definitions, and I could correct it. The mass is conserved, if the sum of the forcing terms is zero, i.e. \sum_i F_i = 0, where i = 0, …, 8 for D2Q9 and F_i is the term you add to the right-hand-side of the lattice Boltzmann equation. For this reason, I assume that there is a mistake in the definition of your forcing term, or even in the lattice Boltzmann equation.
What kind of boundary conditions do you use? Is it possible that you have a leak, because the walls do not cover all of your fluid boundary?

Best,
Timm

#3

Dear Timm,

thanks for your reply! I checked my forcing terms, and they sum up to 0 nicely, so that can not be the cause. Also, I am using very simple boundary conditions. The wall just reflects the populations like f_{6,7,8}(i=0) = f_{2,3,4}(i=1) and similar for the lower wall, and the periodicity is handled like f_{1,2,8}(beginning) = f_{1,2,8}(end) and vice versa - might I eventually lose some population there? How would you handle simple periodic boundaries here?

Regards,
Philippe

#4

Hi,

ok, maybe you should recheck the equilibrium definition, the propagation and collision steps. Make also sure that the wall encloses the total volume, i.e. that there are no holes in your domain. However, I had holes in my simulations once, and the mass loss rate was quite small. Somehow I don’t think that this is the reason.
Try the following:

  1. Run a simulation without force or pressure gradient, with initially constant density and zero velocity. Everything should stay this way forever. If you observe time evolution, you know that the error is somewhere in the equilibrium, propagation or collision step, maybe also when you compute the periodicity of the flow.
  2. Run another simulation without force, but start with a given initial velocity. The system should then tend to zero velocity after some time (due to viscosity). But the mass has to be strictly conserved. If this works, then there will be an error in the force computations somewhere. Be sure that the body force is not too large (try 10[sup]-6[/sup] for example).
    By the way, are you in the small Mach number regime? Could you give some more information on your simulation parameters?
    Regarding the periodicity of the flow: I don’t know how your population indexes are defined, but you cannot do too many mistakes there.

Timm

#5

Hi,

ok, here is some more detail on my simulation parameters: I use a 50x3 cells grid (which doesn’t matter because of my periodic boundaries), and the body force applied is in the range of 10[sup]-3[/sup]-10[sup]-5[/sup] giving me a Mach number in that range as well.
I tried the things you suggested, and my results are that the mass stays constant for several thousand steps in both cases, with and without initial velocity. If I add an initial velocity, it decays very quickly to 0 (within 10 iterations or so), no matter what value it has, the value just scales the graph. Increasing viscosity gives me slower decay, which I think is correct that way. So my boundaries and collision steps should be alright I assume.
The error must be buried somewhere inside my force computation, I do it strictly following the book of Wolf-Gladrow, adding to each population f_i the scalar product of the force vector K and the corresponding lattice velocity c_i within the collision step, is there something I have to take care of?

Philippe

#6

I believe there is something wrong. A velocity decay in ten simulation steps, independent of the initial value, sounds very strange. But maybe you are in a very small Re-domain, when inertia is negligible. What are your Reynolds number? What is the value of your relaxation parameter? Test an initial velocity of 0.1.

#7

I use a relaxation parameter of 0.5 thus my Reynolds number (if I got the units right) should be around 2. I now tried out lower relaxation parameters, they slow down the velocity`s decay but it is still independent of initial values. Plus, I observed that as soon as I add an initial velocity, mass decay starts again, I lose around 0,3% of my mass during the first coupe of steps and then the mass stays there forever. A lower relaxation parameter also slows down this process and gives me slightly less mass loss.

#8

How is your relaxation parameter defined? Is it the number in the denominator or the numerator of your collision step. If it is in your denominator, then it is the relaxation time. And it cannot be smaller than 0.5. Even a value slightly above 0.5 can produce problems. Maybe this is the reason.

#9

It is the numerator, and unfortunately changing it just changes the results I get quantitatively, but they are messy anyway. At the moment, I lowered it to around 0.25 which produces slightly better results. To what degree should mass be conserved in LBM?
I think the problem still lies within the collision step, this is what I do for each cell after propagating:


do k=0,8
  eq = equilibrium(i,j,k)
  grid(i,j)%f(k) = grid(i,j)%f(k) + omega*(eq - grid(i,j)%f(k))
  grid(i,j)%f(k) = grid(i,j)%f(k) + dot_product(c(:,k),lbExtForce)/18.d0
enddo

where omega is the relaxation parameter, k denotes the velocity indices and lbExtForce is a (2D) vector holding the external force. I don’t see any mistake therein, maybe I missed something obvious…

#10

I now wrote a workaround that simply resets the density to 1 in each averaging step. That way, my results are correct but this can’t be the universal solution I think…

#11

Hi all,
I am a new comer in the world of Lattice boltzmann. I need a code written in fortran 90 for the channel flow so that I can understand the working of LBM.

Thanks

#12

Hi, Philippe,

the workaround might work (I use such a function as well, if there is mass increase due to velocity boundary conditions), but I would not recommend using it in general. There is a bug in your code, and you have to find it. Otherwise you can never be sure what else is affected by the bug. Your entire simulation may be garbage.
But I have also spent a few hours today to identify a bug in my code… that is our fate. :smiley:

Timm

#13

It seems like I found it - something in my equilibrium distribution function was wrong. Now I get nice parabolic profiles and mass conservation up to machine precision :slight_smile:

#14

Wonderful! :wink:

#15

Hi Philippe,

I don’t suppose you remember what the error was in the equilibrium dist. function? I’m having a similar problem; I’m trying to test my code with gravity-driven flow in a slit, like in Chap 5. of Sukop, and my density disappears over time – even though the parabolic profile is correct. Hmm…

#16

Today I had a similar problem again. I have rewritten my propagation step. Then I realized that my system was gaining mass, very slowly. I found out that one of the populations was propagated in a wrong way. Now, everything is fine again. My experience is that the most probable causes for mass violations are

  1. wrong definition of body force in collision step, first moment of lattice forces must be zero
  2. “open” boundaries, i.e. the fluid domain is not fully surrounded by obstacles and/or periodic walls
  3. wrong propagation step (finding errors here can be a pain in the ass)
  4. wrong equilibrium definition
#17

Thanks, Timm. I think I narrowed it down to my equilibrium distribution function, because when I supply these conditions

rho: 0.1
u_x: 0.05
u_y: 0

solve for associated feq[9], and recalculate macroscopic observables with this feq[9], I get

equilibrated
rho: 0.099875
u_x: 0.0500626
u_y: 0.

Shouldn’t the equilibrium function preserve rho and u_x, at least within some reasonable approximation? Sorry if this is a really elementary question, as I suspect it is; I just got started last week. I feel like u_x = 0.05 is well below the speed of sound, 1/sqrt(3)…

#18

So you have plugged rho, u_x and u_y into the equilibrium, computed all equilibrium populations, and then you have recomputed the velocity and the density from it? And you skipped collision, propagation etc., right?
Maybe there is also an error in the computation of the density and velocity. Which definition of the velocity do you use? Just post your equilibrium code here, if you like.

#19

Ah, I figured it out. It was a really stupid error; C++ thought my constants were type int, rather than double, causing error in my calculation of feq. Thanks, now it’s back to lurking…

#20

Dear Timm,

I am trying to introduce particle forces (particle suspension) to the LBM scheme via the forcing term. I use the Newton equation to compute force from fluid to particle and vice versa, i.e., two-way coupling. During the simulation I have observed that the density is slowly increased with time (about 1% after 50000 time steps). In the scheme I am using, the particle volume fraction is not incorporated into the LBE yet and I think, from a point of view of a beginner, that the density is increased because of the coupling of the particle forces via the forcing terms.

After I read your answer above that

“The mass is conserved, if the sum of the forcing terms is zero, i.e. \sum_i F_i = 0, where i = 0, …, 8 for D2Q9 and F_i is the term you add to the right-hand-side of the lattice Boltzmann equation”,

I am wondering if this also valid to my case. Could you please explain in more in detail, how could one make the sum of the forcing terms equal to zero ?

Regards,

Radompon