Lattice Boltzmann models based on the vielbein formalism for the simulation of flows in curvilinear geometries

In this paper, we consider the Boltzmann equation with respect to orthonormal vielbein fields in conservative form. This formalism allows the use of arbitrary coordinate systems to describe the space geometry, as well as of an adapted coordinate system in the momentum space, which is linked to the physical space through the use of vielbeins. Taking advantage of the conservative form, we derive the macroscopic equations in a covariant tensor notation, and show that the hydrodynamic limit can be obtained via the Chapman-Enskog expansion in the Bhatnaghar-Gross-Krook approximation for the collision term. We highlight that in this formalism, the component of the momentum which is perpendicular to some curved boundary can be isolated as a separate momentum coordinate, for which the half-range Gauss-Hermite quadrature can be applied. We illustrate the capabilities of this formalism by considering two applications. The first one is the circular Couette flow between rotating coaxial cylinders, for which benchmarking data are available for all degrees of rarefaction, from the hydrodynamic to the ballistic regime. The second application concerns the flow in a gradually expanding channel. We employ finite-difference lattice Boltzmann models based on half-range Gauss-Hermite quadratures for the implementation of diffuse reflection, together with the fifth-order weighted essentially nonoscillatory and third-order total variation diminishing Runge-Kutta numerical methods for the advection and time stepping, respectively.

Figure 1

Circular Couette flow setup. (a) Cartesian grid and momentum space decomposition along the Cartesian axes; (b) cylindrical grid and momentum space decomposition along the Cartesian axes; (c) cylindrical grid and momentum space decomposition adapted to the curvilinear coordinates.

Figure 2

Effect of grid stretching on 16 points between $R_{\mathrm{in}}=1$ and $R_{\mathrm{out}}=2$. (a) The parameter $\delta$ controls the positioning of the stretching center (i.e., the point where the grid is the coarsest). (b) The parameter $A$ contains the amplitude of the stretching, with $A=0$ and 1 corresponding to equidistant and infinitely stretched points, respectively.

Figure 3

Circular Couette flow setup.

Figure 4

Density profile $n\left(R\right)$ for static cylinders having radii $R_{\mathrm{in}}=1$ and $R_{\mathrm{out}}=2$. The numerical results were obtained using the $\mathrm{H}(4;5)\times \mathrm{H}(4;5)$ model with $\tau =\mathrm{Kn}/n$, where $\mathrm{Kn}={10}^{-3}$. In all cases, $N_R=16$ nodes were used and the stretching was performed according to $\delta =0.5$ and $A=0.95$.

Figure 5

Numerical results obtained using the $\mathrm{H}(4;5)\times \mathrm{H}(4;5)$ models and the $\chi$ implementation for various values of the angular velocity $\mathrm{\Omega }$ in the case of the rigid rotation considered in Sec. 5c3. (a) Density profile compared to Eq. (5.27). (b) Azimuthal velocity compared to the analytic solution $u^{\widehat{\phi }}=\mathrm{\Omega }R$.

Figure 6

Comparison of the $\tilde{f}$ and $\chi$ formulations for $\mathrm{\Omega }=0.5$, using the $\mathrm{H}(4;5)\times \mathrm{H}(4;5)$ model.

Figure 7

Azimuthal velocity profile $u^{\widehat{\phi }}\left(R\right)$ for the angular velocity of the inner cylinder of ${\mathrm{\Omega }}_{\mathrm{in}}=0.01$. The curves correspond to various values of $\beta =R_{\mathrm{in}}/R_{\mathrm{out}}$. Our numerical results, obtained using the $\mathrm{H}(2;3)\times \mathrm{H}(2;3)$ model in the $\chi$ implementation, are overlapped with the analytic solution (5.30).

Figure 8

Comparison between the numerical and analytic results for the profiles of (a) $n$ [the analytical curve is obtained numerically by solving Eq. (5.35)]; (b) $u^{\widehat{\phi }}$ [Eq. (5.30)]; (c) $T$ [Eq. (5.32)], together with Eq. (5.37) giving the position of the maximum in the temperature profile; (d) $q^{\widehat{R}}$ [Eq. (5.33)]. The curves correspond to various values of $\beta =R_{\mathrm{in}}/R_{\mathrm{out}}$. The inner cylinder rotates with ${\mathrm{\Omega }}_{\mathrm{in}}=0.5$, while the outer cylinder is kept at rest. The numerical results, obtained with the $\chi$ implementation using the $\mathrm{H}(4;5)\times \mathrm{H}(4;5)$ model, are overlapped with the analytic solutions.

