help on adding body force

Lattice Boltzmann Method LBM: Theory
#1

Dear all,

I’m new to this forum and indeed to LBM. However I’m trying to add body-force (gravity) to an existing code by just adding an additional term to velocity (as explained e.g. in Sukop and Thorne). I’ve carefully read the “howto” on units by J.Latt, but I still have some problems in “howto scale the gravity”. I.e. so far my code works quite satisfying, I even get a readsonable velocity profile for a sample gravity driven 2D poiseuille flow, but I’m struggling with units. I’d like to compare results to analytical solution, and I’m wondering how to get accelleration scaled, in order to get comparable results.

Please, could anyone send me a short code snipplet with some explanations on the necessary scaling/unit conversion, or just some brief explanation on that topic?

Kind Regards, Francois

(Sorry about my “novice level”, but since I have no other resource available for asking, I’m trying here…)

#2

Okay, here goes. What we want is a quantity defined in physical units converted to lattice units. Let me explain it by a few examples and hopefully you then get the principle. I use _p to denote physical units and _lb to denote lattice units.

Let’s say we have the quantity distance, which is defined for example in meters. Let’s take a 2D square box with sides of a size of 4 meters. We want to subdivide this into 100 by 100 cells. Then, we have dx_p = 0.04 meter. So in lattice units, we can now think of space as multiples of dx_p. That is, a space of 2 lattice units corresponds to a space of 2 * 0.04 = 0.08 meters in physical units. Likewise, a space of 2 meters corresponds to 2 / 0.04 = 50 lattice units of space. So, to convert space from physical to lattice units we divide by dx_p.

Although you are freely to choose your time step, it is recommended to choose it such that we have dt_p ~ dx_p^2 (check Jonas Latt’s thesis or units document for that). So, let’s choose dt_p = (4/100)^2 = 0.0016. If our simulation uses seconds as time scale, we define it as 0.0016 seconds. So in lattice units, we can now think of time as multiples of dt_p. That is, two time steps in lattice units corresponds to 2 * 0.0016 = 0.0032 seconds in physical units. Likewise, 2 seconds corresponds to 2 / 0.0016 = 1250 lattice units of time. So, to convert time from physical to lattice units we divide by dt_p.

Now let’s get down to business for our first real conversion. We have a velocity which has the value in physical units of 5 meters/second, or said otherwise 5/1 meters/second. Now, we divide the space by dx_p and the time by dt_p, we obtain 5/dx_p per 1/dt_p, i.e. 5/0.04 per 1/0.0016, i.e. (5/0.04) / (1/0.0016), i.e. 5 * (0.0016/0.04), i.e. 5 * (dt_p/dx_p). So, to convert velocity from physical units to lattice units we multiply the velocity physical units with (dt_p/dx_p).

Now let’s proceed to the gravity. Gravity has the value in physical units of 9.8 meters/second^2, or said otherwise 9.8/1 meters/second^2. Now, we divide the space by dx_p and the time squared by dt_p squared, we obtain 9.8/dx_p per 1/(dt_p^2), i.e. 9.8/0.04 per 1/0.0016^2, i.e. (9.8/0.04) / (1/0.0016^2), i.e. 9.8 * (0.0016^2/0.04), i.e. 9.8 * (dt_p^2/dx_p). So, to convert gravity from physical units to lattice units we multiply the gravity in physical units with (dt_p^2/dx_p).

Summarizing it:


dx_p        = 0.04
dt_p        = 0.0016
velocity_p  = 5
velocity_lb = 5 * (dt_p / dx_p)
gravity_p   = 9.8
gravity_lb  = 9.8 * (dt_p^2 / dx_p)

Hopefully you get it now

One small remark about implementing the gravity force. There exist various ways to implement a gravity force in LBM, but it is in general not recommended to adjust your velocity to implement the gravity force. Perhaps somebody else can explain this well? It is better to add the gravity term in your collision operator, check this thread for a simple way to do this: http://www.lbmethod.org/forum/read.php?3,781

