Search

[LBM Study] 5. Non-Dimensionalization and Choice of Simulation Parameters

Category
Study
Kewords
LBM
Kruger
3 more properties
Table of Contents

5.1 Non-dimensionalization

5.1.1 Unit Scales and Conversion Factors

•
Physical quantities require reference scales for dimensional representation. ⇒ For example, a length ℓ\ell can be reported in multiples of a unit scale with length ℓ0=1\ell_0 = 1 m.
•
Non-dimensionalization is achieved by dividing a dimensional quantity by a chosen reference quantity. ⇒ Resulting in a dimensionless number called the lattice value or the quantity's value in lattice units.
•
The reference quantity is called the conversion factor, denoted by CC. ⇒ Non-dimensionalised quantities are denoted by a star ∗*, giving the relation: ℓ∗=ℓ/Cℓ\ell^* = \ell/C_\ell.
•
Any mechanical quantity qq has a dimension which is a combination of the dimensions of length ℓ\ell, time tt, and mass mm: [q]=[ℓ]qℓ[t]qt[m]qm[q] = [\ell]^{q_\ell}[t]^{q_t}[m]^{q_m} where the exponents are numbers.
◦
Since three fundamental dimensions are sufficient to generate the dimension of any mechanical quantity, one requires exactly three independent conversion factors to define a unique non-dimensionalization scheme.
◦
Three conversion factors CiC_i (i=1,2,3i = 1, 2, 3) are independent if the following relations have no solution for the numbers aia_i and bib_i:
[Ci]=[Cj]ai[Ck]bi(i≠j≠k≠i).[C_i] = [C_j]^{a_i}[C_k]^{b_i} \quad (i ≠ j ≠ k ≠ i).
•
In LB simulations, one usually takes CℓC_\ell, CtC_t (or CuC_u) and CρC_\rho as basic conversion factors because length, time (or velocity) and density are natural quantities in any LB simulation.
◦
For 2D simulations, the system should be treated as a 3D system with thickness of one lattice constant to enable the use of 3D conversion factors without restrictions.

5.1.2 Law of Similarity and Derived Conversion Factors

•
Physics is independent of units which are an arbitrary human construct.
◦
Ratios of physical phenomena are what matter, and the physical outcome should not depend on whether we use dimensional or dimensionless quantities.
◦
The law of similarity in fluid dynamics states that two incompressible flow systems are dynamically similar if they have the same Reynolds number and geometry.
•
The Reynolds number is defined as:
Re=ℓUν=ρℓUη.\text{Re} = \frac{\ell U}{\nu} = \frac{\rho \ell U}{\eta}.
◦
ℓ\ell and UU are typical length and velocity scales.
◦
ρ\rho, ν\nu, and η\eta are density, kinematic viscosity, and dynamic viscosity of the fluid, respectively.
•
The Reynolds number must be identical in both physical and lattice systems: $$\ell^* U^/\nu^ = \ell U/\nu$$, which yields the relation $$C_\ell C_u/C_\nu = 1$$.
•
From this relation, we can derive that $$C_\nu = C_\ell C_u$$, showing how the conversion factor for viscosity is related to those for length and velocity.
•
Any derived conversion factor $$C_q$$ can be constructed directly by writing down a suitable combination of basic conversion factors: $$C_q = C_1^{q_1} C_2^{q_2} C_3^{q_3}$$ without any additional numerical prefactor.
•
This construction is unique and guarantees that the conversion is consistent - the physics of the system and characteristic dimensionless numbers are kept invariant.
•
Unit systems must not be mixed as this causes inconsistencies in the definition of conversion factors and subsequent violation of the law of similarity.

5.2 Parameter Selection

5.2.1 Parameters in the Lattice Boltzmann Method

