Study on lattice Boltzmann method/large eddy simulation and its application at high Reynolds number flow

Abstract

Lattice Boltzmann method combined with large eddy simulation is developed in the article to simulate fluid flow at high Reynolds numbers. A subgrid model is used as a large eddy simulation model in the numerical simulation for high Reynolds flow. The idea of subgrid model is based on an assumption to include the physical effects that the unresolved motion has on the resolved fluid motion. It takes a simple form of eddy viscosity models for the Reynolds stress. Lift and drag evaluation in the lattice Boltzmann equation takes momentum-exchange method for curved body surface. First of all, the present numerical method is validated at low Reynolds numbers. Second, the developed lattice Boltzmann method/large eddy simulation method is performed to solve flow problems at high Reynolds numbers. Some detailed quantitative comparisons are implemented to show the effectiveness of the present method. It is demonstrated that lattice Boltzmann method combined with large eddy simulation model can efficiently simulate high Reynolds numbers’ flows.

Introduction

Recently, the lattice Boltzmann method (LBM) has emerged as a well-known alternative of computational technique in fluid dynamics for modeling fluid flow in a way that is consistent with the Navier–Stokes equation,^1,2 due to its intrinsic advantages over conventional Navier–Stokes schemes. The LBM is an innovative numerical method based on kinetic theory to simulate various hydrodynamic systems; it is a reasonable candidate for simulation of turbulence, flow-induced noise, and sound propagation.

The development of the lattice Boltzmann equation (LBE) was independent of the continuous Boltzmann equation. It was introduced to solve some of the difficulties of the lattice gas automata (LGA). A parameter matching procedure based on the Chapman–Enskog analysis of the LGA needs to construct a set of relaxation equations so that the correct hydrodynamic equations are derived.

Compared to the second-order Navier–Stokes equations, the LBM is derived from a set of first-order partial differential equations the LBM formulation does not include nonlinear convection term while Navier–Stokes solvers have to treat such nonlinear one. Moreover, it is necessary for the incompressible Navier–Stokes equations to solve the Poisson’s equation, and then the pressure can be obtained. However, the pressure in the LBM is determined directly from the state equation; the computing cost can be decreased. Furthermore, it has been shown that complex boundary geometries are treated easily in LBM. As a result, LBM has widely emerged as a computational tool for practical engineering applications in the simulations of fluid flows^3,4 and aerospace industries. ⁵

It is well known that the LBM is often used as a direct numerical simulation tool without any assumptions for the relationship between the turbulence stress tensor and the mean strain tensor. Thus, the smallest captured scale in the LBM is the lattice unit, and the largest scale depends on the characteristic length scale in the simulation. These scales are often determined by the available computer memory. Consequently, the LBM is able to resolve relatively low Reynolds number flows. Numerical studies ⁶ have shown that LBM can result in the numerical instability for simulating high Reynolds number flows if unresolved small-scale effects on large-scale dynamics are not considered. A better option is to combine the LBM and large eddy simulation (LES) model in order to solve the problem at high Reynolds numbers. A subgrid model is often used as an LES model in the numerical simulation for traditional Navier–Stokes equation. The idea of subgrid model ⁷ is based on an assumption to include the physical effects that the unresolved motion has on the resolved fluid motion. The model often takes a simple form of eddy viscosity models for the Reynolds stress.

Smagorinsky model ⁷ is one of all subgrid models, including the standard Smagorinsky model and dynamic subgrid model. The standard Smagorinsky model uses a positive eddy viscosity to represent small-scale energy damping. Nevertheless, the shortcoming of the standard one is that the eddy viscosity could be large and positive at some scales smaller than test grids, while it can be large and negative at others sometimes. ⁷ Compared with the standard one, dynamic subgrid models ⁸ could include the dependence of the subgrid model coefficients on local quantities to overcome these side effects.

LBM

The governing equation of LBM for describing fluid flows is Boltzmann equation modeled by adopting the Bhatnagar, Gross, and Krook ⁹ (BGK) collision model. In contrast to the traditional schemes of solving the discretized Navier–Stokes equations, LBM approach focuses on the microscopic scales via the discrete Boltzmann equation and tracks particle distributions on a lattice. The standard LBE involving BGK model can be written into the following form

$$\begin{array}{c}\multicolumn{1}{c}{\frac{f_{\alpha }(x+e_{\alpha x}\delta t,y+e_{\alpha y}\delta t,t+\delta t)-f_{\alpha }(x,y,t)}{\delta t}}\\ \multicolumn{1}{c}{=\frac{f_{\alpha }^{\mathrm{eq}}(x,y,t)-f_{\alpha }(x,y,t)}{\tau },\alpha =0,1,\dots ,M}\end{array}$$

