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 a 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 : *jcjo@kins.re.kr b 1. Introduction 6 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. 7 ) 2 cos{(α − 1)π / 2 + π / 4}, 2 sin{(α − 1)π / 2 + π / 4} τ (f α − fα ~ fα in which ρ is the fluid density and expressed as ρ = ∑ fα r v n (9) 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 (2) the nodes. Plugging Eq. (9) into Eq. (8) yields [H ]{f&α }+ erα ⋅ [∇H ]{ f α } + 1 [H ]({ f α } − {~fα }) = 0 (10) τ for each finite element. The superimposed dot denotes the temporal derivative. Applying the weighted residual formulation to Eq. (10) gives the following expression ~ ⎞ r 1 (11) ⎛ ∑ ∫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 { }) e and the summation is performed over the total number of elements. Furthermore, {w} is a column vector of the weighting functions. The size (4) of {w} is equal to the number of nodes per element. Rewriting Eq. (11) yields (12) [M ]{F&α }+ [K ]{Fα } + [C ]{Fα } − [C ]{F~α } = 0 is the fluid velocity. They can be (5) where α and r r ρ 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 α + − ⎥ c 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 (3) ~ 1 Ωα = − α ∂t direction, and Ωα denotes the collision operator. The discrete velocity vector is given below and shown in Fig. 1 for the D2Q9. ( (1,-1) In addition, ωα is the weighting parameter for each velocity direction as given below: α =0 ⎞ ⎛ 4/9 (7) ωα = ⎜⎜ 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 ~ 1 (8) + e ⋅ ∇f + ( f − f ) = 0 where fi is the particle velocity distribution function along the α− r direction, t represents time, eα is the discrete velocity vector along the α− α =0 α = 1 to 4 8 4 (0,-1) Figure 1. D2Q9 lattice showing discrete nine velocity vectors (1) (0, 0) c ( cos{(α − 1)π / 2},sin{(α − 1)π / 2}) 1 (1,0) 3 (-1,-1) The lattice Boltzmann equation is expressed as ∂f α r + eα ⋅ ∇f α = Ωα ∂t (1,1) 5 0 (-1,0) 2. Finite Element Based Lattice Boltzmann Method ⎛ r ⎜⎜ eα = ⎜ ⎜c ⎝ (0,1) 2 (-1,1) [M ] = ∑ [m] = ∑ ∫S {w}[H ]dS (13) [K ] = ∑ [k ] = ∑ ∫S {w}(erα ⋅ [∇H ])dS (14) e (6) α e 489 Transactions of the Korean Nuclear Society Autumn Meeting Pyeong Chang, Korea, October xx-xx, 2008 [C ] = ∑ [c] = ∑ ∫S 1 {w}[H ]dS (15) {Fα } = ∑ { f α } (16) 1 ~ ~ Fα = ∑ f α (17) 0.8 τ { } { } 1.2 Norm. Velocity e 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 T (18) [M ] = [m] = [H ] [H ]dS ∑ -0.2 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figure 2. Comparison of steady state velocity profile of flow between two co-axial cylinders using four-node quadrilateral element. (19) T 0 Norm. Dist. e (20) [C ] = ∑ [c ] = ∑ ∫S 1 [H ] [H ]dS e Exact Galerkin Moment 0 Se T 0.4 0.2 ∑∫ [K ] = ∑ [k ] = ∑ ∫S [H ] (erα ⋅ [∇H ])dS 0.6 τ 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 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. (21) Equation (21) is solved for the given initial and boundary conditions. 3. Numerical Results REFERENCES 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. 490

1/--pages