- Generated on Sun Jun 16 2024 23:45:16 for MAIA by
1.9.5
|
MAIA bb96820c
Multiphysics at AIA
|
Discretizing the Boltzmann equation according to Chen1992 and Qian1992 yields the lattice Boltzmann equation
\begin{equation} f_i(\mv{x} + \mv{\xi}_i \Delta t, t + \Delta t) = f_i(\mv{x}, t) + \omega (f_i^{eq}(\mv{x}, t)-f_i(\mv{x}, t)) \,, \label{eq:lattice_boltzmann} \end{equation}
with the discrete PPDF \( f_i \), the time increment \( \Delta t \) and the discrete particle velocity \( \mv{\xi}_i \) in the discrete direction \( i\). For time integration, a simple Euler forward scheme is used.
The discrete Maxwell distribution function is defined by
\begin{equation} f_i^{eq} = t_p \rho \left[ 1 + \frac{\mv{\xi}\cdot\mv{u}}{c_s^2} + \frac{1}{2} \left( \frac{\mv{\xi} \cdot \mv{u}}{2c_s^2} \right)^2 - \frac{\mv{u}^2}{c_s^2} \right] \label{eq:equilibrium} \,, \end{equation}
where \( c_s=1/\sqrt{3} \) is the isothermal speed of sound. The coefficient \( t_p \) is given for the most popular discretization models mentioned by Qian1992 in the table below.
| \( t_p \) | D2Q9 | D3Q19 | D3Q27 |
|---|---|---|---|
| \( \frac{\mv{\xi}_i^2}{\mv{\xi}_0^2}=0 \) | \( \frac{4}{9} \) | \( \frac{1}{3} \) | \( \frac{8}{27} \) |
| \( \frac{\mv{\xi}_i^2}{\mv{\xi}_0^2}=1 \) | \( \frac{1}{9} \) | \( \frac{1}{18} \) | \( \frac{2}{27} \) |
| \( \frac{\mv{\xi}_i^2}{\mv{\xi}_0^2}=2 \) | \( \frac{1}{36} \) | \( \frac{1}{36} \) | \( \frac{1}{54} \) |
| \( \frac{\mv{\xi}_i^2}{\mv{\xi}_0^2}=3 \) | - | - | \( \frac{1}{216} \) |
The denotation DxQy describes the spatial dimension (x) and the number of molecular velocities (y). Molecules are only allowed to move from a node to a neighbouring node ( \( \Delta x \)) in on time step ( \( \Delta t \)), the mean molecular velocity can also be defined as \( \xi_0=\Delta x / \Delta t \).
The algorithm to solve Eq. \(\eqref{eq:lattice_boltzmann}\) cosists of two steps:
Collision step:
\(\tilde{f_i}\left(\mv{x},t\right) = f_i\left(\mv{x},t\right)+ \omega \left(f_i^{eq}\left(\mv{x},t\right)-f_i\left(\mv{x},t\right)\right)\)
Propagation step:
\(f_i\left(\mv{x}+\xi_i\Delta t,t+\Delta t\right)=\tilde{f_i}\left(\mv{x},t\right)\)
Using the Chapman-Enskog expansion, the relation between \( \omega \) and the kinematic shear viscosity \( \nu \) can be derived as
\begin{equation} \nu = c_s^2 \left(\frac{\Delta t}{\omega}-\frac{\Delta t}{2}\right). \label{eq:omega_nu} \end{equation}
Considering the limits of Eq. \(\eqref{eq:omega_nu}\), the stability limits of the LB method are shown. Since the viscosity only has a physical meaning for \( \nu>0 \), the upper limit of the relaxation frequency is \( \omega<2 \). For a small collision frequency, i.e., \( \omega \rightarrow 0 \), Eq. \(\eqref{eq:omega_nu}\) yields a large viscosity \( \nu \rightarrow \inf \). If \( \omega \rightarrow 2 \), the LB method tends to become unstable.
The LB scheme described above is valid for low Mach numbers, and hence the energy equation decouples from the momentum equations. To simulate the temperature distribution, additional sets of PPDFs are solved on the same lattice. Three approaches are implemented for computing the temperature. The basic thermal approach is based on Shan1997 and the temperature distribution function \( g_i \) is calculated by
\begin{equation} g_i(\mv{x} + \mv{\xi}_i \Delta t, t + \Delta t) = g_i(\mv{x}, t) + \omega_g (g_i^{eq}(\mv{x}, t)-g_i(\mv{x}, t)) \,. \label{eq:lattice_boltzmann_thermal} \end{equation}
The thermal relaxation frequency \( \omega_g \) is a function of the thermal diffusivity \( \kappa \)
\begin{equation} \kappa = c_s^2 \left(\frac{\Delta t}{\omega_g}-\frac{\Delta t}{2} \right) \,. \label{eq:omega_g} \end{equation}
The thermal equilibrium distribution function \( g_i^{eq} \) for the basic approach is given by
\begin{equation} g_i^{eq} = t_p T \left[ 1 + \frac{\mv{\xi}\cdot\mv{u}}{c_s^2} + \frac{1}{2} \left( \frac{\mv{\xi} \cdot \mv{u}}{2c_s^2} \right)^2 - \frac{\mv{u}^2}{c_s^2} \right] \label{eq:equilibrium_thermal} \,, \end{equation}
With their internal energy LB method He1998 improved the method of Shan1997. That is, the internal energy \( \epsilon = DRT/2 \) multiplied by \( \rho \) is transported by the internal energy distribution function \( h_i \)
\begin{equation} h_i(\mv{x} + \mv{\xi}_i \Delta t, t + \Delta t) = h_i(\mv{x}, t) + \omega_h (h_i^{eq}(\mv{x}, t)-h_i(\mv{x}, t)) \,. \label{eq:lattice_boltzmann_inner} \end{equation}
The relation between the thermal relaxation frequency \( \omega_h \) to \( \kappa \) differs from Eq. \(\eqref{eq:omega_g}\)
\begin{equation} \kappa = c_s^2 \frac{D+2}{D} \left(\frac{\Delta t}{\omega_h}-\frac{\Delta t}{2} \right) \,. \label{eq:omega_h} \end{equation}
The thermal equilibrium reads
\begin{equation} h_i^{eq} = t_p \rho \epsilon \left[ \frac{\mv{\xi}^2}{Dc_s^2} + \left(\frac{\mv{\xi}^2}{Dc_s^2} - \frac{2}{D} \right) \times \frac{\mv{u} \cdot \mv{\xi}}{c_s^2} + \frac{\mv{u}\cdot\mv{u}}{2c_s^2}\cdot\left(\frac{\mv{\xi} \cdot \mv{\xi}}{c_s^2}-\delta\right)\right] \label{eq:equilibrium_inner} \,, \end{equation}
where \( \delta \) stands for the Kronecker delta.
The third thermal approach considers the total energy \( E = e + 1/2 u^2 \) according to Guo2007, where \( u \) is the velocity magnitude. The total energy distribution function is formulated as
\begin{equation} k_i(\mv{x} + \mv{\xi}_i \Delta t, t + \Delta t) = k_i(\mv{x}, t) + \omega_k (k_i^{eq}(\mv{x}, t)-k_i(\mv{x}, t)) + \\ (\omega_k - \omega) \left(\mv{\xi} \cdot \mv{u} - \frac{\mv{u} \cdot \mv{u}}{2} \right) (f_i(\mv{x}, t)-f_i^{eq}(\mv{x}, t)) \,. \label{eq:lattice_boltzmann_total} \end{equation}
The collision frequency is in this case related to the thermal conductivity by
\begin{equation} \kappa = c_s^2 \left(\frac{\Delta t}{\omega_k}-\frac{\Delta t}{2} \right) \,. \label{eq:omega_k} \end{equation}
The equilibrium distributions are given by
\begin{equation} k_i^{eq} = t_p \rho c_s^2 \left[ \frac{\mv{\xi} \cdot \mv{u}}{c_s^2} + \left( \frac{\mv{\xi}\cdot\mv{u}}{c_s^2} \right)^2 - \frac{\mv{u}^2}{c_s^2} + \frac{1}{2} \left( \frac{\mv{\xi}^2}{c_s^2} - D \right) \right] + E f_i^{eq} \label{eq:equilibrium_total} \,. \end{equation}
The macroscopic variables are calculated by summing over all PPDFs, i.e., the density and velocities by calculating
\begin{equation} \rho(\mv{x},t) = \sum_i f_i(\mv{x},t), \label{eq:lbdensity} \end{equation}
\begin{equation} \rho(\mv{x},t)\mv{u}(\mv{x},t) = \sum_i \xi_i f_i(\mv{x},t). \label{eq:lbvelocity} \end{equation}
For the basic thermal approach, the temperature is computed by
\begin{equation} T(\mv{x},t) = \sum_i g_i(\mv{x},t). \label{eq:lbtemperature} \end{equation}
For the approaches considering the internal and total energy, the temperature is derived from
\begin{equation} \rho(\mv{x},t)\epsilon = \sum_i h_i(\mv{x},t), \label{eq:lbtemperature_inner} \end{equation}
and
\begin{equation} \rho(\mv{x},t)E = \sum_i k_i(\mv{x},t). \label{eq:lbtemperature_total} \end{equation}
The static pressure \( p_s \) is derived from the ideal gas law by
\begin{equation} p_s = \rho RT = \rho c_s^2 = \rho/3 \label{eq:lbpressure_static} \end{equation}
The table below shows the PPDF directions of the D3Q27 model.
|
0 |
1 |
2 |
3 |
4 |
5 |
6 |
7 |
8 |
9 |
10 |
11 |
12 |
13 |
14 |
15 |
16 |
17 |
18 |
19 |
20 |
21 |
22 |
23 |
24 |
25 |
26 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
|
-1 |
1 |
0 |
0 |
0 |
0 |
-1 |
-1 |
1 |
1 |
-1 |
-1 |
1 |
1 |
0 |
0 |
0 |
0 |
-1 |
-1 |
-1 |
-1 |
1 |
1 |
1 |
1 |
0 |
Besides the BGK approach presented above, further collision models are implemented in the LB solver (List). Some of the implemented collision steps are an combination of the normal BGK collision step with a thermal LB extension. Furthermore, the user can choose between incompressible and compressible collision steps, which are still based on the BGK model. The remaining approaches are more complex than the BGK model, but might offers advantages for specific applications. For instance, the MRT model conducts the collision in momentum space. That is, it provides a separate relaxation frequency for each moment. Hence, it is also referred to as multi-relaxation model. Its collision equation reads
\begin{equation} f_i(\mv{x} + \mv{\xi}_i \Delta t, t + \Delta t) = f_i(\mv{x}, t) + \mv{\mv{M}}^{-1} \mv{\mv{S}} (\mv{m}^{eq}(\mv{x}, t)-\mv{m}(\mv{x}, t)) \,. \label{eq:MRT} \end{equation}
Due to the use of one relaxation frequency per moment, the MRT provides increased numerical stability for many simulations. However, the correct adjustment of the relaxation frequencies might be difficult.
The cascaded LB Geier2006 and the cumulant Geier2015 LB are successors of the MRT collision step. The cascased LB performs the collision using the central moments. The cumulant performs the collision in cumulant space. These transformation allows the user to simulate high Reynolds number flows, where the relaxation frequency \( \omega \) is close to the limit of \( 2 \).
In addition to the so far mentioned collision steps, nearly all models are equipped with a Smagorinsky model to perform an LES. To activte the Smagorinsky model the property solverMethod must be adapted.
For the calculation non-dimensional variables and parameters are used. The reference values are the density \(\rho_0\), the mean thermodynamic velocity \(\xi_0\), the grid distance \(\Delta x\) and the temperature \(T_0\). In the code implementation the following non-dimensional variables are used:
| Variable | Equation |
|---|---|
| Velocity | \(u^*=\frac{u}{\xi_0}\) |
| Isothermal speed of sound | \(c_s^*=\frac{c_s}{\xi_0}=\sqrt{\frac{1}{3}}\) |
| Molecular velocity | \(\xi_0^*=\frac{\xi_0}{\xi_0}=1\) |
| Grid distance | \(\Delta x^*=\frac{\Delta x}{\Delta x}=1\) |
| Characteristic Length | \(L^*=\frac{L}{\delta x}\) |
| Time step | \(\Delta t^*=\frac{\Delta t \cdot \xi_0}{\Delta x} = 1\) |
| Density | \(\rho^*=\frac{\rho}{\rho_0}\) |
| Pressure | \(p^*=\frac{p}{\rho_0 \cdot \xi_0^2}\) |
| Collision Operator | \(\Omega^* = \omega \cdot \Delta t = \frac{2}{6 \cdot \nu^* + 1}\) |
| Temperature | \(T^*=\frac{T}{T_0}\) |
| Thermal diffusivity | \(\kappa^*=\frac{\nu^*}{Pr}\) |
| Collision Operator (thermal) | \(\Omega_g^* = \omega_g \cdot \Delta t = \frac{2}{6 \cdot \frac{\nu^*}{Pr} + 1}\) |
The non-dimensional kinematic viscosity \(\nu^*\) can be calculated using the Reynolds similarity ( \(Re = Re^*\)) and the Mach-Number \(Ma^* = \frac{u^*}{c_s^*} = \sqrt{3} \cdot u^*\) (which is given in the properties_run file)
\begin{equation}Re = \frac{L \cdot u}{\nu} = \frac{L^* \cdot u^*}{\nu^*} \Leftrightarrow \nu^* = \frac{u^* \cdot L^*}{Re} = \frac{L}{\sqrt{3} \cdot \Delta x} \cdot \frac{Ma^*}{Re}\end{equation}
Furthermore, it can be shown, that the Mach-Number \(Ma^*\) only affects the temporal accuracy:
\begin{equation}Ma^* = \frac{u^*}{c_s^*} = \sqrt{3} \cdot u^* = \sqrt{3} \cdot \frac{u}{\xi_0} \Leftrightarrow \xi_0 = \frac{\sqrt{3} \cdot u}{Ma^*} \Leftrightarrow \Delta t = Ma^* \cdot \frac{\Delta x}{\sqrt{3} \cdot u}\end{equation}
By determining the time step \(\Delta t\) and afterwards the mean thermodynamic velocity \(\xi_0 = \frac{\Delta x}{\Delta t}\) the velocity \( u \) can be calculated. The density \(\rho\) and the temperature \( T \) can easily be calculated from the non-dimensional variable and the corresponding reference value.
In the following, the most important boundary conditions are explained. A list of all boundary conditions can be found here: List of LB boundary conditions
BC 1000 BC 1060 (Thermal)
Velocity is prescribed, the boundary normal is used as velocity vector, the pressure is extrapolated in opposite direction (von Neumann condition). Define the properties
lbControlInflow = 1
inflowSegmentIds = 1,2,... ;
to prescribe a velocity field in a parabolic manner.
BC 1002 BC 1022 (Thermal)
Velocity is prescribes. Use this for a periodic channel, main direction is X, periodic in Z-direction. channel center (Y-component) is placed at y=0. Profile is given by ($h$ is half of channel height):
\(u(y)=u_{max}\left(h^2-y^2\right)=\frac{3}{2}\frac{\bar{u}}{h^2}\left(h^2-y^2\right)=\frac{3}{2}\frac{\frac{Ma}{\sqrt{3}}}{h^2}\left(h^2-y^2\right)\)
BC 1061 (only thermal)
Prescribes temperature and velocity with a Blasius profile for boundary-layer flows.
BC 1080
Generates a set of boundary conditions for in- and outflows of internal flow configurations such as a pipe. The temporal variation of the volume flow (or the main stream velocity) is reconstructed using Fourier coefficients determined by the FFTW library. The non-dimensional frequency \(\alpha\) (m_alpha) is defined as
\begin{equation}\alpha = \frac{\mathrm{Wo}}{\sqrt{2\pi}} = \frac{\mathtt{m\_diameter}}{\sqrt{\mathtt{m\_nu}\times N}}\end{equation}
where \(N\) is the number of time-steps per one cycle and the Womersley number is \(\mathrm{Wo}=L\sqrt{\omega/\nu}\) at an angular frequency \(\omega\). For a given Reynolds number $\mathrm{Re}=UL/\nu$ and a Mach number \(\mathrm{Ma}=U/\sqrt{RT}\) the number of time steps is
\begin{equation}N = \frac{2\pi \mathrm{Re}}{\mathrm{Wo}^2}\times \mathtt{m\_diameter} \times \frac{\sqrt{3}}{\mathrm{Ma}}\end{equation}
In practice a validation can use the flow configuration in literature (Varghese et al., J. Fluids Mech., Vol. 582, 2007). In September 2019 the simulation using LB (120 and 200 million cells) shows a good agreement with the DNS results of stenotic flows at a Reynolds number 600.
BC 2000
interpolated bounce back with BFL-rule of Bouzidi M, Firdaouss M, Lallemand P. Momentum transfer of a Boltzmann-lattice fluid with boundaries. Physics of Fluids. 2001;13(11):3452-3459.
BC 2001
simple bounce back, no interpolation, wall is located at the cell surface
BC 2002
simple bounce back, no interpolation, wall is located at the cell center
BC 2020
Same as BC 2000, but also sets the temperature to the value specified by the property
For a slip condition specular reflection of the PPDFs is performed at the surface.
BC 302d
Replace d by the direction for which slip should be implemented
+-x +-y +-z
d 0,1 2,3 4,5
BC 3000
Zero gradient in direction of boundary normal
BC 301d
Zero gradient in axes-directions (d: see below)
<b>{=html}Definition of d:</b>{=html} the extrapolation direction is always pointing out of the domain
-x +x -y +y -z +z
d 0 1 2 3 4 5
BC 4000 BC 4110 (thermal)
Pressure is prescribed, the velocity is extrapolated normal to the surface
BC 4071 BC 4081 (thermal)
Pressure is prescribed, the velocity is extrapolated normal to the surface. The prescribed pressure is read from the property <i>{=html}<b>{=html}rho1</b>{=html}</i>{=html}.
BC 4070 BC 4080
Pressure is prescribed according to appearing volume flux and pressure distribution from last iteration (use BC 4080 for thermal simulations)
Using the equation of Saint Vernant and Wanzel, one can obtain the pressure by using
\(p_1 = p_0\left(1 - \frac{\gamma - 1}{2\gamma}\cdot\frac{\rho_0}{p_0}\cdot v_0^2\right)^{\frac{\gamma}{\gamma -1}}\)
A non-dimensionalization with
\(\rho^\ast = \frac{\rho}{\rho_0},\quad c_{s}^{\ast} = \frac{c_s}{\xi_0},\quad v^{\ast} = \frac{v}{\xi_0}\)
and using $p=\rho c_s^2$, $c_s=\frac{1}{\sqrt{3}}\xi_0$ and extending with $\frac{\rho^{\ast^2}}{\rho^{\ast^2}}$ leads to
\(\rho_{i+1}=\left(1 - \frac{\gamma - 1}{2\gamma}\cdot \frac{3}{\rho^{\ast^2}_{i}}\left(\rho^{\ast} v^\ast\right)^2\right)^{\frac{\gamma}{\gamma -1}}\)
So, first all the velocities are extrapolatet to the boundary cells and the new density ratio is calculated.
BC 4072 BC 4082
Pressure is prescribed according to a Reynolds number that should be reached. Define the density via the property
\(\delta\rho = 1.0 - \rho_1,\) which is used for the maximal step size to in- or decrease the density ratio. In each timestep the local Reynolds number is calculated and the desity ratio to be prescribed is adjusted accordingly.
BC 4073
BC 4073 works similar as BC 4072, except that the Reynolds number is measured at a specified location which can be defined in the geometry file by
where P1, P2, P3, N1, N2, and N3 define a plane point and normal.
If you are running a restart define:
To generate a perodic link the corresponding stl files are declared as periodic segment in the geometry file. No bc is needed.
For example:
The property initMethods determines the LBM initalization method. The LB solver provides several initialization conditions. Some of these are quite specific. Others are more general. In the following, only the general ICs are presented. A full list of ICs can be found here: List of initial conditions
1.9.5