`Mathematical Problems in EngineeringVolume 2010, Article ID 724578, 13 pageshttp://dx.doi.org/10.1155/2010/724578`
Research Article

## An Alternative Scheme to Calculate the Strain Rate Tensor for the LES Applications in the LBM

State Key Laboratory of Hydroscience and Engineering, Department of Thermal Engineering, Tsinghua University, Beijing 100084, China

Received 7 September 2010; Revised 17 November 2010; Accepted 18 November 2010

Copyright © 2010 Jun Li and Zhengwei Wang. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

Large eddy simulations (LES) based on the Smagorinsky model can be conveniently used in the lattice Boltzmann method (LBM) because the strain rate tensor, , used to determine the eddy kinematic viscosity can be calculated from the second-order moment of the nonequilibrium distribution function, and the current total nondimensional relaxation time can be determined explicitly. A new method is developed where the distribution function after the relaxation subroutine differs from that after the motion subroutine leading to a similar method to determine , but its application is inconvenient due to the implicit feature. However, the derivation also leads to an alternative explicit scheme for calculating based on physical analysis of the momentum transport process, where the stress tensor, , is calculated first, and then is determined from using the constitutive relationship for Newtonian fluid. The current total nondimensional relaxation time is also given explicitly so that this LES model can be easily used in the LBM.

#### 1. Introduction

The lattice Boltzmann method (LBM) [13] originated from lattice gas (LG) automata [47] can be derived from the Bhatnagar-Gross-Krook (BGK) model equation [8] which is a good approximation to the Boltzmann equation. A general rule for systematically formulating LBM models, which are sufficient for the Navier-Stokes, Burnett fluids, and beyond, was given in [9, 10]. In recent years, many improvements have been made to the LBM with the resulting formulations having computational advantages over traditional continuum methods [1113]. Many successful models have been proposed to extend the scope of LBM applications, including models for incompressible flows [1416] and flows involving with thermal energy exchange [1721]. Some sophisticated solid-fluid boundaries are proposed for regular geometries [22, 23], but there are also general boundary-fitted models [24, 25] available.

Large eddy simulations (LES) are useful for numerical predictions of complex turbulent flows and its application in LBM is referred to as the LBM-LES algorithm. Among the successful subgrid stress models used in LES [26, 27], the Smagorinsky model [28] in LBM-LES is known to be quite convenient because the needed strain rate tensor can be calculated directly from the second-order moment of the nonequilibrium distribution function. In this paper, the subgrid stress model used in LBM-LES is the Smagorinsky model unless otherwise stated. There are many successful calculations [2936] of turbulent flows using LBM-LES, and Yu et al. [31] gave a detailed description of the LBM-LES algorithm.

The calculation of the strain rate tensor is the key step in the implementation of LES in LBM. In the present work, the distribution function after the relaxation subroutine differs from that after the motion subroutine, and then the evolution equation of the LBM algorithm can be written as two different but equivalent formulas based on which two approaches to calculate the strain rate tensor are obtained by mathematical analysis using the Chapman-Enskog expansion [37]. The first approach results in an explicit form which can be easily used in the LBM-LES algorithm. By contrast, the second approach is inconvenient due to its implicit feature. Additionally, an alternative explicit approach for calculating the strain rate tensor is developed based on a physical analysis of the momentum exchange process, to understand which the above distinction is also important. The consistency between the first two approaches is proven theoretically, and numerical results are given to show the consistency between the first and third explicit approaches. Then, predictions of the LBM-LES algorithm-based on the third alternative explicit approach are verified by comparison with a Navier-Stokes equation-based method and the influences of the Smagorinsky constant and cell size on the numerical LBM-LES results are discussed.

#### 2. Two Equivalent Evolution Equations for the LBM Algorithm