Figure 9

Effect of lattice spacing on the behavior of $q^{\widehat{R}}$ near the wall when the $\chi$ implementation on an equidistant grid is employed. The inset shows that the amplitude of the deviation of $q^{\widehat{R}}$ with respect to the expected analytic value (5.33) is proportional to ${\left(\delta R\right)}^{0.58}$.

Figure 10

Comparison of the $\tilde{f}$ (stretched and equidistant) and $\chi$ (stretched and equidistant) radial heat flux using $N_R=32$ grid points.

Figure 11

Trajectory of a free-streaming particle originating from the inner cylinder, which passes through point $\mathcal{P}$ at distance $R$ from the symmetry axis with momentum $\mathit{p}$. Since the particle travels downwards, $\theta <0$ and $\varphi >0$.

Figure 12

Effect of the quadrature order on the oscillations in the stationary profile of the pressure and azimuthal velocity in the ballistic regime. The results were obtained with the models $\mathrm{HH}(4;Q_R)\times \mathrm{HH}(4;10)$ and the curves correspond to various values of $Q_R$. The parameters employed in these simulations are $R_{\mathrm{in}}=1$, $R_{\mathrm{out}}=2$, ${\mathrm{\Omega }}_{\mathrm{in}}=1.0$, ${\mathrm{\Omega }}_{\mathrm{out}}=0$, and $\delta t=2\times {10}^{-5}$. Our simulations reached the stationary state after $500000$ iterations. We used $N_R=16$ grid points together with the stretching given by $A=0.95$ and $\delta =0.75$. The thick continuous curves correspond to the analytic solutions derived in Sec. 5e2.

Figure 13

Comparison between our numerical results and the analytic predictions in the ballistic regime. (a) $P=nT$ [Eqs. (5.50) and (5.61)]; (b) $u^{\widehat{\phi }}/u_{\mathrm{in}}^{\widehat{\phi }}$ [Eq. (5.54)]; (c) $q^{\widehat{R}}$ [Eq. (5.62)]; (d) $q^{\widehat{\phi }}$ [Eq. (5.62)]. The radii of the inner and outer cylinders were $R_{\mathrm{in}}=1$ and $R_{\mathrm{out}}=2$, such that $\beta =0.5$. The quadrature used was $\mathrm{HH}(4;200)\times \mathrm{HH}(4;10)$. The curves correspond to various values of ${\mathrm{\Omega }}_{\mathrm{in}}$. The analytic solution for $u^{\widehat{\phi }}$ reported by Willis [121] and reproduced in Eq. (5.56) is shown in (b) alongside the exact expression (5.54). The time step was set to $\delta t=2\times {10}^{-5}$ and the number of nodes was $N_R=16$.

Figure 14

Same as Fig. 13 in the case of (a) $T^{\widehat{R}\widehat{R}}$ [Eq. (5.60a)]; (b) $T^{\widehat{\phi }\widehat{\phi }}$ [Eq. (5.60b)]; (c) $T^{\widehat{R}\widehat{\phi }}$ [Eq. (5.59)]; (d) $T$ [Eq. (5.61)]. The simulation parameters are the same as in Fig. 13.

Figure 15

Comparison between our simulation results and those reported in Ref. [85]. (a) $n\left(R\right)$; (b) $T\left(R\right)$; (c) $u^{\widehat{\phi }}\left(R\right)/u_{\mathrm{in}}^{\widehat{\phi }}$. The angular velocity of the inner cylinder was set to ${\mathrm{\Omega }}_{\mathrm{in}}=0.5\sqrt{2}$, while the Knudsen number is related to the relaxation time through Eq. (5.65). The models employed in these simulations are summarized in the non-negligible Mach section of Table 3.

Figure 16

Comparison of our simulation results for the constants $Q$ (a) and $T_{\mathrm{in}}^{\widehat{R}\widehat{\phi }}$ (b) against the analytic solutions given in Eqs. (5.34) and (5.31) in the hydrodynamic limit, as well as in Eqs. (5.63) and (5.59) in the free molecular flow limit. The Knudsen number is related to the relaxation time through Eq. (5.65). In (b), the analytic result reported by Willis [121] is represented using a solid black line.

Figure 17

