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

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) [1–3] originated from lattice gas (LG) automata [4–7] 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 [11–13]. Many successful models have been proposed to extend the scope of LBM applications, including models for incompressible flows [14–16] and flows involving with thermal energy exchange [17–21]. 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 [29–36] 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 𝑓_mot,𝛼 after the motion subroutine and to 𝑓_rel,𝛼 after the relaxation subroutine. The evolution algorithms for these two different subroutines can be summarized as follows: 𝑓_mot,𝛼𝑥+𝑒_𝛼Δ𝑡,𝑡+Δ𝑡=𝑓_rel,𝛼(),𝑓𝑥,𝑡_rel,𝛼(𝑥,𝑡)=𝑓_mot,𝛼()+1𝑥,𝑡𝜏𝑓^eq_𝛼(𝑥,𝑡)−𝑓_mot,𝛼(𝑥,𝑡), (2.1)(2.2) where 𝜏 is the nondimensional relaxation time. Here, 𝑓_mot,𝛼 and 𝑓_rel,𝛼 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 𝑓_mot,𝛼𝑥+𝑒_𝛼Δ𝑡,𝑡+Δ𝑡=𝑓_mot,𝛼()+1𝑥,𝑡𝜏𝑓^eq_𝛼(𝑥,𝑡)−𝑓_mot,𝛼(𝑥,𝑡), (2.3) which describes the evolution of 𝑓_mot,𝛼. Making an identical transformation with (2.2) by replacing (𝑥,𝑡) with (𝑥+𝑒_𝛼Δ𝑡,𝑡+Δ𝑡) and then substituting (2.1) into (2.2) gives: 𝑓_rel,𝛼𝑥+𝑒_𝛼Δ𝑡,𝑡+Δ𝑡=𝑓_rel,𝛼()+1𝑥,𝑡𝜏𝑓^eq_𝛼𝑥+𝑒_𝛼Δ𝑡,𝑡+Δ𝑡−𝑓_rel,𝛼(𝑥,𝑡), (2.4) which describes the evolution of 𝑓_rel,𝛼. 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 𝑓_mot,𝛼 is somewhat different from that of 𝑓_rel,𝛼.

3. Traditional Mathematical Analysis for the Calculation of 𝑆_𝑖𝑗

In the present work, the D2Q9 model [3] is used. After making the distinction between 𝑓_mot,𝛼 and 𝑓_rel,𝛼, 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 𝑆_𝑖𝑗≈−32𝜌𝑐²𝜏Δ𝑡_𝛼𝑒_𝛼,𝑖𝑒_𝛼,𝑗𝑓_mot,𝛼−𝑓^eq_𝛼, (3.1) but using (2.4) leads to a new formula 𝑆_𝑖𝑗≈−32𝜌𝑐²𝜏𝜏Δ𝑡𝜏−1_𝛼𝑒_𝛼,𝑖𝑒_𝛼,𝑗𝑓_rel,𝛼−𝑓^eq_𝛼. (3.2) This two formulas have shown the importance of the distinction between 𝑓_mot,𝛼 and 𝑓_rel,𝛼. Both (3.1) using 𝑓_mot,𝛼 and (3.2) using 𝑓_rel,𝛼 can be used to calculate the strain rate tensor and there must be a relationship between them. Notice that 𝑓^eq_𝛼 in (3.1), (3.2), and (2.2) are equivalent to each other and (2.2) can be rewritten as 𝜏𝑓𝜏−1^rel,𝛼−𝑓^eq_𝛼=𝑓_mot,𝛼−𝑓^eq_𝛼. (3.3) 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 𝑓_rel,𝛼, which is omitted here for conciseness, can help reveal the importance of the distinction between 𝑓_mot,𝛼 and 𝑓_rel,𝛼.

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 𝜐_eddy=𝐶_𝑆Δ𝑥²2_𝑖,𝑗𝑆_𝑖𝑗𝑆_𝑖𝑗, (3.4) 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: 1𝜐=3(𝜏−0.5)𝑐²Δ𝑡,𝜐+𝜐_eddy=13𝜏+𝜏_eddy−0.5𝑐²Δ𝑡, (3.5) and (3.1) should be changed to𝑆_𝑖𝑗=−32𝜌𝑐²𝜏+𝜏_eddyΔ𝑡_𝛼𝑒_𝛼,𝑖𝑒_𝛼,𝑗𝑓_mot,𝛼−𝑓^eq_𝛼=−32𝜌𝑐²𝜏+𝜏_eddy𝑄Δ𝑡^𝑖𝑗, (3.6)𝜏_eddy and 𝜏_total=𝜏+𝜏_eddy can be determined from (3.4)–(3.6) 𝜏_eddy=0.5𝜏²+18𝐶_𝑆Δ𝑥²𝜌𝑐⁴Δ𝑡²⁻¹2_𝑖,𝑗𝑄_𝑖𝑗𝑄_𝑖𝑗−𝜏, (3.7) which is an explicit formula because 𝑄_𝑖𝑗=∑_𝛼𝑒_𝛼,𝑖𝑒_𝛼,𝑗(𝑓_mot,𝛼−𝑓^eq_𝛼) is known before calculating 𝜏_eddy. At the end of the current time step, 𝑓_rel,𝛼 can be determined using (2.2), where 𝜏 should be replaced by 𝜏_total, 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, 𝜐_eddy.

