Body force implementation without lattice effects

Lattice Boltzmann Method LBM: Theory
#1

Dear all

According to Guo et al. (2002) paper: “Discrete lattice effects on the forcing term in the lattice Boltzmann method” the correct form of the body force (to be introduced in the collision operation) is:

F(i)=(1-1/(2*tau))w(i)[(e(i,j)-u(k))/cs^2+(e(i,k)*u(k))*e(i,j)/cs.^4]*F(j)

where tau is the relaxation time, w(i) lattice weights, e(i,j) lattice speeds, u(k) fluid speed, cs sound speed, F(j) physical force to be implemented and F(i) lattice force.

For the D2Q9 model the indexes are:

i=1,2,3,…9
j=1,2
k=1,2

Could anyone please tell me how to expand the right hand side of the F(i) equation?

According to my index notation knowledge (from college times) the result should be, if I did not make any mistake writing it:

F(i)=(1-1/(2*tau))w(i)([(e(i,1)-u(1))/cs.^2+(e(i,1)*u(1)+e(i,2)*u(2))*e(i,1)/cs^4]*F(1)+[(e(i,2)-u(2))/cs.^2+(e(i,1)*u(1)+e(i,2)*u(2))*e(i,2)/cs^4]*F(2)

However, after several computations of flow benchmark cases I believe that I did some mistake here. Don´t know what, but I think I did!!

If someone have implemented Guo’s body force model and/or find any mistake in my result please tell me. I am really in a dead end here!

Thanks in advance.

Goncalo

#2

Hi Goncalo,

I had a similar problem. To make things easier I have written a class for physical vectors including operators for scalar product, cross product etc. This way, one can easily code expressions like the Guo forcing term. I have tested the computing speed, and it is not slower as if you would only use double or integer values. The big advantage is that the code is very simple and errors are minimized.
But coming back to your question: There is a mistake in your first equation, it should read


F(i) = (1 - 1 / (2 * tau)) * w(i) *
[(e(i,j) - u(j)) / cs^2 + (e(i,k) * u(k)) * e(i,j) / cs^4]
*F(j)

The expanded code, however, seems to be correct:


F(i) = (1 - 1 / (2 * tau)) * w(i) *
([(e(i,1) - u(1)) / cs^2 + (e(i,1) * u(1) + e(i,2) * u(2)) * e(i,1) / cs^4] * F(1) + 
[(e(i,2) - u(2)) / cs^2 + (e(i,1) * u(1) + e(i,2) * u(2)) * e(i,2) / cs^4] * F(2)

Why do you think that you have a flaw in your code? Have you checked all parentheses?

Timm

#3

Hello Tim

Thank you for helping me once more :slight_smile:

You say that my indexes are mistaken and that instead of k it should be j. For me the indexes that you said make more sense than those I wrote however, reporting to Guo’s article, the way to find a proper expression for F(i) is based on a power series expansion of the kind :

F(i)=w(i)*[A+B(j)*e(i,j)/cs^2+C(j,k)e(i,j)e(i,k)/(2cs^4)-C(j,k)/(2cs^2)]

Using the Chapman-Enskog expansion one finds that A=0, B(j)=(1-1/(2tau))F(j) and C(j,k)=(1-1/(2tau))(u(j)*F(k)+u(k)*F(j)).
Gathering all terms one has:

F(i)=(1-1/(2*tau))w(i)[F(j)*e(i,j)/cs^2+(u(j)*F(k)+u(k)*F(j))*e(i,j)e(i,k)/(2cs^4)-(u(j)*F(k)+u(k)F(j))/(2cs^2)]

Since u(j)*F(k)+u(k)F(j) is symmetric one can write it as: 2u(k)F(j) (or as 2u(j)*F(k) but the former is more useful), so:

F(i)=(1-1/(2*tau))w(i)[F(j)e(i,j)/cs^2+(2u(k)*F(j))e(i,j)e(i,k)/(2cs^4)-(2u(k)F(j))/(2cs^2)]

which yields, putting in evidence F(j):

F(i)=(1-1/(2*tau))w(i)[(e(i,j)-u(k))/cs^2+(u(k)*e(i,k))*e(i,j)/(cs^4)]*F(j)

In practice however you and I use the same expression so the fact that it should be a k or a j is not problematic, regardless results are logical.

The benchmark flow test that I performed was as follows:

Using He and Zou boundary conditions (Pressure inlet and oulet for inlet and outlet boundaries) and (Zero velocity for walls) I simulated a Poiseuille flow of the form:

niud2u/dy2=dp/dx-rhog(x) which yields a velocity profile as: u(y)=-1/(2niu)[dp/dx-rho*g(x)]y(y-LY); where L is the channel width.

So then I defined g(x)=A*(x/LX)^n, where n is the order of the polynomial, LX is the channel length and A is a constant easily determined from a momentum balance.

I used n=0,1,2 and 3 and compared the result provided by Guo’s model against F(i)=w(i)*e(i,j)*F(j)/cs^2. The result is that Guo’s model is better only in the case of n=0, i.e. the acceleration is constant and the force only changes because rho changes.

Then I thought: “Since Guo’s model uses a different definition of the equilibrium and fluid velocity, i.e. v(j)=Sum(f(i)*e(i,j))/rho+g(j)/2, maybe I am messing up here”.

So I used Luo’s definition of the body force which is almost equal to Guo’s i.e.: F(i)=w(i)[(e(i,j)-u(k))/cs^2+(u(k)*e(i,k))*e(i,j)/(cs^4)]*F(j) but the velocity is simply given by v(j)=Sum(f(i)*e(i,j))/rho.
From the Chapman-Enskog expansion one can see (for instance in Guo’s paper) that for tau=1 Luo’s body force and F(i)=w(i)*e(i,j)*F(j)/cs^2 body force yield exactly the same macroscopic equations. So if I compare them both against a Navier-Stokes analytical solution the error should be the same. I did that and the error is not the same :frowning:

That lead me to the conclusion that [(e(i,j)-u(k))/cs^2+(u(k)*e(i,k))*e(i,j)/(cs^4)]*F(j) was not well written.

Did you checked your code with the Guo body force model against other models? Did it always yield better results?
By the way what boundary conditions did you use? I believe that because of the different velocity definition which incorporates the force, the bounceback boundary conditions cannot be used in Guo’s model!!

Thanks for your help Tim.

Regards

Goncalo

#4

Hello Goncalo,

I have not tested the Guo forcing term in every detail, but for a body force driven 3D Poiseuille flow (constant force), the errors were very similar to the results from the conventional forcing term, and I strongly believe that my implementation is working properly.
I assume that your implementation is not correct. Can you post your forcing term code or send it to me as a pm?
At the walls I use bounceback. And you think that this is not working when using Guo’s model? Very interesting! We have to discuss this further. I will come back to this point when I have more time.

Timm

#5

Hello,

I haven’t checked exactly what you have written. But for the steady Poiseuille case the body force is a constant (independent on space and time). Therefore there should be small differences between Guo forcing and the “standard” way to impose it . These differences are due to the fact that you are approaching the steady state with a different numerical scheme. But at the end you should get very close results.

What are these differences you are talking about? Are they really important?

#6

Dear Timm and Malaspin

First of all, thank you for your interest in my problem.

Reporting to my problem. It consists of a fully developed flow driven simultaneously by an inlet and outlet pressure boundary conditions and a body force.
Here it is how I conceived the problem:


niu*d2u/dy2=dp/dx-rho*g(x)

where I allowed the body force, g, to change along the streamwise direction x.

This yields a parabolic velocity profile of the form:


u(y)=-1/(2*niu)*[dp/dx-rho*g(x)]*y*(y-LY)

LY= channel width

I relate dp/dx and rho*g(x) so that in the end this two terms compensate each other and I will obtain a fully developed flow exclusively dependent on my Reynolds number.

Now to model the rho*g(x) term I used three different models:

  1. Guo’s model

F(i)=(1-1/(2*tau))*w(i)*[(e(i,j)-u(k))/cs^2+(e(i,k)*u(k))*e(i,j)/cs.^4]*rho*g(j)
  1. Luo’s model

F(i)=w(i)*[(e(i,j)-u(k))/cs^2+(e(i,k)*u(k))*e(i,j)/cs.^4]*rho*g(j)
  1. Traditional model

F(i)=w(i)*e(i,j)*rho*g(j)/cs.^2

Note: If using tau=1 the macroscopic equations of model 2) and 3) are the same.