Velocity profile in the low Mach number regime ${\mathrm{\Omega }}_{\mathrm{in}}=0.01$. Our simulation results are compared with the LB results reported by Watari in Ref. [43]. The models employed in these simulations are summarized in the low Mach section of Table 3 alongside the corresponding time step. The number of grid points was $N_R=16$, stretched according to Eq. (3.8) with $A=0.95$ and $\delta =0.5$.

Figure 18

Time (in seconds) required to achieve steady state using the models summarized in Table 3.

Figure 19

Number of millions of site updates per second in the context of the circular Couette flow for a system with $N_R=128$ nodes in the radial direction when the half-range and full-range Gauss-Hermite quadratures of orders $Q_R$ and $Q_{\phi }$ are employed on the radial and azimuthal directions, respectively. The total number of velocities is $N_{\mathrm{vel}}=2Q_R{Q}_{\phi }$. The curves correspond to the cases when $Q_R=Q_{\phi }$; when $Q_R=4$ is kept fixed and $Q_{\phi }$ is varied; and when $Q_{\phi }=4$ is kept fixed and $Q_R$ is varied. The simulation data are represented using lines with points, while the solid lines correspond to Eqs. (5.70), (5.71), and (5.72), where the parameters $a$, $b$, and $c$ are obtained using a fitting routine.

Figure 20

(a) Geometry of the gradually expanding channel for ${\mathrm{Re}}_c=100$. The vertical lines correspond to constant values of $\xi$ chosen equidistantly between $-10$ and 40. The horizontal lines are drawn at constant values of $\lambda$ which are stretched towards the boundaries via Eq. (6.20). (b) Fluid domain in $(\lambda ,\xi )$ coordinates, corresponding to the upper half of the channel.

Figure 21

The orientation of the principal axes of the momentum space expressed with respect to the vielbein. The vectors ${\mathit{p}}_{\widehat{\lambda }}=e_{\widehat{\lambda }}$ and ${\mathit{p}}_{\widehat{\xi }}=e_{\widehat{\xi }}$ can be regarded as unit vectors along these axes. The vertical lines correspond to constant values of $\xi$, while the horizontal lines represent lines of constant values of $\lambda$. The geometry corresponds to setting ${\mathrm{Re}}_c=6$ in Eq. (6.24).

Figure 22

Streamlines for the flow through the gradually expanding channel corresponding to ${\mathrm{Re}}_c=100$, obtained for $Q_0=0.1$ and $\mathrm{Kn}=0.001$ ($\mathrm{Re}=100$), corresponding to the incompressible hydrodynamic limit.

Figure 23

Simulation results for the flow through the gradually expanding channel with ${\mathrm{Re}}_c=100$, $Q_0=0.1$, and $\mathrm{Kn}=0.001$. (a) Normalized wall pressure $\mathrm{\Delta }{P}_{\mathrm{w}}$ [Eq. (6.49)] with respect to the normalized coordinate $y/L_c=3y/{\mathrm{Re}}_c$ along the channel, validated against the results reported by Cliffe [126]. (b) Normalized streamwise velocity $u^y/Q_0$ at $\xi =100/12\simeq 8.33$, validated against the results reported by Roache [97].

Figure 24

Streamlines for the flow through the gradually expanding channel corresponding to ${\mathrm{Re}}_c=100$, obtained for $(Q_0,\mathrm{Kn})\in \{(0.1,0.001);(0.5,0.005);(1.2,0.012)\}$, such that $\mathrm{Re}=100$, highlighting the outermost contour of the vortex. Only the region around the vortex in the upper half of the channel is represented.

Figure 25

Numerical results for the gradually expanding channel flow for (a) normalized local particle flow rate $Q\left(x\right)/Q_0$ and (b) temperature $T$ across the channel at $\xi ={\mathrm{Re}}_c/12$, as well as (c) the normalized wall pressure difference $\mathrm{\Delta }{P}_{\mathrm{w}}$ [Eq. (6.49)] against the normalized streamwise coordinate $y/L_c=3y/{\mathrm{Re}}_c$, at various values of $\mathrm{Kn}$. The particle flow rate is varied according to $Q_0=100\mathrm{Kn}$ in order to maintain $\mathrm{Re}=100$ for all simulations.

Figure 26