•
Any LB simulation is characterized by a set of parameters:
◦
The lattice constant Δx\Delta x is the distance between neighbouring lattice nodes in physical units, i.e. [Δx\Delta x]=[m].
◦
The physical length of a time step is denoted Δt\Delta t, thereforce [Δt\Delta t] = [s].
◦
The BGK relaxation parameter τ\tau is a relaxation time with [τ\tau] = [s], where τ∗\tau^* denotes the dimensionless relaxation parameter.
◦
The dimensionless fluid density ρ∗\rho^* has an average value usually set to unity: ρ0∗=1\rho_0^ *= 1. ⇒ This situation is more complicated for multicomponent or multiphase simulations.
◦
The typical simulated velocity U∗U^* is usually part of the simulation output but may need to be specified on boundaries as inlet and outlet velocities. ⇒ It is desired to estimate the magnitude of U∗U^* before the simulation is started in order to avoid unstable situations or very long computing times.
◦
In the standard LBM, the lattice speed of sound cs∗c_\mathrm{s}^* is 1/3≈0.577\sqrt{1/3} \approx 0.577, and all simulated velocities must be significantly smaller: U∗≪cs∗U^* \ll c_s^*. ⇒ In practice this means that the maximum value of U∗U^* should be below 0.2.
•
It is common to set Δx∗=1\Delta x^* = 1, Δt∗=1\Delta t^* = 1, and ρ0∗=1\rho_0^* = 1.
◦
This means that the conversion factors for length, time, and density equal the dimensional values for the lattice constant, time step, and density:
Cℓ=Δx,Ct=Δt,Cρ=ρ.C_\ell = \Delta x, \quad C_t = \Delta t, \quad C_\rho = \rho.
◦
The units defined by Δx∗=1\Delta x^* = 1 and Δt∗=1\Delta t^* = 1 are called lattice units.
•
The parameters τ\tau and τ∗\tau^* are connected through the conversion factor for Δt\Delta t because τ\tau has the dimension of time:
τ=τ∗Ct=τ∗Δt.\tau = \tau^* C_t = \tau^* \Delta t.

Viscosity

•
The kinematic lattice viscosity is related to the relaxation parameter according to ν∗=cs∗2(τ∗−1/2)\nu^* = c_\mathrm{s}^{*2}(\tau^* - 1/2). ⇒ The typical problem is to relate the dimensionless relaxation parameter τ∗\tau^* to the physical kinematic viscosity ν\nu since the latter is usually given by an experiment and the former has to be defined for a simulation.
•
The kinematic viscosity is related to the simulation parameters according to:
v=cs∗2(τ∗−12)Δx2Δt.v = c_\mathrm{s}^{*2} \left( \tau^* - \frac{1}{2} \right) \frac{\Delta x^2}{\Delta t}.
◦
This is a consistency equation showing that τ∗\tau^*, Δx\Delta x, and Δt\Delta t are not independent. → Only two of them can be chosen freely.

Pressure, Stress and Force

