Maximum stream-function comparison for the (r, z) with the proposed results in .
According to the expansion in Eq. (12), every time scale is grouped starting with the terms of order O(1)
followed by terms of order O(ε)
and finally the terms of order O(ε2)
where (ei)α is the component of the velocity vector ei on the α-coordinate direction.
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 [14–17, 24]. The D2Q5 model has discrete velocity directions given by Eq. (16)
Considering the evolution equations as the main purpose of our discretization, the equilibrium functions are defined for uθ and ωθ as
where is the speed of sound in the lattice. These equilibrium functions fulfill the lattice constraints: , and for uθ. The same holds for ωθ replacing fi by gi.
In order to recover the governing equations, zeroth and first moments are taken to Eqs. (14) and (15). Defining for l = 1, 2 and the fact that for k ≥ 1 (see reference ) the zeroth moment of Eq. (14) produces Eq. (19), while the zeroth moment of Eq. (15) produces Eq. (20).
which is replaced into Eq. (20) and produces
Finally, the source terms are redefined in order to eliminate the time derivative of the source within the evolution equation .
The same procedure is applied to ωθ replacing fi by gi obtaining an equation similar to Eq. (26) but in terms of ωθ.
Comparing Eqs. (26) and (27) with Eqs. (2) and (3) it is observed that the terms and have to be neglected in order to match both momentum equations. An order analysis is done for these terms assuming Uc, Lc, and tc as the characteristic velocity, length, and time scales, respectively. Considering Eq. (27) for azimuthal vorticity, the term is the same order of and the term is the same order of . Taking the ratio of the order of the latter terms, we obtain
where M = Uc/cs is the Mach number of the lattice. Eq. (28) shows that the term is very small compared with 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 .
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.
Consequently, the source terms for uθ are defined as
and for ωθ as
where and .
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(ε2) 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 is eliminated in Eqs. (29) and (30) since .
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).
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 . is the weight coefficient for the source terms, and they must satisfy the constrain . In the present study, the source term S is defined in Eq. (35)
and the equilibrium distribution is defined as in Eq. (36)
where . ur and uz are calculated with Eq. (5), employing a central difference scheme in the lattice domain.
With every evolution equation discretized in the LBM framework, an algorithm of the method is finally proposed:
Numerical parameters definition: Re, Ω, R, Ar, N = (Grid size).
Initial and boundary conditions definition
Time loop until steady state is reached:
LBM solver → ψt + Δt until:
d. Computation of ur and uz.
As shown in the algorithm the time loop uses a unified Lattice-Boltzmann solver in which five steps are performed:
Collision step for every particle function distribution.
4. Streaming step for every particle function distribution.
5. Variable recover through the summation of the distribution functions, i.e., , and .
As stated in the algorithm, the time loop is performed until the steady state is reached which numerically is considered when
where the relative error of the velocity field is calculated between 1000 consecutive time steps.
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.
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:
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 Rin = 0.5m and angular velocity Ωin = 0.2rad/s, and outer cylinder’s radius Rout = 1m and angular velocity Ωout = 0. With these parameters, the analytic solution for this flow is given by Eq. (39)
Results were obtained with a lattice resolution of Δx = (Rout − Rin)/(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.
A Taylor-Couette flow consists of a viscous fluid confined between two concentric rotating cylinders of length H = Ar(Rout − Rin) with aspect ratio Ar = 3.8. The Reynolds number is defined as Re = RinΩ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
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. . There is a good agreement between the present formulation and the hybrid scheme demonstrating the versatility of the proposed numerical method.
|Re||ψ max||ψ max |
|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|
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.  are observed.
Cylindrical cavity steady flow have been studied both numerically [9, 18] and experimentally . 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  was able to summarize the flow regimes in the Escudier diagram varying the Reynolds number Re = R2Ω/ν 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.
The parameters in the simulation were set to R = 1 and Ω = 0.1, which makes the characteristic velocity Uc = ΩR = 0.1 ensuring that . The boundary conditions for the primitive variables are defined as
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|
In order to verify the precision of the present method, the maximum axial velocity is compared with experimental data  and to numerical results proposed previously [3–5, 9] in steady state. Three of these numerical results were taken under consideration: Zhou  in which a LBM in p-v is employed, Bhaumik et al.  results in which a 3DLBM using MRT is employed, and the improved model results proposed by Li et al. . Also a DNS solution is considered .
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 (upresent − uexp)/uexp. 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.
The streamlines shown in previous figures are in agreements with those presented in mentioned references.
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.
© 2017 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
786total chapter downloads
Login to your personal dashboard for more detailed statistics on your publications.Access personal reporting
Edited by Hector Perez-De-Tejada
By Ildebrando Pérez‐Reyes, René Osvaldo Vargas‐Aguilar, Eduardo Valente Gómez‐Benítez and Iván Salmerón‐Ochoa
Edited by Hector Perez-De-Tejada
By Nick Verhelst and Jacques Tempere