Study on LBM/LES and Its Application at High Reynolds Number
Flow
a
Si Haiqing a,1, Shi Yan
College of Civil Aviation and Flight, Nanjing University of Aeronautics and Astronautics,
Nanjing, China, 210016
ABSTRACT
Lattice Boltzmann method (LBM) combining with LES is developed in the paper
to simulate fluid flow at high Reynolds numbers. A subgrid model is used as an LES
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 LBM/LES 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 LBM combining with
LES model can efficiently simulate high Reynolds numbers flows.
KEYWORDS: Lattice Boltzmann Method; Subgrid Model; Wall Boundary Condition;
Computational Fluid Dynamics;
1. 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;
1
Corresponding author: Si Haiqing, E-mail address: [email protected] (H.Q. Si)
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 need 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 PDE’s, 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
equation, and then the pressure can be obtained. However, 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 been widely emerged as a computational tool for
practical engineering applications in the simulations of fluid flows[3-4] and aerospace
industries[5].
It is well known that the LBM is often used as a direct numerical simulation tool
without any assumptions for the relation 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[6]
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[7] 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[7] 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[7]. Compared with the standard one, dynamic subgrid models[8]
could include the dependence of the subgrid model coefficients on local quantities to
overcome these side effects.
2. LATTICE-BOLTZMANN METHOD
The governing equation of LBM for describing fluid flows is Boltzmann equation
modeled by adopting the Bhatnagar, Gross, and Krook[9] (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 Lattice Boltzmann equation
involving of BGK model can be written into the following form:
fα ( x + eα xδ t , y + eα yδ t , t + δ t ) − fα ( x, y, t )
δt
=
fαeq ( x, y, t ) − fα ( x, y, t )
τ
, α = 0,1,L , M (1)
Where, fα is the density distribution function along the α -direction at position ( x, y ) ;
fαeq is local equilibrium distribution function; δ t is time-step; eα (eα x , eα y ) is the
discrete velocity of particles at time t ; M is the number of the direction for the
particle velocity; τ is the relaxation time, its definition is related to the kinematics
viscosity and temperature by the following formulation:
τ=
ν
T
+
Δt
2
(2)
The macroscopic quantities, density and momentum density, can be obtained
directly from the distribution function fα , their definitions are in the following form:
M
M
α =0
α =0
ρ ( x, y, t ) = ∑ fα ( x, y, t ) , ρ u( x, y, t ) = ∑ fα ( x, y, t )eα
(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 Lattice Boltzmann equations using a multi-scale analysis.
In the simulation, the particle velocity model D2Q9 [10] is used in the present code.
The discrete velocity set eα , equilibrium distribution function f αeq and weighting
factors wα are given by
⎧(0, 0)
α =0
⎪
eα = ⎨(cos θα ,sin θα )c, θα = (α − 1)π / 2
α = 1, 2,3, 4
⎪
⎩ 2(cos θα ,sin θα )c, θα = (α − 5)π / 2 + π / 4 α = 5, 6, 7,8
⎡ 3(e ⋅ u) 9(eα ⋅ u) 2 3u 2 ⎤
f αeq = wα ρ ⎢1 + α2 +
− 2⎥
4
2
2c ⎦
c
c
⎣
⎧4 / 9, α = 0
⎪
wα = ⎨1/ 9, α = 1, 2,3, 4
⎪1/ 36, α = 5, 6, 7,8
⎩
.
(8)
(9)
(10)
For the treatment of boundary condition, Zou’s nonequilibrium bounce-back
scheme[11] 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 curved wall is handled with second-order accuracy based on Mei’s work[12]. Lift
and drag evaluation in the lattice Boltzmann equation 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.
3. SUBGRID MODEL FOR LES
In the LES, the most common subgrid model is Smagorinsky one, where the
anisotropic part of the Reynolds stress term is modeled as
−2C Δ 2 S Sij
(11)
in which S = 2Sij Sij is the magnitude of the large scale strain rate tensor
1 ⎛ ∂u ∂u ⎞
Sij = ⎜ i + j ⎟ ,
2 ⎜⎝ ∂x j ∂xi ⎟⎠
(12)
C is the Smagorinsky constant, Δ is the filter width.
Next, LBM is modified to simulate the filtered physical variables in the LBE.
Firstly, the filtered particle distribution fα , is denoted as follows:
fα ( x) = ∫ fα ( x)G ( x, x′)dx′
(13)
Then the LBE can be written into a kinetic equation for the filtered particle
distribution function:
fα ( x + eα xδ t , y + eα yδ t , t + δ t ) − fα ( x, y, t )
δt
=
fαeq ( x, y, t ) − fα ( x, y, t )
(14)
τ∗
The above equation for fα has the similar form as LBE, with the only difference in
the relaxation time τ ∗ . Here, the turbulence relaxation time τ 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:
τ ∗ = τ +τt
(15)
Then the total viscosity can be denoted as follows:
ν ∗ = ν +ν t
Where ν =
(16)
τ
2τ − 1
is laminar viscosity, ν t = t is turbulence viscosity. According to
3
6
the present Smagorinsky subgrid model,ν t should be C Δ 2 S .
4. NUMERICAL RESULTS
4.1 Circular cylinder in low Reynolds numbers flow
To validate the program developed in the paper, low Reynolds numbers flow
around a circular cylinder is firstly simulated here. As depicted in Fig.1, fluid with the
velocity U = 0.1 flows around the cylinder from the left to the right. Re =
UD
ν
is
Reynolds number, and D is diameter of the cylinder. Uniform grids in Fig.2 are used
in the simulation; cylinder is immersed in the present grid. Two different cases for
Re = 20 and Re = 40 are shown in the following figures. Fig. 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 found in Fig.4, symmetry vortex pairs attached on the cylinder are
also obviously seen in the figure. Drag evaluation in the lattice Boltzmann equation
takes momentum-exchange method for curved boundary. Fig.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 paper 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 is extracted at non-dimensional time 20.
Next, high Reynolds flows will be tested according to the coupled LBM and LES.
Fig.1 Computational domain
Fig.2 Computational mesh
Re=20
Re=40
Fig.3Stream line for Re=20 and 40
Re=20
Re=40
Fig.4Vortex contours for cylinder
Re=20
Re=40
Fig.5Time history of Cd for cylinder
Table1 Comparisons of lift and drag
Lift
Drag
Re=20
Re=40
Keller and Nieuwstadt [13]
1.786
4.357
Dennis and Chang [14]
1.88
4.69
The present
1.85
4.72
Keller and Nieuwstadt [13]
2.053
1.55
Dennis and Chang [14]
2.045
1.52
The present
2.24
1.62
4.2 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 Fig.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 Fig.6.
Fig.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 the Table 2
in order to show the ability of the present LBM/LES for simulating high Reynolds
flow. Compared with Ghia’s results[15], the calculated results in the paper 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.
Fig.6 Computational domain and grids
Fig.7 Stream lines and vortex contours for cavity flow
Table 2 Comparisons of Vorticity positions inside the cavity
Main vortex
Lower left
Coordinates
Ghia’s result
The present result
X
0.5117
0.5098
Y
0.5322
0.5267
X
0.0645
0.0681
Y
0.1504
0.1571
X
0.7813
0.7953
Y
0.0625
0.0652
vortex
Lower right
vortex
4.2 NACA0012 airfoil at high Reynolds numbers
NACA0012 airfoil at high Reynolds flow is simulated at Re = 5 ×105 and
Re = 6.6 ×105 using the present method, where Re = UC / ν , 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 experimental result
in order to show the ability to simulate high Reynolds numbers flow. The angle of
attack (AOA) in the first calculation is 70, a range of AOA from 00 to 140 is applied in
the second one. The computational domain and mesh are depicted in Fig.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
boundary, respectively. Non-reflective boundary condition is applied at the top and
bottom boundary. In the present study, density distribution function at all the far field
boundaries is estimated using the non-equilibrium bounce-back scheme from Zou and
He[11]. No-slip boundary is implemented on the body wall, where Guo’s
non-equilibrium extrapolation scheme[16] is used to calculate the density distribution
function. Lift and drag evaluation in the lattice Boltzmann method also takes
momentum-exchange method for curved airfoil boundary.
Table 3 shows the comparison of lift and drag coefficients between the present and
the referred results for Re = 5 ×105 and AOA=70. From the table, it is clearly seen
that the calculated result from LBM/LES could obtain good agreements with others
for Cl and Cd from Imamura’s[17].
Fig.8 Computational domain, boundary and grid
Table 3 Detailed comparisons of Cl and Cd at Re = 5 ×105 , AOA=70
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
Detailed quantitative comparisons are performed to compare with experimental
results at Re = 6.6 ×105 , where a range of AOA from 00 to 140 is applied in the
simulation. The relation between Cl, Cd and AOA would be shown in the figure using
the presented method. As depicted in the Figs.9-10, the present result could provide
reasonable agreement with experimental and other referred one from Imamura’s
GILBM[17] 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 Fig.9-10. The similar phenomena are also
found in the other numerical methods[17]. The predicted stall angle is 110 here using
the present method, which approaches the experimental one.
Experimental data
GILBM
Present
1.2
1
Cl
0.8
0.6
0.4
0.2
0
0
5
10
15
AOA
Fig.9 Lift coefficient vs different AOA at Re = 6.6 × 105
0.035
Experimental data
GILBM
Present
0.03
Cd
0.025
0.02
0.015
0.01
0
5
10
15
AOA
Fig.10 Drag coefficient vs different AOA at Re = 6.6 × 105
5. CONCLUSIONS
Compared with the other methods, the present method combining the LBM and
large eddy simulation (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 curved wall is handled using second-order accuracy
scheme. Lift and drag evaluation in the lattice Boltzmann equation takes
momentum-exchange method for curved boundary. The developed program is firstly
to simulate low Reynolds number circular cylinder flow field, it proves that it could
provides 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 combination LBM with LES for simulating high Reynolds number flow.
Some comparisons demonstrate that the present calculated results could obtain good
agreement with experimental ones.
ACKNOWLEDGEMENTS
This work was supported by Special Foundation of China Postdoctoral Science
(Grant No.: 201104565) and National Natural Science Foundation of China (Grant
No.: 11272151).
CONFLICT OF INTERESTS
The authors declare that there is no conflict of interests regarding the publication of
this paper.
REFERRENCES
[1] Shan, X., Yuan X.F., and Chen, H., Kinetic theory representation of hydrodynamics: a way
beyond the Navier-Stokes equation, Journal of Fluid Mechanic, 550: 413-441, 2006.
[2] Chen, S. and Doolen, G., Lattice Boltzmann method for fluid flows, Ann. Rev. Fluid Mech., 30:
329-364, 1998.
[3] Adam, J.-L., Ricot, D., Dubief, F., and Guy, C., Aeroacoustic simulation of automobile
ventilation outlets, JASA. 123(5):3250, 2008.
[4] Balasubramanian, G., Crouse, B., and Freed, D., Numerical simulation of real world effects on
sunroof buffeting of an idealized generic vehicle, AIAA Paper 2009-3348, 2009.
[5] Keating A., Dethioux, P., Satti, R., et al., Computational aeroacoustics validation and analysis
of a nose landing gear, AIAA Paper 2009-3154, 2009.
[6] Sterling J.D. and Chen S., Stability analysis of lattice Boltzmann methods, Journal of
Computational Physics, 123: 196-206, 1996.
[7] Hou S., et al., A lattice Boltzmann subgrid model for high Reynolds number flows, Pattern
Formation and Lattice Gas Automata, 6: 149, 1996.
[8] Germano M., Piomeelli U., Moin P., et al., A dynamic subgrid-scale eddy viscosity model,
Physics of Fluids A, 3(7): 1760-1765,1991.
[9] Bhatnagar, P., Gross, E.P., and Krook, M.K., A model for collision processes in Gases, 1. small
amplitude processes in charged and neutral one-component systems, Physics Review, 94(3):
515-525, 1954.
[10] He X. and Luo L. S., Theory of the lattice Boltzmann method: from the Boltzmann equation
to the lattice Boltzmann equation, Phys. Rev. E, 56: 6811-6817, 1997.
[11] Zou Q. and He X., On pressure and velocity boundary conditions for the lattice Boltzmann
BGK model, Physics of Fluids, 9(6): 1591-1598. 1997.
[12] Mei R., Yu D. and Shyy W., Force evaluation in the lattice Boltzmann method involving
curved geometry, Physical Review E, 65(041203):1-13, 2002.
[13] Nieuwstadt F., Keller H.B., Viscous flow past circular cylinder, Computers and Fluids,
42:471-489, 1973.
[14] Dennis S.C.R., Chang G.Z., Numerical solutions for steady flow past a circular cylinder
at Reynolds number up to 100, Journal of Fluid Mechanics, 42:471-489, 1970.
[15] Ghia U., Ghia K N. and Shin C T. High-Re solutions for incompressible flow using the
Navier-Stokes equations and a multigrid method, Journal of Computational Physics, 48(3):
387-410, 1982.
[16] Guo ZL, Zheng CG, Shi BC, Non-equilibrium extrapolation method for velocity and
boundary conditions in the lattice Boltzmann method, Chinese Physics, 11(4):366-374, 2002.
[17] Imamura T., Suzuki K., Nakamura T., et al., Flow simulation around an airfoil using lattice
Boltzmann method on generalized coordinates, AIAA-2004-244, 2004.