Plaintext
Choice of units in lattice Boltzmann simulations
Jonas Latt
April 2008
This document is part of the LBMethod.org documentation project, and may be used under the
terms of a Creative Commons Attribution-Share Alike 3.0 License. See www.lbmethod.org/legal
for details.
1 Rationale
Lattice Boltzmann (LB) simulations are supposed to represent the physics of an actually existing,
real system. During implementation, the question invariably pops up how to chose the units of the
simulated quantities, the “lattice variables”. Two constraints determine the choice of units. First,
the simulation should be equivalent, in a well defined sense, to the physical system. Second, the
parameters should be fine-tuned in order to reach the required accuracy, i.e. the grid should be
sufficiently resolved, the discrete time step sufficiently small and so on.
The present text assumes that your final aim is to solve a macroscopic fluid equation, such
as the incompressible Navier-Stokes equation. Therefore, discretization of the system is discussed
in terms of macroscopic variables, as it is common in computational fluid dynamics. A different
approach should be chosen if you investigate a fluid from a microscopic point of view, that is, if you
explicitly wish to solve the continuum Boltzmann equation, to represent for example high Knudsen
number flows.
The approach presented here consists of two steps. A physical system is first converted into
a dimensionless one, which is independent of the original physical scales, but also independent of
simulation parameters. In a second step, the dimensionless system is converted into a discrete
simulation. The correspondence between these three systems (the physical one (P), the dimen-
sionless one (D), and the discrete one (LB) ) is made through dimensionless, or scale-independent
numbers. The solutions to the incompressible Navier-Stokes equations for example depend only on
one dimensionless parameter, which is the Reynolds number (Re). Thus, the three systems (P),
(D), and (LB) are defined so as to have the same Reynolds number. The transition from (P) to (D)
is made through the choice of a characteristic length scale l0 and time scale t0 , and the transition
from (D) to (LB) through the choice of a discrete space step δx and time step δt . Here’s a graphical
representation of this relationship:
Re,l0 ,t0 Re,δx ,δt
Physical system (P) ←→ Dimensionless system (D) ←→ Discrete system (LB)
This is not the only possible approach. Obviously, one could directly pass from (P) to (LB). There
are however several reasons for which you are encouraged to take an intermediate step and explicit
the dimensionless system. First of all, the discrete variables δx and δt are important parameters
which have impact on the accuracy and the stability of a simulation. They do not depend on the
scales of the considered physical system, and definitely not on the arbitrary choice of physical units.
Further, you are very likely to get one day into a situation in which you don’t have any physical
system to refer to at all. This is for example the case when you try to reproduce the results of a
benchmark computation which was obtained by another model than LB. And finally, getting rid
1
of the original physical units is an important abstraction that gets you in the right mood for what
you are going to do, namely numerics.
2 Example 1: Incompressible fluid flow
2.1 Governing equations
In an incompressible fluid, the density takes a constant value ρ = ρ0 which does not vary in time
and space. The equations of motion, the Navier-Stokes equations, are governed by the laws for
mass and momentum conservation. The conservation law for mass states that the velocity field is
divergence-less:
∇p · up = 0. (1)
Here, u is the velocity, and the index “p” indicates that variables and derivatives are evaluated in
physical units. The conservation law for momentum leads to the following relation:
1
∂tp up + (up · ∇p )up = − ∇p pp + νp ∇2p up , (2)
ρ0p
where pp is the pressure and νp the kinematic viscosity in physical units.
2.2 Dimensionless formulation
Equations (1) and (2) are now cast into a dimensionless form. For this, a length scale l0 and a time
scale t0 are introduced which are representative for the flow configuration. The length l0 could for
example stand for the size of an obstacle which is immersed in the fluid, and t0 could be the time
needed by a passive scalar in the fluid to travel a distance l0 . The physical variables such as the
time tp and the position vector rp , are replaced by their dimensionless counterpart:
tp rp
td = and rd = (3)
t0,p l0,p
In the same manner, a unit conversion is introduced for other variables, based on a dimensional
analysis:
2
l0,p
l0,p 1 1
up = ud , ∂tp = ∂td , ∇p = ∇d , and pp = ρ0 pd . (4)
t0,p t0,p l0,p t20,p
Plugging this change of variables into Eqs. 1 and 2 leads to the dimensionless version of the Navier-
Stokes equations:
1 2
∂td ud + (ud · ∇d )ud = −∇d pd + ∇ ud and (5)
Re d
∇d · ud = 0, (6)
where the dimensionless Reynolds number has been defined as
l02
Re = , (7)
t0 ν
evaluated in any system of units. Two flows obeying the Navier-Stokes equations are equivalent if
they are embedded in the same geometry (except for a scaling factor) and have the same Reynolds
number.
Note that by expressing reference variables in the dimensionless system, one finds that l0,d = 1
and t0,d = 1. To help intuition, one may therefore consider the dimensionless system as “the system
in which l0 and t0 are unity”. It is also remarked that the viscosity in the dimensionless system is
νd = 1/Re.
2
2.3 Discretization of the dimensionless system
The discrete space interval δx is defined as the reference length divided by the number of cells N
used to discretize this length. In the same way, δt is defined as the reference time divided by the
number of iteration steps Niter needed to reach this time. Recall that both reference variables are
unity in the dimensionless system. Thus, the discretization parameters are
δx = 1/N and (8)
δt = 1/Niter . (9)
Other variables, such as velocity and viscosity, are easily converted between (D) and (LB) through
a dimensionless analysis:
δx
ud = ulb and (10)
δt
δx2
νd ≡ 1/Re = νlb , (11)
δt
and thus,
δt
ulb = ud and (12)
δx
δt 1
νlb = . (13)
δx2 Re
Finally, defining the reference velocity u0 ≡ l0 /t0 , one finds, by definition,
u0,d = 1 and (14)
δt
u0,lb = . (15)
δx
2.4 Simulation setup
Let us implement numerically a 2D flow which is confined within a box of size 3cm × 3cm, and
brought into movement by the top-lid, which moves to the right at a speed of 2cm/min. The
viscosity of the fluid is 5cm2 /min (raspberry jam).
1. Define a characteristic length and time. Given the geometry of the problem, it makes sense
to define a reference length l0 as the length of the box, and a reference speed u0 as the speed
of the top lid. Thus, l0,p = 3cm, u0,p = 2cm/min. The characteristic time is t0,p = l0,p /u0,p =
3/2min.
2. Compute the Reynolds number: Re = u0,p l0,p /νp = 2 · 3/5 = 6/5.
3. Chose discretization parameters. Say that we want a grid which consists of 101×101 nodes. If
the physical location of the boundary is on top of a grid node (as is the case for example with
the Zou/He boundary condition), the size of the system is 100×100 in lattice variables. Thus,
the discrete grid spacing is δx = 1/100. Furthermore, let’s chose a given time resolution, say
δt = 2 · 10−4 .
4. Compute the value of variables in the simulation. The velocity ulb , which is used to implement
the boundary condition on the top lid, is evaluated from Eq. 12. The lattice viscosity νlb is
computed from Eq. 13. If you are using single-relaxation time BGK, the relaxation time is
computed through the relation τ = νlb /c2s +1/2, where the speed of sound is a lattice constant:
c2s = 1/3.
3
2.5 Discussion: How to chose δt ?
There is no straightforward intuitive way to chose δt . In other numerical schemes than LB, δt is
often linked to δx by stability considerations. In explicit time-stepping schemes, it is common to
use the relation δt ∼ δx2 to keep the model numerically stable.
In lattice Boltzmann, δt and δx are linked by a different constraint. From the preceding discus-
sion, it is clear that the velocity, measured in lattice units, is u0,lb ∼ δt /δx . Even if you allow for
the possibility that your fluid is compressible, the value of ulb may not be larger than the speed
of sound√cs , because the model does not support supersonic flows. This leads to the constraint
δt < δx / 3. An even stronger constraint appears for the simulation of incompressible flows, i.e.
when you explicitly want to solve Eq. 2 and nothing else. The LB model is a quasi-compressible
fluid solver. This means that it enters a slightly compressible regime to solve the pressure equation
of the fluid. Compressibility effects do however impact the numerical accuracy. As these effects
scale like the square Mach-number, Ma2 , they are kept under control by keeping the Mach number
low. The Mach number is nothing else than u0,lb /cs , which means, it is proportional to u0,lb . The
compressibility error ε(Ma) therefore scales like ε(Ma) ∼ δt2 /δx2 . Now imagine that you increase the
resolution of your lattice. The reason for doing this is most probably to reduce the numerical error
of the model. As the LB model is second-order accurate, the lattice error is ε(δx ) ∼ δx2 . There’s of
course no meaning in reducing the lattice error if after this the compressibility error takes over. A
sensible thing to do is therefore to keep both error terms at the same order: ε(Ma) ∼ ε(δx ). This
immediately leads to the relation
δt ∼ δx2 , (16)
which, not surprisingly, is the same constraint as the one incurred in other explicit fluid solvers.
In practice, a good value to start with, on a lattice size of N ∼ 100, is ulb ∼ 0.02. If you have
little concern about the accuracy of your simulation, you may want to increase this value a bit to
increase δt , and thus speed up your simulation. But if you start increasing the lattice resolution,
never forget to scale the time resolution as imposed by Eq. 16, or you will get no benefit at all from
the larger lattice.
3 Example 2: Thermal flow with Boussinesq approximation
3.1 Governing equations
A thermal fluid can be approximated by solving the equation of temperature separately through an
advection-diffusion equation. In the Boussinesq approximation, the fluid-temperature interaction
is represented by a linear buoyancy term which acts as a body force on the fluid. The macroscopic
equations are
∂tp up + (up · ∇p )up = −∇pp + νp ∇2p up + α(Tp − Tc,p )g0,p j for the fluid, and (17)
∂t,p T + ∇ · (up Tp ) = dp ∇2p Tp for the temperature. (18)
Here, the continuity equation is omitted, and ρ0 is taken unity without loss of generality, to keep the
discussion shorter. The parameter α is the coefficient of volume expansion, d the thermal diffusivity,
g0 the gravitational acceleration, and j a unit vector pointing in direction of gravitation.
The constant Tc determines an absolute offset of the temperature field. It is enforced in the
choice of the initial condition, by choosing for example the average temperature to be initially equal
to Tc , and in the choice of the temperature boundary condition. In Rayleigh-Benard convection for
example, the boundary condition for the temperature on two opposite boundaries is defined as
1
Tp,hot = Tp,c + ∆Tp and (19)
2
1
Tp,cold = Tp,c − ∆Tp ,
2
4
where ∆T is the temperature difference between the two boundaries. With this initial and boundary
condition, the constant Tc can be chosen arbitrarily and does not affect the flow evolution, as the
buoyancy term only depends on the difference T − Tc .
3.2 Dimensionless formulation
As in the previous example, the system is made dimensionless by introducing a reference length
l0 and a reference time t0 . Furthermore, a reference value T0 is used to make the temperature
dimensionless. Actually, this system possesses only two degrees of freedom instead of three for the
choice of reference units. This is because the occurrence of the constant g0 in Eq. 17 introduces
a natural reference frame. You may recall the example of a lid driven cavity flow in the previous
section, in which the value of t0 was derived from a characteristic value of the velocity, namely the
velocity of the moving lid. It appears meaningful to use a similar approach in the present case, this
time based on the acceleration g0 . A straightforward attempt would be to define a relation of the
type g0 = l0 /t20 , but an even smarter choice exists with which one gets rid of the volume-expansion
coefficient at the same time:
s
l0 l0
g0 = 2 =⇒ t0 = . (20)
t0 αT0 g0 αT0
Finally, the various dimensionless variables are defined through the relations
2
l0,p
1 1 l0,p
Tp = T0,p Td + Tc,p , ∂t,p = ∂t,d , ∇p = ∇d , up = ud , and pp = pd . (21)
t0,p l0,p t0,p t20,p
Plugging Eqs. 20 and 21 into Eqs. 17 and 18 yields the dimensionless equations of motion:
r
Pr 2
∂td ud + (ud · ∇d )ud = −∇pd + ∇ ud + Td j for the fluid, and (22)
Ra d
r
1
∂t,d Td + ∇ · (ud Td ) = ∇2 Td for the temperature, (23)
Ra · Pr d
where the dimensionless Prandtl number Pr and Rayleigh number Ra have been defined as
ν g0 α∆T l02
Pr = and Ra = . (24)
d νd
The temperature boundary condition of Eq. 19 becomes:
1
Td,hot = ∆Td and (25)
2
1
Td,cold = − ∆Td
2
It is interesting at this point to think of an appropriate definition of T0 . As discussed before, the
constant Tc is not eligible as a reference, as its value has no influence on the flow evolution. Instead,
a good choice is T0 ≡ ∆T , in which case one gets ∆Td = 1.
3.3 Discretization of the dimensionless system
Getting a discrete version of the equations of motion is straightforward. Again, the discrete space
and time steps are defined as δx = 1/N and δt = 1/Niter , the inverse resolution of the reference
length and the reference time. Furthermore, the value of T0 in lattice units needs to be picked out.
The following discussion is restricted to the choice of T0 ≡ ∆T . Making a choice for the value of
5
∆Tlb is a less crucial step than fixing the value of δx and δt . Indeed, ∆Tlb does not intervene in
the calculation of lattice variables, except in determining the value of Tlb (for example in the initial
or boundary condition). Even then, the temperature scale imposed through ∆Tlb is not expected
to affect the numerics in a significant way, as the advection-diffusion equation for the temperature
is linear. Finally, the possibility of having a temperature with non-zero offset Tc,lb is considered
for generality. The conversion formula between lattice units and the dimensionless system for the
temperature becomes:
1
Td = (Tlb − Tc,lb ) =⇒ Tlb = Td ∆Tlb + Tc,lb . (26)
∆Tlb
The fluid viscosity and thermal diffusivity are recovered as follows from a dimensional analysis:
r r
δt Pr δt 1
νlb = 2 and dlb = 2 . (27)
δx Ra δx Ra · Pr
Finally, the buoyancy force term f for the simulation is computed as follows:
δx δt2 Tlb − Tc,lb
fd = Td = flb =⇒ flb = . (28)
δt2 δx ∆Tlb
3.4 Simulation setup
The following steps lead to the definition of a simulation for Rayleigh-Benard convection:
1. Chose the discrete steps δx and δt . As before, δt should be chosen in such a way as to keep
the Mach number small. This is however not an easy task, because the Mach number is
not necessarily proportional to δt /δx , as it was in the simple Navier-Stokes case. Instead,
you should use your knowledge of the physics of the flow to estimate a typical value of ulb
during the flow evolution. If you have no such knowledge, run a few simulation to obtain the
information numerically.
2. Chose Tc,lb and ∆Tlb . If you have no specific reason for doing otherwise, pick Tc,lb = 0 and
∆Tlb = 1 for simplicity. This leads to the boundary condition Thot = 1/2 and Tcold = −1/2.
3. Compute the viscosity and thermal diffusivity in lattice units from Eq. 27. If you use a BGK
model for the fluid and for the temperature, you may determine the corresponding relaxation
times as τfluid = νlb /c2s + 1/2 and τtemp = dlb /c2s + 1/2.
4. Compute the value of the force term in lattice units from Eq. 28.
6