If (3.2) is used (note: 𝜏 in (3.2) should also be replaced by 𝜏_total as in (3.6)) to determine the strain rate tensor, 𝑓_rel,𝛼 must be known in advance which means that the current 𝜏_total also must be known, but 𝜏_total is determined in turn by the strain rate tensor. This implicit relationship can be expressed as 𝑆=_𝑖𝑗^Eq()3.2𝐹₁𝑓_rel,𝛼,𝜏=_total^Eq()2.2𝐹₂𝜏=_total^Eq()3.5𝐹₃𝜐=_total^Eq()3.4𝐹₄𝑆_𝑖𝑗, (3.8) 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 𝜏_eddy. 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 𝑇_𝑖𝑗=−𝑝𝛿_𝑖𝑗+2𝜌𝜐𝑆_𝑖𝑗, (4.1) 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 𝑇_𝑖𝑗=−^+∞_−∞𝑓𝑐_𝑖−𝑢_𝑖𝑐_𝑗−𝑢_𝑗𝑑𝑐,⃗ (4.2) 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 𝑓_mot,𝛼 or 𝑓_rel,𝛼 of the LBM.

In the original Boltzmann equation, 𝑓 can be divided into two parts as a function of 𝑐_𝑖−𝑢_𝑖: 𝑓_{𝑐𝑖−𝑢𝑖>0} is related to the transport from one side to the other side through the imagined interface moving at the macrovelocity 𝑢_𝑖 and 𝑓_{𝑐𝑖−𝑢𝑖<0} is related to the backward transport through the same interface. Note that there are two sides but only one interface involved when dividing 𝑓 into 𝑓_{𝑐𝑖−𝑢𝑖>0} and 𝑓_{𝑐𝑖−𝑢𝑖<0}, 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 𝑓_mot,𝛼 can be divided into two parts as a function of 𝑒_𝛼,𝑖 (here the imagined interface is static), (𝑓_mot,𝛼)_{𝑒𝑖𝛼,>0} is related to the transport from the negative direction of the 𝑖-axis to the concerned node and (𝑓_mot,𝛼)_{𝑒𝑖𝛼,<0} 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 𝑓_mot,𝛼 with the severe limitation that 𝑓_mot,𝛼 contains only the gains at the concerned node from its neighboring nodes but not the backward losses. Thus, 𝑓_mot,𝛼 cannot be used to construct an equation for the stress tensor which needs the information for both the gains and the losses. Similarly with 𝑓_rel,𝛼, (𝑓_rel,𝛼)_{𝑒𝑖𝛼,>0} is related to the transport from the concerned node in the positive direction of the 𝑖-axis and (𝑓_rel,𝛼)_{𝑒𝑖𝛼,<0} is related to the transport from the concerned node in the negative direction of the 𝑖-axis. Thus, 𝑓_rel,𝛼 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 𝑓_mot,𝛼 and 𝑓_rel,𝛼 cannot be used independently to construct an equation for the stress tensor, they can be combined to construct an effective distribution function 𝑓^eff_𝛼 which can be used as 𝑓 in the original Boltzmann equation to calculate the stress tensor.