1 Like
#3

Dear Bart,

thank you for your reply and detailed explanation! If I’m right, then for a simple 2D poiseuille flow I’d like to run I will have to do the following…


physical extension 
 x = 0.100m 
 y = 0.010m  (direction of flow)
 u_p = 1m/s
 nu_p = 135.2m^2/s (air)
 g_p = 9.81 m/s^2

with 50 cells resolution in x direction, and 5 cells in y...
 dx_p = 2e-3
 dt_p = dx_p^2 = 4e-6

in lattice units

 u_lb = u_p * (dt_p / dx_p) = 1 * (4e-6 / 2e-3) = 2e-3
 nu_lb = nu_p * (dt_p / dx_p^2) = 135.2e-7 * (4e-6 / 2e-3^2) = 135.2e-7
 g_lb = g_p * (dt_p^2 / dx_p) = 9.81 * (4e-6^2 / 2e-3) = 78.48e-9

So far, I understood well. Referring to J.Latt’s documents, I thought I could use the reynolds number to get viscosity by…


nu_lb = ( u_lb * resolution ) / Re

thus specifying the characteristic velocity for the case, u_lb, implicitly!?

However, I’ve set up a test for that simple 2D poiseuille flow, using the above conversion. I use periodic boundaries in flow direction, thus I do not specify a velocity, only dx, dt, viscosity and acceleration (all as converted above). In my results I receive a reasonable velocity profile, but since I do not specify a velocity (only periodic boundaries and an acceleration), how to interpret the results?

First I thought of using the following…


same extend in x and y, thus having...
 dx_p = 2e-3

with e.g. Re=10, and u_lb = 0.01
 nu_lb = (u_lb * resolution) / Re = (0.01 * 50) / 10 = 0.05
 dt_p = dx_p * u_lb = 20e-6

for gravity, the same as above...
 g_lb = g_p * (dt_p^2 / dx_p) = 9.81 * (20e-6^2 / 2e-3) = 1.962e-6

What, if using this approach? Is it wrong? Results looking good from qualitative point of view, just the maximum magnitude of the velocity is about 0.00122, which is far from u_lb=0.01… I’m really confused!

So, finally, I think I understand what you explained concerning conversion, but I think I do not understand how to interpret the settings according to the results!?

Finally, in my code I’ve just used something like…


 // u_x already computed as sum of pdf's with respect of direction vectors,
 // just add gravity...
 u_x = u_x + tau*g_x

Maybe, I totally misunderstood the principal of how to implement it?

Please, I would appreciate continue dicussion on that topic

Thanks in advance, Francois

#4

The units conversion you specify in the beginning of your post looks correct to me.

I have never implemented a Poiseuille flow, so I can’t help you on that part. Perhaps somebody else here can help you with that. I do have my doubts about this formula you are giving: nu_lb = ( u_lb * resolution ) / Re. Couldn’t find that in the units document nor the thesis of Jonas Latt.

As for problems with implementing the gravity by adding a term to the velocity, a quote of Alex:
“shifting velocity method introduces the viscosity dependent results for macroscopic variables. Even it’s easier and actually is better in terms of stability, but it’s not recommended to be used.” More info on this can be found in this topic: http://www.lbmethod.org/forum/read.php?3,2098

A better way to integrate gravity in your simulation is by adding a term to the distribution functions in the collision operator. More info on this can be found in the topic I also gave in my first post: http://www.lbmethod.org/forum/read.php?3,781 Or you could check out several papers on this subject, especially Guo’s method seems to be a good one: http://www.lbmethod.org/literature:guo_02b

#5

Hi there,

