Simulation of Axisymmetric Flows with Swirl in Vorticity- Stream Function Variables Using the Lattice Boltzmann Method

Abstract

In the present work, a Lattice Boltzmann formulation in vorticity-stream function variables is proposed for axisymmetric flows with swirl. For this purpose, several source terms are proposed and implemented. Although containing velocity gradients, these sources are in the Lattice Boltzmann framework and fulfill the Euler and Navier-Stokes equations in their conservative form. The main characteristics of the proposed method are: First, the momentum equation is solved using a unified Lattice Boltzmann solver; second, the proposed sources are consistent with the nonviscous and viscous momentum equations; and third, the implemented method is second-order accurate in space. Numerical tests on the Taylor-Couette flow with finite aspect ratio of 3.8 and the lid-driven cylindrical cavity flow were carried out showing good agreement with numerical and experimental results found in the literature, evidencing the ability of the implemented method to solve axisymmetric flows with swirl. In the case of the lid-driven cylindrical cavity flow, the implemented method is able to correctly reproduce some qualitative characteristics of this flow such as the vortex breakdown close to the cavity axis at different Reynolds numbers and cavity aspect ratio.

1. Introduction

The Lattice Boltzmann method (LBM) was created in the late 1980s as a derivation of the Lattice Gas Automata (LGA). This method has shown to be an efficient solver not only for the Navier-Stokes (NS) equations but also for some other nonlinear partial differential equations [1]. In several cases, LBM has been used for the solution of axisymmetric flows, modeling the problem in two dimension (2D) [2–7] and three dimension (3D) [8–10].

Although it is well known that axisymmetric flows can be mathematically described as a 2D problem, considering the governing equations in a cylindrical coordinate form, there is an inherent discrepancy between the cylindrical behavior of the flow and the two dimensional spatial discretization (i.e., the type of lattice used). In order to overcome these discrepancies, Halliday et al. [11] included position and time-dependent sources in the Lattice Boltzmann evolution equation to achieve the proper evolution equations through the Chapman-Enskog methodology. Subsequent studies have successfully used this idea to improve the Lattice Boltzmann method in order to solve different axisymmetric flows with or without swirl. Huang et al. [12] proposed a hybrid Lattice Boltzmann finite-difference axisymmetric formulation where the planar velocities were solved through a pressure-velocity (p − v) LBM. The azimuthal component through a finite-difference scheme was solved by inserting source terms into the two-dimensional Lattice Boltzmann equation. In order to avoid the use of different frameworks and solving every momentum equation within the LBM approach, Guo et al. [13] proposed a simple and straightforward incorporation of source terms in a LB p − v formulation. This way, they were able to predict steady and unsteady axisymmetric flows starting from the general Lattice Boltzmann equation. Li et al. [4] proposed an improved axisymmetric Lattice Boltzmann scheme where a multiple relaxation time (MRT) collision model was used to insert source terms that contained no velocity gradients. The same approach was also explained and implemented by Zhou [5]. Regarding LBM formulations in vorticity-stream function, Chen et al. [14] developed an axisymmetric Lattice Boltzmann formulation without swirl considering vorticity and the stream-function as its primitive variables. Improvements were made coupling a thermal LBM [15] and finally proposing an axisymmetric formulation [16, 17].

In the present work, different source terms are proposed for the Lattice Boltzmann implementation of axisymmetric flows with swirl. Although containing velocity gradients, these sources are in the Lattice Boltzmann framework and fulfill the Euler and NS equations in their conservative form. This implementation is tested with some flows such as the lid-driven cylindrical cavity flow: a benchmark case for axisymmetric flow solvers, deeply studied both numerically [18–20] and experimentally [21]. The main characteristic of the lid-driven cylindrical cavity flow is that for certain Reynolds number (Re) and cavity aspect ratio (Ar), it can undergo structural changes such as vortex breakdowns that are triggered by azimuthal vorticity stretching and its interaction with the azimuthal velocity [20, 22].

2. Axisymmetric Lattice Boltzmann implementation

2.1. Governing equations

Let $\overrightarrow{v}\left(r,z,t\right)=u_r{\widehat{e}}_r+u_{\theta }{\widehat{e}}_{\theta }+u_z{\widehat{e}}_z$ be the velocity field of an axisymmetric, viscous flow with swirl with the corresponding vorticity ($\overrightarrow{\omega }$) defined as shown in Eq. (1)

$\overrightarrow{\omega }=-\frac{\partial u_{\theta }}{{\partial }_z}{\widehat{e}}_r+{\omega }_{\theta }\left(r,z,t\right){\widehat{e}}_{\theta }+\frac{1}{r}\frac{\partial \left(r{u}_{\theta }\right)}{\partial r}{\widehat{e}}_z$E1

For such flows, the 3D Navier-Stokes equations are equivalent to the following simplified vorticity-stream function formulation

$\frac{\partial u_{\theta }}{\partial t}+\nabla \cdot \left(\overrightarrow{u}{u}_{\theta }\right)=-2u_r\frac{u_{\theta }}{r}+\nu \left[\frac{{\partial }^2{u}_{\theta }}{\partial r^2}+\frac{{\partial }^2{u}_{\theta }}{\partial z^2}+\frac{1}{r}\frac{\partial u_{\theta }}{\partial r}-\frac{u_{\theta }}{r^2}\right]$E2