(1)

where $f_{\alpha }$ is the density distribution function along the $\alpha$ -direction at position $(x,y)$ , $f_{\alpha }^{\mathrm{eq}}$ is the local equilibrium distribution function, $\delta t$ is the time-step; $e_{\alpha }(e_{\alpha x},e_{\alpha y})$ is the discrete velocity of particles at time $t$ , $M$ is the number of the direction for the particle velocity, and $\tau$ is the relaxation time, and its definition is related to the kinematics viscosity and temperature by the following formulation

$$\tau =\frac{\nu }{T}+\frac{\Delta t}{2}$$

(2)

The macroscopic quantities, density and momentum density, can be obtained directly from the distribution function $f_{\alpha }$ , and their definitions are in the following form

$$\rho (x,y,t)=\sum _{\alpha =0}^M{f}_{\alpha }(x,y,t),\rho \mathrm{u}(x,y,t)=\sum _{\alpha =0}^M{f}_{\alpha }(x,y,t)e_{\alpha }$$

(3)

The pressure $p$ can be obtained by the state equation of an ideal gas. Under some suppositions, the full compressible viscous Navier–Stokes equations can be derived from LBEs using a multi-scale analysis.

In the simulation, the particle velocity model D2Q9 ¹⁰ is used in the present code. The discrete velocity set $e_{\alpha }$ , equilibrium distribution function $f_{{}_{\alpha }}^{\mathrm{eq}}$ , and weighting factors $w_{\alpha }$ are given by