•
The equation of state of the LB fluid is:
p∗=cs∗2ρ∗.p^* = c_\mathrm{s}^{*2}\rho^*.
◦
This is, however, not the entire truth. ⇒ Only the pressure gradient ∇p\nabla p rather than the pressure pp by itself appears in the NSE.
◦
The total pressure pp does appear in the energy equation. → This euation is not relevant for non-thermal LB models.
◦
The reference pressure is thus irrelevant; only pressure changes matter.
•
To connect with the physical pressure, one decomposes the LB density into its constant average ρ0∗\rho_0^* and deviation ρ′∗\rho^{'*} from the average:
ρ∗=ρ0∗+ρ′∗.\rho^* = \rho_0^* + \rho'^*.
•
Generally, the LB density can be converted to the physical pressure for non-thermal models as:
p=p0+p′=p0+p′∗Cp,p′∗=cs∗2ρ′∗.p = p_0 + p' = p_0 + p'^* C_p, \quad p'^* = c_\mathrm{s}^{*2} \rho'^*.
◦
Cp=CρCℓ2/Ct2=CρCu2C_p = C_\rho C_\ell^2/C_t^2 = C_\rho C_u^2 is the conversion factor for pressure.
◦
p0p_0 is the physical reference pressure which can be freely specified by the user.
•
The components of the stress tensor have the same dimension as a pressure, therefore the conversion factor CρCℓ2/Ct2C_\rho C_\ell^2/C_t^2 is always identical: σ=σ∗Cσ\boldsymbol{\sigma} = \boldsymbol{\sigma}^* C_\sigma with Cσ=CpC_\sigma=C_p.
•
The conversion factor for any force (no matter if body force or surface force) is Cf=CρCℓ4/Ct2C_f = C_\rho C_\ell^4 / C_t^2.
•
The conversion factor for a body force density FF is obviously CF=Cf/Cℓ3=CρCℓ/Ct2C_F = C_f / C_\ell^3 = C_\rho C_\ell / C_t^2.

5.2.2 Accuracy, Stability and Efficiency

Accuracy and Parameter Scaling

•
There are several error terms which affect the accuracy of an LB simulation:
◦
The spatial discretization error scales like Δx2\Delta x^2.
◦
The time discretization error scales like Δt2\Delta t^2.
◦
The compressibility error for simulations in the incompressible limit is proportional to Ma2∝U∗2\text{Ma}^2 \propto U^{*2}. ⇒ Since U∗U^* decreases with increasing Cu=Cℓ/Ct=Δx/ΔtC_u = C_\ell/C_t=\Delta x / \Delta t, this error scales like Δt2/Δx2\Delta t^2/\Delta x^2.
◦
The BGK truncation error in space is proportional to (τ−1/2)2(\tau - 1/2)^2, indicating that τ∗\tau^* should not be much larger than unity.
•
One needs to come up with certain relationships between Δx\Delta x and Δt\Delta t to control the error.
•
The diffusive scaling Δt∝Δx2\Delta t \propto \Delta x^2 guarantees that the leading order of the overall error scales like Δx2\Delta x^2, though the LB algorithm becomes effectively first-order accurate in time.
◦
The diffusive scaling leaves τ∗\tau^*, and hence the non-dimensional viscosity ν∗\nu^*, unchanged.
◦
The diffusive scaling is the standard approach to test if an LB algorithm is second-order accurate: one performs a series of simulations, each with a finer resolution Δx\Delta x than the previous.
◦
The overall velocity error should then decrease proportionally to Δx2\Delta x^2.
•
The acoustic scaling Δt∝Δx\Delta t \propto \Delta x keeps the compressibility error unchanged and must be chosen when the speed of sound is a physically relevant parameter.
◦
For incompressible simulations, the numerical solution can only converge to the incompressible NSE when Δt∝Δxγ\Delta t \propto \Delta x^\gamma with γ>1\gamma > 1.

Stability

•
The relaxation parameter τ∗\tau^* should not be too close to 1/2, and the velocity U∗U^* should not be larger than about 0.4 for τ∗≥0.55\tau^* \geq 0.55.
◦
For τ∗<0.55\tau^* < 0.55, the stability criterion can be approximated as:
τ∗>12+αUmax⁡∗.\tau^* > \frac{1}{2} + \alpha U_{\max}^*.
▪
α\alpha is a numerical constant which is of the order of 1/81/8.
•
The grid Reynolds number Reg\text{Re}_g is defined by taking the lattice resolution Δx\Delta x as length scale:
Reg=Umax⁡∗Δx∗ν∗=Umax⁡∗cs∗2(τ∗−12)⟹τ∗=12+Umax⁡∗cs∗2Reg.\mathrm{Re}_g = \frac{U_{\max}^* \Delta x^*}{\nu^*} = \frac{U_{\max}^*}{c_\mathrm{s}^{*2} \left( \tau^* - \frac{1}{2} \right)} \quad \Longrightarrow \quad \tau^* = \frac{1}{2} + \frac{U_{\max}^*}{c_\mathrm{s}^{*2} \mathrm{Re}_g}.
◦
The grid Reynolds number Reg\text{Re}_g should not be much larger than O(10)O(10).
◦
The physical interpretation is that the lattice should always be sufficiently fine to resolve local vortices. ⇒ The simulation usually remains stable as long as all relevant hydrodynamic length scales are resolved.
•
Instabilities are often triggered at boundaries rather than in the bulk, so boundary treatment is crucial for stability.

Efficiency

•
The total number of site updates required is NsNtN_s N_t where NsN_s is the total number of lattice sites and NtN_t is the required number of time steps.
•
The memory requirements are proportional to NsN_s, making LB quite memory-hungry method.
•
The total required runtime TT and memory MM obey:
T∝1ΔxdΔt,M∝1Δxd.T \propto \frac{1}{\Delta x^d \Delta t}, \quad M \propto \frac{1}{\Delta x^d}.
◦
dd is the spatial dimensions.
◦
It is important to choose Δt\Delta t and especially Δx\Delta x as large as possible to reduce the computational requirements while maintaining accuracy and stability.

5.2.3 Strategies for Parameter Selection

Mapping of Dimensionless Physical Parameters

•
Any physical system can be characterized by dimensionless parameters like Reynolds or Mach numbers. ⇒ The first step before setting up a simulation is to identify these parameters and assess their relevance.
◦
Inertia is relevant as long as Re is larger than order unity - Only flows with vanishing small Reynolds numbers (Stokes flow) do not depend on the actual value of Re.
•
LB is mostly used for the simulation of incompressible fluids where the Mach number is small
◦
It is not necessary to map the exact value of Ma\text{Ma}. → It is sufficient to guarantee that Ma\text{Ma} is "small" in the simulation.
•
A lattice Mach number Ma(l)=U/cs∗\text{Ma}^{(l)} = U^/c_s^* is considered small if Ma(l)<0.3\text{Ma}^{(l)} < 0.3.
•
The relation of simulation parameters in terms of Reynolds and Mach numbers can be written in the useful form:
Re(l)=ℓ∗U∗ν∗=ℓ∗U∗cs∗2(τ∗−12)⟹Re(l)Ma(l)=ℓcs∗(τ∗−12)Δx.\mathrm{Re}^{(l)} = \frac{\ell^* U^*}{\nu^*} = \frac{\ell^* U^*}{c_\mathrm{s}^{*2} \left( \tau^* - \frac{1}{2} \right)} \quad \Longrightarrow \quad \frac{\mathrm{Re}^{(l)}}{\mathrm{Ma}^{(l)}} = \frac{\ell}{c_\mathrm{s}^* \left( \tau^* - \frac{1}{2} \right) \Delta x}.
◦
ℓ∗=ℓ/Δx\ell^* = \ell / \Delta x is a typical system length scale in lattice units.
◦
It should be noted that Ma(l)/Re(l)\mathrm{Ma}^{(l)} / \mathrm{Re}^{(l)} is the lattice Knudsen number Kn\text{Kn}. ⇒ Hydrodynamic behaviour is only expected for sufficiently small Knudsen numbers. → This sets an upper bound for the lattice constant Δx\Delta x.

Parameter Selection Strategies

•
The first step is to set the lattice density ρ0∗\rho_0^*, typically to unity. ⇒ It is a pure scaling parameter without significant effect on accuracy, stability, or efficiency.
•
A typical scenario is having a maximum lattice size ℓ∗\ell^* that can be handled by the computer, suggesting the lattice constant Δx\Delta x should be set next.
•
Choose the lattice Mach number (or the non-dimensional velocity U∗U^*) reasonably, e.g., Ma(l)=0.1\text{Ma}^{(l)} = 0.1 or U∗=0.1U^* = 0.1.
•
For given system size ℓ∗\ell^*, velocity U∗U^*, and Reynolds number, τ∗\tau^* can be computed from the following relation:
Reg=Umax⁡∗Δx∗ν∗=Umax⁡∗cs∗2(τ∗−12)⟹τ∗=12+Umax⁡∗cs∗2Reg.\mathrm{Re}_g = \frac{U_{\max}^* \Delta x^*}{\nu^*} = \frac{U_{\max}^*}{c_\mathrm{s}^{*2} \left( \tau^* - \frac{1}{2} \right)} \quad \Longrightarrow \quad \tau^* = \frac{1}{2} + \frac{U_{\max}^*}{c_\mathrm{s}^{*2} \mathrm{Re}_g}.
•
Check via stability criteria whether the chosen values for U∗=Ma(l)cs∗U^* = \text{Ma}^{(l)} c_\mathrm{s}^* and τ∗\tau^* provide stable simulations.
τ∗>12+αUmax⁡∗.\tau^* > \frac{1}{2} + \alpha U_{\max}^*.
◦
If τ∗\tau^* is too small, either decrease Δx\Delta x (increase ℓ∗\ell^*, more expensive) or increase Ma(l)\text{Ma}^{(l)} (increase U∗U^*, less accurate/stable).
•
Once the parameters U∗U^* and Δx\Delta x are fixed, we can calculate Δt\Delta t.
This scenario reveals some problems arising when large Reynolds numbers are simulated: It requires either large lattices, small relaxation parameters or large Mach numbers.
⇒ It is possible only to a limited extent to reach large Reynolds numbers by increasing the Mach number and decreasing the relaxation parameter due to accuracy and stability issues.
•
The maximum achievable Reynolds number for a given lattice size, assuming τ∗\tau^* near the stability limit (τ∗=1/2+U∗/4\tau^* = 1/2 + U^*/4):
Re(l)=ℓ∗U∗cs∗2(τ∗−12)=4ℓ∗cs∗2.\mathrm{Re}^{(l)} = \frac{\ell^* U^*}{c_\mathrm{s}^{*2} \left( \tau^* - \frac{1}{2} \right)} = \frac{4\ell^*}{c_\mathrm{s}^{*2}}.
◦
This shows that the achievable Reynolds number is limited by O(10)ℓ∗O(10) \ell^*.
•
Alternative methods involve setting Mach number and viscosity first, or setting resolution and viscosity first, then verifying the validity of all parameters.

Small Reynolds Numbers

•
Small Reynolds numbers can be reached by choosing a large Δx\Delta x, a large relaxation parameter τ∗\tau^*, or small lattice Mach number Ma(l)\text{Ma}^{(l)}.
•
The resolution cannot be arbitrarily decreased.
◦
At some point the lattice domain is so smal that the details of the flow are finer than the spacing between lattice nodes.
•
Using τ∗≫1\tau^* \gg 1 is not advisable due to increased numerical errors.
◦
This can be avoided by using advanced collision operators such as TRT or MRT).
•
The only unbounded way to reduce Re(l)\text{Re}^{(l)} is to decrease Ma(l)\text{Ma}^{(l)} and therefore the flow velocity U∗U^*. ⇒ But this means Δt\Delta t becomes very small since Δt∝Ma(l)\Delta t \propto \text{Ma}^{(l)}.
•
In some situations, when Reynolds number is not important (e.g., capillary flows, microfluidics), we can use a numerical Reynolds number which is larger than the physical Reynolds number to accelerate the simulations: Re(l)>Re(p)\text{Re}^{(l)} > \text{Re}^{(p)}.
◦
One should always check if the simulation results are still valid.
•
LB is particularly useful for flow problems with intermediate Reynolds numbers, especially between O(1)O(1) and O(100)O(100).

5.3 Examples

5.3.1 Poiseuille Flow I

•
Consider force-driven 2D Poiseuille flow with physical parameters: channel diameter $$w = 10^{-3}$$ m, kinematic viscosity $$\nu = 10^{-6}$$ m²/s, density $$\rho = 10^3$$ kg/m³, gravity $$g = 10$$ m/s².