Multiple-relaxation-time lattice Boltzmann model for the axisymmetric convection diffusion equation

International Journal of Heat and Mass Transfer 67 (2013) 338–351

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

Multiple-relaxation-time lattice Boltzmann model for the axisymmetric convection diffusion equation Like Li ⇑, Renwei Mei, James F. Klausner Department of Mechanical and Aerospace Engineering, University of Florida, Gainesville, FL 32611-6250, United States

a r t i c l e

i n f o

Article history: Received 27 April 2013 Received in revised form 11 August 2013 Accepted 13 August 2013 Available online 8 September 2013 Keywords: Lattice Boltzmann model Multiple-relaxation-time (MRT) Axisymmetric convection–diffusion equation Non-equilibrium distribution functions

a b s t r a c t A multiple-relaxation-time (MRT) lattice Boltzmann (LB) model with five discrete velocities in two dimensions (D2Q5) for simulating the axisymmetric convection–diffusion equation (CDE) is proposed. The axisymmetric CDE in the cylindrical coordinate system is rearranged as an ordinary CDE in the Cartesian coordinate system with the extra terms considered as additional sources. The gradient of the macroscopic scalar in the source terms is directly obtained from the non-equilibrium components of the distribution functions; thus no finite-difference calculation is required. The proposed model is consistent with the LB models for the axisymmetric hydrodynamic equations and they can be directly coupled. A series of numerical tests are conducted to validate the applicability and accuracy of the proposed model, including: (i) steady convection–diffusion in an annulus with mass injection on the boundaries, (ii) steady and transient convection–diffusion in a circular pipe, (iii) steady heat conduction inside a sphere, (iv) steady natural convection in an annulus between two coaxial vertical cylinders, and (v) steady swirling flows in a vertical cylindrical container. Analytical solutions are available for tests (i)–(iii), and the numerical results show that the proposed model is second-order accurate in space and first-order accurate in time. The surface-averaged Nusselt numbers for test (iv), the flow pattern and the velocity and temperature distributions for test (v) agree well with published results. Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction Heat and mass convection–diffusion in axisymmetric systems such as circular cylinders and annular cavities formed by coaxial cylinders are widely encountered in engineering practices. Examples of applications include heat exchangers, desalination systems, solar energy collectors, electronic cooling devices and chemical reactors. Fluid flows and heat and mass transfer in cylindrical systems are essentially three-dimensional (3-D) problems, while under axisymmetric conditions they can be reduced to quasi-2-D problems, which greatly reduces the computational effort. The lattice Boltzmann equation (LBE) method has become an effective numerical method for the convection–diffusion equation (CDE) [1–9]. Especially, the multiple-relaxation-time (MRT) lattice Boltzmann (LB) model [5,8,9] is robust and retains the spatial second-order accuracy for curved boundaries. Although 3-D LB models can be directly applied to simulate the axisymmetric CDE when the curved boundary condition is properly handled, quasi2-D LB models, which require far less computational effort and bypass the curved boundary treatment for straight pipes, are

attractive alternatives for solving the axisymmetric CDE as demonstrated by the LB models developed for axisymmetric athermal flows [10–20] and thermal problems [21–26]. The standard LB models break the computation into a collision step and a streaming step. In order for axisymmetric LB models to work using the existing LB models developed in 2-D Cartesian coordinate systems in terms of the standard collision-streaming procedure, either the governing equations for momentum and energy transport in the cylindrical coordinate system are represented in a pseudo-Cartesian coordinate system, or the effect of axisymmetry is modeled at the particle distribution function level [17,18,25,26]. The first approach results in additional terms, including the velocity and temperature gradients, in the transformed governing equations, and they can be treated as additional source terms in the respective macroscopic transport equations. This treatment is straightforward since it is carried out at the macroscopic level. For example, the typical axisymmetric CDE is

    @/ @ @ @ @/ @ @/ þ þ G; rDrr Dzz þ ðrur /Þ þ ðuz /Þ ¼ @t r@r @z r@r @r @z @z ð1Þ

⇑ Corresponding author. Tel.: +1 352 328 2963; fax: +1 352 392 1071. E-mail addresses: likelichina@ufl.edu (L. Li), rwmei@ufl.edu (R. Mei), klaus@ ufl.edu (J.F. Klausner). 0017-9310/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijheatmasstransfer.2013.08.039

where G is the typical source term representing possibly the effects of viscous dissipation, pressure work, volumetric heating, heat

L. Li et al. / International Journal of Heat and Mass Transfer 67 (2013) 338–351

released or mass created from chemical reaction for the energy or species transport equation. In Eq. (1), ur and uz are the velocity components in the radial r- and axial z-directions, respectively, Drr and Dzz are the diagonal diffusion coefficients of the diffusion tensor Dij. It can be rewritten as

    @/ @ @ @ @/ @ @/ þ þ Sa1 þ Sa2 þ G; Drr Dzz þ ður /Þ þ ðuz /Þ ¼ @t @r @z @r @r @z @z ur / Drr @/ with Sa1   ; Sa2  : ð2Þ r r @r Eq. (2) becomes an ordinary CDE in the pseudo-Cartesian coordinate system (r, z) with three source terms, Sa1, Sa2, and G, of different origins. In the second approach, new LB evolution equations for the particle distribution functions must be carefully constructed so that the correct macroscopic transport, Eq. (1), can be recovered through the Chapman–Enskog expansion for the new evolution equation. When MRT-LB models [27,5] are desired, such modifications at the microscopic level for the evolution equations of the distribution functions may involve redefining the transformation matrix and the equilibrium distribution functions. The first approach is thus preferred if the gradients can be efficiently obtained from the distribution functions instead of using finite-difference schemes. A review of the various LB models for axisymmetric fluid momentum equations can be found in the recent work of Zhou [20], which also introduced an improved axisymmetric LB model based on his original version in [14]. The main advantages of Zhou’s revised model compared with the other axisymmetric LB models are as follows: (i) the source terms represent the same additional terms in the transformed governing equations and no extra source terms are needed in the evolution equation, (ii) it preserves the original computational procedure for both the distribution functions and the macroscopic variables including density and velocities as those in the conventional LBE method, and (iii) no calculation for velocity gradients is required. It is well known that the MRT-LB models [27,28] have better numerical stability and accuracy than BGK-LB models, i.e., the single-relaxation-time (SRT) collision models based on the work of Bhatnagar et al. [29]. Thus when the fluid momentum equations are involved, an MRT version of Zhou’s axisymmetric LB model [20] is used in the present work. Several axisymmetric thermal LB models, which simulate the axisymmetric CDE for temperature, have also been proposed in the literature [21–26]. Peng et al. [21] and Huang et al. [22] proposed hybrid approaches in which the fluid velocities in the radial and axial directions are simulated using the athermal axisymmetric LB model by Halliday et al. [10], while the azimuthal velocity and the temperature field are solved using finite-difference schemes. Chen et al. [23] proposed an LB model for axisymmtric buoyancy-driven thermal flows in which the vorticity-streamfunction formulation was used for the velocity field so that the governing equations for vorticity and temperature became CDEs with source terms. The vorticity and temperature gradients in the source terms were obtained from finite-difference calculations. In addition, a Poisson equation for the stream function must be solved during each time step in [23]. To avoid the finite-difference computation of the temperature gradient, i.e., Sa2 in Eq. (2), in LB models, Li et al. [24] proposed an implicit evolution equation for the distribution functions to recover the axisymmetric energy equation. Additional body forces can also be modeled by including a forcing term in this model. The implicitness is then eliminated by introducing a new distribution function, which includes not only the original distribution function but also the source term. As a result, the modified evolution equation is more difficult to implement than that in standard LB models. Furthermore, the

339