Only the change that occurs between 𝑓_mot,𝛼 and its previous 𝑓_rel,𝛼 in the motion subroutine represents the transport effect, while the change between 𝑓_rel,𝛼 and its previous 𝑓_mot,𝛼 in the relaxation subroutine is unrelated to the transport process. The components 𝑇_𝑖𝑗(𝑛,𝑙₁,𝑙₂,𝑙₃), 𝑖=1 of the stress tensor for node (𝑙₁,𝑙₂,𝑙₃) in the three-dimensional case at the 𝑛th time step can be calculated using 𝑓_rel,𝛼(𝑛−1,𝑙₁,𝑙₂,𝑙₃) and 𝑓_mot,𝛼(𝑛,𝑙₁,𝑙₂,𝑙₃) to construct 𝑓^eff_𝛼(𝑛,𝑙₁+0.5,𝑙₂,𝑙₃) and 𝑓^eff_𝛼(𝑛,𝑙₁−0.5,𝑙₂,𝑙₃)𝑓^eff_𝛼𝑛,𝑙₁+0.5,𝑙₂,𝑙₃=𝑓_rel,𝛼𝑛−1,𝑙₁,𝑙₂,𝑙₃,𝑒_𝛼,1𝑓>0,_mot,𝛼𝑛,𝑙₁,𝑙₂,𝑙₃,𝑒_𝛼,1𝑓<0,_rel,𝛼𝑛−1,𝑙₁,𝑙₂,𝑙₃+𝑓_mot,𝛼𝑛,𝑙₁,𝑙₂,𝑙₃2,𝑒_𝛼,1𝑓=0,_eff𝛼𝑛,𝑙₁−0.5,𝑙₂,𝑙₃=𝑓_mot,𝛼𝑛,𝑙₁,𝑙₂,𝑙₃,𝑒_𝛼,1𝑓>0,_rel,𝛼𝑛−1,𝑙₁,𝑙₂,𝑙₃,𝑒_𝛼,1𝑓<0,_rel,𝛼𝑛−1,𝑙₁,𝑙₂,𝑙₃+𝑓_mot,𝛼𝑛,𝑙₁,𝑙₂,𝑙₃2,𝑒_𝛼,1=0, (4.3) where 𝑓^eff_𝛼(𝑛,𝑙₁+0.5,𝑙₂,𝑙₃)(or 𝑓^eff_𝛼(𝑛,𝑙₁−0.5,𝑙₂,𝑙₃)) contains both the gains and losses of node (𝑙₁,𝑙₂,𝑙₃) due to transport from its neighboring node in the positive (or negative) direction of the 𝑖=1 axis through the interface in the middle. Taking the average of 𝑓^eff_𝛼(𝑛,𝑙₁+0.5,𝑙₂,𝑙₃) and 𝑓^eff_𝛼(𝑛,𝑙₁−0.5,𝑙₂,𝑙₃) as 𝑓^eff_𝛼(𝑛,𝑙₁,𝑙₂,𝑙₃) to calculate 𝑇_𝑖𝑗(𝑛,𝑙₁,𝑙₂,𝑙₃), 𝑖=1 gives a simple formula:𝑓^eff_𝛼𝑛,𝑙₁,𝑙₂,𝑙₃=𝑓_rel,𝛼𝑛−1,𝑙₁,𝑙₂,𝑙₃+𝑓_mot,𝛼𝑛,𝑙₁,𝑙₂,𝑙₃2, (4.4) which can also be used to calculate 𝑇_𝑖𝑗(𝑛,𝑙₁,𝑙₂,𝑙₃), 𝑖=2,3. 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 𝑓^eff_𝛼 as follows: 𝑝^eff=13𝑐²𝜌^eff=13𝑐²_𝛼𝑓^eff_𝛼,𝑢^eff_𝑖=1𝜌^eff_𝛼𝑓^eff_𝛼𝑒_𝛼,𝑖. (4.5) Now, 𝑓^eff_𝛼 defined by (4.4) can be used as 𝑓 in (4.2) to calculate 𝑇_𝑖𝑗𝑇_𝑖𝑗=−_𝛼𝑓^eff_𝛼𝑒_𝛼,𝑖−𝑢^eff_𝑖𝑒_𝛼,j−𝑢_j^eff, (4.6) with 𝑆_𝑖𝑗 then determined from (4.1) and (3.5) 𝑆_𝑖𝑗=12𝜐𝜌^eff𝑇_𝑖𝑗+𝑝^eff𝛿_𝑖𝑗=32(𝜏−0.5)𝑐²Δ𝑡𝜌^eff𝑇_𝑖𝑗+𝑝^eff𝛿_𝑖𝑗. (4.7) Since 𝑓_rel,𝛼(𝑛−1,𝑙₁,𝑙₂,𝑙₃) and 𝑓_mot,𝛼(𝑛,𝑙₁,𝑙₂,𝑙₃) are known for the 𝑛th time step, 𝑓^eff_𝛼, 𝜌^eff, 𝑝^eff, 𝑢^eff_𝑖, 𝑇_𝑖𝑗 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, 𝜏_total=𝜏+𝜏_eddy should be used to replace 𝜏 in (4.7). Since (3.4)–(3.5) are also valid here, the explicit formula for 𝜏_eddy is 𝜏_eddy=0.5(𝜏−0.5)²+18𝐶_𝑆Δ𝑥²𝜌^eff𝑐⁴Δ𝑡²⁻¹2_𝑖,𝑗𝑄^eff_𝑖𝑗𝑄^eff_𝑖𝑗−(𝜏−0.5), (4.8) where Q^eff_𝑖𝑗=𝑇_𝑖𝑗+𝑝^eff𝛿_𝑖𝑗 and 𝜌^eff are already known before calculating 𝜏_eddy. At the end of the 𝑛th time step, 𝑓_rel,𝛼 can be determined by (2.2) where 𝜏 should be replaced by 𝜏_total.

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 Re=1000 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 𝑦=1 m moves in the positive 𝑥-axis at a velocity 1 m/s and the kinematic viscosity is 0.001 m²/s making Re=1000 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.