$$e_{\alpha }=\{\begin{array}{c}\hfill (0,0)\alpha =0\hfill \\ \multicolumn{1}{c}{(\cos {\theta }_{\alpha },\sin {\theta }_{\alpha })c,{\theta }_{\alpha }=(\alpha -1)\pi /2\alpha =1,2,3,4}\\ \multicolumn{1}{c}{\sqrt{2}(\cos {\theta }_{\alpha },\sin {\theta }_{\alpha })c,{\theta }_{\alpha }=(\alpha -5)\pi /2+\pi /4\alpha =5,6,7,8}\end{array}$$

(4)

$$f_{{}_{\alpha }}^{\mathrm{eq}}=w_{\alpha }\rho [1+\frac{3({\mathbf{e}}_{\mathbf{\alpha }}·\mathbf{u})}{c^2}+\frac{9{({\mathbf{e}}_{\mathbf{\alpha }}·\mathbf{u})}^2}{2c^4}-\frac{3{\mathbf{u}}^2}{2c^2}]$$

(5)

$$w_{\alpha }=\{\begin{array}{c}\hfill 4/9,\alpha =0\hfill \\ \multicolumn{1}{c}{1/9,\alpha =1,2,3,4}\\ \multicolumn{1}{c}{1/36,\alpha =5,6,7,8}\end{array}$$

(6)

For the treatment of boundary condition, Zou and He’s ¹¹ non-equilibrium bounce-back scheme is applied on the far-field boundary in the computational domain due to its second-order accuracy. The boundary condition for the particle distribution functions on the curved wall is handled with second-order accuracy based on Mei et al.’s ¹² work. Lift and drag evaluation in the LBE takes momentum-exchange method for curved boundary. Mei has shown that it is reliable, accurate, and easy to implement for both two-dimensional and three-dimensional flows.

Subgrid model for LES

In the LES, the most common subgrid model is the Smagorinsky one, where the anisotropic part of the Reynolds stress term is modeled as

$$-2C{\Delta }^2\left|\overline{S}\right|{\overline{S}}_{\mathrm{ij}}$$

(7)

in which $\left|\overline{S}\right|=\sqrt{2{\overline{S}}_{\mathrm{ij}}{\overline{S}}_{\mathrm{ij}}}$ is the magnitude of the large-scale strain rate tensor

$${\overline{S}}_{\mathrm{ij}}=\frac{1}{2}(\frac{\partial {\overline{u}}_i}{\partial x_j}+\frac{\partial {\overline{u}}_j}{\partial x_i})$$

(8)

$C$ is the Smagorinsky constant, and $\Delta$ is the filter width.

Next, LBM is modified to simulate the filtered physical variables in the LBE. First, the filtered particle distribution ${\overline{f}}_{\alpha }$ is denoted as follows

$${\overline{f}}_{\alpha }(x)=\int f_{\alpha }(x)G(x,x\prime )\mathrm{dx}\prime$$

(9)

Then the LBE can be written into a kinetic equation for the filtered particle distribution function

$$\begin{array}{cc}\multicolumn{1}{c}{}& \frac{{\overline{f}}_{\alpha }(x+e_{\alpha x}\delta t,y+e_{\alpha y}\delta t,t+\delta t)-{\overline{f}}_{\alpha }(x,y,t)}{\delta t}\hfill \\ \multicolumn{1}{c}{}& =\frac{{\overline{f}}_{\alpha }^{\mathrm{eq}}(x,y,t)-{\overline{f}}_{\alpha }(x,y,t)}{{\tau }^{*}}\hfill \end{array}$$

(10)

The above equation for ${\overline{f}}_{\alpha }$ has the similar form as LBE, with the only difference in the relaxation time ${\tau }^{*}$ . Here, the turbulence relaxation time ${\tau }_t$ is introduced into the present relaxation time to consider the effects of small-scale fluid motion. Therefore, the efficient relaxation time can be written into the following

$${\tau }^{*}=\tau +{\tau }_t$$

(11)

Then the total viscosity can be denoted as follows

$${\nu }^{*}=\nu +{\nu }_t$$

(12)

where $\nu =(2\tau -1)/6$ is the laminar viscosity, and ${\nu }_t={\tau }_t/3$ is the turbulence viscosity. According to the present Smagorinsky subgrid model, ${\nu }_t$ should be $C{\Delta }^2\left|\overline{S}\right|$ .

Numerical results

Circular cylinder in low Reynolds numbers’ flow

To validate the program developed in the article, low Reynolds numbers’ flow around a circular cylinder is first simulated here. As depicted in Figure 1, the fluid with the velocity $U=0.1$ flows around the cylinder from the left to the right. $\mathrm{Re}=(\mathrm{UD})/\nu$ is the Reynolds number, and $D$ is diameter of the cylinder. Uniform grids in Figure 2 are used in the simulation; cylinder is immersed in the present grid. Two different cases for $\mathrm{Re}=20$ and $\mathrm{Re}=40$ are shown in the following figures. Figure 3 reveals stream line in the flow field for the different cases; it is clearly seen that their stream lines between upper half-plane and lower half-plane are symmetry. Vorticity contours for the different Re can be found in Figure 4; symmetry vortex pairs attached on the cylinder are also obviously seen in the figure. Drag evaluation in the LBE takes momentum-exchange method for curved boundary. Figure 5 shows the time history of drag for these two cases. Comparisons between the present result and the referred one are given in Table 1; it is further demonstrated from the table that the program developed in the article could be able to simulate low Reynolds flow field, which provides a possibility to model high Reynolds number flow problems using the presented method. Here, the data in Table 1 are extracted at non-dimensional time 20. Next, high Reynolds flows will be tested according to the coupled LBM and LES.

Figure 1.

Computational domain.

Figure 2.

Computational mesh.

Figure 3.

Stream line for Re = 20 and 40.

Figure 4.

Vortex contours for cylinder.

Figure 5.

Time history of Cd for cylinder.

Table 1.

Comparisons of lift and drag.

		Re = 20	Re = 40
Lift	Nieuwstadt and Keller 13	1.786	4.357
	Dennis and Chang 14	1.88	4.69
	The present model	1.85	4.72
Drag	Nieuwstadt and Keller 13	2.053	1.55
	Dennis and Chang 14	2.045	1.52
	The present model	2.24	1.62

Lid-driven cavity flow at high Reynolds numbers

Lid-driven cavity flow is considered here to validate the simulation for high Reynolds flow using the present developed LBM/LES model. It is a classic problem for studying the complex flow in a simple geometry. As shown in Figure 6, flow is driven by the uniform velocity U on the top boundary of square cavity; vortex will be formed inside the cavity depending on the different Reynolds number Re. A uniform grid is used in the present simulation in Figure 6.

Figure 6.

Computational domain and grids.

Figure 7 reveals the stream line and vortex contours at Reynolds number 7500 and U = 0.1. From this figure, the complex vortex structures are found, main vortex lies in the center of the cavity, small vortex also exists in the corner of the cavity around the main vortex, and a second vortex is also clearly seen in the lower right of the cavity. Detailed quantitative comparisons of the vortex position can be found in Table 2 in order to show the ability of the present LBM/LES for simulating high Reynolds flow. Compared with Ghia et al.’s ¹⁵ results, the calculated results in the article can keep consistent with them for the different vortex positions. It is expected that LBM/LES is used to simulate high Reynolds number flow around the complex bodies.

Figure 7.

Stream lines and vortex contours for cavity flow.

Table 2.

Comparisons of vorticity positions inside the cavity.

	Coordinates	Ghia’s result	The present result
Main vortex	X	0.5117	0.5098
	Y	0.5322	0.5267
Lower left vortex	X	0.0645	0.0681
	Y	0.1504	0.1571
Lower right vortex	X	0.7813	0.7953
	Y	0.0625	0.0652

NACA0012 airfoil at high Reynolds numbers

NACA0012 airfoil at high Reynolds flow is simulated at $\mathrm{Re}=5\times {10}^5$ and $\mathrm{Re}=6.6\times {10}^5$ using the present method, where $\mathrm{Re}=\mathrm{UC}/\nu$ , $C$ is the chord of NACA0012 airfoil. $U$ is set to be 0.2 in the numerical simulations. The calculations at two different Reynolds numbers are performed to compare with the experimental result in order to show the ability to simulate high Reynolds numbers’ flow. The angle of attack (AOA) in the first calculation is 7°, and a range of AOA from 0° to 14° is applied in the second one. The computational domain and mesh are depicted in Figure 8, where the different boundary conditions are imposed at the inlet, outlet, top boundary, bottom boundary, and wall. Fixed velocity and pressure are imposed at the inlet and outlet boundaries, respectively. Non-reflective boundary condition is applied at the top and bottom boundaries. In this study, density distribution function at all the far-field boundaries is estimated using the non-equilibrium bounce-back scheme from Zou and He, ¹¹ No-slip boundary is implemented on the body wall, where Guo et al.’s ¹⁶ non-equilibrium extrapolation scheme is used to calculate the density distribution function. Lift and drag evaluation in the LBM also takes momentum-exchange method for curved airfoil boundary.

Figure 8.

Computational domain, boundary, and grid.

Table 3 shows the comparison of lift and drag coefficients between the present and the referred results for $\mathrm{Re}=5\times {10}^5$ and AOA = 7°. From the table, it is clearly seen that the calculated result from LBM/LES could obtain good agreements with others for C_l and C_d from Imamura et al. ¹⁷

Table 3.

Detailed comparisons of C_l and C_d at $\mathrm{Re}=5\times {10}^5$ , AOA = 7°.

Method	PowerFlow	T. Imamura (B-L model)	T. Imamura (laminar)	CFL3D (S-A model)	Present (LBM/LES)
Cl	0.63	0.6979	0.55 ± 0.11	0.7449	0.714
Cd	0.028	0.0177	0.03 ± 0.015	0.0157	0.0162

B-L model: Baldwin-Lomax model; LBM: lattice Boltzmann method; LES: large eddy simulation; S-A model: Spalart-Allmaras model.

Detailed quantitative comparisons are performed to compare with the experimental results at $\mathrm{Re}=6.6\times {10}^5$ , where a range of AOA from 0° to 14° is applied in the simulation. The relationship between C_l , C_d , and AOA is shown in the figure using the presented method. As depicted in Figures 9 and 10, the present result could provide reasonable agreement with the experimental and other referred one from Imamura et al.’s ¹⁷ GILBM before fluid separation happens on the upper surface of the airfoil. However, it is also obviously seen that there exists some differences with the experimental result after separation in Figures 9 and 10. The similar phenomena are also found in the other numerical methods. ¹⁷ The predicted stall angle is 11° here using the present method, which approaches the experimental one.

Figure 9.

Lift coefficient versus different AOA at $\mathrm{Re}=6.6\times {10}^5$ .

Figure 10.

Drag coefficient versus different AOA at $\mathrm{Re}=6.6\times {10}^5$ .

Conclusion

Compared with the other methods, the present method combining the LBM and LES model could solve the problem at high Reynolds numbers’ flows. As one of all subgrid models, Smagorinsky model can provide a positive eddy viscosity to represent small-scale energy damping. For the treatment of boundary condition, non-equilibrium bounce-back scheme has been applied on the far-field boundary in the computational domain. The boundary condition for the particle distribution functions on the curved wall is handled using second-order accuracy scheme. Lift and drag evaluation in the LBE takes momentum-exchange method for curved boundary. The developed program is first to simulate low Reynolds number circular cylinder flow field; it proves that it could provide a possibility to model high Reynolds problems using the present method. Next, lid-driven cavity flow and NACA0012 airfoil flow are calculated to show the ability of combining LBM with LES for simulating high Reynolds number flow. Some comparisons demonstrate that the present calculated results could obtain good agreement with the experimental ones.