$\frac{\partial {\omega }_{\theta }}{\partial t}+\nabla \cdot \left(\overrightarrow{u}{\omega }_{\theta }\right)=2\frac{u_{\theta }}{r}\frac{\partial u_{\theta }}{\partial z}+\nu \left[\frac{{\partial }^2{\omega }_{\theta }}{\partial r^2}+\frac{{\partial }^2{\omega }_{\theta }}{\partial z^2}+\frac{1}{r}\frac{\partial {\omega }_{\theta }}{\partial r}-\frac{{\omega }_{\theta }}{r^2}\right]$E3

${\omega }_{\theta }=-\frac{1}{r}\frac{{\partial }^2\psi }{\partial z^2}-\frac{\partial }{\partial r}\left(\frac{1}{r}\frac{\partial \psi }{\partial r}\right)$E4

$u_r=-\frac{1}{r}\frac{\partial \psi }{\partial z}{u}_z=\frac{1}{r}\frac{\partial \psi }{\partial r}$E5

where $\overrightarrow{u}=\left(u_r,u_z\right)$. The definition of the velocity in terms of the stream function ψ fulfills the continuity equation $\nabla \cdot \overrightarrow{u}=0$. Note that the convective terms of these governing equations are written in a conservative form in order to match the operators achieved by the numerical discretization based on LBM.

2.2. Numerical method

The discrete Lattice Boltzmann equation (LBE) is given by

$f_i\left(\mathbf{x}+{\mathbf{e}}_{\mathbf{i}}\Delta t,t+\Delta t\right)-f_i\left(\mathbf{x},t\right)={\Omega }_i\left(f_i\left(\mathbf{x},t\right)\right)$E6

where f_i is the particle distribution function along the ith direction, e_i is a vector in the direction of the microscopic velocities and Ω_i(f_i(x, t)) is the collision operator. Δx and Δt are space and time increments, and Δx/Δt = |e_i| = c is the magnitude of the microscopic velocity. Employing a second-order Taylor expansion on the convective part (LHS of Eq. (6)) and using the BGK approximation of the collision operator Ω_i(f_i(x, t)), Eq. (7) is achieved.

$\left({\partial }_t+\nabla \cdot {\mathbf{e}}_{\mathbf{i}}\right)f_i+\frac{1}{2}\left({\partial }_t^2+2{\partial }_t\nabla \cdot {\mathbf{e}}_{\mathbf{i}}+\nabla \nabla :{\mathbf{e}}_{\mathbf{i}}{\mathbf{e}}_{\mathbf{i}}\right)f_i=-\frac{1}{\tau }\left[f_i-f_i^{\mathit{eq}}\right]$E7

where τ is the dimensionless relaxation time of the distribution function f_i and $f_i^{\mathit{eq}}$ is the equilibrium function distribution. Using the Chapman-Enskog expansion, the distribution function f_i is expanded as

$f_i=f_i^{\left(0\right)}+\epsilon f_i^{\left(1\right)}+{\epsilon }^2{f}_i^{\left(2\right)}+\dots$E8

where ε is a formal parameter in the expansion that allows to keep track of different orders of magnitude. It will be considered only as a label and will be dropped out of the final results setting ε = 1 [23]. The time and space derivatives are also expanded in terms of ε as shown in Eqs. (9) and (10).

${\partial }_t=\epsilon {\partial }_t^{\left(1\right)}+{\epsilon }^2{\partial }_t^{\left(2\right)}$E9

${\partial }_x=\epsilon {\partial }_x^{\left(1\right)}$E10

According to Wolf-Gladrow [23], the reasons behind the different expansions in time and space lie in the fact that different macroscopic processes such as convection and diffusion can be distinguished by their time scales but act on similar spatial scales. Replacing Eq. (8) through Eq. (10) in Eq. (7), Eq. (11) is obtained

$D_i\left(f_i^{\left(0\right)}+\epsilon f_i^{\left(1\right)}\right)+\frac{1}{2}{D}_i^2\left(f_i^{\left(0\right)}+\epsilon f_i^{\left(1\right)}\right)=-\frac{1}{\tau }\left[f_i^{\left(0\right)}+\epsilon f_i^{\left(1\right)}-f_i^{\mathit{eq}}\right]$E11

where $D_i=\epsilon {\partial }_t^{\left(1\right)}+{\epsilon }^2{\partial }_t^{\left(2\right)}+\epsilon {\partial }_x^{\left(1\right)}\cdot {\mathbf{e}}_{\mathbf{i}}$ is the total derivative operator expanded through ε. Finally, source terms (h_i) are included in the RHS of Eq. (11) in order to fulfill the momentum equations constraints where $h_i^{\left(1\right)}$ and $h_i^{\left(2\right)}$ are the expansion of h_i in terms of ε.

$D_i\left(f_i^{\left(0\right)}+\epsilon f_i^{\left(1\right)}\right)+\frac{1}{2}{D}_i^2\left(f_i^{\left(0\right)}+\epsilon f_i^{\left(1\right)}\right)=-\frac{1}{\tau }\left[f_i^{\left(0\right)}+\epsilon f_i^{\left(1\right)}-f_i^{\mathit{eq}}\right]+\epsilon h_i^{\left(1\right)}+{\epsilon }^2{h}_i^{\left(2\right)}$E12

According to the expansion in Eq. (12), every time scale is grouped starting with the terms of order O(1)

$f_i^{\left(0\right)}=f_i^{\mathit{eq}}$E13

followed by terms of order O(ε)