temperature calculation in [24] requires the radial velocity component and thus it deviates from the standard calculation of macroscopic variables in the LBE method. Zheng et al. [25] extended their athermal axisymmetric LB model [17] to simulate thermal flows by including a D2Q4 LB model for the temperature field. The radial coordinate and the temperature distribution functions were lumped together and the source term in the evolution equation of the temperature distribution functions was properly constructed so that the scalar source term, Sa1, and the gradient term, Sa2 given in Eq. (2), do not show up explicitly in the thermal LB model. It was realized that, however, due to its special construction of the source term, the extension of Zheng et al.’s model [25] to include additional sources, such as viscous dissipation and pressure work, is not straightforward. To overcome this limitation, Zheng et al. [26] also proposed a coupled LB model for axisymmetric thermal flows with viscous dissipation and pressure work starting from the continuous axisymmetric Boltzmann equation. It is noticed that the velocity distribution functions are explicitly coupled in the evolution equation of the energy distribution functions in [26]. It is emphasized that all the existing LB models for axisymmetric thermal flows in [23–26] used the BGK collision operator. A robust and straightforward MRT-LB model that can simulate the general axisymmetric CDE is highly desired. Furthermore, we notice that numerical validation regarding the order-of-accuracy of the axisymmetric thermal LB models has not been well established. In this work, we propose an MRT-LB model for the general axisymmetric CDE for scalar transport such as temperature, mass concentration and the azimuthal velocity component in rotational flows. When the CDE is coupled with the hydrodynamic equations, the radial-axial velocity field is solved using Zhou’s axisymmetric LB model [20] with the BGK collision model replaced by an MRTbased model. The axisymmetric CDE in the cylindrical coordinate system is transformed to an ordinary CDE in the Cartesian coordinate system, and the extra terms due to the coordinate transformation are considered as additional sources. The scalar source term, Sa1 in Eq. (2), is obtained from the summation of the distribution functions, and the gradient term, Sa2 in Eq. (2), is calculated from the non-equilibrium parts of the distribution functions; thus no finite-difference calculations are required. The transformed CDE in the Cartesian coordinate system is then simulated with the D2Q5 MRT-LB model proposed in [5]. Compared with the existing LB models [23–26] for axisymmetric thermal flows, the proposed model for the axisymmetric CDE has the following features: (i) the macroscopic variable and its gradient terms in the CDE, such as Sa1 and Sa2 in Eq. (2), are naturally treated as additional source terms, and they are obtained from the distribution functions; thus no special treatment for the LB evolution equation or the equilibrium distribution functions is required; (ii) it is easier to implement when compared with Li et al.’s model [24] and the standard calculation of the macroscopic variables is preserved; (iii) the forms of the source terms can be very general and they can be directly included in the present model, (iv) viscous dissipation and pressure work can be conveniently treated as source terms in the present model; thus the explicit coupling of the velocity distribution function in the evolution equation for the energy distribution function in [26] is avoided, and (iv) rather than using BGK collision operators as in [23–26], an MRT collision operator is applied to the present model. To validate the applicability and accuracy of the proposed LB model for the axisymmetric CDE, several numerical tests are conducted in this study. Analytical solutions are available for the first three tests, and the results show that the proposed model possesses second-order accuracy in space and first-order accuracy in time. The axisymmetric LB model for the radial-axial velocity field by Zhou

340

L. Li et al. / International Journal of Heat and Mass Transfer 67 (2013) 338–351

[20] and the present LB model for the temperature CDE are coupled in the fourth and fifth tests, where natural convection in an annulus between two coaxial vertical cylinders and confined swirling flows in a vertical cylindrical container are studied, respectively. The numerical results from these two tests are compared with published data, and good agreement is obtained for each case. The rest of this paper is organized as follows. In Section 2, the D2Q5 MRT-LB model for the CDE is introduced; and the treatment of the macroscopic scalar variable and its gradient in the source terms in the transformed CDE is also provided. Numerical validation with five tests is presented in Section 3, and Section 4 concludes the paper.

3 x2 0

2

x1

1

4 Fig. 1. Discrete velocity set {ea} for the D2Q5 LB model. ea = (0, 0), (±1, 0), and (0, ±1).

  m Drr @/ 1 1 dx X ¼ 1 ear g ðneqÞ a r 2srr dt a¼1 r @r   X 1 1 dx m ¼ 1 ear ½g a  g ðeqÞ a ; r 2srr dt a¼1

Sa2  2. Lattice Boltzmann equation method The treatment of the additional terms, Sa1, and Sa2, in Eq. (2) due to the polar-to-Cartesian coordinate transformation as additional source terms is straightforward in the LBE method. While it is preferable to obtain these macroscopic quantities, especially the gradient term Sa2, from the microscopic distribution functions rather than using finite-difference schemes on the solved macroscopic field since the standard collision-streaming procedure is carried out on the distribution function level. Yoshida and Nagaoka [5] gave a formula in their MRT-LB model to obtain the scalar gradient based on the distribution functions. The accuracy of the gradient formula was examined and verified by Li et al. [8] with numerous computational examples. Li et al. [9] also provided an evaluation formula to obtain the gradient from the moment of the local non-equilibrium components of the distribution functions and showed that these two formulas are essentially equivalent. The non-equilibrium approach [9] for the scalar gradient evaluation will be used in this work. The details about the present MRT-LB model for Eq. (2) based on the model proposed in [5] and its implementation are described in the following. When the CDE (2) is coupled with the momentum transport equation, the axisymmetric LB model by Zhou [20] is applied for the radial-axial velocity field and the BGK collision operator is extended to an MRT version for consistency. The extension of Zhou’s BGK model [20] to the MRT model is given in Appendix A. 2.1. Lattice Boltzmann models The LBE governs the evolution of the scalar distribution function, g(x, n, t), which is defined in the discrete phase space (x, n) and time t, where x = rer + zez is the spatial vector and n is the particle velocity usually discretized to a small set of discrete velocities {na|a = 1, 2, . . . , m}. With the objective to recover the CDE (2), we propose the following evolution equation for the distribution functions

g a ðx þ ea dt; t þ dtÞ  g a ðx; tÞ ¼ ½Lgðx; tÞa þ xa ½Sa1 ðx; tÞ þ Sa2 ðx; tÞ þ Gðx; tÞdt;

ð3Þ

where ga(x, t)  g(x, na, t), ea is the ath discrete velocity vector (the discrete velocity set {ea} for the two-dimensional five-velocity (D2Q5) LB model is shown in Fig. 1), dt is the time step, L is the standard collision operator, and xa is the weight coefficient. The leading order solution of the CDE is obtained from m X /¼ ga ;

ð4Þ

a¼0

and the additional source terms are obtained from the distribution functions as

Sa1

m m ur ur X ur X  /¼ ga ¼  g ðeqÞ ; r r a¼0 r a¼0 a

ð5Þ

ðeqÞ

ð6Þ

ðneqÞ

where g a (defined below in Eq. (8)) and g a are the equilibrium and non-equilibrium components of the distribution function ga, ðneqÞ ðeqÞ respectively, with g a ¼ g a  g a , srr is the relaxation coefficient (see Eq. (12) below) and e ar is the projection of the discrete velocity ea in the radial direction. The derivation and numerical validation of the second-order accuracy of the gradient evaluation in Eq. (6) in ðneqÞ terms of g a can be found in [9]. The development of the axisymmetric LB models in this study follows the procedure in [5] and thus the present LB models are able to preserve their second-order accuracy in space and first-order accuracy in time. The D2Q5 discrete lattice model is applied for the CDE (2) rather than the D2Q9 model that is widely used for simulating the velocity field for athermal problems since (i) the CDE is considered as a scalar transport equation so that only the conservation of the leading-order moment of the distribution functions is required as shown in Eq. (4). While both the leading-order and first-order moments of the distribution functions must be conserved [8,27] in athermal LB models for velocity field simulations, which necessitates more discrete lattice velocities than the D2Q5 model; (ii) the D2Q5 model is more efficient to implement and requires much less effort than the D2Q9 model in the boundary condition treatment, especially for curved boundaries; and (iii) its robustness and accuracy have been verified with detailed analysis and numerical tests in [5,8,9]. The details about of the present axisymmetric LB model formulation are given in the following and the validation of the spatial and temporal accuracy is provided in Section 3. 2.1.1. BGK model The BGK collision operator with a single-relaxation-time coefficient s is

1 ½Lgðx; tÞa ¼  ½g a ðx; tÞ  g ðeqÞ a ðx; tÞ:

ð7Þ

s

ðeqÞ

Following [5], the equilibrium distribution function g a ðx; tÞ is defined as

g ðeqÞ a ðx; tÞ ¼



xa þ

 dtuj eaj xa /ðx; tÞ; dxeD

ð8Þ

where dx is the square lattice grid spacing and eD is a constant related to the dimensions. In the D2Q5 model, the weight coefficients are chosen to be [5]

xa ¼



1=3; ða ¼ 0Þ; 1=6; ða ¼ 1; 2; 3; 4Þ;

and the constant eD in Eq. (8) is eD = 1/3.

ð9Þ

341

L. Li et al. / International Journal of Heat and Mass Transfer 67 (2013) 338–351

2.1.2. Multiple-relaxation-time (MRT) model The MRT-LB models, which allow the various distribution functions and their moments to evolve toward their equilibria at different relaxation rates, usually have better numerical stability and accuracy than the BGK-LB model as demonstrated in [5,8,9]. In addition, anisotropic diffusion can be effectively simulated in MRT-LB models [5]. The MRT collision operator is thus applied in the present LB model and it can be written as

Lgðx; tÞ ¼ M1 S½mðx; tÞ  mðeqÞ ðx; tÞ;

ð10Þ

where M is a matrix to transform the distribution functions g to their moments m by

2

1

1

1

1

1

3

7 6 0 7 6 0 1 1 0 7 6 0 1 1 7 m ¼ Mg ¼ 6 7g; 60 0 7 6 4 1 1 1 1 5 4 0 1 1 1 1

2.2.1. Dirichlet boundary condition Given / = Ud at the boundary node xw, the boundary condition for the distribution function can be expressed as [8]

ð11Þ

1

S

s0

6 60 6 ¼6 60 6 40 0

0

3

0

0

0

srr srz srz szz

0