The results that I got for a D2Q9 model with He and Zou boundary conditions in all boundaries are:

  1. Using the incompressible model of Luo and g(j)=const the three models yield machine accuracy. Perfect!

  2. Using the incompressible model of Luo and g(j)=const*x the Guo’s model yield the better accuracy, which is logical. However, Luo’s and traditional model yield different accuracies. Luo’s accuracy was sightly better. Not so perfect!

  3. Using the standard model and g(j)=const the Guo’s model yields better accuracy and Luo’s and traditional model yield exactly the same accuracy. This is as I would expect… Perfect!

  4. Using the standard model and g(j)=const.*x the Guo’s model yields the worst accuracy followed by Luo’s model. The traditional model yield the better accuracy by far. …Totally wrong!

Below are data quantifying my point 3) and 4) respectively:

Constant g(j)
Re=10 tau=1
Guo model Luo model Traditional model
Ma L2 Ma L2 Ma L2
0.1 0.010257 0.1 0.010271 0.1 0.010271
0.09 0.0082138 0.09 0.0082225 0.09 0.0082225
0.07 0.0050043 0.07 0.0050075 0.07 0.0050075
0.05 0.002608 0.05 0.0026089 0.05 0.0026089

Linearly changing g(j)
Re=10 tau=1
Guo model Luo model Traditional model
Ma L2 Ma L2 Ma L2
0.1 0.016795 0.1 0.016282 0.1 0.0051127
0.09 0.013415 0.09 0.01299 0.09 0.0042043
0.07 0.0082803 0.07 0.0080116 0.07 0.0024848
0.05 0.0043711 0.05 0.004226 0.05 0.0012543