The evolution of the distribution function in LBM is split into a series of time steps, and within each time step, is updated to after the motion subroutine and to after the relaxation subroutine. The evolution algorithms for these two different subroutines can be summarized as follows: where is the nondimensional relaxation time. Here, and are used to make the distinction between the values of after the motion and relaxation subroutines during the evolution process to clarify which should be used to calculate the low-order moments needed in the LBM algorithm. Although this distinction is needless and can be neglected in the usual LBM algorithm where only collision-invariant moments, like density and velocity, are involved, it becomes useful when discussing the LBM-LES algorithm where the second-order moment of the nonequilibrium distribution function, which is not collision-invariant, is used to determine the strain rate tensor. The importance of this distinction will be shown in the following discussion.

Substituting (2.2) into (2.1) gives which describes the evolution of . Making an identical transformation with (2.2) by replacing with and then substituting (2.1) into (2.2) gives: which describes the evolution of . Equations (2.3) and (2.4) are equivalent forms of the evolution equation for the LBM algorithm with the small difference between them meaning that the evolution algorithm of is somewhat different from that of .

#### 3. Traditional Mathematical Analysis for the Calculation of

In the present work, the D2Q9 model [3] is used. After making the distinction between and , both (2.3) and (2.4) can be used to get the control equations for the macroquantities by mathematical analysis using Chapman-Enskog and Taylor expansion but the calculational formulas for the strain rate tensor are different. Particularly, using (2.3) gets but using (2.4) leads to a new formula This two formulas have shown the importance of the distinction between and . Both (3.1) using and (3.2) using can be used to calculate the strain rate tensor and there must be a relationship between them. Notice that in (3.1), (3.2), and (2.2) are equivalent to each other and (2.2) can be rewritten as Thus, (3.1) and (3.2) are consistent. Although (3.2) can be obtained directly from (3.1) and (3.3), the independent derivation based on (2.4) using , which is omitted here for conciseness, can help reveal the importance of the distinction between and .

Although (3.1) and (3.2) are equivalent, the application of (3.1) in the LBM-LES algorithm is more conveniently. For convenience, use the notation and to denote their filtered variables of the resolved scale in the LBM-LES algorithm. The eddy kinematic viscosity can be calculated according to Smagorinsky model where is the Smagorinsky constant and the cutoff length is equal to the cell spacing. In the LBM-LES algorithm, the relationship between the nondimensional relaxation time and the kinematic viscosity is: and (3.1) should be changed to and can be determined from (3.4)–(3.6) which is an explicit formula because is known before calculating . At the end of the current time step, can be determined using (2.2), where should be replaced by , which makes the current effective kinematic viscosity implemented in the LBM-LES algorithm equal the sum of the physical kinematic viscosity, , and the current eddy kinematic viscosity, .

If (3.2) is used (note: in (3.2) should also be replaced by as in (3.6)) to determine the strain rate tensor, must be known in advance which means that the current also must be known, but is determined in turn by the strain rate tensor. This implicit relationship can be expressed as which is more complicated than the system of (3.4)–(3.6) used to obtain (3.7) because (2.2) is involved in (3.8) which makes it difficult to obtain an explicit formula similar to (3.7) for the calculation of . Thus, the application of (3.2) in the LBM-LES algorithm is much more difficult than that of (3.1).

#### 4. An Alternative Scheme to Calculate Based on Physical Analysis

The key step to implement LES in LBM is to determine the current so that the current eddy kinematic viscosity and total nondimensional relaxation time can be determined. Section 3 presented two formulas for based on mathematical analysis using the Chapman-Enskog expansion. In addition, can first be determined by a physical analysis of the momentum exchange rate and then can be calculated according to the constitutive relationship for Newtonian fluid which leads to an alternative scheme for calculating .

Generally, is related to the exchange rate of momentum along the -axis through a unit area vertical to the -axis per unit time. Using the original Boltzmann equation theory which is continuous in the phase space where is the molecular velocity similar to in the LBM and is the macrovelocity already used in the LBM. However, of the original Boltzmann equation cannot be simply replaced by or of the LBM.

In the original Boltzmann equation, can be divided into two parts as a function of : is related to the transport from one side to the other side through the imagined interface moving at the macrovelocity and is related to the backward transport through the same interface. Note that there are two sides but only one interface involved when dividing into and , so the net transport of is related to the net momentum exchange through just one interface which can be used to define the stress tensor at the interface by (4.2).

