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

Authors: Jun Li; Zhengwei Wang
Journal: Mathematical Problems in Engineering
 (2010)

---

## 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,Sij, 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 Sij,  but its application is inconvenient due to the implicit feature. However, the derivation also leads to an alternative explicit scheme for calculating Sij based on physical analysis of the momentum transport process, where the stress tensor, Tij, is calculated first, and then Sij is determined from Tij 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.


---

## Body

## 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 functionfα in LBM is split into a series of time steps, and within each time step, fα is updated to fmot,α after the motion subroutine and to frel,α after the relaxation subroutine. The evolution algorithms for these two different subroutines can be summarized as follows: (2.1)fmot,α(x⃗+e⃗αΔt,t+Δt)=frel,α(x⃗,t),(2.2)frel,α(x⃗,t)=fmot,α(x⃗,t)+1τ[fαeq(x⃗,t)-fmot,α(x⃗,t)],
where τ is the nondimensional relaxation time. Here, fmot,α and frel,α are used to make the distinction between the values of fα 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 (2.3)fmot,α(x⃗+e⃗αΔt,t+Δt)=fmot,α(x⃗,t)+1τ[fαeq(x⃗,t)-fmot,α(x⃗,t)],
which describes the evolution of fmot,α. Making an identical transformation with (2.2) by replacing (x⃗,t) with (x⃗+e⃗αΔt,t+Δt) and then substituting (2.1) into (2.2) gives: (2.4)frel,α(x⃗+e⃗αΔt,t+Δt)=frel,α(x⃗,t)+1τ[fαeq(x⃗+e⃗αΔt,t+Δt)-frel,α(x⃗,t)],
which describes the evolution of frel,α. 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 fmot,α is somewhat different from that of frel,α.


## 3. Traditional Mathematical Analysis for the Calculation ofSij

In the present work, the D2Q9 model [3] is used. After making the distinction between fmot,α and frel,α, 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 Sij are different. Particularly, using (2.3) gets (3.1)Sij≈-32ρc2τΔt∑αeα,ieα,j(fmot,α-fαeq),
but using (2.4) leads to a new formula (3.2)Sij≈-32ρc2τΔtττ-1∑αeα,ieα,j(frel,α-fαeq).
This two formulas have shown the importance of the distinction between  fmot,α and frel,α. Both (3.1) using fmot,α and (3.2) using frel,α can be used to calculate the strain rate tensor and there must be a relationship between them. Notice that fαeq in (3.1), (3.2), and (2.2) are equivalent to each other and (2.2) can be rewritten as (3.3)ττ-1(frel,α-fαeq)=fmot,α-fαeq.
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 frel,α, which is omitted here for conciseness, can help reveal the importance of the distinction between fmot,α and frel,α.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 f and Sij 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 (3.4)υeddy=(CSΔx)22∑i,jSijSij,
where CS 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: (3.5)υ=13(τ-0.5)c2Δt,υ+υeddy=13[(τ+τeddy)-0.5]c2Δt,
and (3.1) should be changed to(3.6)Sij=-32ρc2(τ+τeddy)Δt∑αeα,ieα,j(fmot,α-fαeq)=-32ρc2(τ+τeddy)ΔtQij,τeddy and τtotal=τ+τeddy can be determined from (3.4)–(3.6) (3.7)τeddy=0.5[τ2+18(CSΔx)2(ρc4Δt2)-12∑i,jQijQij-τ],
which is an explicit formula because Qij=∑αeα,ieα,j(fmot,α-fαeq) is known before calculating τeddy. At the end of the current time step, frel,α 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, frel,α 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 (3.8)Sij=Eq(3.2)F1(frel,α,τtotal)=Eq(2.2)F2(τtotal)=Eq(3.5)F3(υtotal)=Eq(3.4)F4(Sij),
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 CalculateSij Based on Physical Analysis