(a)

(b)

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 𝜏_eddy 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.

(a) 𝑆𝑥𝑥

(b) 𝑆𝑥𝑦

(c) 𝑆𝑦𝑦

(d) 𝜏eddy

Figure 2 

Results for Re=1000 at the 200th time step, traditional scheme: solid line, current scheme: dashed line, 𝐶𝑠=0.17.

(a) 𝑆𝑥𝑥

(b) 𝑆𝑥𝑦

(c) 𝑆𝑦𝑦

(d) 𝜏eddy

Figure 3 

Results for Re=1000 after convergence, traditional scheme: solid line, current scheme: dashed line, 𝐶𝑠=0.17.

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 Re=1000. 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 𝐶_𝑠=0.17 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 𝐶_𝑠=0.1 was found to yield better energy spectra results than the typical value of 𝐶_𝑠=0.17 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.

(a)

(b)

Figure 4 

Steady-state velocity profiles for Re=1000 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 𝐶_𝑠=0.17. The velocity profiles given in Figure 5 show that the refined cells improve the accuracy of the velocity profiles as expected.

(a)

(b)

Figure 5 

Steady-state velocity profiles for Re=1000 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, 𝐶_𝑠=0.17 was found to yield better velocity profiles than 𝐶_𝑠=0.1, 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.