${\partial }_t^{\left(1\right)}{f}_i^{\left(0\right)}+\left({\partial }_x^{\left(1\right)}\cdot {\mathbf{e}}_{\mathbf{i}}\right)f_i^{\left(0\right)}=-\frac{1}{\tau }{f}_i^{\left(1\right)}+h_i^{\left(1\right)}$E14

and finally the terms of order O(ε²)

$\begin{array}{l}{\partial }_t^{\left(2\right)}{f}_i^{\left(0\right)}+\frac{1}{2}{\left({\partial }_t^{\left(1\right)}\right)}^2{f}_i^{\left(0\right)}+{\partial }_t^{\left(1\right)}\left({\partial }_x^{\left(1\right)}\cdot {\mathbf{e}}_{\mathbf{i}}\right)f_i^{\left(0\right)}+\frac{1}{2}\left({\partial }_x^{\left(1\right)}{\partial }_x^{\left(1\right)}:{\left({\mathbf{e}}_{\mathbf{i}}\right)}_{\alpha }{\left({\mathbf{e}}_{\mathbf{i}}\right)}_{\beta }\right)f_i^{\left(0\right)}\\ +{\partial }_t^{\left(1\right)}{f}_i^{\left(1\right)}+\left({\partial }_x^{\left(1\right)}\cdot {\mathbf{e}}_{\mathbf{i}}\right)f_i^{\left(1\right)}=\frac{1}{\tau }{f}_i^{\left(2\right)}+h_i^{\left(2\right)}\text{.}\end{array}$E15

where (e_i)_α is the component of the velocity vector e_i on the α-coordinate direction.

2.3. Lattice

The D2Q5 lattice model (two dimensions and five directions) has shown to be adequate for advection-diffusion problems based on its easy implementation and its inherent orthogonality $\left({\mathbf{e}}_{\mathbf{i}}{\mathbf{e}}_{\mathbf{j}}={\delta }_j^i\right)$ [14–17, 24]. The D2Q5 model has discrete velocity directions given by Eq. (16)