However, the LBM is not continuous in the phase space and although can be divided into two parts as a function of (here the imagined interface is static), is related to the transport from the negative direction of the -axis to the concerned node and is related to the transport from the positive direction of the -axis to the concerned node. Thus, there are three (not two as in the division of of the Boltzmann equation) terms and two (not just one) interfaces involved in the transport information contained in with the severe limitation that contains only the gains at the concerned node from its neighboring nodes but not the backward losses. Thus, cannot be used to construct an equation for the stress tensor which needs the information for both the gains and the losses. Similarly with , is related to the transport from the concerned node in the positive direction of the -axis and is related to the transport from the concerned node in the negative direction of the -axis. Thus, also cannot be used to construct an equation for the stress tensor because it contains only the losses from the concerned node to its neighboring nodes but not the backward gains. Although both and cannot be used independently to construct an equation for the stress tensor, they can be combined to construct an effective distribution function which can be used as in the original Boltzmann equation to calculate the stress tensor.

Only the change that occurs between and its previous in the motion subroutine represents the transport effect, while the change between and its previous in the relaxation subroutine is unrelated to the transport process. The components , of the stress tensor for node in the three-dimensional case at the th time step can be calculated using and to construct and where (or ) contains both the gains and losses of node due to transport from its neighboring node in the positive (or negative) direction of the axis through the interface in the middle. Taking the average of and as to calculate , gives a simple formula: which can also be used to calculate , . In the following discussion, the notation will be omitted for conciseness with the variables defined at point in the phase space unless otherwise stated. The effective density, pressure and macrovelocity are defined from as follows: Now, defined by (4.4) can be used as in (4.2) to calculate with then determined from (4.1) and (3.5) Since and are known for the th time step, , , , , and then can all be calculated explicitly using (4.4)–(4.7). This approach is based on a physical analysis rather than a strict mathematical proof, but its validity will be proven in the following Section using numerical results.

In the LBM-LES algorithm, should be used to replace in (4.7). Since (3.4)–(3.5) are also valid here, the explicit formula for is where and are already known before calculating . At the end of the th time step, can be determined by (2.2) where should be replaced by .

#### 5. Numerical Results

The proposed explicit scheme for the LBM-LES algorithm was validated using the benchmark problem of lid-driven flow in a cubic cavity as shown in Figure 1 where the reliable steady solutions for are given in [38]. Simulations of this problem by LBM were first reported by Hou [39] who did not use the LES algorithm but had a very refined mesh in the LBM method. In the present model, the size of cubic cavity is 1 m, the wall boundary at  m moves in the positive -axis at a velocity 1 m/s and the kinematic viscosity is 0.001 m2/s making so that reliable literature results can be used for comparison. The computational domain was divided into 51-51-51 uniform cells and the time step was 0.001155 s.

Figure 1: Schematic model (a) and velocity profiles on the symmetry plane cited from [38] (b).
##### 5.1. Comparison between the Two Explicit LBM-LES Schemes

Noted earlier, and can be calculated either by (3.6)–(3.7) or by (4.7)–(4.8). The two explicit schemes should have consistent results. Numerical results on the symmetry plane are given in Figure 2 for intermediate results at the 200th time step and in Figure 3 for the final steady results. The results for the two explicit schemes at the initial stage of the evolution process in Figure 2 agree well with each other with only small differences observed between the solid lines and the dashed lines. After convergence as shown in Figure 3, the differences between them almost completely disappear and the solid line overlaps the dashed line. Thus, the current alternative scheme based on (4.7)–(4.8) can be used as equivalent to the traditional scheme based on (3.6)–(3.7). The current alternative scheme has about the same computational cost as the traditional scheme.

Figure 2: Results for at the 200th time step, traditional scheme: solid line, current scheme: dashed line, .
Figure 3: Results for after convergence, traditional scheme: solid line, current scheme: dashed line, .
##### 5.2. Comparison of the Current Explicit Scheme with an N-S Equation-Based Method

