Hasan Thesis
Hasan Thesis
Hasan Thesis
School of Naval Architecture and Ocean Engineering University of Ulsan January, 2012
Abstract
Free surface has long been one of the most popular topics in the field of Computational Fluid Dynamics (CFD) because of its importance in many engineering applications. On the other hand, Lattice Boltzmann Method (LBM) has also evolved the reputation in CFD world due to its simplicity in coding and ease of implementation.
In this thesis, an algorithm about how to simulate two dimensional dam breaking with lattice Boltzmann method (LBM) is presented. LBM is a relatively new simulation technique for complex fluid systems and has attracted interest from researchers in computational physics. Unlike the traditional CFD methods, which solve the conservation equations of macroscopic properties (i.e., mass, momentum, and energy) numerically. LBM considers a typical volume element of fluid to be composed of a collection of particles that is represented by a particle velocity distribution function for each fluid component at each grid. We use the modified Lattice Boltzmann Method for incompressible fluid and surely there is a difference in the LBM between the compressible and incompressible fluid. We will represent detailed information on single phase flow which considers only the water instead of both air and water. Interface treatment and conservation of mass are the most important things in simulating free surface where the interface is treated by mass exchange within the water region. We consider free surface boundary condition on the interface and bounce back boundary condition for the treatment of solid obstacles. Finally, the results of the simulation will be compared with some other methods like Volume of Fluid, Level Set Method etc and some examples related with free surface will also be provided in the thesis.
Acknowledgements
Foremost, I would like to thank my Allah (God) for creating me in this world and for giving me such a wonderful life. Then I wish to express my sincere gratitude to my advisor Prof. Rho-Taek Jung, for the continuous support of my Masters study and research, for his patience, motivation, enthusiasm, and immense knowledge. His guidance helped me all the time of my research and writing of this thesis. I could not have imagined having a better advisor and mentor for my research.
Besides my advisor, I would like to thank the member of my thesis committee, Prof. Hyung Taek Ahn and Prof. Weoncheol Koo, for their encouragement and insightful comments.
My sincere thanks also go to all the professors of the school of NAOE for giving the guidance through their lectures in the class to carry out my research. I am really very much thankful to my lab mates specially, Md. Sayeemuzzaman, for the stimulating discussions, for the sleepless nights we were working together before deadlines, and for all the fun we have had during the last two years. I also thank my other colleagues in the department.
I would like to express the deepest appreciation to my parents Md. Mosharraf Hossaian and Mrs. Kamrun Nesa, to give me their endless love and my brother Md. Samiul Hasana and cousine Abu Lahel, to assist my in various occasions since my childhood. Finally, I would like to give warm thanks to my wife, Tanziha Tasnim, for her support, concern and always kind & friendly encouragements throughout my study.
ii
Contents
Abstract Acknowledgement Contents List of Figures List of Tables CHAPTER 1 Introduction CHAPTER 2 Numerical Methods in Fluid Study CHAPTER 3 Lattice Boltzmann Method 3.1 Lattice Gas Automata 3.2 Boltzmann Gas concept 3.3 LBM Framework 3.4 Advantage & Disadvantage of LBM CHAPTER 4 Free Surface Treatments 4.1 Governing Equations 4.2 Gravitational Effect in Free Surface 4.3 Free Surface 4.4 Interface Advections 4.5 Boundary Conditions 4.6 Handling of excess mass 4.7 Force on the boundary CHAPTER 5 Implementation and Result 5.1 Algorithm of the Process to Simulate Breaking Dam 5.2 Numerical Results CHAPTER 6 Concluding Remarks References i ii iii iv vi 1 3 6 6 9 11 21 23 23 23 24 25 28 32 33 35 35 38 51 53
iii
List of Figures
Figure 3.1: FHP unit velocity vectors Figure 3.2: Zero net momentum, head-on, 2- and 3- particle collisions Figure 3.3: D2Q9 lattice & velocity directions Figure 4.1: Mass exchange between two interface cells Figure 4.2: Reconstruction of distribution function coming from empty cells Figure 4.3: Modified bounceback Figure 5.1: Computational Domain Figure 5.2(a): time 0.00212s Figure 5.2(b): time 0.00637s Figure 5.2(c): time 0.0106s Figure 5.2(d): time 0.0148s Figure 5.2(e): time 0.0212s Figure 5.2(f): time 0.0255s Figure 5.2(g): time 0.0318s Figure 5.2(h): time 0.0425s Figure 5.2(i): time 0.0531s Figure 5.2(j): time 0.0637s Figure 5.2(k): time 0.0743s Figure 5.2(l): time 0.0828s Figure 5.3(a): time 0.00094s Figure 5.3(b): time 0.0094s Figure 5.3(c): time 0.0141s Figure 5.3(d): time 0.0188s Figure 5.3(e): time 0.0236s Figure 5.3(f): time 0.0283s Figure 5.3(g): time 0.033s Figure 5.3(h): time 0.0424s Figure 5.3(i): time 0.0576s Figure 5.3(j): time 0.0708s Figure 5.3(k): time 0.0755s Figure 5.3(l): time 0.085s Figure 5.4(a): time 0.00053s 8 9 12 27 29 30 38 39 39 39 39 39 39 40 40 40 40 40 40 42 42 42 42 42 42 43 43 43 43 43 43 44
iv
Figure 5.4(b): time 0.008s Figure 5.4(c): time 0.016s Figure 5.4(d): time 0.0217s Figure 5.4(e): time 0.024s Figure 5.4(f): time 0.027s Figure 5.4(g): time 0.0324s Figure 5.4(h): time 0.0398s Figure 5.4(i): time 0.0478s Figure 5.4(j): time 0.053s Figure 5.4(k): time 0.0637s Figure 5.4(l): time 0.0743s Figure 5.5(a): Horizontal position of free surface at various time steps. Figure 5.5(b): Conservation of mass at various time steps. Figure 5.5(c): Computational forces on right vertical wall at various time steps Figure 5.5(d): Physical forces on right vertical wall at various time steps
44 44 44 44 44 45 45 45 45 45 45 47 48 49 50
List of Tables
Table 4.1: velocity directions and weight functions of D2Q9 13
vi
CHAPTER 1 Introduction
One of the important applications of free surface is breaking dam as it is directly related to free surface flow. A free surface is the surface of a fluid that is subject to constant perpendicular normal stress and zero parallel shear stress, such as the boundary between two homogenous fluids, for example liquid water and the air in the Earth's atmosphere. This means free surface flow is defined as a multiphase flow. But in case of particular application such as dam breaking, the air part can be neglected and so the multiphase is converted to single phase. The reason for the term "free" comes from the huge difference in the densities of the gas and liquid. For example, water has a density of 1000 kgm-3 whereas air has a density of 1.2041 kgm-3. A low gas density means that its inertia can generally be ignored compared to that of the liquid. In this sense the liquid moves independently, or freely, with respect to the gas. The only influence of the gas is the pressure it exerts on the liquid surface. In other words, the gas-liquid surface is not constrained, but free. In the numerical study we use the Lattice Boltzmann Method for incompressible flow. The entire framework is constituted from a particle perspective where collisions, streaming and particle-particle interactions. LBM is quite handy to simulate single as well as multiphase flow. A rich variety of behaviors, including unsteady flows, phase separation, evaporation, condensation, cavitations, solute and heat transport, buoyancy and fluid structure interactions can easily be simulated. Besides, due to the ease of extensibility and very good computational efficiency, Lattice Boltzmann Method has achieved a promising alternative over the conventional Navier-Stokes equations. Free surface treatment, as one of the heart touching topics in computational fluid dynamics, has been carried out already by many fluid dynamists. Interactive free surface fluids [1], Stable free surface flows on adaptively coarsened grids [2], Implementation of virtual river [3], Interfacial wave modeling [4] are some of major works done with Lattice Boltzmann Method. In the thesis, a computational study is developed to simulate the breaking dam. In
1
order to simulate the free surface one has to be very careful in handling the interface as the proper treatment of interface gives good simulation result. We will provide mathematical presentation of interface movement by mass exchange between interface-interface and interface-fluid region. Interface advection is done with the help of the fluid fraction, which is updated by recording the inflow and outflow of mass via distribution functions. Maintaining the conservation of mass is one of the vital tasks in free surface simulation. We apply the second order bounce back boundary conditions for the treatment of solid and for the free surface boundary condition we give a reconstruction method in the interface between air and water.
modeled correctly, the system should behave as a fluid [3]. Different situations can be modeled by changing the average energy of the molecules and their separation.
The main disadvantage with such an approach is that large computer resources are required, many simulations taking hours to evolve a fraction of a second [4]. The system must be updated in small time steps, the new position and velocity of all particles being calculated, at every time step, from a knowledge of their previous position and velocity, taking into account any external forces which are acting on them. Any particles which collided during the previous time step have to be identified and their new trajectories calculated. This can be restrictively time consuming when considering even a very small volume of fluid. Even when a gas is being considered where there are fewer molecules and a larger time step can be used, because of the longer mean free path of the molecules, the number of molecules which can be considered is severely limited.
lattice sites and only at discrete times can run much faster on a computer than a molecular dynamics simulation. The lattice gas model has another big advantage over molecular simulation since all the collisions occur at the same time. This is a particular advantage if the simulation is being run on a parallel computer. These two time saving advantages of the lattice gas model allow simulations of a significantly large scale to be performed.
function. The distribution function is a statistical parameter from which the macroscopic properties of the fluid can be found. The simulation of the Boltzmann equation is performed on a regular lattice and the forms of the collision function is taken to be the BGK [6] collision operator which was first considered to represent collisions in the nondiscrete Boltzmann equation. The computation reflects the evolution from the lattice gas model. The model is updated in the same manner as the lattice gas model except that now, instead of considering individual particles to be travelling along the links, it is the distribution function which is evolved. The lattice Boltzmann model has the advantages associated with the lattice gas model but all the lattice gas difficulties have been overcome. Thus the lattice Boltzmann model is an ideal tool in fluid simulation.
Typically, the rule for updating the state of cells is the same for each cell and does not change over time, and is applied to the whole grid simultaneously, though exceptions are known. The main motivation for the transition from LGA to LBM was the desire to remove the statistical noise by replacing the Boolean particle number in a lattice direction with its ensemble average, the so-called density distribution function. Accompanying this replacement, the discrete collision rule is also replaced by a continuous function known as the collision operator. In the LBM development, an important simplification is to approximate the collision operator with the Bhatnagar-Gross-Krook (BGK) relaxation term [6]. This lattice BGK (LBGK) model makes simulations more efficient and allows flexibility of the transport coefficients. On the other hand, it has been shown that the LBM scheme can also be considered as a special discretized form of the continuous Boltzmann equation. From Chapman-Enskog theory, one can recover the governing continuity and Navier-Stokes equations from the LBM algorithm. In addition, the pressure field is also directly available from the density distributions and hence there is no extra Poisson equation to be solved as in traditional CFD methods.
indicating the presence (1) or absence (0) of particles moving from a lattice gas site at x to a neighboring site at x + ei.
The evolution of a lattice gas model proceeds in two steps that take place during each time step. The first step is a propagation or streaming step in which the particles move to new sites according the their previous positions and their velocities. Next, the particles collide and scatter according to collision rules.
Collision Rules:
There are several possible types of collisions on the hexagonal lattice. Only two types are considered in the simplest FHP model; two-body collisions involve two particles while three-body collisions involve three. Two critical features of the lattice gas that allow it to simulate the Navier-Stokes equations are mass conservation and momentum conservation. Thus, it is essential that the microscopic-scale collisions honor mass and momentum conservation. In the lattice gas, all particles have the same mass and speed so that the momentum conservation reduces to the conservation of the vector sum of the velocities. Head-on collisions between two particles (or three particles approaching one another from 120 separation) have no net momentum. Hence the results of these collisions must have zero net momentum.
The vectors shown in these figures represent velocity vectors attributable to particles at the centre of each hexagon just prior to and just after the collision step.
This is the streaming process. There are however collisions that result in some phase points starting at (x,p) not arriving at (x + p/mdt,p + Fdt) = (x + dx, p + dp) and some not starting at (x,p) arriving there too. We set (-) dxdpdt equal to the number of molecules that do not arrive in the expected portion of phase space due to collisions during time dt. Similarly, we set (+) dxdpdt equal to the number of molecules that start somewhere other than (x,p) and arrive in that portion of phase space due to collisions during time dt. If we start with the above equation and add tha changes in f(1) due to these collisions we obtain f(1)( x + dx, p + dp, t+dt) dx dp = f(1)(x,p,t) dxdp + [(+) - (-)] dxdpdt (3.2)
10
The first order terms of a taylor series expansion of L.H.S of Eq. (3.2), f(1)( x + dx, p + dp, t+dt) dx dp = f(1)(x,p,t) + dx x f(1) + dp p f(1) + gives the Boltzmann equation [ f(1 )(x,p,t) +dx x f(1) + dp p f(1) + + [(+) - (-)]dxdpdt From which we get the following: v x f(1) + F p f(1) +
) ) )
dt +.. ,
(3.3)
(3.4)
In its complete form with the collision operator written more explicitly, the Boltzmann equation is a nonlinear integral differential equation and is particularly complicated. According to [7], 50 years elapsed from the time that Boltzmann derived the equation before an approximate solution was found. With Lattice Boltzmann methods, we approximately solve the equation from the particle perspective and focus on an equation strongly akin to Eq. (3.2); it explicitly contains the collide and stream notion central to LBM.
The above figure shows the directional specification and the velocities ei where i = 0,1,..,8 is a direction index and e0 indicates particles at rest. This model is known as D2Q9 as it is two dimensional and contains 9 velocities.. The lattice unit (lu) is the fundamental measure of length in the Lattice Boltzmann models and time steps (ts) are the time unit. Besides this model there are some other models too. In case of two dimensional models there some models like D2Q5, D2Q6. D2Q5 is a two dimensional square lattice which has five velocity directional velocities and D2Q6 is a hexagonal lattice having six directions. On the other hand, there are 3-D models in Lattice Boltzmann study like D3Q15, D3Q27. There are 15 and 27 velocity directions in the case of D3Q15 and D3Q27 respectively. In our study, we concentrate on D2Q9 lattice. The velocity magnitude of e1 through e4 is 1 lattice unit per time step and the velocity magnitude of e5 through e8 is 2 lattice unit per time step. These velocities are exceptionally convenient in that all of their x- and ycomponents are either 0 or 1.
12
The discrete velocity directions are obtained from the following equations. ei = (0,0) for i = 0 ei = [ cos{ 0.5(i-1)} , sin{ 0.5(i-1)} ] for i = 1-4 ei = [ 2cos{ 0.25(2i-9)} , 2sin{ 0.25(2i-9)} ] for i = 5-8 These equations are used to find out the velocity directions for D2Q9 lattice. Actually there are different types of equations in order to get the velocity directions for each of the model. Below a table is given which shows the velocity directions for D2Q9 LB model MODEL Directions Discrete velocity Vectors ey 0 0 1 0 -1 1 1 -1 -1 ex 0 0 1 1 2 0 3 -1 D2Q9 4 0 5 1 6 -1 7 -1 8 1
D2Q9 4 0 5 1 6 -1 7 -1 8 1
-1
-1
-1
Table 4.1: velocity directions and weight functions of D2Q9 The next step is to incorporate the single-particle distribution function f, which is essentially the one that appears in equation (3.2), except that it has only nine discrete bins instead of being a continuous function. These distribution functions along nine directions (for D2Q9) give us the hydrodynamic variables. The macroscopic fluid density is the summation of all distribution functions of various directions. = (3.5)
13
The macroscopic velocity u is the average of the microscopic velocities ea weighted by the directional density functions. u= (3.6)
This simple equation allows us to pass from the discrete microscopic velocities that comprise the LBM back to a continuum of microscopic velocities representing the fluids motion.
In the following paragraphs, we outline the major steps of the brief derivation of the relevant equations and relations of the lattice Boltzmann method. Neglecting external forces, the transport equation for f(x, , t) can be expressed by the Boltzmann equation as
= (f,f)
(3.7)
The collision term (f,f) is quadratic in f, consisting of a complex integra-differential expression. A suitable simplification of the collision integral for the near equilibrium state of low Mach number hydrodynamics is the single relaxation time approximation, the socalled Bhatnagar Gross-Krook (BGK) model [6]. (f,f) = - ( f -f(0) ) where f
(0)
(3.8)
relaxation time which controls the rate of approaching equilibrium, or in other words the
14
viscosity of the fluid. The BGK relaxation still fulfills Boltzmanns H-theorem and locally conserves mass and momentum.
In order to solve for distribution function numerically, (3.7) is first discretized in the velocity space using a finite set of velocity vectors discrete Boltzmann equation, (a = 0, 1. . .,N) leading to the velocity
+ ei
=-
i = 0,1,.,N
(3.9)
Where N is the total number of velocity directions. It should be noted that equivalent to f (x, ei , t).
(x, t) is
To obtain the main equation of the lattice Boltzmann approach, (13.3) is discretized numerically in a very special manner. The discretization of space and time is accomplished by an explicit finite difference approximation. By scaling the lattice spacing, the time step and the discrete velocities appropriately, the discretized equations take the following explicit form: ) , ) = , ) , )
(3.10)
Where = /t is the dimensionless relaxation time and x is a point in the discretized physical space.
For simulating two-dimensional flows, the 9-velocity D2Q9 model (i = 0, . . . , 8), and for three-dimensional simulations, both, the 15-velocity D3Q15 (i = 0, . . . , 14) and the 19velocity D3Q19 model (i = 0, . . . , 18) are widely used. Surprisingly, the discretization with such a low number of collocation points is sufficient to describe the fluid in the nearequilibrium state of low Mach number hydrodynamics. Collisions of the fluid particles are considered as a relaxation towards a local equilibrium. In Lattice Boltzmann method, for all of these models there is a suitable equilibrium function.
15
= With c =
(3.11) is
the computational time step. In the Eq. (3.11) we can see the weighting factors which depend on the lattice model [8] and are given in the following table with three models (2D and 3D).
The above mentioned equilibrium function is usually used for compressible fluid flow but there is also another equilibrium function for incompressible fluid and surely there is difference between these two functions. As we know the main difference between incompressible fluid and compressible fluid lies in the density. Similarly, in order to distinguish the two functions in the case of lattice Boltzmann equilibrium function the difference is in density. We will describe briefly in the next section about it.
Now some explanation of the Eq. (3.10) will be given. The right hand side of (3.10) is usually called collision step and the left hand side is called streaming step. In streaming, we move the direction-specific densities to the nearest neighbor lattice nodes.
Collision is the last key element of a Lattice Boltzmann Model. Traverse the domain and implement of the governing equation at each node is the main task at collision process. For the collision step, the equilibrium distribution function has to be calculated at each cell and at each time step from the local density using (3.5) and the local macroscopic flow velocity u using (3.6). In particular when dealing with complicated formulations of the boundary conditions, it can be advantageous or even necessary to split the update process into the following two equations. , )= , ) , ) , )
(3.12)
16
)=
, )
(3.13)
Where and
, ) denotes the distribution values after collision (but before propagation), , ) are the values after collision and propagation, thus the values entering the
neighboring cell as data for the next time step. The Navier-Stokes equations (up to second order accuracy in space and time) can be derived formally from the lattice Boltzmann equation through the Chapman-Enskog expansion by a standard multi-scale expansion with time and space rescaled and the distribution function fi expanded up to second order [11, 12, 13].
The relation between the relaxation time and the kinematic shear viscosity , including a correction for the truncation error due to the discretization, can be obtained from the result of the Chapman-Enskog expansion. As the discretization error is known a priori from this analysis, it can be corrected and thus, the lattice Boltzmann method does not suffer from numerical diffusion as many other finite difference methods do. The final relation for the kinematic viscosity for D2Q9 is
(3.14)
(physical) viscosity. Numerical difficulties arise when the relaxation parameter approaches is the speed of sound of the fluid.
Now we can draw a flow chart on how to implement the lattice Boltzmann Method. A straightforward implementation would consist of three nested loops over the three spatial dimensions and treat the collision step separately from the propagation process. This algorithm would read the values of the current time step from the local cell, execute the relaxation and write the results back to a temporary array, as this can be done independently for all cells. In a separate nested loop which would only contain copy operations, these
17
values would then be propagated back to adjacent cells in the original array.
Propagation: Streaming of the distribution functions according to their direction to the next nodes
The lattice Boltzmann (LB) method has emerged as a new and effective numerical technique of computational fluid dynamics. Modeling of the incompressible Navier-Stokes equation is among many of its wide applications. Indeed, the lattice Boltzmann equation (LBE) was first proposed to simulate the incompressible Navier-Stokes equations. The incompressible Navier-Stokes equations can be derived from the lattice Boltzmann equation through the Chapman-Enskog procedure if the density fluctuation is assumed to be negligible. Unfortunately, this is not always the case in numerical simulations by using LBE method. The compressible effect in the existing LBE models may produce some serious errors in numerical simulations. There have been efforts to reduce or to eliminate the compressible effect in the LBE method (6-8). However, the results are not entirely
18
satisfactory mainly for the reason that the existing models for incompressible flow are only valid for steady flows in theory. Therefore, it is still necessary to improve the LBE method for simulations of the incompressible Navier-Stokes equations in general, especially for unsteady flows.
Actually, the incompressibility can be achieved only when the mass density becomes a constant. However, it is practically impossible to maintain a constant density in lattice Boltzmann models. Theoretically the lattice Boltzmann equation always simulates the compressible Navier-Stokes equation instead of the incompressible one, because the spatial density variation is not zero in LBE simulations. In order to correctly simulate the incompressible Navier-Stokes equation in practice, one must ensure that the Mach number, M, and the density variation, p, are of the order () and order (), respectively, where is the Knudsen number [16]. However, in numerical simulations such as flow through porous media, a pressure gradient is applied to drive the system, and the pressure gradient is established by maintaining a density gradient in the system. Moreover, this is the only way to implement the boundary condition of a pressure gradient in the system by using the method of the LBE or the lattice gas automata (LGA), because of the simple ideal gas equation of state of the system, and because pressure is not an independent dynamical variable in the LBE or LGA methods. Under this circumstance, the assumption of constant density is no longer valid and the magnitude of the density variation may be rather significant. This would inevitably bring in an error in the LBE simulations of the incompressible Navier-Stokes equations. It is well known that the lattice Boltzmann method is only applicable to the low Mach number hydrodynamics, because a small velocity expansion is (implicitly) used in derivation of the Navier-Stokes equation from the lattice Boltzmann equation. It should be noted that the small Mach number limit is equivalent to the incompressible limit. Thus, it should be possible to develop a lattice Boltzmann model to simulate the incompressible Navier-Stokes equation properly. In our thesis we used the equilibrium distribution function for incompressible fluid flow [16].
19
The equilibrium distribution function of incompressible fluid for Lattice Boltzmann method is defined as
(3.15)
If we compare the equilibrium function for incompressible fluid and compressible fluid then we can see the difference actually in the density which also satisfies the real difference between compressible and incompressible fluid. In the case of compressible fluid, the density is used to multiply the bracket which consists of velocity and lattice speed whereas in the case of incompressible fluid, the density is added with the term inside the bracket.
20
needs to designate a particular node as a solid obstacle and no special programming treatment is required. It should be noted that bounceback boundaries come in several variants and do not work perfectly. Nevertheless, with proper considerations of the effective boundary location and for =1, quite satisfactory results can be obtained.
Besides these two there are some other boundary conditions. Von Neumann (flux) boundary condition is one of them where a velocity vector is specified from which density/ pressure is computed on the basis of conditions inside the domain. Macroscopic density/ pressure is the only part needs to be computed. Moreover, there is another one named Dirichlet (pressure) boundary condition. It constrains the pressure/density at the boundaries.
[1] The Navier-Stokes equations are second-order partial differential equations (PDEs); the discrete velocity Boltzmann equation from which the lattice Boltzmann model is derived, consists of a set of first order PDEs.
[2] Navier-Stokes solvers inevitably need to treat the nonlinear convective term uu. The lattice Boltzmann method totally avoids the nonlinear convective term, because the convection becomes a simple advection (uniform data shift). The non-linearity of the Navier-Stokes equations is hidden in the quadratic velocity terms of the equilibrium distribution function.
[3] CFD solvers for the incompressible Navier-Stokes equations need to solve the Poisson equation for the pressure. This involves global data communication, while in the lattice Boltzmann method data communication is always local and the pressure is obtained through an equation of state.
21
[4] In the lattice Boltzmann method, the Courant-Friedrichs-Lewy (CFL) number is proportional to , in other words, the grid CFL number is equal = = 1. Consequently, the time
dependent lattice Boltzmann method is inefficient for solving steady-state problems, because its speed of convergence is dictated by acoustic propagation, which is very slow.
[5] Boundary conditions involving complicated geometries require a careful treatment in both Navier-Stokes and lattice Boltzmann solvers. In NavierStokes solvers, normal and shear stress components require appropriate handling of geometric estimates of normals and tangents, as well as onesided extrapolations. In lattice Boltzmann solvers, the boundary condition issue arises because the continuum framework, such as the no-slip condition at the wall, does not have a direct counterpart.
[6] Since the Boltzmann equation is kinetic-based, the physics associated with molecular level interactions can be incorporated more easily. Hence, the lattice Boltzmann model might be fruitfully applied to micro-scale fluid flow problems.
[7] The spatial discretization in the lattice Boltzmann method is dictated by the discretization of the particle velocity space. This coupling between discretized velocity space and configuration space leads to regular square grids. This is a limitation of the lattice Boltzmann methods, especially for aerodynamic applications where both the far field boundary condition and the near wall boundary layer need to be carefully implemented. Local grid refinement techniques cannot completely solve this issue. However, the algorithm of the LBM is similar to the structure of Jacobis method for the iterative solution of linear systems on structured meshes. In particular, the time loop in the LBM corresponds to the iteration loop in Jacobis method.
22
(4.1)
(4.2)
The density and the velocity can be calculated according to the Eq. (3.5) and Eq. (3.6) respectively. The process of getting the weight functions and discrete velocity directions have already been described in Table (3.1).
So the governing equation of the Lattice Boltzmann method is very simple. Moreover, the implementation of LBM is also very easy as we will discuss about it in the next chapter 5.
density variation increases and the constant density approximation becomes less valid. So there are some other methods also considered for application when there is a non-negligible density variation. Among these methods only one method will be considered for the free surface treatment that is considered in this thesis.
(4.3)
Where
as in Eq. (3.6) and is relaxation parameter. Here the body forces F = mg and g is the gravitational acceleration. So the equilibrium function is now not a function of & u. Instead, it is a function of & , where contains the gravity.
24
interface cells, which also store a fluid fraction value, similar to volume-of-fluid methods for conventional NS solvers [20]. Furthermore cell type conversions need to be handled if interface cells become completely filled or empty. The boundary conditions presented here do not compute the gas phase as a separate fluid, but assume a viscosity difference between gas and fluid phase that is high enough to approximate the gas velocity near the interface with the fluid velocity. This is especially suitable for simulations with large gas regions, since these do not require any computations. Empty cells that contain no fluid need not be considered in the algorithm until they are eventually converted to interface cells as described below. In our simulation we will also face the solid obstacles and the solid boundary is treated by modified bounceback method.
, )=
Where
(4.4)
and v indicates the mass, density and volume of each cell at position x and at time t.
During the streaming the interface advection is done by the mass exchange between interface cells and fluid cells. It is done through two ways: mass exchange between interface-interface and another is mass exchange between interface-fluid. As it has been already mentioned that Gas cells are separated from liquid cells by a layer of interface cells
25
the interface cells form a completely closed boundary, so it is assumed that no distribution function is directly advected from gas to fluid cells and also to interface cells. In this case we will need a reconstruction step that will be described later. This is a crucial point to assure mass conservation since mass coming from the liquid or mass transferred to the liquid always passes through the interface cells where the total mass is balanced. Hence, global conservation laws are fulfilled if mass and momentum conservation is ensured for interface cells. Per definition, the volume fraction of fluid and gas cells is 1 and 0,
respectively. The fluid mass content of a cell is denoted with m = m(x, t). The mass content is a function of the volume fraction and the density. For a gas cell the fluid mass content m is zero whereas that of a fluid cell is given by its density that means density and mass are equivalent to each other. Fluid cell gains and losses mass due to distribution functions. I f we consider the interface cell, density and mass are not equivalent and we have to account for the partially filled state by introducing a second parameter, as it has already been defined in Eq.(4.4) the volume fraction = (x, t). The fluid mass content m, the volume
fraction and the density are related by m(x, t) = (x, t) (x, t) for x interface cells.
All cells are able to change their state. It is important to notice that direct state changes from fluid to gas and vice versa are not possible. Hence, fluid and gas cells are only allowed to transform into interface cells whereas interface cells can be transformed into both gas and fluid cells. A fluid cell is transformed into an interface cell if a direct neighbor is transformed into a gas cell. At the moment of transformation the fluid cell contains a certain amount of fluid mass m which is stored. During further development the interface cell may gain mass from or lose mass to the neighboring cells. These mass currents are calculated and lead to a temporal change of m. If m drops below zero, the interface cell is transformed into a gas cell and this procedure will be described in the next section.
As it is mentioned that the free surface movement is done by the mass exchange [19], this mass exchange is carried out by only two factors; fluid fraction and distribution function. At each time step the mass exchange is calculated and is added to the old mass in order to get the new mass. This new mass will be used as old mass in the next time step. Let us
26
assume that the position of a cell is x and the position of another cell is x + ei. Then the exchanged mass is calculated by the following equation. , )= , ) , )
(4.5)
In the above equation = -I, this means the two distribution functions are opposite to each other. So the mass exchange is calculated from the difference of the two opposite distribution functions. The first term inside the bracket indicates the distribution function coming from the opposite direction of the distribution function which is indicated by the second term. This is quite understandable from the position of the two cells; one is at the other is at . Figure (4.1) also illustrates the phenomena very clearly. and
The Eq. (4.5) is used to calculate the mass exchange between a fluid cell and an interface cell. But in the case of mass exchange between an interface and interface, an additional term is multiplied with this equation. The interchange between an interface and a fluid cell should be the same as that of two fluid cells since the cell boundary is completely covered with liquid. In this case, the mass exchange can be directly calculated from the particle distribution functions. The interchange between two interface cells is approximated by assuming that the mass current is weighted by the mean occupied volume fraction. This equation looks like the following equation:
27
)=
, )
, )
, )
, )
(4.6)
In the Eq. (4.6) we can see an extra term related to fluid fraction. As the mass exchange between the two opposite interface cells is needed to calculate so the average fluid fraction of the two cells are considered. Both equations are completely symmetric, as of course the amount of fluid leaving one cell, has to enter the other one and vice versa, which means )= e t). This also corresponds the conservation of mass. It should be
noted that there is no mass transfer between interface cells and gas cells. Also no mass transfer is computed between the fluid cells.
After getting calculating the mass exchange for all the interface cells, new mass is calculated by adding the exchanged mass with the old mass. , , )= , ) , )
(4.7)
So
) is new mass at the new time step. We need to remember that these mass
reconstruction procedure in order to find out this distribution function. From Figure (4.2) we can get a quite clear understanding.
In Fig. (4.2) a dotted line is seen. This is distribution function that we do not have in our computations. So this distribution function is reconstructed by the following way [19]. ,t t) = , ) , ,
(4.8)
Equilibrium distribution function is a function of density and velocity. In the Eq. (4.8) we can see there is a difference in density. In the previous cases we have used the density. But here in the reconstruction term, has been used. Actually, for the
is the density of
air. We know that . So in the interface for the distribution functions that are coming from empty cells is used to satisfy free surface physical boundary condition. It should be
noted that in the Lattice Boltzmann method the pressure is derives directly from the density (distribution function).
In the Eq. ((4.8) the velocity of the interface is used. This also satisfies the free surface kinematic boundary condition as according to this boundary condition the air moves the same way as the interface does. So air is moving with velocity similar to the velocity of the interface.
29
The lattice Boltzmann method has traditionally been used with the bounce back boundary condition for the particles at the slip interface boundary. Consider, for example, the particles at the bottom boundary. From here and on, boundary conditions are considered for the representative case of the bottom boundary. However, arguments given below can be also applied for slip boundaries other than the bottom one.
Let us denote a unit cell at the bottom boundary in Fig. (4.3). Let us first assume that the fluid-solid interface boundary is at the line A. Then, the lattice sites 6, 2, 5 above the boundary belong to the fluid and the lattice sites 7, 4, 8 below the boundary belong to the solid.
The lattice sites 3, 0, 1 are right on the boundary. At the time t, the particles are at the node 0 of which the position is x. The particles at (x, t) moving in the directions of e4, e8, e7 into the solid region reflect back into the position x. turning to the opposite directions e2, e5, e6, respectively, at the next time t + 1.
30
This can be explicitly written as, for the particle distributions f2, f5, f6. , , , 1) = 1) = 1) = , ) , ) , ) , ) , ) , ) , ) , ) , )
The bounce back method that used Eqs. (4.9) (4.11) has traditionally been used in conjunction with the boundary at the regular lattice sites, as denoted by the line A in Fig. (4.3). This bounce back method is known as being of the first order in error, implying that the numerical error of the LBGK solution from the exact solution reduces linearly with the mesh size. We refer to this boundary method as the traditional bounce back or TBB method.
A variant of the TBB method is the one obtained by locating the boundary off the regular lattice sites, halfway between the nodes. Let us denote this boundary by the line B in Fig (4.3). In this boundary condition, the lattice sites 3, 0, 1 as well as the lattice sites 6, 2, 5 belong to the fluid and the lattice sites 7, 4, 8 belong to the solid. The bounce back method that uses Eqs. (4.9) - (4.11) has been frequently referred to as the shifted bounce back or SBB method, since the boundary is shifted from the regular lattice sites. The SBB method is known to have better error characteristics than the TBB method. Indeed, the SBB method is found to be of the second order in error while the TBB method is of the first order. The specular reflection exactly at the regular lattice sites, as denoted by the line A in Fig. (4.3), can be realized by letting it take place for a unit time interval from t to t + 1, instead of the instant reflection at t + . A particle, for example, moving in the direction of e4 finds that the destination is not allowed and then spends a unit time to turn to the opposite direction into the fluid, keeping its speed unchanged. For the paticles moving in the directions of e4, e7, e8 into the solid at (x. t) at the node 0 in the bottom boundary, the particle distributions at the new time t +1can be simply written as following: , , , 1) = 1) = 1) = , ) , ) , )
31
This is the standard or alternative bounceback slip boundary conditions (ABB). With the ABB method of the delayed reflection (4.12) - (4014) at hand, the bounce back method can be generalized. This generalized bounce back method can be used for boundaries that lie at the arbitrary locations between the neighboring regular lattice sites. This generalization is done by linearly combining the above two different boundary methods, namely, the SBB and ABB methods. For the bottom boundary, if an interpolation parameter to is used for the particles moving in the directions of e4, e7, e8, the particle distributions , 1), , 1) are then computed by , 1) ,
, , ,
1) = 1 1) = 1 1) = 1
) ) )
, ) , ) , )
, ) , ) , )
, ) , ) , )
, ) , ) , )
This bounceback method is called the modified bounceback method and in the thesis we have used this method with the value of = 1.
When we get the new fluid fraction for each of the interface cells, sometimes it is found that the fluid fraction is greater than one or less than zero for the interface cells. This means the interface has some positive or negative excess mass that are needed to be distributed to
32
the surrounding cells in order to keep the conservation of mass. If filled cell and if calculated as = )
< 0 it is called emptied cells. For the filled cells the excess mass is
(4.18)
As we are having some excess mass in the interface cell, so this excess mass needs to be distributed to the surrounding cells. Also we have to keep the fact in mind that the layer of the interface cells has to be closed. The excess masses for the filled cells are distributed to the surrounding interface and empty cells. In order to do that first we have to make the empty cell into interface and then the excess mass is distributed evenly to the surrounding interface and empty cells. To do that the distribution function of the new interface cell, which was an empty neighbor of the filled cell, is initialized by the equilibrium distribution function which is a function of average density and average velocity of the surrounding fluid and interface cells and the flag of the filled cell is re-initialized to fluid cell. In case of emptied cells, the excess mass is distributed evenly among the interface and fluid cells and the flag of the emptied cell is re-initialized to empty cell.
33
Where
and are opposite to each other. The fluid fraction, volume and velocity are . So the force is found out by the above equation which
represented by , and
illustrates the momentum exchange. After getting the force in all directions, the total force is calculated by summing up all the directional forces. = (4.21)
34
= =
= 1000
Computational density of water, = 1 Fluid fraction for interface cell initially, = . and fluid cell, = 1.
Computational Gravitational acceleration = Computational Kinematic viscosity, = CFL Condition, From the equation (5.2) by using the value of Physical time step physical time step (5.1) (5.2) , we can easily get the
. Then using this physical time step in the equation (5.1) along with
the value of physical Kinematic viscosity we get the computational Kinematic viscosity, .
35
From this equation we can easily calculate the relaxation parameter, sound speed of the fluid and the physical time step.
Step [2]: Initialization of equilibrium function: Equilibrium function is a function of density and velocity. From the equation (4.3) the gravity is incorporate in the velocity term. At the point (i, j) and in the direction n the equilibrium function is written as follows , , )= , ) , ) , ) , ) , )
Step [3]: Initialization of distribution function: The distribution function is initialized by the equilibrium distribution function. At point (i, j) and in the direction n the distribution function is written as follows , , )= Step [4]: Streaming & collision are carried out: After defining the equilibrium distribution function and initial distribution function along with the relaxation parameter, the streaming and collision are carried out that means the governing equations is applied. At point (i,j) and in the direction of n it looks like the following: , , )= Here the step t. Step [10]: Applying the boundary conditions: The free surface boundary conditions are applied in the interface. To treat the solid obstacles we use the modified bounceback method in our implementation. , , ) , , ) , , ) , , ) is at time , , )
Step [6]: calculation of density and velocity: Having all the distribution function at the current time step the new density and velocity are calculated according to the Eq.(3.5) and (3.6) respectively.
36
Mass exchange is performed according to equation (4.5) & (4.6).So mass exchange between two interface cell is given by the equation below: , , 6) = 1, 1, ) , , 6) 1, 1) 2 , )
Where (i,j) are the grid points. After calculating the exchanged mass, this mass is added to the old mass to get the new mass. , )= , )
Step [8]: Calculation of fluid fraction: After having the new density and new mass of an interface cell, the fluid fraction is calculated according to the formula given in the Eq. (4.4). Again it is being mentioned that there will be no computation of the fluid fraction for the fluid cell and empty or gas cell. Fluid fraction of fluid cell is assumed always equal to 1 whereas for the case of gas cell it is equal to 0. Step [9]: Distribution of excess mass: When we will have the new fluid fraction for the interface cell, we may have filled cell and emptied cell according to the section (4.6). If we have either positive or negative excess mass, then these masses should be distributed to the surrounding cells according to the section (4.6).If we do not have any excess mass, then there will be nothing to be distributed. Step [10]: Re-calculation of fluid fraction: We need this step if we have step 8. That means if there is filled cells or emptied cells, then we need this step. The neighbors of filled cell have to be converted into interface cells. So the fluid fraction of new interface cell is calculated which was a gas cell in previous time step. Similarly, the neighbors of emptied cell have to be converted into interface cells. So the fluid fraction of new interface cell is calculated which was a fluid cell in previous time step.
37
50
40
30
j
20 10
10
20
30
40
50
38
100 90 80 70 60
j
50 40 30 20 10 50 100
50 40 30 20 10 50 100
50 40 30 20 10 50 100
50 40 30 20 10 50 100
j
j
50 40 30 20 10 50 100
50 40 30 20 10 50 100
39
100 90 80 70 60
100 90 80 70 60
50 40 30 20 10 50 100
50 40 30 20 10 50 100
50 40 30 20 10 50 100
50 40 30 20 10 50 100
100 90 80 70 60
100 90 80 70 60
50 40 30 20 10 50 100
50 40 30 20 10 50 100
40
In the above figures we can see the advection of interface as time goes by. Results are shown for various computational time steps. We can also observe the velocity vectors and their directions. Due to the gravitational effect the water keeps going downward and shown from the fig: 5.2(a) to fig: 5.2(c). After some time it hits the right hand side solid boundary. We can see the water moves up after hitting the wall from the fig: 5.2(d) to fig: 5.2(f) which is quite physical. Then the water comes down because of the gravity and goes to the left side which is observed from fig: 5.2(g) to fig: 5.2(j). Then the simulation results of fig: 5.2(k) to fig: 5.2(l) shows us that the fluid is coming from the left side and going to the right side.
41
42
50
100
150
50
100
150
50
100
150
43
50
100
150
200
50
100
150
200
175
150
150 125 100 75 50 25
125 100 75 50 25
50
100
150
200
50
100
150
200
50
100
150
200
50
100
150
200
50
100
150
200
50
100
150
200
50
100
150
200
50
100
150
200
50
100
150
200
50
100
150
200
45
The movement of interface has been shown by the figures given above. The illustration of the figures is same as for the domain size 101x101. Due to the gravity the water comes down and moves in the same way as the fluid moves naturally. Fig 5.3(a)-5.3(l) gives the result for the 151x151 domain size whereas the results for the 201x201 are given by the figures from 5.4(a)-5.4(l).
46
, with n = 1
Fig 5.4(a) shows us the position of free surface at the horizontal solid boundary. As time goes by the value of the position increases that means the free surface goes forward. The results of different grid size have been provided here. We can see the numerical results show good characteristics compared with the computational results.
Z
1.4 1.3 1.2 1.1 1 0
0.5
47
Fluid Fraction
0.5 0
-0.5 -1
48
1.5E-05
1E-05
F
5E-06 0 1000
2000
3000
T
Figure 5.5(c): Computational forces on right vertical wall at various time steps.
49
2.2E-06
2E-06
1.8E-06
1.6E-06
1.4E-06
1.2E-06
F
1E-06 8E-07 6E-07 4E-07 2E-07 0 0.03 0.04 0.05
0.06
0.07
T
Figure 5.5(d): Physical forces on right vertical wall at various time steps.
50
B. The LBM actually works for the compressible fluid flow. But as we are considering the water, which is an incompressible fluid, we need some change in the compressible equation of LBM. So in the case of incompressible flow we have made a change in the density of the equilibrium function of the Lattice Boltzmann equation.
C. In the case of boundary condition, we have used two kinds of conditions. One is free surface boundary condition another is the solid boundary condition. For the free surface condition we made a special treatment of the distribution functions which are coming from the gas cells. And for the solid boundary conditions we used the modified slip bounceback boundary condition, not the standard bounceback condition.
D. The viscous effect of the fluid is shown in the results. We have presented results for three values of relaxation parameter which shows that how the viscosity works in the fluid. The results also show us the velocity contour and also directions that shows the physical significance of the free surface.
51
The lattice Boltzmann method (LBM) has evolved to a promising alternative to the well established methods based on finite elements/volumes for computational fluid dynamics simulations. Ease of implementation, extensibility, and computational efficiency are the major reasons for LBMs growing field of application and increasing popularity. Since established finite volume/difference/element methods are the result of an evolution over many decades, one might expect that the simple LBM cannot compete. However, the works done by various authors so far and extensive comparison of Lattice Boltzmann method with the other methods demonstrate that LBM is competitive although there is still room for improvements. From the current point of view, LBM based solvers are in particular suited for problems involving complex geometries, complex physics or weakly compressible transient flows with Mach numbers up to about 0.2.
52
References
[1] P. Roach. Computational Fluid Dynamics. Hermosa Publishing, 1972. [2] J. Connor and C. Brebbia. Finite Element Techniques for Fluid Flow. NewnwsButterworths, 1976. [3] J. Duppy and A. Dianoux. Microscopic Structure and Dynamics of Liquids. Plenum Press, 1978. [4] J. Goodfellow. Molecular Dynamics. Macmillan Press, 1991. [5] U. Frisch, B. Hasslacher, and Y. Pomeau. Lattice-gas automata for the Navier-Stokes equation. Physical Revies Letters, 56:1505-1508G, 1986. [6] P.L Bhatnagar, E.P. Gross, and M Krook. A model for collision process in gases. I. Small Amplitude Processes in Charged and Neutral One-Component Systems. Phys. Rev., 94(3):511-525, 1954. [7] Physical Revies, 94(3): 511-525, 1954. [8] An introduction to the theory of the Boltzmann equation. Holt, Reinhart and Winston, Inc., New York, 1971. [9] Qian, YH, dHumieres D, Lallemand P (1992) Lattice BGK models for NavierStokes equation. Europhys Lett 17:479-484 [10] S. Wolfram. Cellular Automaton Fluids 1: Basic Theory. J. Stat. Phys., 3/4:471-526, 1986 [11] D.H Rothman and S. Zaleski. Lattice Gas Cellular Automata. Simple models of complex Hydrodynamics. Cambridge University Press, 1997. [12] S. Chapman and T.G. Cowling. The Mathematical Theory of Non-Uniform Gases. Cambridge University Press, 1995. [13] U. Frisch, D. dHumieres, B. Hasslacher, P. Lallemand, Y. Pomeau, and J.P. Rivet. Lattice Gas Hydrodynamics in Two and Three Dimensions. Comples Systems. 1:649-707, 1987. [14] X.He and L.-S. Luo. A Priori Derivation of the Lattice Boltzmann Equation. Phys. Rev. E, 55(6):R6333-R6336, 1997. [15] Carlin Krner, Thomas Pohl, Ulrich Rde, Nils Threy, and Thomas Zieser. Parallel Lattice Boltzmann Methods for CFD Applications.
53
[16] Xiaoyi He and Li-Shi Luo. Lattice Boltzmann Model for the Incompressible NavierStokes Equation. Journal of Statistical Physics, Vol. 88, Nos. 3/4, 1997. [17] Li-Shi Luo, PhD thesis, Georgia Institute of Technology, 1993. [18] Yu, R. Mei, L.-S Luo, and W. Shyy. Visous Flow Computaions with the Method of Lattice Boltzmann Equation. Progr. Sci., 39:329-367, 2003. [19] J. M. Buick and C. A. Greated. Gravity in a Lattice Boltzmann model. Physical review E, Volume 61, Number 5, May 2000. [20] Nils Threy and Ulrich Rde. Stable free suface flows with the lattice Boltzmann method on adaptively coarsened grids. Computing and Visualization in Science manuscript No. DOI: 10.1007/s00791-008-0090-4. [21] Hirt, C.W., Nichols, B.D.: Volume of Fluid (VOF) Method for the Dynamics of Free Boundaries. J. Comp. Phys. 39, 201-225(1981). [22] Michael C. Sukop, Daniel T. Thorne, Jr. Lattice Boltzmann Modeling, An Introduction for Geoscientists and Engineers. [23] Huabing Li, Xiaoyan Lu, Haiping Fang, Yuehong Qian. Force evaluations in Lattice Boltzmann simulations with moving boundaries in two dimensions. Physical Revies E 70, 026701 (2004).
54
55