${\mathbf{e}}_{\mathbf{i}}=c\times \left(\begin{array}{ll}\overrightarrow{0}\hfill & i=0,\hfill \\ \left[\mathit{cos}\left(i-1\right)\pi /2,\mathit{sin}\left(i-1\right)\pi /2\right]\hfill & i\ne 0\text{.}\hfill \end{array}\right.$E16

Considering the evolution equations as the main purpose of our discretization, the equilibrium functions are defined for u_θ and ω_θ as

$f_i^{\mathit{eq}}=\frac{u_{\theta }}{5}\left[1+\frac{\overrightarrow{u}\cdot {\mathbf{e}}_{\mathbf{i}}}{c_s^2}\right]$E17

$g_i^{\mathit{eq}}=\frac{{\omega }_{\theta }}{5}\left[1+\frac{\overrightarrow{u}\cdot {\mathbf{e}}_{\mathbf{i}}}{c_s^2}\right]$E18

where $c_s=\left|c\right|\sqrt{2/5}$ is the speed of sound in the lattice. These equilibrium functions fulfill the lattice constraints: ${\displaystyle {\sum }_{i=0}^4}{f}_i=u_{\theta }$, ${\displaystyle {\sum }_{i=0}^4}{f}_i{\mathbf{e}}_{\mathbf{i}}=u_{\theta }\overrightarrow{u}$ and ${\displaystyle {\sum }_{i=0}^4}{\mathbf{e}}_{\mathbf{i}}{\mathbf{e}}_{\mathbf{i}}{f}_i=c_s^2{u}_{\theta }$ for u_θ. The same holds for ω_θ replacing f_i by g_i.

2.4. Recovery of the governing equations

In order to recover the governing equations, zeroth and first moments are taken to Eqs. (14) and (15). Defining $H^{\left(l\right)}={\displaystyle {\sum }_{i=0}^4}{h}_i^{\left(l\right)}$ for l = 1, 2 and the fact that ${\displaystyle {\sum }_{i=0}^4}{f}_i^{\left(k\right)}=0$ for k ≥ 1 (see reference [24]) the zeroth moment of Eq. (14) produces Eq. (19), while the zeroth moment of Eq. (15) produces Eq. (20).

${\partial }_t^{\left(1\right)}{u}_{\theta }+{\partial }_x^{\left(1\right)}\cdot \left(\overrightarrow{u}{u}_{\theta }\right)=H^{\left(1\right)}$E19

$\begin{array}{l}{\partial }_x^{\left(1\right)}\cdot \left({\displaystyle \sum _{i=0}^4}{f}_i^{\left(1\right)}{\mathbf{e}}_{\mathbf{i}}\right)+\frac{1}{2}\left(2{\partial }_t^{\left(2\right)}{u}_{\theta }+{\partial }_t^{\left(1\right)}{\partial }_t^{\left(1\right)}{u}_{\theta }\right.\\ \left.+2{\partial }_t^{\left(1\right)}{\partial }_x^{\left(1\right)}\cdot \left(\overrightarrow{u}{u}_{\theta }\right)+c_s^2{\partial }_x^{\left(1\right)}{\partial }_x^{\left(1\right)}{u}_{\theta }\right)=H^{\left(2\right)}\end{array}$E20

The first term on the LHS of Eq. (20) is rewritten with the first moment of Eq. (14)

${\displaystyle \sum _{i=0}^4}{f}_i^{\left(1\right)}{\mathbf{e}}_{\mathbf{i}}=\tau \left({\displaystyle \sum _{i=0}^4}{h}_i^{\left(1\right)}{\mathbf{e}}_{\mathbf{i}}-{\partial }_t^{\left(1\right)}\left(\overrightarrow{u}{u}_{\theta }\right)-c_s^2{\partial }_x^{\left(1\right)}{u}_{\theta }\right)$E21

which is replaced into Eq. (20) and produces

$\begin{array}{l}\tau {\partial }_x^{\left(1\right)}\cdot \left({\displaystyle \sum _{i=0}^4}{h}_i^{\left(1\right)}{\mathbf{e}}_{\mathbf{i}}-{\partial }_t^{\left(1\right)}\left(\overrightarrow{u}{u}_{\theta }\right)-c_s^2{\partial }_x^{\left(1\right)}{u}_{\theta }\right)+\frac{1}{2}\left(2{\partial }_t^{\left(2\right)}{u}_{\theta }\right.\\ \left.+{\partial }_t^{\left(1\right)}{\partial }_t^{\left(1\right)}{u}_{\theta }+2{\partial }_t^{\left(1\right)}{\partial }_x^{\left(1\right)}\left(\overrightarrow{u}{u}_{\theta }\right)+c_s^2{\partial }_x^{\left(1\right)}{\partial }_x^{\left(1\right)}{u}_{\theta }\right)=H^{\left(2\right)}\end{array}$E22

The second time derivative over u_θ on the LHS of Eq. (22), i.e., ${\partial }_t^{\left(1\right)}{\partial }_t^{\left(1\right)}{u}_{\theta }$, is rewritten by taking ∂_t on Eq. (19) as follows:

${\partial }_t^{\left(1\right)}\left({\partial }_t^{\left(1\right)}{u}_{\theta }\right)={\displaystyle \sum _{i=0}^4}{\partial }_t^{\left(1\right)}{h}_i^{\left(1\right)}-{\partial }_x^{\left(1\right)}\cdot \left(\overrightarrow{u}{u}_{\theta }\right)$E23

Replacing Eq. (23) into Eq. (22) leads to

$\begin{array}{l}\tau {\partial }_x^{\left(1\right)}\cdot \left({\displaystyle \sum _{i=0}^4}{h}_i^{\left(1\right)}{\mathbf{e}}_{\mathbf{i}}-{\partial }_t^{\left(1\right)}\left(\overrightarrow{u}{u}_{\theta }\right)-c_s^2{\partial }_x^{\left(1\right)}{u}_{\theta }\right)+\frac{1}{2}\left(2{\partial }_t^{\left(2\right)}{u}_{\theta }+{\displaystyle \sum _{i=0}^4}{\partial }_t^{\left(1\right)}{h}_i^{\left(1\right)}\right.\\ \left.-{\partial }_t^{\left(1\right)}{\partial }_x^{\left(1\right)}\left(\overrightarrow{u}{u}_{\theta }\right)+2{\partial }_t^{\left(1\right)}{\partial }_x^{\left(1\right)}\left(v{u}_{\theta }\right)+c_s^2{\partial }_x^{\left(1\right)}{\partial }_x^{\left(1\right)}{u}_{\theta }\right)=H^{\left(2\right)}\end{array}$E24

Finally, the source terms are redefined in order to eliminate the time derivative of the source within the evolution equation [25].

$h_i=\epsilon h_i^{\left(1\right)}+{\epsilon }^2\left(h_i^{\left(2\right)}+\frac{1}{2}{\partial }_{t1}{h}_i^{\left(1\right)}\right)$E25

Equation (25) combined with Eqs. (19) and (24) produces the momentum equation for u_θ

$\begin{array}{l}{\partial }_t{u}_{\theta }+\epsilon {\partial }_x^{\left(1\right)}\cdot \left(\overrightarrow{u}{u}_{\theta }\right)+{\epsilon }^2\left(\tau {\partial }_x^{\left(1\right)}\cdot {\displaystyle \sum _{i=0}^4}{h}_i^{\left(1\right)}{\mathbf{e}}_{\mathbf{i}}+\right.\\ \left(\frac{1}{2}-\tau \right)\left(\left.{\partial }_t^{\left(1\right)}{\partial }_x^{\left(1\right)}\cdot \left(\overrightarrow{u}{u}_{\theta }\right)+c_s^2{\partial }_x^{\left(1\right)}{\partial }_x^{\left(1\right)}{u}_{\theta }\right)\right)=\epsilon H^{\left(1\right)}+{\epsilon }^2{H}^{\left(2\right)}\end{array}$E26

The same procedure is applied to ω_θ replacing f_i by g_i obtaining an equation similar to Eq. (26) but in terms of ω_θ.

$\begin{array}{l}{\partial }_t{\omega }_{\theta }+\epsilon {\partial }_x^{\left(1\right)}\cdot \left(\overrightarrow{u}{\omega }_{\theta }\right)+{\epsilon }^2\left(\tau {\partial }_x^{\left(1\right)}\cdot {\displaystyle \sum _{i=0}^4}{h}_i^{\left(1\right)}{\mathbf{e}}_{\mathbf{i}}+\right.\\ \left(\frac{1}{2}-\tau \right)\left(\left.{\partial }_t^{\left(1\right)}{\partial }_x^{\left(1\right)}\cdot \left(\overrightarrow{u}{\omega }_{\theta }\right)+c_s^2{\partial }_x^{\left(1\right)}{\partial }_x^{\left(1\right)}{\omega }_{\theta }\right)\right)=\epsilon H^{\left(1\right)}+{\epsilon }^2{H}^{\left(2\right)}\end{array}$E27

Comparing Eqs. (26) and (27) with Eqs. (2) and (3) it is observed that the terms ${\partial }_t^{\left(1\right)}{\partial }_x^{\left(1\right)}\cdot \left(\overrightarrow{u}{u}_{\theta }\right)$ and ${\partial }_t^{\left(1\right)}{\partial }_x^{\left(1\right)}\cdot \left(\overrightarrow{u}{\omega }_{\theta }\right)$ have to be neglected in order to match both momentum equations. An order analysis is done for these terms assuming U_c, L_c, and t_c as the characteristic velocity, length, and time scales, respectively. Considering Eq. (27) for azimuthal vorticity, the term ${\partial }_t^{\left(1\right)}{\partial }_x^{\left(1\right)}\cdot \left(\overrightarrow{u}{\omega }_{\theta }\right)$ is the same order of $U_c/t_c^2{L}_C$ and the term $c_s^2{\partial }_x^{\left(1\right)}{\partial }_x^{\left(1\right)}{\omega }_{\theta }$ is the same order of $c_s^2/L_c^2{t}_c$. Taking the ratio of the order of the latter terms, we obtain

$O\left(\frac{{\partial }_t^{\left(1\right)}{\partial }_x^{\left(1\right)}\cdot \left(\overrightarrow{u}{\omega }_{\theta }\right)}{c_s^2{\partial }_x^{\left(1\right)}{\partial }_x^{\left(1\right)}{\omega }_{\theta }}\right)=O\left(\frac{U_c/t_c^2{L}_c}{c_s^2/L_c^2{t}_c}\right)=O\left(\frac{U_c^2}{c_s^2}\right)=O\left(M^2\right)$E28

where M = U_c/c_s is the Mach number of the lattice. Eq. (28) shows that the term ${\partial }_t^{\left(1\right)}{\partial }_x^{\left(1\right)}\cdot \left(\overrightarrow{u}{\omega }_{\theta }\right)$ is very small compared with $c_s^2{\partial }_x^{\left(1\right)}{\partial }_x^{\left(1\right)}{\omega }_{\theta }$ and it can be neglected if M < < 1, according to the LBM dynamics. This procedure is also valid for the azimuthal velocity leading to neglect the term ${\partial }_t^{\left(1\right)}{\partial }_x^{\left(1\right)}\cdot \left(\overrightarrow{u}{u}_{\theta }\right)$.

After the order analysis is performed, Eqs. (26) and (27) are rewritten as

$\begin{array}{l}{\partial }_t{u}_{\theta }+\epsilon {\partial }_x^{\left(1\right)}\cdot \left(\overrightarrow{u}{u}_{\theta }\right)+{\epsilon }^2\left(\tau {\partial }_x^{\left(1\right)}\cdot {\displaystyle \sum _{i=0}^4}{h}_i^{\left(1\right)}{\mathbf{e}}_{\mathbf{i}}+\right.\\ \left(\frac{1}{2}-\tau \right)\left.c_s^2{\partial }_x^{\left(1\right)}{\partial }_x^{\left(1\right)}{u}_{\theta }\right)=\epsilon H^{\left(1\right)}+{\epsilon }^2{H}^{\left(2\right)}\end{array}$E29

$\begin{array}{l}{\partial }_t{\omega }_{\theta }+\epsilon {\partial }_x^{\left(1\right)}\cdot \left(\overrightarrow{u}{\omega }_{\theta }\right)+{\epsilon }^2\left(\tau {\partial }_x^{\left(1\right)}\cdot {\displaystyle \sum _{i=0}^4}{h}_i^{\left(1\right)}{\mathbf{e}}_{\mathbf{i}}+\right.\\ \left(\frac{1}{2}-\tau \right)\left.c_s^2{\partial }_x^{\left(1\right)}{\partial }_x^{\left(1\right)}{\omega }_{\theta }\right)=\epsilon H^{\left(1\right)}+{\epsilon }^2{H}^{\left(2\right)}\end{array}$E30

2.5. Source definition

As it was stated in the introduction, there exists a discrepancy between the lattice dimension and the dimensional nature of the flow. The discrepancy arises in the operators achieved through the multiscale analysis (Cartesian) and those that are natural to the momentum equations (cylindrical). As shown in the RHS of Eq. (31), there is an additional derivative that contains the swirl of the flow which is not captured by the operators achieved through the multiscale analysis.

${\frac{D}{Dt}}_{\mathit{car}}=\frac{\partial }{\partial t}+u_r\frac{\partial }{\partial r}+u_z\frac{\partial }{\partial z}\ne \frac{\partial }{\partial t}+u_r\frac{\partial }{\partial r}+u_z\frac{\partial }{\partial z}+u_{\theta }\frac{\partial }{\partial \theta }={\frac{D}{Dt}}_{cyl}$E31

In order to overcome this problem, the inclusion of a source terms is needed and therefore defined to match the governing equations, Eqs. (2, 3) with Eqs. (29, 30).

Consequently, the source terms for u_θ are defined as

$h_i^{\left(1\right)}\text{≡}-t_{i}2u_r\frac{u_{\theta }}{r}\mathrm{and}{h}_i^{\left(2\right)}\text{≡}{t}_i\left(\tau -\frac{1}{2}\right)c_s^2\left[\frac{1}{r}\frac{\partial u_{\theta }}{\partial r}-\frac{u_{\theta }}{r^2}\right]$E32

and for ω_θ as

$h_i^{\left(1\right)}\text{≡}{t}_{i}2\frac{u_{\theta }}{r}\frac{\partial u_{\theta }}{\partial z}\mathrm{and}{h}_i^{\left(2\right)}\text{≡}{t}_i\left(\tau -\frac{1}{2}\right)c_s^2\left[\frac{1}{r}\frac{\partial {\omega }_{\theta }}{\partial r}-\frac{{\omega }_{\theta }}{r^2}\right]$E33

where ${\displaystyle {\sum }_{i=0}^4}{t}_i=1$ and $\tau =\nu /c_s^2+1/2$.

The source terms of O(ε) were defined, both for u_θ and ω_θ, in order to reproduce the Euler equations in their conservative form. Then, the terms of O(ε²) were defined in order to reproduce the cylindrical terms of the Laplacian operator that appears in Eqs. (2) and (3). Finally, due to the symmetries of the lattice, the term ${\partial }_x^{\left(1\right)}\cdot {\displaystyle {\sum }_{i=0}^4}{h}_i^{\left(1\right)}{\mathbf{e}}_{\mathbf{i}}$ is eliminated in Eqs. (29) and (30) since ${\displaystyle {\sum }_{i=0}^4}{h}_i^{\left(1\right)}{\mathbf{e}}_{\mathbf{i}}=0$.

3. LBM algorithm

This section describes the algorithm used to solve every Lattice Boltzmann evolution equation within the same framework producing a LBM solver able to solve axisymmetric flows. The implementation of the algorithm is based on the key steps in LBM: streaming and collision that are given by Eq (6).

3.1. Poisson equation solver

In order to solve Eq. (4), which is a Poisson equation for ψ, the model proposed by Chai et al. [26] is employed. The evolution equation is given by

$L_i\left(t+\Delta t,x+\Delta x\right)-L_i\left(t,x\right)=-\frac{1}{{\tau }_{\psi }}\left[L_i\left(t,x\right)-L_i^{\mathit{eq}}\left(t,x\right)\right]+S\Delta t{\overline{\psi }}_i{c}^2\left(0.5-{\tau }_{\psi }\right)/2$E34

where L is the distribution function associated with ψ and S accounts for the source terms. τ_ψ is the dimensionless relaxation time that is set with accurate results to 1 [26]. ${\overline{\psi }}_i$ is the weight coefficient for the source terms, and they must satisfy the constrain ${\displaystyle {\sum }_{i=1}^4}{\overline{\psi }}_i=1$. In the present study, the source term S is defined in Eq. (35)

$S\text{≡}-r{\omega }_{\theta }+uz$E35

and the equilibrium distribution is defined as in Eq. (36)

$L_i^{\mathit{eq}}\left(x,t\right)=\left\{\begin{array}{lll}-\psi \left(x,t\right),\hfill & :\hfill & i=0\hfill \\ \frac{1}{4}\psi \left(x,t\right),\hfill & :\hfill & i=1,2,3,4\hfill \end{array}\right.$E36

where $\psi \left(x,t\right)={\displaystyle {\sum }_{i=1}^4}{L}_i\left(x,t\right)$. u_r and u_z are calculated with Eq. (5), employing a central difference scheme in the lattice domain.

3.2. Algorithm

With every evolution equation discretized in the LBM framework, an algorithm of the method is finally proposed:

1. Numerical parameters definition: Re, Ω, R, Ar, N = (Grid size).

2. Initial and boundary conditions definition

3. Time loop until steady state is reached:

   1. $u_{\theta }^t\to$ LBM solver $\to u_{\theta }^{t+\Delta t}$

   2. $u_{\theta }^{t+\Delta t}\to$ LBM solver $\to {\omega }_{\theta }^{t+\Delta t}$

   3. ${\omega }_{\theta }^{t+\Delta t}\to$ LBM solver → ψ^t + Δt until:

$\frac{D\psi }{Dt}={\nabla }^2\psi +{\omega }_{\theta }\mathrm{tends}\mathrm{t}\mathrm{o}\frac{D\psi }{Dt}=0\to {\nabla }^2\psi =-{\omega }_{\theta }$E37

  d. Computation of u_r and u_z.

3.2.1. LBM solver

As shown in the algorithm the time loop uses a unified Lattice-Boltzmann solver in which five steps are performed:

1. Equilibrium function calculation through Eqs. (17, 18, and 36).

2. Source term calculations through Eqs. (32, 33, and 35) using the information obtained at time t

3. Collision step for every particle function distribution.

${\Omega }_i\left(f_i\left(\mathbf{x},t\right)\right)=-\frac{1}{\tau }\left[f_i-f_i^{\mathit{eq}}\right]$ E38

 4. Streaming step for every particle function distribution.

$f_i\left(\mathbf{x}+{\mathbf{e}}_{\mathbf{i}}\Delta t,t+\Delta t\right)={\Omega }_i\left(f_i\left(\mathbf{x},t\right)\right)+f_i\left(\mathbf{x},t\right)$ E390

 5. Variable recover through the summation of the distribution functions, i.e., ${\displaystyle {\sum }_{i=0}^4}{f}_i=u_{\theta }$, ${\displaystyle {\sum }_{i=0}^4}{g}_i={\omega }_{\theta }$ and ${\displaystyle {\sum }_{i=1}^4}{L}_i=\psi$.

As stated in the algorithm, the time loop is performed until the steady state is reached which numerically is considered when

$\frac{\left|\left|u\left(t\right)-u\left(t-1000\delta t\right)\right|\right|}{\left|\left|u\left(t\right)\right|\right|}<{10}^{-7}$E37

where the relative error of the velocity field is calculated between 1000 consecutive time steps.

4. Numerical results

In the present section, the proposed source terms are validated using some well-known benchmarks flows, including the circular Couette flow, the Taylor-Couette flow, and the swirling flow within the lid-driven cylindrical cavity. All cases were validated for a laminar regime. For each case, the boundary and initial conditions will be discussed.

4.1. Circular Couette flow

In this case, the flow between two infinitely long concentric cylinders is simulated. The inner cylinder rotates at constant speed Ω, while the outer is stationary (see Figure 1). The analytic solution of this flow is used to prove that the proposed method is second order. The boundary conditions for the fluid variables are as follows:

Figure 1.

Configuration of the Circular Couette and Taylor-Couette Flow.

$\begin{array}{l}r=R_{\mathit{in}}:\left(u_r,u_{\theta },u_z\right)=\left(0,\Omega R_{\mathit{in}},0\right),{\omega }_{\theta }={\nabla }^2\psi ,\psi =0,\\ r=R_{\mathit{out}}:u_r=u_{\theta }=u_z=0,{\omega }_{\theta }={\nabla }^2\psi ,\psi =0\text{.}\end{array}$E38

where the Laplacian of ψ on the boundaries is calculated using a second-order Taylor approximation employing the inner nodes values. Symmetry is imposed on the top and bottom boundaries.

Figure 2 compares the analytic solution with the numerical results for a laminar flow with Re < 10, inner cylinder’s radius R_in = 0.5m and angular velocity Ω_in = 0.2rad/s, and outer cylinder’s radius R_out = 1m and angular velocity Ω_out = 0. With these parameters, the analytic solution for this flow is given by Eq. (39)

Figure 2.

Azimuthal velocity comparison of the Circular Couette flow. (−): Analytical solution, (o): present LBM solution.

$u_{\theta }=Ar+\frac{B}{r}$E39

where

$A=\frac{{\Omega }_{\mathit{in}}{R}_{\mathit{in}}^2-{\Omega }_{\mathit{out}}{R}_{\mathit{out}}^2}{R_{\mathit{in}}^2-R_{\mathit{out}}^2}=-0.0667\mathrm{and}=R_{\mathit{in}}^2{R}_{\mathit{out}}^2\frac{{\Omega }_{\mathit{in}}-{\Omega }_{\mathit{out}}}{R_{\mathit{out}}^2-R_{\mathit{in}}^2}=0.0667\text{.}$E400

Results were obtained with a lattice resolution of Δx = (R_out − R_in)/(N − 1), i.e., N = 50. It is clear that the numerical results are in good agreement with the analytical solution. The relative global error, defined by Eq. (40), is presented in Figure 3 for different mesh sizes.

$E=\left|\right|u_a-u_{LBM}{\left|\right|}_2$E40

In Eq. (40) u_LBM is the azimuthal velocity predicted by the present method. The slope of the fitting in Figure 3 is 2.04, which shows that the proposed method is second-order accurate in space.

Figure 3.

Convergence analysis. (−): least-square fitting with slope 2.04.

4.2. Taylor-Couette flow with finite aspect ratio

A Taylor-Couette flow consists of a viscous fluid confined between two concentric rotating cylinders of length H = Ar(R_out − R_in) with aspect ratio Ar = 3.8. The Reynolds number is defined as Re = R_inΩ_inA/ν, where Ω_in is the angular velocity of the inner cylinder and A is the gap of the annulus. In this case, the boundary conditions for z, due to the finite length, have to be specified, besides the boundary conditions for r used in the Circular Couette flow, as

$\begin{array}{l}z=0:\psi =u_z=0,{\omega }_{\theta }={\nabla }^2\psi ,\\ z=H:\psi =u_z=0,{\omega }_{\theta }={\nabla }^2\psi \text{.}\end{array}$E41

Three different Re were simulated and analyzed; Re = 85, 100 and 150. The maximum stream-function values in the r − z plane are listed in Table 1 and compared to those presented by Huang et al. [12]. There is a good agreement between the present formulation and the hybrid scheme demonstrating the versatility of the proposed numerical method.

Re	ψ max	ψ max [12]
85	4.32 × 10^− 2	4.810 × 10^− 2
	5.252 × 10^− 2	5.501 × 10^− 2
	6.38 × 10^− 2	6.427 × 10^− 2

Table 1.

Maximum stream-function comparison for the (r, z) with the proposed results in [12].

In order to validate the qualitative characteristic of the flow in terms of axisymmetric toroidal vortices, Figure 4 shows the contours of stream-function and vorticity for Re = 150 and Ar = 3.8. Similar flow patterns consistent with those reported by Huang et al. [12] are observed.

Figure 4.

Contours of stream-function (left) and vorticity (right).

4.3. Lid-driven cylindrical cavity flow

Cylindrical cavity steady flow have been studied both numerically [9, 18] and experimentally [21]. One of the interesting features of this flow is that vortex breakdowns takes place within the cavity producing recirculating zones located in the cavity axis. In 1984, Escudier [21] was able to summarize the flow regimes in the Escudier diagram varying the Reynolds number Re = R²Ω/ν and the aspect ratio Ar = H/R of the cavity. The flow problem consists of a cylinder with top and bottom walls, where the top wall rotates at a constant angular velocity (see Figure 5). Four cases were chosen from the Escudier diagram: Case 1 (Ar = 1.5, Re = 990), Case 2 (Ar = 1.5, Re = 1290), Case 3 (Ar = 2.5, Re = 1010), and Case 4 (Ar = 2.5, Re = 2200) in order to demonstrate the quantitative and qualitative accuracy of the present formulation.

Figure 5.

Configuration of the cylindrical cavity flow.

The parameters in the simulation were set to R = 1 and Ω = 0.1, which makes the characteristic velocity U_c = ΩR = 0.1 ensuring that $M=U_c/c_s^2=0.25<<1$. The boundary conditions for the primitive variables are defined as

$\begin{array}{l}r=0:u_r=u_{\theta }=0\mathrm{and}\frac{\partial u_z}{\partial r}=0,{\omega }_{\theta }=0,\psi =0\\ r=R:u_r=u_{\theta }=u_z=0,{\omega }_{\theta }={\nabla }^2\psi ,\psi =0\\ z=0:u_r=u_{\theta }=u_z=0,{\omega }_{\theta }={\nabla }^2\psi ,\psi =0\\ z=H:u_r=u_z=0,u_{\theta }=\Omega r,{\omega }_{\theta }={\nabla }^2\psi ,\psi =0\end{array}$E42

Figure 6.

Normalized axial velocity u_z of Case 1 depending on mesh size: (×): 200 × 300, (−): 150 × 225, (⋅): 100 × 150, (⋅ − ⋅ −): 70 × 105, (− − −): 50 × 75, (⊲): 30 × 45.

In order to establish mesh independence in the solution, different grid sizes were used on the r − z domain. Case 1 was simulated using the following meshes: 30 × 45, 50 × 75, 70 × 105, 100 × 150, 150 × 225, 200 × 300, and the results are shown in Figure 6. It is observed that the solution is independent if the grid size used is larger or equal to 150 × 225. Based on this fact, it is believed that a lattice size of Δx = R/150 = 0.0067 will produce accurate solutions and will be used to simulate the other cases.

	Re = 990		Re = 1290		Re = 1010
Reference	u z,max/u0	x max /H	u z,max/u0	x max /H	u z,max/u0	x max /H
Present	0.0901	0.2098	0.0727	0.1696	0.1109	0.498
Expt.[27]	0.097	0.21	0.068	0.14	0.103	0.46
Zhou [5]	0.0992	0.207	0.0706	0.147	0.105	0.448
DLBM [9]	0.093	0.22	0.072	0.16	0.102	0.52
DNS [9]	0.099	0.19	0.0665	0.125	0.106	0.44
Li [4]	0.0987	0.213	0.0716	0.147

Table 2.

Comparisons of maximum axial velocities.

In order to verify the precision of the present method, the maximum axial velocity is compared with experimental data [27] and to numerical results proposed previously [3–5, 9] in steady state. Three of these numerical results were taken under consideration: Zhou [5] in which a LBM in p-v is employed, Bhaumik et al. [9] results in which a 3DLBM using MRT is employed, and the improved model results proposed by Li et al. [4]. Also a DNS solution is considered [9].

Figure 7.

Stream function contours for Case 1: Re = 990 and Ar = 1.5 (y-coordinate corresponds to Z and x-coordinate to R).

Figure 8.

Stream function contours for Case 2: Re = 1290 and Ar = 1.5 (y-coordinate corresponds to Z and x-coordinate to R).

The maximum axial velocities for Cases 1, 2, and 3 are shown in Table 2, comparing the solution of the present model with previous reported results. The relative error is calculated as (u_present − u_exp)/u_exp. Case 1 presents a relative error of 7.1%, Case 2 of 6.9%, and Case 3 of 7.6%. Finally, in order to establish the present formulation proficiency of solving complex flows we present in Figures 7–9 stream function contours for three different cases. Figures 8 and 9 clearly show the formation of recirculating bubbles close to the axis of symmetry, known as vortex breakdowns.

Figure 9.

Stream function contours for Case 3: Re = 2200 and Ar = 2.5 (y-coordinate corresponds to Z and x-coordinate to R).

The streamlines shown in previous figures are in agreements with those presented in mentioned references.

5. Conclusion

In the present work, a LBM in vorticity-stream function variables for axisymmetric flow with swirl is presented and implemented. By considering the convective term of the evolution equations in their conservative form, proper sources are carefully proposed being able to reproduce both nonviscous and viscous momentum equations. Through a multiscale analysis, performed in Cartesian coordinates, it was found that a discrepancy between the operators and the governing equations was achieved. This difference is overcome by the definition and inclusion of source terms in the proposed LB formulation. The Chapman-Enskog analysis was used to achieve conservative operators that arise naturally in the implemented LBM. As a result a unified LBM algorithm was built in which every evolution equations is solved with the same algorithm. The numerical method proved to be second-order accurate in space. Finally, the method was able to reproduce complex flows, such as the Taylor-Couette flow where toroidal vortices were observed and good agreement was found with qualitative results proposed in the literature. Furthermore, the proficiency of the method to solve the lid-driven cylindrical cavity flow quantitatively showed an error below 8% when compared with experimental data. Qualitatively, the method solved the flow through many flow regimes observing one and two vortex breakdown located in the cavity axis.