0 0

s3

7 07 7 07 7: 7 05

0 0

0 0

ð12Þ

s4

The equilibrium moments of the distribution functions in Eq. (10) are explicitly defined as [8]

mðeqÞ ¼ ð0; ur /; uz /; a/; 0ÞT ;

ð13Þ

where a is a constant related to the weight coefficient as a = (5x0  1) = 2/3. The relaxation coefficients, sij, in the relaxation matrix are related to the diffusion coefficients, Dij, as

1 2

sij ¼ dij þ

dt

eD ðdxÞ2

Dij ;

g a ðxf ; t þ dtÞ ¼ cd1 g^a ðxf ; tÞ þ cd2 g^a ðxff ; tÞ þ cd3 g^a ðxf ; tÞ þ cd4 eD Ud ;

and S is a matrix of multiple-relaxation-time parameters [5]

2

that for the distribution functions. As shown in Fig. 2, a boundary condition is required during the streaming step in Eq. (16) as the distribution functions at the exterior nodes (xe in Fig. 2) are not defined. Considering the LBE method is second-order accurate in the interior nodes, second-order accurate boundary conditions are usually pursued. The following second-order accurate Dirichlet and Neumann boundary conditions by Li et al. [8] based on the ‘‘bounce-back’’ idea and spatial interpolation are used in the present work. The fraction of the intersected link inside the computational field is denoted by D = ||xf  xw||/||xf  xe|| with 0 6 D 6 1, where xf is the first interior lattice node inside the field (see Fig. 2).

where xff is the second lattice node inside the computational field along the lattice velocity direction ea (see Fig. 2), i.e., xff ¼ xf þ ea dt, and cd1–cd4 are the coefficients related to the local D value along ea direction. The asymptotic analysis in [8] showed that second-order accuracy is preserved for the boundary condition treatment given by Eq. (17) when the following relationship is maintained, where cd1 (–1) is an adjustable variable

cd2 ¼ 

2Dcd1 þ 1 ; 2D þ 1

cd1 þ 2D ; 2D þ 1

and cd4 ¼

cd1 þ 1 : 2D þ 1

ð18Þ