Only a half of the flow domain is used in [38] when solving the N-S equation by the finite element method with the symmetry boundary condition at the symmetry plane using both a 51-51-25 uniform mesh and a 51-51-25 nonuniform mesh for . Their results showed that the nonuniform mesh gave more accurate solutions but the iterative solution with the uniform mesh converged faster. The LBM-LES algorithm is most easily used with uniform cells and so the 51-51-51 uniform cells were used in the whole flow domain as in Section 5.1.

The velocity profiles for and along the central axes of the symmetry plane are used to verify the accuracy of the LBM-LES algorithm-based on the current explicit scheme given by (4.7)–(4.8) by comparison with the predictions in [38] using the 51-51-25 uniform mesh. As shown in Section 5.1, the current scheme based on (4.7)–(4.8) is equivalent to the traditional scheme based on (3.6)–(3.7), so the results with the traditional scheme will not be shown here.

The steady-state velocity profiles in Figure 4 show that the LBM-LES algorithm with 51-51-51 uniform cells and agrees well with the N-S equation-based results with the absolute values of the velocity components and predicted by the LBM-LES algorithm being somewhat smaller than the N-S equation results in [38] near the vertical and bottom walls. In addition, this difference increases when is decreased to 0.1, which differs from the conclusion of Yu et al. [31], where was found to yield better energy spectra results than the typical value of in simulations of decaying homogeneous isotropic turbulence. Thus, the influence of on the numerical results may depend on the specific problem and also on which kind of results are considered whether the distribution of macroquantities or the spectra distributions. Therefore, this conclusion is not yet definitive and more work is needed to evaluate these conclusions for other problems.

Figure 4: Steady-state velocity profiles for on the symmetry plane for different .

The influence of the cell size on the result accuracy was evaluated using 65-65-65 uniform cells with . The velocity profiles given in Figure 5 show that the refined cells improve the accuracy of the velocity profiles as expected.

Figure 5: Steady-state velocity profiles for on the symmetry plane for different numbers of cells.

#### 6. Conclusions

An alternative explicit scheme using the average of the two successive values of the distribution functions before and after the motion subroutine is proposed here to calculate the strain rate tensor for the application of LES in the LBM. It was proven to be equivalent to the traditional explicit scheme using the nonequilibrium distribution function and requires about the same computational time as the traditional method. This scheme was not derived using the traditional mathematical analysis but using a physical analysis of the momentum transport process of the LBM algorithm with ideas from the original Boltzmann theory, which provides an alternative perspective on the extension of the LBM algorithm. This physical analysis can also be used to construct LBM models for other problems and its further application is ongoing.

The LBM-LES algorithm-based on the current scheme was validated using the benchmark problem of the lid-driven flow in a cubic cavity by comparison with reliable velocity profiles obtained by the N-S equation-based method. For this problems, was found to yield better velocity profiles than , which differs from the conclusions of Yu et al. [31] obtained by comparing energy spectra results for the simulated decaying homogeneous isotropic turbulence. Thus, the influence of on the numerical results seems to depend on the specific problem and also on the kind of results. These conclusions are then not yet complete and need further study.

#### Acknowledgments

Special thanks are due to the National Natural Science Foundation of China (no. 50879036 and no. 50979044) National High Technology Research and Development Program of China (863 Program: 2009AA05Z424), and State Key Laboratory of Hydroscience and Engineering (no. 2009T3) for supporting the present work. J. Li would like to thank Professor Guixiang Cui and Dr. Lan Xu for helpful communications.

#### References