Isn’t this strange?

Tim:
Regarding what I said about the inaccuracy of using bounceback boundaries when u(j) is defined as it is in Guo’s model I may be incorrect.
I said that based on a shallow interpretation of Ladd and Verbeg paper: Lattice Boltzmann simulations of particle-fluid suspensions. 2001. Journal of Statistical Physics.

The authors on pag. 1221, when talking about bounceback boundary conditions, say and I quote:

“We note that boundary conditions that have been tuned to give u=0 at the solid-fluid interface are no longer exact if we take the velocity field: rhou’=rhou+1/2fdt”

However, in this sentence they reference the paper: Analytical Solutions of simple flows and anaysis of nonslip boundary conditions for the lattice boltzmann bgk model. 1997 by He, Zou, Luo and Dembo.

From reading that paper the boundary conditions that are tuned to give u=0 are those from Noble and Inamuro. Bounceback always yield a slip velocity at wall. So maybe rhou’=rhou+1/2fdt is inaccurate only for Noble’s and Inamuro’s boundary conditions, not bounceback,

Nevertheless, I used He and Zou boundary conditions and I had to change them a little bit to account for the new velocity definition. And this approach of mine is right! I am sure :slight_smile:

Does anyone ever tested several body forces in a benchmark flow situation similar to the one that I am exposing? Are my results strange?

Thanks for the help.

Regards

P.S. Timm in the spreadsheet that I sent you it is written constant and linear body force. Instead of body force I meant to say acceleration…

#7

Timm you are right regarding the indexes that I used. It is j and not k in the velocity of Guo’s body force expression. I did some mistake there!! I still don’t know where but I sure did!
The good news are that my code considers that index as being in fact j (because k was not logical in practice) so my problem remains… :confused:

#8

Dear Goncalo,

I am pretty new in LB simulation. Actually, I was trying to add forcing term to a LB code and confused a little bit.
By adding forcing term like you proposed recently (Physica A 390 (2011) 1085-1095), could we say that our approach is 2nd order accurate in space an time?

Thanks and regards,

Saeid

#9

Dear Saeid

First of all thank you for your interest in my work.

Yes, the body force expression discussed in that Physica A paper is 2nd order accurate in space, but not in time!
You should note that, the work you referred to focused on the recovering of the correct macroscopic body force behavior in a steady-state process.

Nevertheless, the extension of that work for both time-independent and time-dependent processes has been published recently:
http://journals.cambridge.org/action/displayAbstract?fromPage=online&aid=8539110&fulltextType=RA&fileId=S0022112012000833

The point you should be aware of is that both my studies apply to the incompressible version of the LBM (this point was implicit in the Physica A paper, but now is fully explicit in the JFM paper). The incompressible LBM recovers the Navier-Stokes equations in artificial compressibility form, which in the steady-state coincide with the incompressible Navier-Stokes eqs.

