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