The key step to implement LES in LBM is to determine the currentSij so that the current eddy kinematic viscosity and total nondimensional relaxation time can be determined. Section 3 presented two formulas for Sij based on mathematical analysis using the Chapman-Enskog expansion. In addition, Tij can first be determined by a physical analysis of the momentum exchange rate and then Sij can be calculated according to the constitutive relationship for Newtonian fluid (4.1)Tij=-pδij+2ρυSij,
which leads to an alternative scheme for calculating Sij.Generally,Tij is related to the exchange rate of momentum along the j-axis through a unit area vertical to the i-axis per unit time. Using the original Boltzmann equation theory which is continuous in the phase space (4.2)Tij=-∫-∞+∞f(ci-ui)(cj-uj)dc⃗,
where ci is the molecular velocity similar to eα,i in the LBM and ui is the macrovelocity already used in the LBM. However, f of the original Boltzmann equation cannot be simply replaced by fmot,α or frel,α of the LBM.In the original Boltzmann equation,f can be divided into two parts as a function of ci-ui: fci-ui>0 is related to the transport from one side to the other side through the imagined interface moving at the macrovelocity ui and fci-ui<0 is related to the backward transport through the same interface. Note that there are two sides but only one interface involved when dividing f into fci-ui>0 and fci-ui<0, so the net transport of f 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 althoughfmot,α can be divided into two parts as a function of eα,i (here the imagined interface is static), (fmot,α)eα,i>0 is related to the transport from the negative direction of the i-axis to the concerned node and (fmot,α)eα,i<0 is related to the transport from the positive direction of the i-axis to the concerned node. Thus, there are three (not two as in the division of f of the Boltzmann equation) terms and two (not just one) interfaces involved in the transport information contained in fmot,α with the severe limitation that fmot,α contains only the gains at the concerned node from its neighboring nodes but not the backward losses. Thus, fmot,α cannot be used to construct an equation for the stress tensor which needs the information for both the gains and the losses. Similarly with frel,α, (frel,α)eα,i>0 is related to the transport from the concerned node in the positive direction of the i-axis and (frel,α)eα,i<0 is related to the transport from the concerned node in the negative direction of the i-axis. Thus, frel,α 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 fmot,α and frel,α cannot be used independently to construct an equation for the stress tensor, they can be combined to construct an effective distribution function fαeff which can be used as f in the original Boltzmann equation to calculate the stress tensor.Only the change that occurs betweenfmot,α and its previous frel,α in the motion subroutine represents the transport effect, while the change between frel,α and its previous fmot,α in the relaxation subroutine is unrelated to the transport process. The components Tij(n,l1,l2,l3), i=1 of the stress tensor for node (l1,l2,l3) in the three-dimensional case at the nth time step can be calculated using frel,α(n-1,l1,l2,l3) and fmot,α(n,l1,l2,l3) to construct fαeff(n,l1+0.5,l2,l3) and fαeff(n,l1-0.5,l2,l3)(4.3)fαeff(n,l1+0.5,l2,l3)={frel,α(n-1,l1,l2,l3),eα,1>0,fmot,α(n,l1,l2,l3),eα,1<0,[frel,α(n-1,l1,l2,l3)+fmot,α(n,l1,l2,l3)]2,eα,1=0,fαeff(n,l1-0.5,l2,l3)={fmot,α(n,l1,l2,l3),eα,1>0,frel,α(n-1,l1,l2,l3),eα,1<0,[frel,α(n-1,l1,l2,l3)+fmot,α(n,l1,l2,l3)]2,eα,1=0,
where fαeff(n,l1+0.5,l2,l3)(or fαeff(n,l1-0.5,l2,l3)) contains both the gains and losses of node (l1,l2,l3) due to transport from its neighboring node in the positive (or negative) direction of the i=1 axis through the interface in the middle. Taking the average of fαeff(n,l1+0.5,l2,l3) and fαeff(n,l1-0.5,l2,l3) as fαeff(n,l1,l2,l3) to calculate Tij(n,l1,l2,l3), i=1 gives a simple formula:(4.4)fαeff(n,l1,l2,l3)=[frel,α(n-1,l1,l2,l3)+fmot,α(n,l1,l2,l3)]2,
which can also be used to calculate Tij(n,l1,l2,l3), i=2,3. In the following discussion, the notation (n,l1,l2,l3) will be omitted for conciseness with the variables defined at point (n,l1,l2,l3) in the phase space unless otherwise stated. The effective density, pressure and macrovelocity are defined from fαeff as follows: (4.5)peff=13c2ρeff=13c2∑αfαeff,uieff=1ρeff∑αfαeffeα,i.
Now, fαeff defined by (4.4) can be used as f in (4.2) to calculate Tij(4.6)Tij=-∑αfαeff(eα,i-uieff)(eα,j-ujeff),
with Sij then determined from (4.1) and (3.5) (4.7)Sij=12υρeff[Tij+peffδij]=32(τ-0.5)c2Δtρeff[Tij+peffδij].
Since frel,α(n-1,l1,l2,l3) and fmot,α(n,l1,l2,l3) are known for the nth time step, fαeff, ρeff, peff, uieff, Tij and then Sij 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 (4.8)τeddy=0.5[(τ-0.5)2+18(CSΔx)2(ρeffc4Δt2)-12∑i,jQijeffQijeff-(τ-0.5)],
where Qijeff=Tij+peffδij and ρeff are already known before calculating τeddy. At the end of the nth time step, frel,α 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 Figure1 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 y=1 m moves in the positive x-axis at a velocity 1 m/s and the kinematic viscosity is 0.001 m2/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.Schematic model (a) and velocity profiles on the symmetry plane cited from [38] (b).
(a)(b)


### 5.1. Comparison between the Two Explicit LBM-LES Schemes

Noted earlier,Sij 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.Results forRe=1000 at the 200th time step, traditional scheme: solid line, current scheme: dashed line, Cs=0.17.
(a)
Sxx(b)
Sxy(c)
Syy(d)
τeddyResults forRe=1000 after convergence, traditional scheme: solid line, current scheme: dashed line, Cs=0.17.
(a)
Sxx(b)
Sxy(c)
Syy(d)
τeddy


### 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 foru and v 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 Figure4 show that the LBM-LES algorithm with 51-51-51 uniform cells and Cs=0.17 agrees well with the N-S equation-based results with the absolute values of the velocity components u and v 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 Cs is decreased to 0.1, which differs from the conclusion of Yu et al. [31], where Cs=0.1 was found to yield better energy spectra results than the typical value of Cs=0.17 in simulations of decaying homogeneous isotropic turbulence. Thus, the influence of Cs 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.Steady-state velocity profiles forRe=1000 on the symmetry plane for different Cs.
(a)(b)The influence of the cell size on the result accuracy was evaluated using 65-65-65 uniform cells withCs=0.17. The velocity profiles given in Figure 5 show that the refined cells improve the accuracy of the velocity profiles as expected.Steady-state velocity profiles forRe=1000 on the symmetry plane for different numbers of cells.
(a)(b)


## 5.1. Comparison between the Two Explicit LBM-LES Schemes

Noted earlier,Sij 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.Results forRe=1000 at the 200th time step, traditional scheme: solid line, current scheme: dashed line, Cs=0.17.
(a)
Sxx(b)
Sxy(c)
Syy(d)
τeddyResults forRe=1000 after convergence, traditional scheme: solid line, current scheme: dashed line, Cs=0.17.
(a)
Sxx(b)
Sxy(c)
Syy(d)
τeddy


## 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 foru and v 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 Figure4 show that the LBM-LES algorithm with 51-51-51 uniform cells and Cs=0.17 agrees well with the N-S equation-based results with the absolute values of the velocity components u and v 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 Cs is decreased to 0.1, which differs from the conclusion of Yu et al. [31], where Cs=0.1 was found to yield better energy spectra results than the typical value of Cs=0.17 in simulations of decaying homogeneous isotropic turbulence. Thus, the influence of Cs 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.Steady-state velocity profiles forRe=1000 on the symmetry plane for different Cs.
(a)(b)The influence of the cell size on the result accuracy was evaluated using 65-65-65 uniform cells withCs=0.17. The velocity profiles given in Figure 5 show that the refined cells improve the accuracy of the velocity profiles as expected.Steady-state velocity profiles forRe=1000 on the symmetry plane for different numbers of cells.
(a)(b)


## 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,Cs=0.17 was found to yield better velocity profiles than Cs=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 Cs 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.



---

*Source: 724578-2010-12-20.xml*