Three particular schemes were examined in [8] and they are provided below for completeness. Scheme 1: ( ð2DÞg^a ðxf ;tÞ þ ð2D  1Þg^a ðxff ;tÞ þ eD Ud ; ð0 6 D 6 0:5Þ     g a ðxf ; t þ dtÞ ¼  1   2D g^a ðxf ;tÞ þ 1  21D g^a ðxf ; tÞ þ 21D eD Ud ; ðD > 0:5Þ ð19aÞ

Scheme 2:

! ð2D  1Þ2 g^a ðxff ; tÞ 2D þ 1     2D  1 3  2D g^a ðxf ; tÞ þ þ2 eD Ud ; and 2D þ 1 2D þ 1

g a ðxf ; t þ dtÞ ¼ 2ðD  1Þg^a ðxf ; tÞ 

g^a ðx; tÞ ¼ g a ðx; tÞ  ½M1 Sðmðx; tÞ  mðeqÞ ðx; tÞÞa ð15Þ

streaming step:

g a ðx þ ea dt; t þ dtÞ ¼ g^a ðx; tÞ;

cd3 ¼

ð14Þ

in order to obtain the leading-order solution of the CDE [5]. For efficient computations, the evolution equation (3) with the above MRT collision operator is solved in two steps collision step:

þ xa ½Sa1 ðx; tÞ þ Sa2 ðx; tÞ þ Gðx; tÞdt; and

ð17Þ

ð16Þ

where g^a represents the post-collision state. 2.2. Boundary conditions Boundary condition treatment is essential in the LBE method to convert given boundary conditions for the macroscopic variable to

Scheme 3:

  2D  1 g^a ðxff ; tÞ g a ðxf ; t þ dtÞ ¼ g^a ðxf ; tÞ þ 2D þ 1     2D  1 2 g^a ðxf ; tÞ þ þ eD Ud : 2D þ 1 2D þ 1

ð19cÞ

2.2.2. Neumann boundary condition Given the flux Dn @/ ¼ Un in the normal direction on the @n boundary (here Dn represents the diffusion coefficient in the specific normal direction), the second-order accurate boundary condition treatment is unique and it is given in [8] as

  2D  1 g^a ðxff ; tÞ g a ðxf ; t þ dtÞ ¼ g^a ðxf ; tÞ  2D þ 1     2D  1 2 dt g^a ðxf ; tÞ þ Una ; þ 2D þ 1 2D þ 1 dx

Fig. 2. Intersection on the boundary wall by the lattice. Closed circles: field nodes inside the computational domain; closed square: boundary node; open circle: exterior node.

ð19bÞ

ð20Þ

where Una is the flux in the discrete velocity direction ea . When the local boundary normal n is aligned with ea , Una ¼ Un and treatment (20) can be directly applied. When the local boundary normal is not in the ea direction, which is usually encountered on inclined or curved boundaries, Una is not equal to Un and it also depends on

342

L. Li et al. / International Journal of Heat and Mass Transfer 67 (2013) 338–351

the unknown tangential flux. The details about the Neumann boundary condition treatment and its extension to mixed boundary conditions for curved boundaries can be found in [8].

3.2. Convection–diffusion in a circular pipe

3. Numerical validation To validate the accuracy of the proposed LB model for the axisymmetric CDE, we present five numerical tests in this section. The velocity field is specified in each of the first three tests and the analytical solutions for / are available. The accuracy and convergence order of the numerical solutions are thus investigated in detail. The fourth and fifth tests consider axisymmetric thermo-hydrodynamic problems, where the radial-axial velocity fields are simulated with an MRT version of the athermal LB model by Zhou [20], and the presently proposed LB model is applied for the thermal field as well as the azimuthal velocity field in the fifth test. These two LB models are coupled to reach steady state and the numerical results are compared with published results.

3.1. Steady convection–diffusion in an annulus In this example, we consider the axisymmetric convection–diffusion in an infinitely long annulus defined in Ri 6 r 6 Ro, where the variations in both the azimuthal and axial directions are neglected. With a constant volumetric flow rate Q supplied on the inner surface, the radial velocity in the annulus is ur = Q/r. The convection–diffusion inside the annulus is characterized by the Péclet number, Pe ¼ uRi Ri =Dr ¼ Q =Dr . For given Dirichlet conditions /(r = Ri,o) = /i,o, the exact solution for / inside the annulus is

/ex ðRi 6 r 6 Ro ; zÞ ¼

Pe Pe Pe /i ðRPe o  r Þ þ /o ðr  Ri Þ Pe RPe o  Ri

:

ð21Þ

With no axial variation only five lattice nodes are used in the zdirection, and periodic boundary conditions are applied at the ends. Each of the inner and outer boundaries is placed at the center of the lattice links (D = 0.5) so that the standard Dirichlet ‘‘bounce-back’’ scheme, which is recovered by all three schemes in Eqs. (19a)– (19c), is implemented. To assess the accuracy of the proposed LB model, the following relative L-2 norm error is defined in the annulus

E2 ¼

( X r;z

ð/ÞLBE  ð/Þex

2

To recapitulate, the present case verifies that the proposed axisymmtric LB model is second-order accurate in space and the MRT model is more robust and accurate than the BGK model.

, )1=2 X 2 ð/Þex :

ð22Þ

r;z

Fig. 3(a) and (b) shows the error, E2, versus the grid resolution, 1/Ri, at various sD values with the present MRT and the BGK models applied, respectively, for Pe = 20, Ro/Ri = 2, /i = 0.5, and /o = 1.0. In the MRT model, the relaxation coefficients are set to be srr = szz = sD and sp = 1 (p = 0, 3, 4) throughout Section 3. When srr = szz = sp = sD, the MRT model reduces to the BGK model. Second-order accuracy with respect to the grid resolution is observed for both models in Fig. 3. It is also noticed that the error E2 increases as sD increases for both models, while the error increase in the MRT model is less sensitive to sD compared to that in the BGK model. As shown in Fig. 3, for sD 6 1.0, the numerical errors in both models are close to each other while for sD P 1.5 the MRT model leads to much smaller errors. Those observations are consistent with the numerical results reported in [5]. Fig. 4 shows the error, E2, versus 1/Ri, for different Péclet numbers at sD = 0.75. Second-order accuracy is obtained for each Pe number tested and E2 increases as Pe increases. As evidenced in Fig. 3, it is preferable to use smaller relaxation coefficient sD and higher grid resolution in order to simulate large Pe number cases.

In order to obtain analytical solutions, the fluid flow inside a circular pipe of radius R0 is assumed to be a plug flow with a constant axial velocity U and zero radial velocity [8,9]. The transient axisymmetric CDE is thus expressed as

    @/ @/ @ @/ @ @/ Drr @/ þ þ Drr Dzz þU ¼ : @t @z @r @r @z @z r @r

ð23Þ

The last term on the right side of Eq. (23) is treated as a source term in the present model as in Eq. (6). Only isotropic diffusion is considered for simplicity so that Drr = Dzz = D and the Péclet number is Pe = (2R0)U/D. The computational domain and the square lattice in the r–z plane are schematically shown in Fig. 5. A periodic boundary condition in the z-direction is imposed so that f(x + L) = f(x) is valid for both / and the distribution function ga. The centerline, r = 0, is always located half-way between the lattice nodes with D = 0.5; and the Neumann condition o//or|r=0 = 0 is imposed due to axisymmetry. Both the steady-state and transient solutions of Eq. (23) will be investigated to assess the spatial and temporal accuracy of the present axisymmetric LB model in the presence of strong convection.

3.2.1. Steady-state solutions At steady-state, the transient term, @//@t, in the governing equation (23) disappears. This steady problem has been studied in [8,9] where the D3Q7 LB model was used in the 3-D domain and curved boundary condition treatments were required to preserve the exact geometry of the circular pipe. In the present work, in contrast, as the 2-D axisymmetric LB model is used, the effect of the curved boundary can be eliminated. Dirichlet and Neumann boundary conditions on the surface of r = R0 are studied separately in this section. For the Dirichlet condition on the boundary with specific D values, Scheme 2 in Eq. (19b) is implemented as it showed smaller relative errors than the other two schemes in Eqs. (19a) and (19c) for most cases studied in [8]. The Péclet number Pe = 20 and the relaxation coefficient sD = 0.75 are used in all simulations.

3.2.1.1. Dirichlet boundary condition on the wall. When a sinusoidal Dirichlet condition is imposed on the wall

/ðr ¼ R0 ; zÞ ¼ cosðkzÞ;

with k ¼ 2p=L;

ð24Þ

the analytical solution for / is [8]



I0 ðkrÞ ; /ex ðr; zÞ ¼ Real eikz I0 ðkR0 Þ

ð25Þ

where ‘‘Real’’ means taking the real part of a complex number, I0(kr) is the modified Bessel’s function of the first kind of order 0 and qffiffiffiffiffiffiffiffiffiffiffiffiffi iU . k ¼ k 1 þ Dk For illustration purposes, Fig. 6 shows the contours of / in the r– z plane from the numerical simulation using Nr = 42, Nz = 83 (Nr and Nz denote the numbers of lattice nodes in the r- and z-directions, respectively) and D = 0.5 for the pipe wall. To assess the accuracy of the present LB model, the relative L-2 norm error, E2, defined in Eq. (22) is computed. In addition, with the analytical solution given by Eq. (25), the exact solutions for the boundary wall flux, the interior gradients @//@r and @//@z in both directions can be obtained. The numerical boundary flux is obtained

343

L. Li et al. / International Journal of Heat and Mass Transfer 67 (2013) 338–351

10

10

−1

10

τ D = 0.75 τ D = 1.0 τ D = 1.5 τ D = 2.5 τ D = 4.5

10

−2

2

E2

10

0

E

10

−3

10

10

0

τD = 0.75 τ D = 1.0 τ D = 1.5 τ D = 2.5 τ D = 4.5

−1

−2

−3

slope = 2

slope = 2 10

−4 −2

10

10

−1

10

−4

10

−2

10

−1

1/Ri

1/R

i

(a)

(b)

Fig. 3. Relative L-2 norm error E2 versus the grid resolution, 1/Ri, at various sD values, for the 1-D annulus convection–diffusion problem with Pe = 20. (a) MRT model, and (b) BGK model.

10

−1

slope = 2 −2

E

2

10

10

10

−3

Pe = 10.0 Pe = 20.0 Pe = 50.0 Pe = 100.0

−4

−2

10

10

−1

1/R

i

Fig. 4. Relative L-2 norm error E2 versus the grid resolution, 1/Ri, at different Pe numbers, for the 1-D annulus convection–diffusion problem.

Fig. 6. Contours of / in the r–z plane at Nr = 42, Nz = 83, and D = 0.5 for the steady convection–diffusion in a pipe with a Dirichlet boundary condition on the wall.

( E2

gr

X

¼



interior nodes

   2 , X  2 )1=2 @/ @/ @/  ; @r LBE @r ex @r ex interior nodes ð27Þ

( E2

gz

¼

X

interior nodes



   2 , X  2 )1=2 @/ @/ @/  : @z LBE @z ex @z ex interior nodes ð28Þ

Fig. 5. Schematic of the square lattice in the r–z plane of a 3-D circular pipe. The axis of the pipe, r = 0, is fixed at the center of the link with D = 0.5, while the pipe surface, r = R0, is placed at an arbitrary location in the link with 0 6 D 6 1.



using the technique proposed in [8] with Una ¼ 1c2 d1 dx dt h



2 2 ^a ðxf ; tÞ ð2D þ 1Þg^a ðxf ; tÞ þ 2D  1c g^a ðxff ; tÞ þ 1c  1 g d1

d1

þeD Ud ; and the interior gradient is evaluated using the similar idea shown in Eq. (6). The following errors are further defined to assess their numerical accuracy

( E2

qw

¼

  2 ,X  2 )1=2 @/ @/  ; @r LBE @r ex @r ex r¼R

X @/

r¼R0

0

ð26Þ

The results of E2, E2_qw, E2_gr and E2_gz, defined in Eqs. (22) and (26)– (28)versus the grid resolution, 1/R0, at different D values are shown in Fig. 7(a)–(d), respectively. For each case in Fig. 7, second-order convergence is observed. It should be noted that only first-order convergence was obtained for the boundary flux and boundary heat transfer in [8] and [9], respectively, due to the curved boundary effect when the D3Q7 LB model was used. In the present axisymmetric case, the curved boundary effect is eliminated and thus secondorder accuracy is preserved for the interior distribution of /, the boundary flux, and the interior gradients in both directions. 3.2.1.2. Neumann boundary condition on the wall. When the Dirichlet condition in Eq. (24) is replaced with the Neumann condition

344

L. Li et al. / International Journal of Heat and Mass Transfer 67 (2013) 338–351

−1

10 10

−2

slope = 2

slope = 2

2

−2

qw

10

|

E

2

E

10

−3

10

Δ = 0.01 Δ = 0.25 Δ = 0.50 Δ = 0.75 Δ = 0.99

−4

Δ = 0.01 Δ = 0.25 Δ = 0.50 Δ = 0.75 Δ = 0.99

−3

10

−4

−2

10

−1

10

10

−2

10

1/R

1/R

0

10

−1

0

(b)

(a) −1

10 −2

10

slope = 2

slope = 2 −2

|

−3

|

E

2

10

2

E

gz

gr

10

Δ = 0.01 Δ = 0.25 Δ = 0.50 Δ = 0.75 Δ = 0.99

−4

10

Δ = 0.01 Δ = 0.25 Δ = 0.50 Δ = 0.75 Δ = 0.99

−3

10

−4

10

−2

10

10

−1

10

−2

−1

10

1/R

0

1/R0

(c)

(d)

Fig. 7. Relative L-2 norm errors versus the grid resolution, 1/R0, for the steady convection–diffusion in a pipe with Pe = 20 and a Dirichlet boundary condition on the wall: (a) E2 for the interior field of /, (b) E2_qw for the boundary flux, (c) E2_gr for the interior radial gradient o//or, and (d) E2_gz for the interior axial gradient o//oz.

 @/ Un ¼ D  ¼ D cosðkzÞ=R0 ; @r r¼R0

ð29Þ

the analytical solution for / becomes [8]

/exðr;zÞ ¼ Real eikz

I0 ðkrÞ ; kR0 I1 ðkR0 Þ

ð30Þ

where I1 (kr) is the modified Bessel’s function of the first kind of order 1. The contours of / in this case are very similar to that with a Dirichlet boundary condition except for the magnitude; thus the contours are not shown for brevity. Following [8], the boundary values of / are also computed based on the local distribution functions at the adjacent lattice nodes and the given boundary flux information, i.e., Ud ¼ e1D h





^a ðxff ;tÞ þ 1  1c2 g^a ðxf ; tÞ þ 1c2 ð2D þ 1Þg^a ðxf ;tÞ þ 1c2  2D g d1

dt dx

d1

d1

/ðr ¼ R0 ; z; tÞ ¼ cosðkz þ xtÞ;

where x = 2p/C is the frequency and C is the period of the imposed boundary temperature. In addition to the Péclet number, another qffiffiffiffiffiffi ffi qffiffiffiffi2ffi R2 x R characteristic number, the Stokes number, St  20pD ¼ CD0 , is defined for this unsteady problem [8,9]. The analytical solution for / in Eq. (23) subject to the boundary condition Eq. (32) is



I0 ðrrÞ ; /ex ðr; z; tÞ ¼ Real eiðkzþxtÞ I0 ðrR0 Þ sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x þ Uk ¼k 1þi : 2 Dk

( E2

tw

¼

X

r¼R0

,

"

½/LBE  /ex 2

X

E2 ¼

)1=2 /2ex

:

ð31Þ

r¼R0

Fig. 8(a)–(d) shows the errors E2, E2_tw, E2_gr and E2_gz defined in Eqs. (22), (31), (27), and (28), respectively, versus the grid resolution, 1/ R0, at different D values. Second-order accuracy is obtained for all cases in Fig. 8. The proposed LB model for the axisymmetric CDE with Neumann boundary conditions is also verified to be second-order accurate. 3.2.2. Unsteady problem For the unsteady case, the Dirichlet boundary condition in Eq. (24) is replaced with the following

with

r ð33Þ

To assess the spatial and temporal accuracy of the LB results, the time-averaged error for / in the interior domain on the r–z plane is defined as

Una : The following L-2 norm error is further defined to assess the

LB results

ð32Þ

1

C

Z 0

C

1 X ð/  /ex Þ2 dt Nr Nz r;z LBE

#1=2 :

ð34Þ

In the present LB simulations, Pe = 20 and St = 1 are used. For a specific relaxation coefficient sD, the diffusion coefficient D is first obtained from Eq. (14); then the axial flow velocity U and the time period C are determined from the definitions of Pe and St, respectively, for different spatial resolution, 1/R0. It should be noted that the spatial and temporal intervals, ex  dx/R0 and et  dt/C, in LB models for the convection diffusion equation, are directly related to each other. Yoshida and Nagaoka [5] showed in their diffusive scaling analysis that the relationship et/(ex)2 = const. must be maintained to obtain second-order accuracy in space and first-order accuracy in time. In the present test for a given Stokes number St, et =ðex Þ2 ¼ R20 =C ¼ St 2 D ¼ const. is satisfied as dx = dt = 1 in the LB model.

345

L. Li et al. / International Journal of Heat and Mass Transfer 67 (2013) 338–351

−1

−1

10

10

slope = 2

slope = 2 −2

−2

2

tw

10

10

E

2

E

|

Δ = 0.01 Δ = 0.25 Δ = 0.50 Δ = 0.75 Δ = 0.99

−3

10

Δ = 0.01 Δ = 0.25 Δ = 0.50 Δ = 0.75 Δ = 0.99

−3

10

−4

−4

10

−2

10

−1

10

10

10

1/R

−2

10

−1

1/R

0

0

(b)

(a) −1

10

10

−1

slope = 2 10

−2

slope = 2 −2

gz

gr

10 |

10

−3

2

10

E

E2

|

Δ = 0.01 Δ = 0.25 Δ = 0.50 Δ = 0.75 Δ = 0.99

−4

10

−2

−1

10

1/R

Δ = 0.01 Δ = 0.25 Δ = 0.50 Δ = 0.75 Δ = 0.99

−3

10

10

−2

−1

10

1/R

0

0

(c)

(d)

Fig. 8. Relative L-2 norm errors versus the grid resolution, 1/R0, for the steady convection–diffusion in a pipe with Pe = 20 and a Neumann boundary condition on the wall: (a) E2 for the interior field of /, (b) E2_tw for the boundary values of /, (c) E2_gr for the interior radial gradient o//or, and (d) E2_gz for the interior axial gradient o//oz.

−2

2

10

E

z

τ D = 0.51 τ D = 0.55 τ D = 0.65 τ D = 0.75 τ D = 1.00 τ D = 1.50 τ D = 2.00

−3

θ R 0

10

Δy P

Δx

r

−4

10

slope = 2 −2

−1

10

10

1/R

0

Fig. 9. Time-averaged error E2 for the interior field of / versus the grid resolution, 1/R0, for the unsteady convection–diffusion in a pipe with Pe = 20 and St = 1.0.

The initial condition /0 = /(x, t = 0) is obtained from Eq. (33) and the initial treatment

g a ðx; t ¼ 0Þ ¼



 dtu @/ xa þ j eaj xa /0  dx 0 ðM1 S1 Mej xÞa ; dxeD @xj

as proposed in [5] is applied here. The results of E2 defined in Eq. (34) versus the grid resolution, dx/R0 = 1/R0, for different relaxation coefficients of sD are shown in Fig. 9. Both r = 0 and r = R0 are located halfway between the lattice nodes (D = 0.5). Second-order convergence is observed for each case in Fig. 9 as R0/dx increases. It should be noted that the error in Fig. 9 contains both contributions from the temporal error and spatial error. When the spatial ffiffiffiffi2ffi R0 is doubled, in order to keep the same resolution is doubled,qi.e., R Stokes number, St ¼ CD0 , the period of the wall temperature variation, C, must be quadrupled in the LB simulation which equivalently leads to a reduction in the dimensionless time step, dt/C, by a factor of four. Thus although the temporal accuracy is only first

Fig. 10. Schematic representation of the computational domain and the lattice layout for the heat conduction inside a sphere.

order, the doubling in R0 leads to a fourfold decrease in both the spatial error (due to the second-order spatial accuracy) and temporal error (due to quadrupling in the ratio C/dt). As a result, the total error is decreased by a factor of four, as clearly shown in Fig. 9 by the slope of 2 in the log–log scale. This test demonstrates that the proposed LB model can be effectively applied to simulate timedependent axisymmetric convection–diffusion problems, and it verifies that the spatial accuracy is second order and the temporal accuracy is first order. 3.3. Steady heat conduction inside a sphere When there is no temperature variation in the azimuthal direction, the heat conduction inside a sphere is axisymmetric and it can

L. Li et al. / International Journal of Heat and Mass Transfer 67 (2013) 338–351

τ D = 0.75

slope = 1 tw

be represented in the radial-axial cylindrical coordinate system. Thus the axisymmetric LB model can be applied to simulate the heat conduction in a sphere. The computational domain and the lattice layout are schematically shown in Fig. 10. The curved geometry of the sphere is preserved by calculating the intersected link fractions (Dx and Dy) for those field lattice nodes, such as point P in Fig. 10, adjacent to the boundary. After the local D value along the lattice velocity direction ea is obtained for each boundary node, the Dirichlet and Neumann boundary condition treatments in Eqs. (19) and (20) can be implemented, respectively.

|

−2

10

E2, E2

346

coupled Scheme 1, E2 coupled Scheme 2, E2 coupled Scheme 3, E2 coupled Scheme 1, E2 _ tw coupled Scheme 2, E2 _ tw coupled Scheme 3, E2 _ tw −3

10

3.3.1. Dirichlet boundary condition on the wall A Dirichlet boundary condition of a Legendre polynomial Pn(l = cos h) on the sphere surface, r = R0, is considered and the analytical solution for the interior temperature / is [31]

/ex ðr; lÞ ¼ ðr=R0 Þn P n ðlÞ:

−2

−1

10

10

1/R0

Fig. 12. Relative L-2 norm errors E2 and E2_tw for the interior field and boundary values of /, respectively, versus the grid resolution, 1/R0, for the steady heat conduction inside a sphere with a Neumann boundary condition.

ð35Þ 4

2

A representative case of P4(l) = (35l  30l + 3)/8 is studied in this work. The axis of the sphere is placed at the center of the lattice link (see Fig. 10) so that the symmetry condition, o//or|r=0 = 0, which is a Neumann boundary condition, is used in all simulations. Fig. 11 shows the relative L-2 norm errors of the interior temperature field and the boundary heat flux E2 and E2_qw defined in Eqs. (22) and (26), respectively, versus the grid resolution 1/R0, with the Dirichlet Schemes 1–3 in Eqs. (19a)–(19c) implemented. It should be noted that the boundary heat flux in this test is evaluated along the discrete lattice velocity directions using the formula proposed in [8]. For those boundary nodes intersected by the lattice vectors along z-direction, @/=@rjr¼R0 in Eq. (26) is replaced by @/=@zjr¼R0 . For each case, quadratic and linear convergence orders are noticed in Fig. 11 for E2 and E2_qw, respectively. These observations are consistent with the results in [8] for curved boundaries. Second-order accuracy of the temperature field is thus preserved with the present axisymmetric LB model, and the firstorder accurate boundary heat flux is due to the irregularly distributed link fractions in both the radial and axial directions caused by the curved boundary [8].

directions. Those coupled schemes for Neumann problems with curved boundaries [8] are also used in the present case. The L-2 norm errors E2 and E2_tw defined in Eqs. (22) and (31), respectively, are computed for this test and their convergence results are shown in Fig. 12, where on average a first-order convergence is obtained for both E2 and E2_tw. As demonstrated earlier in Fig. 11 for Dirichlet problems with curved boundaries, the evaluation of the boundary flux Una is only first-order accurate even with a second-order accurate temperature field obtained. The same formula for evaluating the boundary flux Una is directly used in the coupled schemes for Neumann problems, thus some first-order error terms are introduced in the boundary condition treatment. The numerical tests in this section serve as an excellent example to show that even with a second-order LB model implemented, the overall accuracy of the LB solutions is directly affected by the geometry of the boundary and the relevant boundary condition treatment.

3.3.2. Neumann boundary condition on the wall When a Neumann flux condition Un ¼ D@/=@rjr¼R0 ¼ DnP n ðlÞ=R0 is imposed on the sphere surface the analytical solution for / is the same as that in Eq. (35) and n = 4 is again studied. To implement the Neumann condition treatment in Eq. (20), the boundary flux Una in each lattice velocity direction is required. Li et al. [8] have developed a Cartesian decomposition method to obtain Una based on the given normal flux Un on curved boundaries by coupling the boundary flux components in both lattice velocity

In the previous tests, the velocity fields are specified. With the objective to show that the proposed LB model for the axisymmetric CDE of / can be effectively coupled with the athermal LB models for axisymmetric hydrodynamic equations, we consider in this test the natural convection in an annulus formed by two coaxial vertical cylinders. The velocity field is solved using an MRT version of Zhou’s athermal axisymmetric LB model [20] and the temperature field is simulated with the present LB model at each time step. The two models are iterated to reach steady state.

τ D = 0.75

3.4. Steady natural convection in an annulus between two coaxial vertical cylinders

z

slope = 1

10−2

∂zT = 0

2 qw

Ri

2

E, E

|

10

slope = 2

−3

g

Scheme 1, E2 Scheme 3, E2 Scheme 1, E

10

To

Ti

Scheme 2, E2

H

Ro

2 _ qw

Scheme 2, E2 _ qw

−4

Scheme 3, E2 _ qw

10

−2

1/R

10−1

∂zT = 0

0

Fig. 11. Relative L-2 norm errors E2 and E2_qw for the interior field of / and the boundary flux, respectively, versus the grid resolution, 1/R0, for the steady heat conduction inside a sphere with a Dirichlet boundary condition.

Fig. 13. Schematic depiction of the computational domain (Ro/Ri = 2.0, and H/ (Ro  Ri) = 2.0) and the thermal boundary conditions on the walls all placed at the center of the lattice links (D = 0.5).

L. Li et al. / International Journal of Heat and Mass Transfer 67 (2013) 338–351

347

Fig. 14. Streamlines in the annulus between two coaxial vertical cylinders at (a) Ra = 103, (b) Ra = 104, and (c) Ra = 105.

The computational domain, the thermal boundary conditions, and the schematic layout of the square lattice are sketched in Fig. 13. Both the radius ratio Ro/Ri and the aspect ratio H/(Ro – Ri) are equal to 2.0 in the present simulations. The vertical cylinder surfaces are maintained at constant temperatures Ti and To, respectively with Ti > To, while the top and bottom walls are insulated. Invoking the Boussinesq approximation, the fluid properties in the annulus are assumed to be constant, except for the density in the buoyancy term in the momentum equation. The Prandtl number (Pr  m/a) is fixed at Pr = 0.7, and the characteristic Rayleigh number is defined as Ra = gb(Ti – To)(Ro – Ri)3/am, where g is the gravitational acceleration, b is the thermal-expansion coefficient, a is the thermal diffusivity and m is the kinematic viscosity of the fluid.

The effect of the temperature distribution on the velocity field is taken into account by adding a body force term qaz = –qgb(T – Tm) in the momentum equation, where Tm = (Ti + To)/2 is the average temperature. It should also be noted that the macroscopic velocities ur and uz are used for calculating the equilibrium moments in Eq. (13) while the velocity distribution functions are not used in the evolution equation (3). Thus the coupling between the athermal and thermal LB models is only on the macroscopic level, and the direct coupling of the distribution functions for the velocity and temperature as in [26] is eliminated. This makes the implementation of the LBE models more straightforward and convenient. The convergence is checked by monitoring the variations of the velocity and temperature values at some representative nodes, and the relative differences between successive iterations

Fig. 15. Isotherms in the annulus between two coaxial vertical cylinders at (a) Ra = 103, (b) Ra = 104, and (c) Ra = 105.

348

L. Li et al. / International Journal of Heat and Mass Transfer 67 (2013) 338–351

heat transfer on the walls is obtained using the technique proposed in [9]. Both the heat flux and heat transfer evaluations are secondorder accurate as verified in [8,9]. A grid-refinement study is first carried out and the Nusselt numbers Nui and Nuo defined in Eq. (36) obtained with different mesh sizes are listed in Table 1. The numerical results for Nui and Nuo evaluated on the inner and outer surfaces are almost identical to each other for all cases, demonstrating the excellent conservation property of the LBE scheme. It is emphasized that such high consistency between Nui and Nuo can only be obtained after the numerical solutions are fully converged. In the present test, all the boundaries are placed at the center of the lattice links (D = 0.5), thus the second-order accuracy is preserved for the LBE method in the interior, for the ‘‘bounce-back’’ hydrodynamic and thermal boundary conditions and the heat transfer evaluation [9] on the straight boundaries. Hence more accurate results can be obtained using the Richardson extrapolation of the raw second-order results before they are rounded in Table 1. In the present case, Nuextrapolated = (cnNufine  Nucoarse)/(cn  1) where c = (Nr  2)fine/ (Nr  2)coarse is the grid resolution ratio, and n = 2 is the order of convergence of the LBE solutions. Table 2 compares the extrapolated Nusselt number with the published data in [24] and [34]. Both the results for Nui and Nuo are presented in [24], and the results from [34] are the average Nusselt numbers on the two cylinder surfaces, i.e., Nu = (Nui + Nuo)/2. The present results in Table 2 are in good agreement with those from [24] and [34], demonstrating that the present LB model can be effectively coupled with the athermal axisymmetric LB models and the interaction between the velocity field and the thermal field is accurately captured.

Table 1 Surface-averaged Nusselt numbers at different grid resolution. Mesh

Average Nusselt numbers Nui, Nuo Ra = 103

102  202 127  252 152  302 202  402

1.69200, 1.69206, 1.69210, 1.69213,

Ra = 104 1.69200 1.69206 1.69210 1.69213

3.21673, 3.21607, 3.21572, 3.21531,

Ra = 105 3.21673 3.21607 3.21572 3.21531

5.83096, 5.81486, 5.80636, 5.79802,

5.83095 5.81486 5.80636 5.79802

Table 2 Comparison of the present surface-averaged Nusselt numbers with published results. Ra

Average Nusselt number Nu Present (extrapolated)

Ref. [24] Nui

1.6922 3.2148 5.7873

 3.216 5.782

Ref. [24] Nuo

Ref. [34]

 3.218 5.787

 3.1629 5.8818

(101  201) 103 104 105

(41  81)

are at least less than 108 for both the velocity and temperature fields. The typical number of time steps to reach steady state is in the order of 105 and the physical computational time is within a few hours on a typical personal computer (processor speed 3.41 GHz). The exact number of time steps depends on the number of lattices used, the magnitude of Ra, and the stopping criterion for steady state convergence. The streamlines and isotherms at steady state are shown in Figs. 14 and 15, respectively, for Ra = 103, 104 and 105 with Nr  Nz = 202  402. The present results agree well with those reported in [24,32–34]. To quantify the accuracy of the present numerical results, the following surface-averaged Nusselt numbers on the cylinder walls are defined and compared with the published results in [24,34]

 @T   dz; and 0 @r i Z H  Ro @T  Nuo ¼  dz: HðT i  T o Þ 0 @r o

Nui ¼ 

Ri HðT i  T o Þ

Z

3.5. Swirling thermal flows in a vertical cylinder Confined axisymmetric flows driven by a rotational velocity is considered in this test to further demonstrate the utility of the present LB model when it is coupled with the axisymmetric LB model for the velocity field. The fluid inside the cylindrical container is assumed to be incompressible and the Boussinesq assumption is valid for the present flow. The top disk rotates at a constant angular speed X, and the rest of the solid boundaries are stationary. The top and bottom disks are maintained at constant temperatures, Th and Tc, respectively (Th > Tc), while the side walls are insulated. The swirling flows are characterized by the Reynolds number Re, and the Richardson number Ri, defined as

H

ð36Þ

It is worth mentioning that in this work the boundary temperature gradients in Eq. (36) are evaluated on the basis of the temperature distribution functions as proposed in [8] rather than using any finite-difference schemes on the temperature field. And the overall

Re ¼ R2 X=v ;

1

0.8

0.8

0.6

0.6

ð37Þ

z/H

z z/ /HH

1

3

and Ri ¼ gbDTh =ðRX2 Þ;

0.4

0.4

0.2

0.2

0

0

0.2

0.4

0.6

r/R

(a)

0.8

1

0

0.2

0.4

R rr // R

0.6

0.8

1

(b)

Fig. 16. Contours of the azimuthal velocity uu (a) and isotherms (b) of the confined swirling flow at Re = 2000 and Ri = 0.

349

L. Li et al. / International Journal of Heat and Mass Transfer 67 (2013) 338–351

1

0.8

0.8

0.6

0.6

z/H

z/H

1

0.4

0.4

0.2

0.2

0

0

0.2

0.4

0.6

0.8

0

1

0.2

0.4

r/R

0.6

0.8

1

r/R

(a)

(b)

Fig. 17. Contours of the azimuthal velocity uu (a) and isotherms (b) of the confined swirling flow at Re = 2000 and Ri = 1.0.

1 x x

0.8

1

x x x xx x x xx x x x x x xx

x

0.8

xx x x x x x x x

0.6

x

x

z/H 0.4

z/H

x x x x x x x x

0.6

x

0.2

0 -0.2

x x

-0.1

x

x

x x x xx xx x x

0

Ri = 0.0, LBE Ri = 0.0, [35] Ri = 1.0, LBE Ri = 1.0, [35]

0.1

ur

x

0

Re = 2000

x x x x x x

0.2

0.2

xx

x

x

x

xx

x

0

0.2

0.4

uϕ

z/H

x x x xx xx x

0.2

-0.05

0

xx xx x x x x x x x x

Re = 2000

0.8

Ri = 0.0, LBE Ri = 0.0, [35] Ri = 1.0, LBE Ri = 1.0, [35]

x

0.6

1

x

Ri = 0.0, LBE Ri = 0.0, [35] Ri = 1.0, LBE Ri = 1.0, [35]

0.05

x x

x x x x x

0.2

0.1

(c)

x x

x

0.4

Re = 2000

xx x

uz

0.8

z/H

x x x x x x x x x x x x x

0.4

0 -0.1

1 x x x

x x

0.6

0.6

(b)

xx

0.8

Ri = 0.0, LBE Ri = 0.0, [35] Ri = 1.0, LBE Ri = 1.0, [35]

x

(a) 1

x

x

x x

0.4

Re = 2000

x

x

x x x x x x x

0

x

-0.4

x

x

x

-0.2

x x x

x xx x x x x

0

T

*

xx x

0.2

x x x xx

0.4

(d)

Fig. 18. Velocity components and temperature profiles along the vertical line at r/R = 0.8 for Re = 2000: (a) radial velocity ur, (b) azimuthal velocity uu, (c) axial velocity uz, and (d) temperature.

where R is the radius of the cylinder, DT = Th – Tc, and h is the aspect ratio of the height to the radius of the cylinder, h = H/R. The flow pattern and heat transfer of the swirling flows in the range of 102 6 Re 6 3  103 and 0 6 Ri 6 1.0 at fixed Prandtl number Pr = 1.0 and cylinder aspect ratio h = 1.0 were investigated in detail by Iwatsu [35] using the vorticity-stream function formulation. Zheng et al. [25] also simulated the same swirling flow prob-

lem and compared their results at Re = 2000, and Ri = 0 and 1.0 with that in [35] to validate their axisymmetric thermal LB model. In the present simulations, the radial and axial velocity components ur and uz are solved using Zhou’s axisymmetric LB model [20] with an MRT collision operator (see Appendix A). While it is noted that the governing equation for the azimuthal velocity, uu, is also a CDE [20,30]

350

L. Li et al. / International Journal of Heat and Mass Transfer 67 (2013) 338–351

! @uu @ @ 2 uu @ 2 uu @ þ ður uu Þ þ ðuz uu Þ ¼ m þ @r @z @t @r2 @z2   2ur m m @uu þ 2 uu þ ;  r r r @r

ð38Þ

where m is the kinematic viscosity of the fluid. Eq. (38) is very similar to Eq. (2) when uu is considered as a scalar variable. The source   terms Sa1 and Sa2 in Eq. (2) can be modified as Sa1 ¼  2ur r þ rm2 uu @u u and Sa2 ¼ mr @r , to account for the scalar variable and its gradient terms in Eq. (38), respectively. Thus the azimuthal velocity component in rotational flows can also be solved using the present LB model for axisymmetric CDEs. In this study, Zhou’s LB model [20] for the radial and axial velocities and the present LB model for the azimuthal velocity and the temperature are coupled at each time step to reach the steady state. In the present case, the steady state is considered to be obtained when the relative differences between 50 successive iterations are less than 108 for both the velocity and temperature fields. Again, the typical number of time steps to reach steady state is in the order of 105, and the physical computational time is within a few hours on a typical personal computer (processor speed 3.41 GHz). The following results are obtained with mesh sizes of 150  150 for (Re, Ri) = (2000, 0), and 200  200 for (Re, Ri) = (2000, 1.0) after a grid-refinement study. The contours of the azimuthal velocity component uu and the isotherms are shown in Figs. 16 and 17 for the extreme values of Ri = 0 and 1.0, respectively, at Re = 2000. The swirling flows driven by the rotating top disk are clearly seen, and thin boundary layers are observed near the top disk in the contours of uu. At Ri = 1.0, the effect of the buoyancy force becomes significant and the isotherms in Fig. 17(b) become more horizontal than that for Ri = 0, especially in the half region near the bottom. It is also noticed that the magnitude of uu is very small and thus not shown in the lower portion of the cylinder at Ri = 1.0 in Fig. 17(a). The present flow pattern and isotherm distributions are very similar to those reported in [25,35]. To further quantify and validate the present results, the profiles of the velocity components and the temperature at r/R = 0.8 are compared with those from [35] in Fig. 18(a)–(d). It should be noted that all the velocity components are normalized by RX and the dimensionless temperature is T⁄ = (T – Tm)/(Th – Tc) with Tm = (Th + Tc)/2. Good agreement is obtained for each case in Fig. 18. The rapid changes of the velocity and temperature profiles across the boundary layers are fully captured in the present simulations. The present LB model can thus be applied to axisymmetric rotational flows where both the temperature and azimuthal velocity fields can be simulated with the proposed LB model for the axisymmetric CDEs.

4. Conclusions A robust multi-relaxation-time (MRT) lattice Boltzmann model for the axisymmetric convection–diffusion equation (CDE) is proposed. The axsiymmetric CDE is transformed to a standard CDE in the Cartesian coordinate system and the additional terms due to the coordinate transformation are treated as source terms at the macroscopic level. The scalar gradient in the source terms is formulated using the local moment of the non-equilibrium components of the distribution functions without the need for finite-difference type of calculations. When the CDE for the scalar variable / is coupled with the hydrodynamic equations, the radial-axial velocity field is solved with an MRT version of the athermal axisymmetric LB model [20]. The macroscopic velocity components, instead of the velocity distribution functions, are coupled in the evolution equation of the distribution functions for /. In this manner, there is no explicit

coupling between the two sets of distribution functions, making the present model much easier to implement and more robust. The second-order spatial accuracy and first-order temporal accuracy of the proposed LB model is validated with a series of numerical tests. The effect of the curved boundary on the convergence order of the LB solutions is also investigated and verified. The agreement between the present LB results and published numerical computations in the literature for the natural convection problem and the confined swirling flows also demonstrates that the present model can be effectively coupled with athermal LB models to yield accurate results for axisymmetric thermohydrodynamic problems. Acknowledgements This paper was prepared with the support of the U.S. Department of Energy under Award No. DE-FE0001321. Appendix A. Extension of Zhou’s BGK-LB model [20] to an MRTLB model The governing equations for the athermal incompressible axisymmetric flows in the cylindrical coordinate system can be written as

@uj ur ¼ ; @xj r

ðA:1Þ

@ui @ui 1 @p @ 2 ui m @ui mui þ uj ¼ þm 2 þ  2 dir ; @t @xj q @xi r @r r @xj

ðA:2Þ

where i, j are the indices denoting the radial r- and axial z-directions, xi and ui are the components of the spatial vector x(r, z) and the velocity vector u(ur, uz), respectively, q is the density, p is the pressure, and dir is the Kronecker’s delta function. Substitution of the continuity equation (A.1) into Eq. (A.2) gives the following momentum equation [19,20]

  @ui @ 1 @p @ @ui @uj þ ðui uj Þ ¼  þm þ @xj @xj @xi @t @xj q @xi   m @ui @ur ui ur 2mui  þ  2 dir ; þ r @r @xi r r

ðA:3Þ

For incompressible flows in the Cartesian coordinate system, the continuity and momentum equations recovered by the standard LBE read

@uj ¼ 0; @xj

ðA:4Þ

  @ui @ 1 @p @ @ui @uj : þ ðui uj Þ ¼  þm þ @xj @xj @xi @t @xj q @xi

ðA:5Þ

As implied in Zhou’s LB model [20], Eq. (A.3) can be considered

as a @ur i and particular form of Eq. (A.5) with explicit force terms mr @u þ @r @xi  uirur  2rm2ui dir , and another implicit source term accounting for the continuity Eq. (A.1). The LBE for the velocity distribution function fa in Zhou’s BGKLB model [20] is formulated as

fa ðx þ ea dt; t þ dtÞ  fa ðx; tÞ ¼ 



ð2s  1Þear dt 1þ 2r s

1

 ðfa  faðeqÞ Þ þ xa hdt þ

dt

je2

eai F i ;

ðA:6Þ

in order to recover the axisymmetric flow equations. For brevity, the special case of r = 0 discussed in [20] is not considered here. The source or sink term h is

L. Li et al. / International Journal of Heat and Mass Transfer 67 (2013) 338–351

h ¼ qur =r;

ðA:7Þ

and Fi is a force term defined by

Fi ¼ 

qui ur r



2qmui dir : r2

ðA:8Þ

For rotational flows with azimuthal velocity component uu, the force term in Eq. (A.8) should be modified to

Fi ¼ 

qui ur r



qu2/ 2qmui d þ dir : ir r2 r

ðA:9Þ

Here we rearrange Eq. (A.6) as

1 1 fa ðx þ ea dt; t þ dtÞ  fa ðx; tÞ ¼  ðfa  faðeqÞ Þ  r s 2s  1 ear dtðfa  faðeqÞ Þ  2s dt þ xa hdt þ 2 eai F i : je

ðA:10Þ

It is noted in Eq. (A.10) that the first bracketed term on the right ðeqÞ side,  1s ðfa  fa Þ, is the standard BGK collision operator in the ðeqÞ LBE method. The second bracketed term,  1r 2s21 s ear dtðfa  fa Þ, dt and the last term,

ai F i , correspond to the velocity gradient je2 e i r and the force term  uirur  2rm2ui dir in Eq. (stress) term mr @u þ @u @r @xi (A.3), respectively. As discussed earlier, the source term, h, included in Eq. (A.10) is necessary to satisfy the continuity equation (A.1). The proof for the correspondence of each term in the LBE (A.10) to their relevant terms recovered in the momentum Eq. (A.3) can be found in [20]. In this manner, the velocity gradient can be considered as a source term and obtained from the non-equilibrium part of the velocity distribution functions as in Eq. (A.10). This is analogous to the treatment of the scalar gradient, Drrr @/ , in Eq. (2) as a source @r term obtained

Pfrom the non-equilibrium distribution functions, ðeqÞ m  1r 1  2s1rr dx a¼1 ear ½g a  g a , in Eq. (6) in the present LB model dt for the axisymmetric CDE. The BGK model by Zhou [20] can then be conveniently extended to an MRT model with the following evolution equation for fa

fa ðx þ ea dt; t þ dtÞ  fa ðx; tÞ ¼ ½M1 Sðm  mðeqÞ Þa þ Gdt;

ðA:11Þ

where M is a transformation matrix converting the velocity distribution functions f to their moment m by m = Mf, m(eq) is the equilibria of the moments, and S is a matrix of relaxation coefficients. According to Eq. (A.10), the combined source term G in Eq. (A.11) is determined from

G¼

1 2s  1 1 ear ðfa  faðeqÞ Þ þ xa h þ 2 eai F i : r 2s je

ðA:12Þ

The standard D2Q9 MRT-LB model in [27] and [28] can be applied to simulate Eq. (A.11) and it is used in the present work. References [1] E.G. Flekky, Lattice Bhatnagar–Gross–Krook models for miscible fluids, Phys. Rev. E 47 (1993) 4247–4257. [2] R.G.M. van der Sman, M.H. Ernst, Convection–diffusion lattice Boltzmann scheme for irregular lattices, J. Comput. Phys. 160 (2000) 766–782. [3] B. Servan-Camas, F.T.-C. Tsai, Lattice Boltzmann method with two relaxation times for advection–diffusion equation: third order analysis and stability analysis, Adv. Water Resour. 31 (2008) 1113–1126.

351

[4] B.C. Shi, Z. Guo, Lattice Boltzmann model for nonlinear convection–diffusion equations, Phys. Rev. E 79 (2009). 016701. [5] H. Yoshida, M. Nagaoka, Multiple-relaxation-time lattice Boltzmann model for the convection and anisotropic diffusion equation, J. Comput. Phys. 229 (2010) 7774–7795. [6] I. Ginzburg, D. d’Humières, A. Kuzmin, Optimal stability of advection–diffusion lattice Boltzmann models with two relaxation times for positive/negative equilibrium, J. Stat. Phys. 139 (2010) 1090–1143. [7] I. Ginzburg, Multiple anisotropic collisions for advection–diffusion Lattice Boltzmann schemes, Adv. Water Resour. 51 (2013) 381–404. [8] L. Li, R. Mei, J.F. Klausner, Boundary conditions for thermal lattice Boltzmann equation method, J. Comput. Phys. 237 (2013) 366–395. [9] L. Li, R. Mei, J.F. Klausner, Heat transfer evaluation on curved boundaries in thermal lattice Boltzmann equation method, ASME J. Heat Transfer (2013), http://dx.doi.org/10.1115/1.4025046. [10] I. Halliday, L.A. Hammond, C.M. Care, K. Good, A. Stevens, Lattice Boltzmann equation hydrodynamics, Phys. Rev. E 64 (2001) 011208. [11] T.S. Lee, H. Huang, C. Shu, An axisymmetric incompressible lattice Boltzmann model for pipe flow, Int. J. Mod. Phys. C 17 (2006) 645–661. [12] T. Reis, T.N. Phillips, Modified lattice Boltzmann model for axisymmetric flows, Phys. Rev. E 75 (2007) 056703. [13] T. Reis, T.N. Phillips, Numerical validation of a consistent axisymmetric lattice Boltzmann model, Phys. Rev. E 77 (2008) 026703. [14] J.G. Zhou, Axisymmetric lattice Boltzmann method, Phys. Rev. E 78 (2008) 036701. [15] S. Chen, J. Tolke, S. Geller, M. Krafczyk, Lattice Boltzmann model for imcompressible axisymmetric flows, Phys. Rev. E 78 (2008) 046703. [16] H. Huang, X.-Y. Lu, Theoretical and numerical study of axisymmetric lattice Boltzmann models, Phys. Rev. E 80 (2009) 016701. [17] Z. Guo, H. Han, B. Shi, C. Zheng, Theory of the lattice Boltzmann equation: lattice Boltzmann model for axisymmetric flows, Phys. Rev. E 79 (2009) 046708. [18] L. Wang, Z. Guo, C. Zheng, Multi-relaxation-time lattice Boltzmann model for axisymmetric flows, Comput. Fluids 39 (2010) 1542–1548. [19] Q. Li, Y.L. He, G.H. Tang, W.Q. Tao, Improved axisymmetric lattice Boltzmann scheme, Phys. Rev. E 81 (2010) 056707. [20] J.G. Zhou, Axisymmetric lattice Boltzmann method revised, Phys. Rev. E 84 (2011) 036704. [21] Y. Peng, C. Shu, Y.T. Chew, J. Qiu, Numerical investigation of flows in Czochralski crystal growth by an axisymmetric lattice Boltzmann method, J. Comput. Phys. 186 (2003) 295–307. [22] H. Huang, T.S. Lee, C. Shu, Hybrid lattice Boltzmann finite-difference simulation of axisymmetric swirling and rotating flows, Int. J. Numer. Meth. Fluids 53 (2007) 1707–1726. [23] S. Chen, J. Tolke, M. Krafczyk, Simulation of buoyancy-driven flows in a vertical cylinder using a simple lattice Boltzmann model, Phys. Rev. E 79 (2009) 016704. [24] Q. Li, Y.L. He, G.H. Tang, W.Q. Tao, Lattice Boltzmann model for axisymmetric thermal flows, Phys. Rev. E 80 (2009) 037702. [25] L. Zheng, B. Shi, Z. Guo, C. Zheng, Lattice Boltzmann equation for axisymmetric thermal flows, Comput. Fluids 39 (2010) 945–952. [26] L. Zheng, Z. Guo, B. Shi, C. Zheng, Kinetic theory based LBE with viscous dissipation and pressure work for axisymmetric thermal flows, J. Comput. Phys. 229 (2010) 5843–5856. [27] P. Lallemand, L.-S. Luo, Theory of the lattice Boltzmann method: dispersion, dissipation, isotropy, Galilean invariance, and stability, Phys. Rev. E 61 (2000) 6546–6562. [28] D. Yu, R. Mei, L.-S. Luo, W. Shyy, Viscous flow computations with the method of lattice Boltzmann equation, Prog. Aerosp. Sci. 39 (2003) 329–367. [29] P.L. Bhatnagar, E.P. Gross, M. Krook, A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems, Phys. Rev. 94 (1954) 511–525. [30] W.P. Graebel, Engineering Fluid Mechanics, Taylor & Francis, New York, 2001. [31] D.W. Hahn, M.N. Ozisik, Heat Conduction, second ed., John Wiley & Sons, New York, 2012. [32] G. de Vahl Davis, R.W. Thomas, Natural convection between concentric vertical cylinders, Phys. Fluids 12 (Suppl. 2) (1969) 198–207. [33] R. Kumar, M.A. Kalam, Laminar thermal convection between vertical coaxial isothermal cylinders, Int. J. Heat Mass Transfer 34 (1991) 513–524. [34] M. Venkatachalappa, M. Sankar, A.A. Natarajan, Natural convection in an annulus between two rotating vertical cylinders, Acta Mech. 147 (2001) 173– 196. [35] R. Iwatsu, Flow pattern and heat transfer of swirling flows in cylindrical container with rotating top and stable temperature gradient, Int. J. Heat Mass Transfer 47 (2004) 2755–2767.