Gradually expanding channel flow results for $\mathrm{\Delta }{P}_{\mathrm{w}}$ [Eq. (6.49)] in the hydrodynamic (a) and slip-flow (b) regimes with respect to the normalized downstream coordinate $y/L_c=3y/{\mathrm{Re}}_c$, as well as (c) for $-dP/dy$ in the upstream ($y/L_c<0$) and downstream ($y/L_c>0.5$) regions. The hydrodynamic limit curves $-dP/dy\simeq 0.317\mathrm{Kn}$ (upstream, $\varphi \simeq -0.018$) and $-dP/dy\simeq 0.0385\mathrm{Kn}$ (downstream, $\varphi \simeq 0.982$) are obtained from Eq. (6.50). The linearized Boltzmann-BGK results for the pressure-driven Poiseuille flow are represented with red dotted lines and are computed using Eq. (6.51) using the values for $G_{\mathrm{P}}^{*}$ reported in Refs. [130, 132]. The half-channel particle flow rate is taken as $Q_0=0.1$.

Figure 27

Numerical results for the normalized local particle flow rate $Q\left(x\right)/Q_0$ across the channel at $\xi \simeq 8.33$ in the hydrodynamic (a) and slip-flow (b) regimes for the gradually expanding channel flow corresponding to ${\mathrm{Re}}_c=100$ in Eq. (6.24). The half-channel particle flow rate is taken as $Q_0=0.1$ for various values of $\mathrm{Kn}$, such that the resulting Reynolds number decreases as $\mathrm{Kn}$ is increased.

Figure 28

Numerical results for the normalized vorticity $-\omega /Q_0$ as a function of (a) $x$; and (b) $d$ [Eq. (6.58)], taken at $y={\mathrm{Re}}_c/16\simeq 8.33$, where ${\mathrm{Re}}_c=100$ defines the channel geometry through Eq. (6.24). The half-channel particle flow rate is taken as $Q_0=0.1$ for various values of $\mathrm{Kn}$, such that the resulting Reynolds number of the flow decreases as $\mathrm{Kn}$ is increased.

Figure 29

Numerical results for the temperature $T$(a) and normalized vorticity $-\omega /Q_0$(b) across the channel at $\xi \simeq 33$ for the gradually expanding channel flow corresponding to ${\mathrm{Re}}_c=100$ in Eq. (6.24). The half-channel particle flow rate is taken as $Q_0=0.1$ for various values of $\mathrm{Kn}$, such that the resulting Reynolds number decreases as $\mathrm{Kn}$ is increased.

Figure 30

Comparison of the simulation results for the normalized wall pressure $\mathrm{\Delta }{P}_{\mathrm{w}}$ [Eq. (6.49)] obtained using the VLB model $\mathrm{HH}(4;6)\times \mathrm{H}(4;5)$ (solid line) and the CLB model $\mathrm{HH}(4;6)\times \mathrm{HH}(4;6)$ (line and points). The results are overlapped.

Figure 31

Comparison of (a) the normalized mass flow rate $Q/Q_0$ and (b) the normalized vorticity $-\omega /Q_0$ at $\lambda =\{0.997,0.982,0.963,0.939,0.910\}$ obtained using the VLB model $\mathrm{HH}(4;6)\times \mathrm{H}(4;5)$ (solid lines) and the CLB model $\mathrm{HH}(4;6)\times \mathrm{HH}(4;6)$ (lines and points).

Figure 32

Convergence study of the normalized vorticity with respect to quadrature orders $Q_{\lambda }$ (a) and $Q_{\xi }$ (b) for the VLB implementation at $y/L_c\simeq 0.25$.

Figure 33

Comparison of the VLB and CLB implementations. Convergence study of the normalized vorticity for the CLB implementation with respect to quadrature order by (a) steadily increasing the quadrature order on both axes at $y/L_c\simeq 0.041$ and (b) keeping $Q_x$ fixed and varying $Q_y$ at $y/L_c\simeq 0.154$.

Figure 34

Number of millions of site updates per second in the context of the gradually expanding channel flow for a system with $N_{\lambda }\times N_{\xi }=30\times 200$ nodes when the VLB $\mathrm{HHLB}(4;Q)\times \mathrm{HLB}(4;5)$ (lines and squares) and CLB $\mathrm{HHLB}(4;Q)\times \mathrm{HHLB}(4;Q)$ (lines and circles) models are employed. The solid lines correspond to Eqs. (6.65) and (6.66), where the parameters $a$, $b$, and $c$ are obtained using a fitting routine.