Just a quick follow-on comment. If you are implementing a body force, such as gravity, directly into the standard lattice Boltzmann equation then I think you certainly should use the method explained in Guo’s article. If you don’t then you won’t get the correct macroscopic equations after applying the Chapman-Enskog expansion. Indeed, something similar is true for the continuous analogue: if you add a term S_i to recover a body force F (a vector, of course), the the zeroth moment of S_i should be zero (conservation of mass), the first moment should give you F (the additional force/acceleration) and the second moment should be of the form Fu+uF (where u is the velocity). If the second moment is not of this form then you will modify the viscous stress tensor. The Guo article is a discrete lattice implementation of this idea.

#6

hi bart,
How are you? I have two question for the LBGK with force.
First,
The gravity is minus, or the direction is downward, why the gravity_p is equal to -9.8, and the gravity_lb is equal to -9.8 * (dt_p^2 / dx_p)?
Thank you in advanced.
Second,
For the LBGK with gravity, whether the macro velocity and density are following:
p[x][y] =f[x][y][0]+f[x][y][1]+…f[x][y][8];
u[x][y] =(f[x][y][0]*e[0][0]+f[x][y][1]*e[1][0]+…f[x][y][8]e[8][0])/p[x][y]+0.5G[0];
v[x][y] =(f[x][y][0]*e[0][1]+f[x][y][1]*e[1][1]+…f[x][y][8]e[8][1])/p[x][y]+0.5G[1];
G[0] stands for the horizontal external force term,
G[1] stands for the vertical external force term.
If the G[0] and G[1] are the gravity, so the G[0]=0 and G[1]=-9.8 * (dt_p^2 / dx_p) (i have no idea the G[1] should be positive or negative)

#7

Hello Bruce,

Thank you for pointing out the mistake. It’s indeed -9.8 and not 9.8, otherwise the fluid would move up instead of down :slight_smile:

[quote=Bruce]p[x][y] =f[x][y][0]+f[x][y][1]+…f[x][y][8];
u[x][y] =(f[x][y][0]*e[0][0]+f[x][y][1]*e[1][0]+…f[x][y][8]e[8][0])/p[x][y]+0.5G[0];
v[x][y] =(f[x][y][0]*e[0][1]+f[x][y][1]*e[1][1]+…f[x][y][8]e[8][1])/p[x][y]+0.5G[1];
G[0] stands for the horizontal external force term,
G[1] stands for the vertical external force term.[/quote]
This looks correct to me. G[1] should be equal to G[1] = -9.8 * (dt_p^2 / dx_p)

You can easily check if everything is okay by initializing the distribution functions of your simulation with the equilibrium distribution functions (without considering gravity). You then run your simulation with gravity enabled, then the fluid should experience a shock downward. So you can simply check the values of v[x][y] of some point in your grid (e.g. center point), the values should be negative at least the first few iterations.

#8

Dear all,

First of all about nondimensionalization, there are quite a lot of posts related to it. The easy way for me to look at it is to match the nondimensional parameters. Important! - that means for some quantities there is no physical connections with LBM quantities. In the example above it is certainly Reynolds number:

Re_phys=U_phys * D/nu_phys = 5*4/nu_phys

Re_phys=Re_lbm = U_lbm * N/nu_lbm

N is usually chosen, nu_lbm is usually taken as to have relaxation parameter to be stable for the LBM (tau=1 is a good choice :)). By equating Re_phys=Re_lbm you can obtain U_lbm.

U_lbm = U_phys * dt/dx
dx=D/N

From here and only here you can obtain dt in your simulation:
dt=U_lbm * dx/ U_phys=D/nu_phys* nu_lbm/N *dx=nu_lbm/nu_phys * dx^2

@Francois - don’t think it’s possible to have different scalings in x and y directions as you do in the example. At least for the Navier-Stokes equation which you simulate. It’s possible with the advection-diffusion equation.

With Guo formulation of force plus BB or Zou-He everywhere one can obtain machine accurate results for the Poiseuillle flow.

Hopefully it will help,
Alex