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:
ℓ∗U∗ν∗=ℓUν.\frac{\ell^* U^*}{\nu^*}=\frac{\ell U}{\nu}.
◦
Plugging in the definition of the conversion factors yields the relation CℓCuCν=1\frac{C_\ell C_u}{C_\nu} = 1 or Cν=CℓCuC_\nu=C_\ell C_u.
•
Any derived conversion factor CqC_q can be constructed directly by writing down a suitable combination of basic conversion factors without any additional numerical prefactor:
Cq=C1q1C2q2C3q3.C_q = C_1^{q_1} C_2^{q_2} C_3^{q_3}.
◦
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}^{(\mathrm{l})} = U^/c_s^* is considered small if Ma(l)<0.3\text{Ma}^{(\mathrm{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}^{(\mathrm{l})} = \frac{\ell^* U^*}{\nu^*} = \frac{\ell^* U^*}{c_\mathrm{s}^{*2} \left( \tau^* - \frac{1}{2} \right)} \quad \Longrightarrow \quad \frac{\mathrm{Re}^{(\mathrm{l})}}{\mathrm{Ma}^{(\mathrm{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}^{(\mathrm{l})} / \mathrm{Re}^{(\mathrm{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}^{(\mathrm{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}^{(\mathrm{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}^{(\mathrm{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}^{(\mathrm{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}^{(\mathrm{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}^{(\mathrm{l})} is to decrease Ma(l)\text{Ma}^{(\mathrm{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}^{(\mathrm{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}^{(\mathrm{l})} > \text{Re}^{(\mathrm{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

A Force-Driven 2D Poiseuille Flow
Channel diameter w=10−3w = 10^{-3} m, kinematic viscosity ν=10−6\nu = 10^{-6} m²/s, density ρ=103\rho = 10^3 kg/m³, gravity g=10g = 10 m/s²
•
The expected centre velocity from analytical solution:
u^=gw28ν=1.25 m/s.\hat{u} = \frac{gw^2}{8\nu} = 1.25 \, \mathrm{m/s}.
•
Reynolds number:
Re=u^wν=1250.\mathrm{Re} = \frac{\hat{u}w}{\nu} = 1250.
•
Initial parameters: Δx=5×10−5\Delta x = 5 \times 10^{-5} m (corresponding to w∗=20w^* = 20), τ∗=0.6\tau^* = 0.6, ρ0∗=1\rho_0^* = 1 (Cρ=103C_\rho = 10^3 kg/m³). ⇒ The system is now fully determined (the user may start with a different set of initial values).
•
From viscosity relation, we can obtain the conversion factor:
Δt=cs∗2(τ∗−12)Δx2v=8.33⋅10−5 s.\Delta t = c_\mathrm{s}^{*2} \left( \tau^* - \frac{1}{2} \right) \frac{\Delta x^2}{v} = 8.33 \cdot 10^{-5} \, \mathrm{s}.
◦
This means that 12,000 time steps are required for 1 s physical time.
•
Expected lattice velocity:
a.
Using the law of similarity and the relation: Re(p)=Re(l)=u^∗w∗/ν∗.\mathrm{Re}^{(\mathrm{p})} = \mathrm{Re}^{(\mathrm{l})} = \hat{u}^* w^* / \nu^*.
b.
Using the conversion factor for the velocity: Cu=Δx/Δt=0.6C_u = \Delta x / \Delta t=0.6 m/s.
◦
Either case: u^∗=2.08\hat{u}^*=2.08 → It is an invalid value because it is much larger than the speed of sound cs∗=1/3c^*_\mathrm{s} = \sqrt{1/3}.
•
Parameter adjustment needed due to Reynolds number being too large for the chosen resolution. ⇒ This means that we have to increase the spatial resolution.
•
Refined parameters: Δx=1×10−5\Delta x = 1 \times 10^{-5} m, τ∗=0.55\tau^* = 0.55 → Resulting in u^∗=0.208\hat{u}^* = 0.208 (acceptable).
•
In order to run the simulation, the lattice value for the force density F=ρgF=\rho g has to be obtained.
◦
Computing the conversion factor for gg: Cg=Δx/Δt2=3.6⋅106C_g = \Delta x / \Delta t^2=3.6 \cdot 10^6 m/s² where we have used the time conversion factor Δt=1.667⋅10−6\Delta t =1.667 \cdot 10^{-6} s.
◦
Resulting in g∗=2.78⋅10−6g^* = 2.78 \cdot 10^{-6} and F∗=ρ0∗g∗=2.78⋅10−6F^* = \rho_0^* g^* = 2.78 \cdot 10^{-6}.
•
Alternative routes exist for finding simulation parameters, but all must satisfy the Reynolds number constraint while meeting stability and accuracy requirements.
◦
Alternative route 1:
◦
Alternative route 2:
•
If the initially guessed simulation parameters give invalid or otherwise unacceptable results, one has to modify one or two parameters (which are not necessarily the initially chosen parameters) while updating the dependent parameters until the desired level of accuracy, stability, and efficiency is obtained.

5.3.2 Poiseuille Flow II

•
Simulations can be set up using only dimensionless parameters without any given dimensional parameters.
◦
For Poiseuille flow, knowing only the Reynolds number (e.g., Re = 1250) is sufficient.
•
Using relation from equation for Reynolds number with initial guess w∗=100w^* = 100 and u^∗=0.1\hat{u}^* = 0.1.
◦
Direct calculation gives τ∗=0.524\tau^* = 0.524, which should yield stable simulation.
•
The required force density F∗=ρ0∗g∗F^* = \rho_0^* g^* can then be obtained from:
u^∗=g∗w∗28ν∗⟹g∗=6.4⋅10−7.\hat{u}^* = \frac{g^* w^{*2}}{8\nu^*} \quad \Longrightarrow \quad g^* = 6.4 \cdot 10^{-7}.
◦
ρ0∗\rho_0^* has been chosen (ρ0∗=1\rho_0^*=1 in this case).
•
All required simulation parameters obtained without any conversion factors.
◦
The simulation is valid for all similar physical systems with the same Reynolds number.
◦
Conversion factors can be obtained a posteriori by setting physical scales after successful simulation.
◦
Full set of conversion factors only obtained when three independent physical parameters are given (e.g. ww, u^\hat{u}, ρ\rho or FF, ρ\rho, ww).
•
In physics it is often the case that systems are only characterized by the relevant dimensionless parameters without giving the scales like length or velocity themselves. ⇒ This is sufficient to set up and run simulations, at least as long as all required dimensionless parameters are known.

5.3.3 Poiseuille Flow III

•
It is very common to drive a Poiseuille flow by a pressure gradient rather than by gravity or a force density.
•
The pressure gradient p′=Δp/ℓ p' = \Delta p/\ell  (where Δp=pin−pout>0\Delta p = p_{in} - p_{out} >0 is the pressure difference between inlet and outlet and ℓ\ell is the length of the channel).
◦
The force density FF is simply p′=Δp/ℓ=Fp' = \Delta p/\ell = F.
•
For 2D Poiseuille flow:
u^∗=Δp∗ρ0∗w∗28ℓ∗ν∗⟹Δp∗ρ0∗=8ℓ∗ν∗u^∗cs∗2w∗2.\hat{u}^* = \frac{\Delta p^*}{\rho_0^*} \frac{w^{*2}}{8\ell^* \nu^*} \quad \Longrightarrow \quad \frac{\Delta p^*}{\rho_0^*} = \frac{8\ell^* \nu^* \hat{u}^*}{c_\mathrm{s}^{*2} w^{*2}}.
◦
We have introduced the average density ρ0∗\rho_0^* to appreciate the fact that the density is not constatn.
◦
We replaced the pressure difference Δp∗\Delta p^* by the density difference via cs∗2Δρ∗c_\mathrm{s}^{*2} \Delta \rho^*.
◦
This equation provides a direct expression for the required relative difference between the inlet and outlet densities. ⇒ The right-hand-side has to be a small quantity. → This is another restriction for an LB simulation and sets additional bounds, especially for the system length ℓ∗\ell ^*.
•
In the incompressible limit, we require the density difference between any two points in the simulation, Δρ∗\Delta \rho^*, to be small compared to the average density ρ0∗\rho_0^*.
A Pressure-Gradient-Driven 2D Poiseuille Flow
u^∗=0.1\hat{u}^* = 0.1, ν∗=0.1\nu^* = 0.1, w∗=50w^* = 50, cs2=1/3c_\mathrm{s}^{2} = 1/3
Allowing density variations of up to 5% (Δρ∗/ρ0∗=0.05\Delta \rho^* / \rho_0^*=0.05)
The maximum channel length becomes ℓ∗≈520\ell^* \approx 520 (only ten times more than the channel diameter).
•
If Δρ∗\Delta \rho^* and the average density ρ0∗\rho_0^* are known, the inlet and outlet densities are set to:
ρin∗=ρ0∗+Δρ∗2andρout∗=ρ0∗−Δρ∗2.\rho_{\mathrm{in}}^* = \rho_0^* + \frac{\Delta \rho^*}{2} \quad \text{and} \quad \rho_{\mathrm{out}}^* = \rho_0^* - \frac{\Delta \rho^*}{2}.
◦
Δρ∗=ρin∗−ρout∗\Delta \rho^* = \rho_{in}^* - \rho_{out}^* is satisfied.
◦
The local pressure p∗p^* can be obtained from the local density ρ∗\rho^* according to p∗=p0∗+cs2(ρ∗−ρ0∗)p^* = p_0^* + c_\mathrm{s}^{2}(\rho^* - \rho_0^*).

5.3.4 Womersley Flow

•
Womersley flow is Poiseuille-like flow (channel width ww, viscosity ν\nu) with an oscillating pressure drop along the length ℓ\ell of the channel.
◦
Δp(t)=Δp0cos⁡ωt\Delta p(t) = \Delta p_0 \cos \omega t with angular frequency ω\omega.
◦
Since there exists an analytical solution for this unsteady flow, it is an ideal benchmark for Navier-Stokes solvers.
•
The Reynolds number for an oscillatory flow is usually defined through the velocity at ω=0\omega = 0: Re=u^0w/ν\text{Re} = \hat{u}_0 w/\nu.
◦
For a 2D channel, u^0\hat{u}_0 is related to Δp0\Delta p_0 according to:
u^0=Δp0w28ℓρν.\hat{u}_0 = \frac{\Delta p_0 w^2}{8 \ell \rho \nu}.
•
Due to the presence of the frequency ω\omega, an additional dimensionless number is required to characterise the flow:
◦
The Womersley number is defined as:
α=ωνw.\alpha = \sqrt{\frac{\omega}{\nu}}w.
•
The Buckingham π theorem: for QQ independent quantities whose physical dimension can be constructed from DD independent dimensions, there are N=Q−DN = Q - D independent dimensionless parameters.
◦
A simple Poiseuille flow:
▪
Q=4Q=4 independent quantities: channel width, flow velocity, viscosity, and density.
▪
D=3D=3 independent dimensions: length, time, and mass.
▪
This gives N=1N=1 parameter which is the Reynolds number.
◦
Womersley flow:
▪
Q=5Q=5 independent quantities: channel width, flow velocity, viscosity, density, and frequency.
▪
D=3D=3 independent dimensions: length, time, and mass.
▪
This gives N=2N=2 parameters which are the Reynolds number and Womersley number.
◦
Characteristic time scales: the advection time ta∼w/u^t_a \sim w/\hat{u}, the diffusion or viscous time tν∼w2/νt_\nu \sim w^2/\nu, and the oscillation time tω∼1/ωt_\omega \sim 1/\omega.
•
All physical time scales which shall be resolved must be sufficiently long compared to the intrinsic sound wave (or acoustic) time scale ts∗∼ℓ∗/cs∗t_\mathrm{s}^* \sim \ell^*/c_\mathrm{s}^* for incompressible flows.
◦
ta∗,tν∗,tω∗≫ts∗t_{\mathrm{a}}^*, t_{\nu}^*, t_{\omega}^* \gg t_{\mathrm{s}}^*.
◦
This defines, for a given lattice size, a lower bound for the oscillation time scale and therefore an upper bound for the frequency (ℓ∗>w∗\ell^* > w^*): IMaximum frequency constraint:
ω∗≪cs∗ℓ∗.\omega^* \ll \frac{c_\mathrm{s}^*}{\ell^*}.
◦
Any unsteady incompressible flow violating this condition is non-physical, and the results will not be a good approximation of the incompressible Navier-Stokes solution.
•
A suggested procedure to set up Womersley flow with known Re\text{Re} and α\alpha:
1.
Select initial values for ℓ∗\ell^* and w∗w^* which are compatible with the Reynolds number, e.g. w∗>0.1Rew^* > 0.1 \text{Re} (stability).
2.
Estimate the sound propagation time scale ts∗∼ℓ∗/cs∗t_\mathrm{s}^* \sim \ell^*/c_\mathrm{s}^* and therefore the recommended maximum for ω∗\omega^* via 2π/ω∗=T∗≫ts∗2 \pi / \omega^*=T^* \gg t_\mathrm{s}^*. ⇒ It requires some trial and error to find a sufficient minimum ratio T∗/ts∗T^* / t_\mathrm{s}^*, but one should at least use a factor of ten.
3.
Balance ω∗\omega^* and ν∗\nu^* to match α\alpha via α=ωνw\alpha = \sqrt{\frac{\omega}{\nu}}w by taking into account the law of similarity for the Womersley number: α(l)=α(p)\alpha^{(\mathrm{l})} = \alpha^{(\mathrm{p})}. ⇒ There is no unique way to choose ω∗\omega^* and ν∗\nu^*, but keep in mind already now that the relaxation parameter τ∗\tau^* should not be too close to 1/21/2 or much larger than unity.
4.
Estimate the velocity u^0∗\hat{u}_0^* from Re\text{Re} using the selected values for w∗w^* and ν∗\nu^*.
5.
Check the validity of u^0∗\hat{u}_0^*:
•
If its value is not in the desired range (too large for reasonable stability or accuracy or too small for a feasible time step), the parameters have to be re-balanced while keeping Re\text{Re} and α\alpha invariant. ⇒ This process is more complicated than for the Poiseuille flow because Re\text{Re} and α\alpha have to be considered simultaneously.
6.
The required pressure difference Δp0∗\Delta p_0^* finally can be obtained from:
u^0=Δp0w28ℓρν.\hat{u}_0 = \frac{\Delta p_0 w^2}{8 \ell \rho \nu}.
Example:

5.3.5 Surface Tension and Gravity

•
Considering a droplet of a liquid in vapour (or the other way around: a vapour bubble in a liquid).
◦
The surface tension of the liquid-gas interface is γ\gamma.
◦
The droplet/bubble may be put on a flat substrate or into a narrow capillary.
•
The density difference of the liquid and vapour phases is defined as Δρ=ρl−ρv>0\Delta \rho = \rho_{\mathrm{l}} - \rho_{\mathrm{v}}>0.
•
In the presence of the gravitational acceleration with magnitude gg, we can define the dimensionless Bond number:
Bo=Δρgr2γ.\text{Bo} = \frac{\Delta\rho gr^2}{ \gamma}.
◦
It characterizes the relative strength of gravity and surface tension effects where rr is a length scale typical for the droplet/bubble.
◦
A common definition for rr is the radius of a sphere with the same volume VV as the droplet/bubble:
r=(3V4π)1/3.r = \left( \frac{3V}{4\pi} \right)^{1/3}.
•
Small Bond number: gravity negligible, surface tension dominates (spherical shapes).
•
Large Bond number: gravity important, droplet/bubble deformed.
•
To find the simulation parameters in lattice units for a given Bond number, first set the lattice densities:
ρv∗=1λρl∗=:1λρ∗.\rho_{\mathrm{v}}^* = \frac{1}{\lambda} \rho_{\mathrm{l}}^* =: \frac{1}{\lambda} \rho^*.
◦
Where λ>1\lambda > 1 is the achievable density ratio of the model.
•
The next step is to choose the radius r∗r^* properly.
◦
The droplet radius must be significantly larger than interface width: r∗≫d∗r^* \gg d^*, typically r∗≥8r^* \geq 8.
•
From Bond number scaling:
γ∗g∗≈ρ∗r∗2Bo(l).\frac{\gamma^*}{g^*} \approx \frac{\rho^* r^{*2}}{\mathrm{Bo}^{(\mathrm{l})}}.
◦
Where we have approximated Δρ≈ρl=ρ\Delta \rho \approx \rho_\mathrm{l}=\rho, which is valid for λ≫1\lambda \gg 1 & Bo(l)=Bo(p)\text{Bo}^{(\mathrm{l})}=\text{Bo}^{(\mathrm{p})}. → All quantities on the right-hand-side are known.
◦
One degree of freedom remains to choose γ∗\gamma ^* and g∗g^* satisfying the constraint.
•
Must consider model limitations on surface tension range and stability constraints on gravity/velocity.
•
Time step effectively defined by g∗g^* via:
g=g∗ΔxΔt2,Δx=rr∗.g = g^* \frac{\Delta x}{\Delta t^2}, \quad \Delta x = \frac{r}{r^*}.
◦
The user has to consider the total runtime of the simulation as well.
◦
Reducing γ∗\gamma^* and g∗g^* will decrease the time step and result in a longer computing time.
Example:

5.4 Summary

•
Purely mechanical systems require exactly three conversion factors.
◦
First define three independent (basic) ones, then derive all others according to Cq=C1q1C2q2C3q3C_q = C_1^{q_1} C_2^{q_2} C_3^{q_3}.
•
The law of similarity is the physical basis for consistent non-dimensionalization. → Systems with same dimensionless parameters are similar.
◦
Unit systems must not be mixed to avoid inconsistencies and wrong physical results.
•
Typical basic conversion factors for LB: length, time (or velocity), and density. → CℓC_\ell, CtC_t (or CuC_u), and CρC_\rho.
•
Key viscosity relation: v=cs∗2(τ∗−12)Δx2Δtv = c_\mathrm{s}^{*2} \left( \tau^* - \frac{1}{2} \right) \frac{\Delta x^2}{\Delta t} constrains simulation parameters.
•
Choose parameters considering intrinsic LB restrictions for accuracy, stability, and efficiency.
◦
Diffusive scaling Δt∝Δx2\Delta t \propto \Delta x^2 preferred for incompressible flows (second-order spatial accuracy).
◦
Acoustic scaling Δt∝Δx\Delta t \propto \Delta x when speed of sound is physically relevant.
◦
For incompressible simulations, Mach number scaling unnecessary. → Increase for efficiency.
•
The relation of simulation parameters in terms of Reynolds and Mach numbers:
Re(l)=ℓ∗U∗ν∗=ℓ∗U∗cs∗2(τ∗−12)⟹Re(l)Ma(l)=ℓcs∗(τ∗−12)Δx.\mathrm{Re}^{(\mathrm{l})} = \frac{\ell^* U^*}{\nu^*} = \frac{\ell^* U^*}{c_\mathrm{s}^{*2} \left( \tau^* - \frac{1}{2} \right)} \quad \Longrightarrow \quad \frac{\mathrm{Re}^{(\mathrm{l})}}{\mathrm{Ma}^{(\mathrm{l})}} = \frac{\ell}{c_\mathrm{s}^* \left( \tau^* - \frac{1}{2} \right) \Delta x}.
•
Set average lattice density ρ0∗=1\rho_0^* = 1 as default.
•
LB simulations possible without specifying conversion factors. → Map to physical system later.
•
Ensure local density variation ρ′∗\rho'^* small compared to ρ0∗\rho_0^* for incompressibility.
•
Number of independent parameters reduced by each imposed dimensionless constraint (Buckingham π theorem).
•
All relevant time scales must exceed acoustic time scale ts∗∼ℓ∗/cs∗t_\mathrm{s}^* \sim \ell^*/c_\mathrm{s}^* for proper physics.
•
LB most suitable for intermediate Reynolds numbers: O(1)O(1) to O(100)O(100).