Log in using OpenID

Transactions of the Korean Nuclear Society Spring Meeting
Gyeongju, Korea, May 29-30, 2008
Finite Element Based Formulation of Lattice Boltzmann Equation
Jong Chull Joa*, Kyung Wan Roha and Young W. Kwonb
Safety Issue Research Dept., Korea Institute of Nuclear Safety, 19 Kusung-dong, Yusung-gu, Daejon 305-338, Korea
Dept. of Mechanical & Astronautical Engineering, Naval Postgraduate School, Monterey, California, 93943, USA
*Corresponding author : *
1. Introduction
The Lattice Boltzmann Method (LBM) has been developed for
application to thermal-fluid problems. Recently, the technique was also
applied to fluid-structure interaction problems [1]. Most of those studies
considered a regular shape of lattice or mesh like square and cubic grids.
In order to apply the LBM to more practical cases, it is necessary to be
able to solve complex or irregular shapes of problem domains. There
have been different kinds of approaches to address the problems. The
most common technique was using the finite volume formulation of the
lattice Boltzmann equation.[2,3] Another approach was a point-wise
interpolation technique for irregular grids.[4] Other techniques were
based on the finite element method. [5,6]
Generally, the finite element method is very powerful for solving twoor three-dimensional complex or irregular shapes of domains using the
isoparametric element formulation [7] which is based on a mathematical
mapping from a regular shape of element in an imaginary domain to a
more general and irregular shape of element in the physical domain. In
addition, there are variety of choices of finite elements such as triangular
or quadrilateral shapes in 2-D, or tetrahedral, triangular prism, or general
six-sided solids in 3-D. As a result, the present study presents a new finite
element formulation for the lattice Boltzmann equation using the general
weighted residual technique. Among the weighted residual formulations,
the collocation method, Galerkin method or method of moments are used
to develop the finite element based LBM.
2 cos{(α − 1)π / 2 + π / 4}, 2 sin{(α − 1)π / 2 + π / 4}
− fα
in which ρ is the fluid density and
expressed as
ρ = ∑ fα
i =1
in which Hi is the spatial interpolation function for the nodal variable
f αi
at the i-th node of the finite element, and n is the number of nodes per
element. In addition, [H ] is a row vector consisting of the interpolation
functions, and { fα } is a column vector containing unknown solutions at
the nodes. Plugging Eq. (9) into Eq. (8) yields
[H ]{f&α }+ erα ⋅ [∇H ]{ f α } + 1 [H ]({ f α } − {~fα }) = 0
for each finite element. The superimposed dot denotes the temporal
Applying the weighted residual formulation to Eq. (10) gives the
following expression
~ ⎞
∑ ∫S {w}⎜⎝ [H ]{f&α }+ eα ⋅ [∇H ]{ f α } + τ [H ] { f α } − f α ⎟⎠dS = 0
where the integration is conducted over each finite element domain S e
α = 5 to 8
denotes the local equilibrium
{ })
and the summation is performed over the total number of elements.
Furthermore, {w} is a column vector of the weighting functions. The size
of {w} is equal to the number of nodes per element. Rewriting Eq. (11)
[M ]{F&α }+ [K ]{Fα } + [C ]{Fα } − [C ]{F~α } = 0
is the fluid velocity. They can be
ρ v = ∑ fα eα
f α = ∑ H i f αi = [H ]{ f α }
distribution of f α . The local equilibrium distribution is
⎡ 3vr ⋅ er 9 ( vr ⋅ erα )2 3vr ⋅ vr ⎤
f%α = ρωα ⎢1 + 2 α +
2c 4
2c 2 ⎥⎦
interpolation functions and nodal variables as given below:
where τ is the relaxation constant and
In order to derive the Finite Element Lattice Boltzmann Method
(FELBM) from Eq. (8), the problem domain is discretized into a number
of finite elements. Then, the variable f α is expressed in terms of the
where c is the lattice speed.
When a single relaxation time technique is used for the collision
operator, the collision operator can be written as
Ωα = −
direction, and Ωα denotes the collision operator. The discrete velocity
vector is given below and shown in Fig. 1 for the D2Q9.
In addition, ωα is the weighting parameter for each velocity direction as
given below:
α =0 ⎞
⎛ 4/9
ωα = ⎜⎜ 1/ 9 α = 1, 2,3, 4 ⎟⎟
⎜ 1/ 36 α = 5, 6, 7,8 ⎟
for D2Q9.
Substitution of Eq. (3) into Eq. (1) results in
∂f α r
+ e ⋅ ∇f + ( f − f ) = 0
where fi is the particle velocity distribution function along the α−
direction, t represents time, eα is the discrete velocity vector along the α−
α =0
α = 1 to 4
Figure 1. D2Q9 lattice showing discrete nine velocity vectors
(0, 0)
c ( cos{(α − 1)π / 2},sin{(α − 1)π / 2})
1 (1,0)
The lattice Boltzmann equation is expressed as
∂f α r
+ eα ⋅ ∇f α = Ωα
2. Finite Element Based Lattice Boltzmann Method
r ⎜⎜
eα =
[M ] = ∑ [m] = ∑ ∫S {w}[H ]dS
[K ] = ∑ [k ] = ∑ ∫S {w}(erα ⋅ [∇H ])dS
Transactions of the Korean Nuclear Society Autumn Meeting
Pyeong Chang, Korea, October xx-xx, 2008
[C ] = ∑ [c] = ∑ ∫S 1 {w}[H ]dS
{Fα } = ∑ { f α }
Fα = ∑ f α
{ }
{ }
Norm. Velocity
Depending on the choice of the weighting functions, the subsequent
technique can be called the Galerkin method, collocation method, method
of moments, least-square method, or sub-domain method. In this study,
the first three techniques are presented.
For the Galerkin method, the weighting function is selected to be the
interpolation functions as used in Eq. (11), i.e. {w} = [H ]T . In this case,
Eqs. (13) through (15) can be expressed as
[M ] = [m] =
[H ] [H ]dS
Figure 2. Comparison of steady state velocity profile of flow between two
co-axial cylinders using four-node quadrilateral element.
Norm. Dist.
[C ] = ∑ [c ] = ∑ ∫S 1 [H ] [H ]dS
[K ] = ∑ [k ] = ∑ ∫S [H ] (erα ⋅ [∇H ])dS
For the collocation method, the weighting functions are selected to be
Dirac delta functions. In the present formulation, the Dirac delta functions
are defined at the nodal points of each element. Therefore, for a 2-D case,
w( x, y ) = δ ( x − xi )δ ( y − y j ) where xi and yj are the nodal coordinate
values. On the other hand, for the method of moment the weighting
functions are chosen to be monomial terms such as xpyq (where p, and q
are non-negative integers) starting from the lowest order.
Once the matrix equation of the first order derivative in time, as given
in Eq. (12), is developed from the weighted residual finite element
formulation, a numerical time integration scheme is applied to the
expression. There are many different numerical techniques for time
integration. Those include, but not limited to, the forward difference
technique, backward difference technique, Crank-Nicolson technique,
Runge-Kutta technique, and predictor-corrector technique. Because of
computational efficiency, the forward difference technique is selected in
this study. When we make the matrix [M] a diagonal matrix, the forward
difference technique becomes an explicit method. As a result, even if a
small time step size ∆t is used for the forward difference scheme because
of conditional stability, the overall computation is efficient. If an
unconditionally stable method is preferred, the Crank-Nicolson technique
may be applied.
Using the forward difference scheme for time integration, Eq. (12) is
expressed as
{Fα }t + ∆t = {Fα }t + ∆t [M ]−1 [C ]{F~α } − [C ]{Fα }t − [K ]{Fα }t
Figure 3. Mesh in converging duct
Figure 4. Comparison of fluid velocity at the center of the duct between
the LBM and CFD results
4. Conclusions
A lattice Boltzmann technique based on the weighted residual finite
element formulation was developed so that the technique could be applied
to a complex shape of domain, which is essential for fluid-structure
interaction problems in commercial nuclear power systems. Some
example problems were solved to demonstrate the developed technique.
Equation (21) is solved for the given initial and boundary conditions.
3. Numerical Results
The first example was a flow between two co-axial circular cylinders.
The ratio of the radii between the outer and inner cylinders was three.
The inner cylinder was kept at rest while the outer cylinder started to
rotate at a constant angular speed. The FELBM solution is plotted in Fig.
2 against the analytical solution. In the figure, the velocity was
normalized with respect to the outer cylinder velocity. In addition, the
distance was normalized such that the position of the inner cylinder was
zero while the position of the outer cylinder was unity. Both the Galerkin
method and the method of moment were used, respectively, with fournode quadrilateral finite elements. The numerical results compare well
with the exact solution as shown in the figure.
The second example was a flow in a linearly converging duct. Its
geometry and mesh is shown in Fig. 3. The pressure difference was
applied between the left inlet and right outlet. Transient flow analyses
were conducted using FELBM and the traditional CFD, respectively.
After 10,000 time steps, the two solutions were compared in Fig. 4, and
they agree well.
[1] Y. W. Kwon, “Development of coupling technique for LBM and FEM
for FSI application”, Engineering Computations, 23(8), 2006, 860-875.
[2] H. Chen, “Volumetric formulation of the lattice Boltzmann method
for fluid dynamics: Basic concept”, Phys. Rev. E., 1998, 3955.
[3] G. Peng, H. Xi, C. Duncan, and S.-H. Chou, “Lattice Boltzmann
method on irregular meshes”, Phys. Rev. E 58, 1998, R4124-R4127.
[4] X. He, L.-S. Luo, and M. Dembo, “Some progress in lattice
Boltzmann method. Part I. Nonuniform mesh grids”, Journal of
Computational Physics, 129, 1996, 357-363.
[5] T. Lee and C.-L. Lin, “A characteristic Galerkin method for discrete
Boltzmann equation”, Journal of Computational Physics, 171, 2001, 336356.
[6] Y. Li, E. J. LeBoeuf, and P. K. Basu, “Least-squares finite-element
lattice Boltzmann method”, Phys. Rev. E 69, 2004, 065701-1 – 065701-4.
[7] Y. W. Kwon and H. Bang, The Finite Element Method Using Matlab,
2nd ed., CRC Press, Boca Raton, Florida 2000.
Report inappropriate content