1. G. R. McNamara and G. Zanetti, “Use of the Boltzmann equation to simulate lattice-gas automata,” Physical Review Letters, vol. 61, no. 20, pp. 2332–2335, 1988.
2. H. Chen, S. Chen, and W. H. Matthaeus, “Recovery of the Navier-Stokes equations using a lattice-gas Boltzmann method,” Physical Review A, vol. 45, no. 8, pp. R5339–R5342, 1992.
3. Y. H. Qian, D. d’Humieres, and P. Lallemand, “Lattice BGK models for Navier-Stokes equation,” Europhysics Letters, vol. 17, no. 6, pp. 479–484, 1992.
4. U. Frisch, B. Hasslacher, and Y. Pomeau, “Lattice-gas automata for the Navier-Stokes equation,” Physical Review Letters, vol. 56, no. 14, pp. 1505–1508, 1986.
5. S. Wolfram, “Cellular automaton fluids 1: basic theory,” Journal of Statistical Physics, vol. 45, no. 3-4, pp. 471–526, 1986.
6. U. Frisch, D. d’Humieres, B. Hasslacher, P. Lallemand, Y. Pomeau, and J.-P. Rivet, “Lattice gas hydrodynamics in two and three dimensions,” Complex Systems, vol. 1, no. 4, pp. 649–707, 1987.
7. D. Montgomery and G. D. Doolen, “Magnetohydrodynamic cellular automata,” Physics Letters A, vol. 120, no. 5, pp. 229–231, 1987.
8. P. L. Bhatnagar, E. P. Gross, and M. Krook, “A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems,” Physical Review, vol. 94, no. 3, pp. 511–525, 1954.
9. X. Shan, X. F. Yuan, and H. Chen, “Kinetic theory representation of hydrodynamics: a way beyond the Navier-Stokes equation,” Journal of Fluid Mechanics, vol. 550, pp. 413–441, 2006.
10. H. Chen and X. Shan, “Fundamental conditions for $N$-th-order accurate lattice Boltzmann models,” Physica D, vol. 237, no. 14–17, pp. 2003–2008, 2008.
11. R. Benzi, S. Succi, and M. Vergassola, “The lattice Boltzmann equation: theory and applications,” Physics Report, vol. 222, no. 3, pp. 145–197, 1992.
12. S. Chen and G. D. Doolen, “Lattice Boltzmann method for fluid flows,” Annual Review of Fluid Mechanics, vol. 30, pp. 329–364, 1998.
13. D. Yu, R. Mei, L. S. Luo, and W. Shyy, “Viscous flow computations with the method of lattice Boltzmann equation,” Progress in Aerospace Sciences, vol. 39, no. 5, pp. 329–367, 2003.
14. Q. Zou, S. Hou, S. Chen, and G. D. Doolen, “A improved incompressible lattice Boltzmann model for time-independent flows,” Journal of Statistical Physics, vol. 81, no. 1-2, pp. 35–48, 1995.
15. X. He and L. S. Luo, “Lattice Boltzmann model for the incompressible Navier-Stokes equation,” Journal of Statistical Physics, vol. 88, no. 3-4, pp. 927–944, 1997.
16. Z. Guo, B. Shi, and N. Wang, “Lattice BGK model for incompressible Navier-Stokes equation,” Journal of Computational Physics, vol. 165, no. 1, pp. 288–306, 2000.
17. X. Shan and H. Chen, “Lattice Boltzmann model for simulating flows with multiple phases and components,” Physical Review E, vol. 47, no. 3, pp. 1815–1819, 1993.
18. X. He, S. Chen, and G. D. Doolen, “A novel thermal model for the lattice Boltzmann method in incompressible limit,” Journal of Computational Physics, vol. 146, no. 1, pp. 282–300, 1998.
19. Z. Guo, B. Shi, and C. Zheng, “A coupled lattice BGK model for the Boussinesq equations,” International Journal for Numerical Methods in Fluids, vol. 39, no. 4, pp. 325–342, 2002.
20. Z. Guo, C. Zheng, B. Shi, and T. S. Zhao, “Thermal lattice Boltzmann equation for low Mach number flows: decoupling model,” Physical Review E, vol. 75, no. 3, Article ID 036704, 2007.
21. Q. Li, Y. L. He, Y. Wang, and G. H. Tang, “An improved thermal lattice Boltzmann model for flows without viscous heat dissipation and compression work,” International Journal of Modern Physics C, vol. 19, no. 1, pp. 125–150, 2008.
22. P. A. Skordos, “Initial and boundary conditions for the lattice Boltzmann method,” Physical Review E, vol. 48, no. 6, pp. 4823–4842, 1993.
23. D. R. Noble, S. Chen, J. G. Georgiadis, and R. O. Buckius, “A consistent hydrodynamic boundary condition for the lattice Boltzmann method,” Physics of Fluids, vol. 7, no. 1, pp. 203–209, 1995.
24. O. Filippova and D. Hänel, “Boundary-fitting and local grid refinement for lattice-BGK models,” International Journal of Modern Physics C, vol. 9, no. 8, pp. 1271–1279, 1998.
25. H. Chen, C. Teixeira, and K. Molvig, “Realization of fluid boundary conditions via discrete Boltzmann dynamics,” International Journal of Modern Physics C, vol. 9, no. 8, pp. 1281–1292, 1998.
26. G. X. Cui, C. X. Xu, L. Fang, L. Shao, and Z. S. Zhang, “A new subgrid eddy-viscosity model for large-eddy simulation of anisotropic turbulence,” Journal of Fluid Mechanics, vol. 582, pp. 377–397, 2007.
27. Z. S. Zhang, G. X. Cui, and C. X. Xu, The Theory and Application of Large-Eddy Simulation of Turbulence, Tsinghua University Press, Beijing, China, 2008.
28. J. Smagorinsky, “General circulation experiments with primitive equation,” Monthly Weather Review, vol. 91, pp. 99–164, 1963.
29. S. Hou, J. Sterling, S. Chen, and G. D. Doolen, “A lattice Boltzmann subgrid model for high Reynolds number flows,” in Pattern Formation and Lattice Gas Automata, vol. 6 of Fields Institute Communications, pp. 151–166, 1996.
30. M. Krafczyk, J. Tölke, and L. S. Luo, “Large-eddy simulations with a multiple-relaxation-time LBE model,” International Journal of Modern Physics B, vol. 17, no. 1-2, pp. 33–39, 2003.
31. H. Yu, S. S. Girimaji, and L. S. Luo, “DNS and LES of decaying isotropic turbulence with and without frame rotation using lattice Boltzmann method,” Journal of Computational Physics, vol. 209, no. 2, pp. 599–616, 2005.
32. Y. H. Dong and P. Sagaut, “A study of time correlations in lattice Boltzmann-based large-eddy simulation of isotropic turbulence,” Physics of Fluids, vol. 20, no. 3, Article ID 035105, 2008.
33. K. N. Premnath, M. J. Pattison, and S. Banerjee, “Generalized lattice Boltzmann equation with forcing term for computation of wall-bounded turbulent flows,” Physical Review E, vol. 79, no. 2, Article ID 026703, 2009.
34. G. Vahala, B. Keating, M. Soe, J. Yepez, L. Vahala, and S. Ziegeler, “Entropic, LES and boundary conditions in lattice Boltzmann simulations of turbulence,” European Physical Journal, vol. 171, no. 1, pp. 167–171, 2009.
35. P. Sagaut, “Toward advanced subgrid models for Lattice-Boltzmann-based Large-eddy simulation: theoretical formulations,” Computers and Mathematics with Applications, vol. 59, no. 7, pp. 2194–2199, 2010.
36. K. N. Premnath, M. J. Pattison, and S. Banerjee, “Dynamic subgrid scale modeling of turbulent flows using lattice-Boltzmann method,” Physica A, vol. 388, no. 13, pp. 2640–2658, 2009.
37. S. Chapman and T. G. Cowling, The Mathematical Theory of Non-Uniform Gases, Cambridge University Press, Cambridge, UK, 3rd edition, 1970.
38. L. Q. Tang, T. Cheng, and T. T. H. Tsang, “Transient solutions for three-dimensional lid-driven cavity flows by a least-squares finite element method,” International Journal for Numerical Methods in Fluids, vol. 21, no. 5, pp. 413–432, 1995.
39. S. Hou, Lattice Boltzmann method for incompressible viscous flow, Ph.D. thesis, Kansas State University, Manhattan, Kan, USA, 1995.