A reversed implication of that is that when the standard LBM is used, where the recovered macroscopic setup is of weakly-compressible type, the difference between the Physica A body force and that from Guo (to name the most well known case) is a compressible error. Therefore, if you intend to reproduce incompressible Navier-Stokes solutions in the standard LBM, it is a bit indifferent whether you use Physica A or Guo body force.
Perhaps, the best way to understand this is by start solving some of the benchmark flows I presented (or some other flows you aware of, which you find useful to test this problem) and see for yourself the differences. Or to reproduce Guo study (Discrete lattice effects on the forcing term in the LBM). By the way the body force model in the Physica A paper corresponds to method 3 in the Guo paper. One thing I think would be interesting for you to do is to reproduce the steady-state test presented in Guo paper. If you implement the problem properly (note Guo did some mistake in solving that flow) you will find that both method 3 and Guo body forces are exactly equal!

Said that, just one more point. In the diffusive scaling, i.e. the scaling you should adopt in order to reproduce incompressible hydrodynamics, the LBM is a 1st order time accurate method. (For a proof of this statement, please check, for instance, the paper from Junk et al. 2005, Asymptotic analysis of the lattice Boltzmann equation).

Regards,

Goncalo

#10

Dear Goncalo,

Thank you so much for your valuable comments,

Saeid

#11

Dear Goncalo,

I have some difficulties with indices. would you please help me.
I can figure out the expansion of the RHS of F(i) equation:

F(i) = (1 - 1 / (2 * tau)) * w(i) *[(e(i,j) - u(j)) / cs^2 + (e(i,k) * u(k)) * e(i,j) / cs^4]*F(j) (**)

But what I can not understand clearly is one step before.

Refering to Guo et al. 2002, Eq(20) is like:
F(i) = (1-1/(2*tau)) * w(i) * [(e(i)-v)/cs^2 + (e(i) * v) * e(i) / cs^4] * F

If I want to expand this I do like this:
F(i) = (1 - 1 / (2 * tau)) * w(i) *[(e(i,j) - u(j)) / cs^2 + (e(i,j) * u(j)) * e(i,j) / cs^4]*F(j)

How did you get eq(**) form this? I confused j and k indices.

Thanks and regards,
saeid

#12

GoncaloSilva
Hello!
Can you send me your paper by email? I meet over the problem about the pressure boundary condtion applied on the inlet and outlet. I use the period boundary condition on the inlet and outlet and the standard bounceback method on the top and bottom, and gain good solutions compared with the analytical solutions, but can not gain bad good solutions with the pressure boundary condtion of Guozhaoli’s method, I don’t know where was the problem. i wish your help.thank you!

#13

Dear Saeid,

You did some confusion with the indices. Guo body force model should read:
F(i) = (1-1/(2*tau)) * w(i) * [(e(i)-v)/cs^2 + (e(i) * v) * e(j) / cs^4] * F --------> One of the particle velocity e(i) should read e(j). Below, is the collision step routine with Guo body force. Although not strictly necessary, in order to avoid screwing up, I split it in 3 steps.

for y=1:LY
    for x=1:LX
        for a=1:9

            rho_col=rho(x,y);                
            u_col=u(x,y);
            v_col=v(x,y);

            F_bracket_1=3.*(g(1,x,y).*c(1,a)+g(2,x,y).*c(2,a));
            F_bracket_2=9.*(u_col.*c(1,a)+v_col.*c(2,a)).*...
                (g(1,x,y).*c(1,a)+g(2,x,y).*c(2,a));
            F_bracket_3=3.*(u_col+v_col).*(g(1,x,y)+g(2,x,y));

            F=(1-1/(2*tau)).*w(a).*rho_col.*...
                (F_bracket_1+F_bracket_2-F_bracket_3);

            fprop(x,y,a)=(1-(1/tau))*f(x,y,a)+(1/tau)*feq(x,y,a)+F;
        end
    end
end

Just one (provocative) suggestion. After putting your code running correctly. It would be interesting for you to compare the accuracy of your results with F_bracket_1, F_bracket_2 and F_bracket_3, i.e. the Guo model, against the forcing model with F_bracket_1 only. Please note that I don’t known what are you modeling but I would say that, if using the incompressible LBM, there is a 50% probability that the terms F_bracket_2 and F_bracket_3 in the forcing model, i.e. Guo model, decrease your solution accuracy. On the other hand, if you use the standard LBM the inclusion of the terms F_bracket_2 and F_bracket_3 are irrelevant because they act at the level of the compressibility errors…For the incompressible case please see:

http://journals.cambridge.org/action/displayAbstract?fromPage=online&aid=8539110&fulltextType=RA&fileId=S0022112012000833

Dear Ydt,

Above there is a link for the paper. If you are not able to access the journal please send a PM with your email and I will send it to you. Yet, if I understood well what you are asking me, I guess that you are in the wrong topic. Here, we are talking about Guo’s forcing model, not Guo’s boundary condition! Regarding this latter I usually adopt ZouHe method to set pressure boundary conditions, so I am not familiar with Guo’s pressure boundary condition, sorry!

Regards,

Goncalo

#14

Dear Goncalo,

Thanks for your nice elaboration. Actually I am studying petroleum engineering and that is why my questions may look silly a little bit. I am using incompressible LBM for simulating fluid flow in porous media.
Would you please tell me what do you mean by standard LBMs. If you mean compressible LBM by that, the inclusion of the terms F_bracket_2 and F_bracket_3 should be relevant because as you said they act at the level of the compressibility errors.

Regards,

saeid

#15

Dear Goncalo:
Thank you very much! My emailt is tywoyangda@126.com. I used Zou/He’s pressure bounday conditions, but I felt it deal with the corner nodes, which is difficult.

#16

Dear Saeidnorouzi,

Would you please tell me what do you mean by standard LBMs.

Roughly speakly, the incompressible LBM recovers the incompressible Navier-Stokes eqs. in artificial compressibility form (up to the second order, to be correct). What I mean by standard LBM, is the LBM where the equilibrium involves \rho and \rho{u}, i.e. mass and momentum density. The transport equations governing these conserved fields are the compressible Navier-Stokes. Generally, since the system runs in the weakly-compressible regime (because of the equilibrium second order approximation) it is said that solutions of this system approximate incompressible Navier-Stokes solutions (the terms responsible for the difference called compressibility errors).

If you mean compressible LBM by that, the inclusion of the terms F_bracket_2 and F_bracket_3 should be relevant because as you said they act at the level of the compressibility errors.
That is correct. Or at least partially. To be precise, the second order Chapman-Enskog expansion tells us that. Yet, we should not forget that what Chapman-Enskog tells us is of approximate character, because of its truncated nature. At this moment I am really really busy, because I urgently need to submit my PhD thesis as soon as possible…(it is a really tough situation), but whenever I have time I would like to come back to this topic and prove what happens exactly! I have already some drafts and suspicion on this…so, I guess, that will not take me too long…

Dear Ydt,

I see your point. Please check your email

Regards

Goncalo

#17

Dear Goncalo,

Thank you again for your time.
wish the best.

Saeid

#18

Hi Goncalo

Thanks for the useful info here. I have a quick question about body force implementation I was hoping you could resolve. In equation 1 of your Physica paper on the RHS we have delta_t*F. My question is this, is delta_t here the lattice delta_t or the physical delta_t?

Up until now I’ve assumed that its the lattice delta_t, but I’m having some issues with one of my models which has prompted this question

Regards,
Bruce

#19

Dear Bruce,

Thank you very much for your interest in my work.

The delta_t that you are referring to comes from the (second-order) discretization of the time/space variables in the discrete velocity Boltzmann equation to obtain the LBE :slight_smile:
As such, that delta_t in eq. (1) must have to be in lattice units.

Having said that and since you are experiencing some problems I might risk saying that your problem may be in the way you scale the physical force rather than the forcing model itself. (recall that to reproduce incompressible hydrodynamics everything will have to respect diffusive scaling.)

Best wishes,

Goncalo

#20

Thanks for the quick response Goncalo, I am fairly certain everything is correct (famous last words).

I’m working with tensor permeability of rock microstructures, and am seeing a discrepancy between permeability tensors computed using zou/he pressure bc’s and body force bc’s. I’m slowly coming to the conclusion that the discrepancy comes from the limitation of non-normal velocity = 0 with zou/he, and the lack of periodicity in the zou/he formulation… It is interesting to note that I have observed little difference when dealing with scalar permeability of simple idealized geometries.