Project Topics
Get Complete Project Material File(s) Now! »
Numerical methodology
Computational Fluid Dynamics (CFD) is an active branch of the uid mechanics. It studies uid ows and their e ects by solving equations with numerical methods. It has become a powerful and crucial tool in all the researches related to the uid mechanics. It allows an access to the instantaneous physical details of the uid ow (such as veloc-ity, pressure and concentration) with a signi cantly low cost comparing to experiments, especially for cases where experiments become di cult, even impossible to realise.
According to the scale of the computing target, CFD numerical methods can be divided into three branches (Tao [2000]): mental model is the Navier{Stokes equations (Ferziger and Peric [2002]).
Meso-scale methods. It simulates the kinetics of particle-groups and obtains values of macro physical parameters by some averaging techniques. The Direct Simulation of Monte Carlo method (DSMC) (Alexander and Garcia [1997]) and the Lattice-Micro-scale methods. It simulates the dynamics of molecules, atoms, etc. and calculates the macro physical parameters by the statistic method. The Molecular Dynamics simulation (MD) (Haile [1997]) is included.
The principle of a CFD numerical method can be listed as blow:
Modelling of the physical problem and establishing the corresponding governing equations (such as a simpli ed or modi ed form of the Navier-Stokes equations and other necessary equations).
Establishment of the physical geometry of the problem. The optimized computa-tional area is de ned as a compromise between the real computational cost and the physical precision.
De nition of the initial and boundary conditions. A correct de nition of conditions is essential to a CFD problem. The conditions generally consist of informations of velocity, pressure, temperature or concentration.
Temporal and spatial discretisation of the computational area and the mesh gen-eration. The spatial discretisation consists of the Finite Volume Method (FVM) (Versteeg and Malalasekera [2007]), the Finite Element Method (FEM) (Zienkiewicz et al. [2005]), the Finite Di erence Method (FDM) (Smith [1985]), etc.
Iterative process to solve the governing equations and convergence criteria. For a system of linear equations, Jacobi method, Gauss{Seidel method, Successive over-relaxation method (Tao [2001]) and the Krylov subspace method (Watkins [2007]) can be employed. A hierarchical discretization method – the Multi-grid method (Hackbusch [2003]) is an e cient iterative method and a large number of algorithms are inspired by this idea.
Analysis and visualization of the simulation results to validate the CFD models and obtain the physical explanation from the mathematical numerical solutions.
Multiphase CFD problems are of great interest for the industry and the academy. Particularly, as a simple case of multiphase ows, the diphasic ow is widely studied. The diphasic ow can consist of a same uid in two di erent phases, two di erent uids in a same phase or two di erent uids in two di erent phases. The precise tracking of the interface between two di erent phases remains di cult, which makes the establishment of an e cient numerical scheme for the diphasic ow a great challenge and also a highly active research subject.
Front Tracking method: it uses a point-group called markers at the interface, where the interface reconstruction is done by connecting the markers with a curve (in two dimensions) or a surface element (in three dimensions) (Chern et al. [1986], Shyy et al. [1996], Tryggvason et al. [2001]).
Level Set method: it uses a level function of which the level 0 is the interface, it provides a high-precision tracking of the interface but the local mass conservation needs to be carefully guaranteed (Granier et al. [1996], Osher and Fedkiw [2001], Osher and Sethian [1988]).
Volume-of-Fluid method (VOF) : it reconstructs the interface from a volume-fraction occupied by one phase in a mesh cell, it has a good local mass conservation and the high-precision tracking of the interface can be achieved by a local mesh re nement (Gretar Tryggvason and Zaleski [2011], Guey er et al. [1999], Hirt and Nichols [1981]).
In this chapter, we are interested in di erent schemes such as VOF, CSF (Continuum Surface Force) and HF (Height-Function) used in a diphasic uid simulation code called GERRIS, which will be employed in the coming chapters. GERRIS is an Open Source, Free Software program (The GERRIS Flow Solver http://gfs.sf.net) and is generally a partial di erential equation solver for uid ows. It allows to solve the time-dependent in-compressible variable-density ows by Euler, Stokes, Navier-Stokes equations or linear and non-linear shallow-water equations. GERRIS is based on the nite-volume discretisation, uses an adaptive quadtree (octree in three dimensions) mesh generation technique and tracks the interface by a Volume-of-Fluid/Piecewise Linear Interface Calculation method (VOF/PLIC1) [Li, 1995] with an accurate surface tension calculation. It supports advect-ed/di used passive tracers to track physical variables. GERRIS supplies portable parallel computations via the MPI library and it has developed special modules for geographic (GfsRiver, GfsOcean) and electro-hydrodynamics (GfsElectroHydro) issues. GERRIS is written in C with an object-oriented style. It is created by Stephane Popinet and is sup-ported by National Institute of Water and Atmospheric research (NIWA) of New Zealand and Institute Jean Le Rond D’Alembert (DALEMBERT) of France.
Numerical schemes
The detail of the numerical schemes used in GERRIS is described in the articles of Popinet [2003, 2009]. Here, we just talk about the main idea of the numerical schemes.
Quadtree/Octree Adaptive Mesh Re nement (AMR)
GERRIS employs a hierarchical adaptive mesh re nement through a fully-threaded tree structure. The adaptation process has two strategies: The rst strategy is the re ning process. Only small cells which satisfy a user-de ned criterion will be re ned. Several criteria can be introduced such as the adaptation of the gradient or a de ned function for a variable. The second strategy is the coarsening process. After each re nement, all cells will be checked. If a cell does not satisfy the re nement criterion (it is over-re ned), the cell will be coarsened.
Then, values of variables in the cell center will be updated. For parent cells, values are calculated by the mass-volumetric average of the children cells in order to make a good local conservation of physical variables (velocity for instance). For children cells, values are calculated by a linear interpolation from their parent cells, which ensures a good local conversation of velocity. A projection step is implemented to avoid the numerical oscillation of vorticity and guarantee the incompressible assumption. Besides, a fractional time-step is introduced to reduce the computing cost.
Temporal discretisation
The diphasic, variable-density Navier-Stokes equations for an incompressible ow with the surface tension can be written as:where (x; t) is the uid density, u = (u; v; w) the uid velocity, p the pressure,
(x; t) the uid viscosity, D the deformation tensor de ned as Dij (@iuj + @jui)=2,
the surface tension, the mean curvature of the interface and n the normal vector to the interface. The distribution function s is a delta function which insures the surface tension is only on the interface. The force Fext represents the external forces applied on the system (for example, the gravitational force g).
A volume-fraction c(x; t) is introduced and the density and viscosity of the system are de ned as :
where 1, 2 and 1, 2 are the density and viscosity of the uid 1 and the uid 2 respectively.
where the right-hand-side term depends only on the values at the step n and n + 1=2, which is a Helmholtz-type equation and is solved using a variant of the multilevel Poisson solver. The Crank{Nicholson discretisation [Crank and Nicolson, 1947] of the viscous terms is implemented and it is formally second-order accurate and unconditionally stable. The velocity advection term un+ 12 run+ 12 is estimated by the Bell{Colella{Glaz second-order unsplit upwind scheme [Bell et al., 1989] and this scheme is stable for a CFL number 2 smaller than 1.
Spatial discretisation
The computational area is discretised in an orthogonal way using square (cubic in three dimensions) nite-volume organised hierarchically by a tree-structure called quadtree in two dimensions (octree in three dimensions).
Each nite volume is de ned as a cell. Each cell could be the parent of up to four children (eight in three dimensions). The cell level is de ned by starting from zero for the root cell and by adding one each time when a group of four children is added. Each cell has a direct neighbour at the same level in each direction (four in two dimensions, six in three dimensions). Each of these neighbours is accessible through a cell face. A mixed cell is also de ned which has two uid phases or is cut by a solid boundary. An example of spatial discretisation and the corresponding tree-structure representation is given in Figure 1.1.
To simplify the calculations (for instance, the gradient and the ux calculation) at the cell boundaries, several constraints are introduced in GERRIS:
The di erence of level between direct-neighbour cells can not be bigger than 1;
The di erence of level between diagonal-neighbour cells can not be bigger than 1;
All the cells directly neighbouring a mixed cell must be at the same level.
All the physical variables (for example velocity, pressure and passive tracers) are col-located in the cell center (square center in two dimensions and cube center in three di-2The CFL condition proposed by Courant et al. [1928, 1967] is a necessary but not su cient condition of stability for solving certain partial di erential equations. It can be described that the dependent physical domain must be contained in the numerical domain in a certain time.
Figure 1.1: Example of two-dimensional quadtree spatial discretisation (left) and the corresponding tree-structure representation (right) [Popinet, 2003].
mensions) and they are interpreted as volume-averaged values in the corresponding nite volume. A collocated variable arrangement makes momentum conservation simpler on dealing with adaptive mesh. It also simpli es the implementation of the Crank{Nicholson discretisation of the viscous terms.
An approximate projection method of spatial discretisation [Almgren et al., 2000] is implemented for the pressure correction equation (1.11) and the associated divergence in the Poisson equation (1.13) in order to avoid the classical decoupling between the pressure and velocity eld caused by a collocated mesh.
At the rst step, an auxiliary cell-centred velocity eld uc is computed by equation (1.14). Then, an auxiliary face-centred velocity eld uf is calculated by averaging the cell-centred values on all the faces of the orthogonal discretisation volumes. When the faces are on the boundary between di erent re nement levels of the quadtree (octree in three dimensions) mesh, an averaging process is performed to guarantee the consistency of the corresponding volume uxes.
where m is the local normal to the interface and x is the position vector.
As m and the local volume fraction c are known, is determined by ensuring that the volume of uid in the cell is equal to c. This volume can be calculated by considering the relative positioning of the line (the plane in three dimensions) in the cell which leads to the matched linear and quadratic (cubic in three dimensions) functions of . A Mixed-Youngs-Centred (MYC) scheme [Aulisa et al., 2007] is employed for interface normal estimation where a 3 3 (3 3 3 in three dimensions) regular Cartesian stencil is used.
Figure 1.2: Left: Cells cut by the interface and its associated volume fraction eld. Right:
Reconstruction with piecewise linear interface [Agbaglah, 2011].
of this volume occupied by the phase 1 is indicated by the dark grey triangle. The area of this triangle is an estimation of the volume ux i+ 12 ;j. For this advection, GERRIS uses a piecewise linear interface calculation scheme [Li, 1995], also known as least-square t and split Lagrangian-Eulerian advection scheme [Scardovelli and Zaleski, 2003]. This geometrical approach is e cient and simple to implement with a Cartesian (or orthogonal) discretisation elements using a quadtree (octree in three dimensions) spatial discretisation.
Figure 1.3: Geometrical ux estimation [Popinet, 2009].
A general form of this advection scheme with a quadtree (octree in three dimensions) discretisation is illustrated in Figure 1.4. In this case, cell C is a coarse cell, so the uxes are computed independently for the top- (Ct) and bottom-halves (Cb). Once this is done, ux calculation in each half-cell is identical to the standard case (Figure 1.3).
Balanced-force surface tension calculation
To get an accurate estimation of the surface tension term ( s n)n+ 12 in the dis-cretised momentum equation (1.9) is one of the most di cult part of a VOF method in surface-tension-driven ows. A Continuum-Surface-Force (CSF) approach [Brackbill r n~ with n~ jrc~j (1.20) where r is a nite di erence operator, c~ is a spatially- ltered volume fraction eld and j:j is a norm operator. Harvie et al. [2006], Lafaurie et al. [1994], Popinet and Zaleski [1999] reported that problematic parasitic currents are created by this approximation in the case of the theoretical equilibrium of a stationary droplet. To solve this problem, it is reported by Francois et al. [2006], Renardy and Renardy [2002] that in this case, the discretised momentum equation can be reduced to then the exact discrete equilibrium between surface tension and pressure gradient is regained by pn+ 1 = cn+ 1 + arbitrary constant (1.23).
Conclusion
In this chapter, we rst give a brief introduction of the numerical schemes used in the code GERRIS. GERRIS is in general a solver for incompressible multiphase ows, it pro-vides solutions for Navier-Stokes-type equations and it has developed special modules for geographic (GfsRiver, GfsOcean) and electro-hydrodynamics (GfsElectroHydro) issues. A staggered temporal discretization is used and it ensures a formally second-order accuracy. The velocity is discretized by a Hodge decomposition and the pressure is calculated by solving a multi-level Poisson equation. A collocated nite-volume spatial discretization is employed with an orthogonal tree-structure mesh arrangement. As one of the highlights in GERRIS, the adaptive quadtree (octree) mesh ensures an e cient computation with a high precision and a low computing cost. The interface advection is calculated by a piece-wise linear Volume-of-Fluid (VOF) geometrical scheme. A height-function scheme is used to calculate the curvature and a balanced-force scheme ensures an accurate computation of the surface tension. GERRIS also provides the parallel computation via MPI and has successful applications in clusters and supercomputers in practice. A visualisation tool – GfsView is developed to realise the on-the- y or o ine visualisation for serial or parallel computations.
In the second part, we introduce the nondimensionalization process in GERRIS with an example of the drop impact. Since GERRIS has a coherent unit system, all coherent unit-systems are acceptable which is a convenience for industrial computation.
Table of contents :
Introduction
1 Numerical methodology
1.1 Numerical schemes
1.1.1 Quadtree/Octree Adaptive Mesh Renement (AMR)
1.1.2 Temporal discretisation
1.1.3 Spatial discretisation
1.1.4 Volume-of-Fluid (VOF) advection scheme
1.1.5 Balanced-force surface tension calculation
1.1.6 Height-function curvature calculation
1.1.7 Parallelisation
1.2 Nondim en sionalization
1.3 Conclusion
2 Validation of the numerical method: meniscus case
2.1 Physical description
2.2 Results and discussions
2.2.1 Dynamics of the contact point
2.2.2 Equilibrium interface
2.3 Conclusion
3 Gas elect in drop impact on a solid substrate
3.1 Physical description
3.2 Convergence test: analysis of the gas layer
3.3 Early stage of the impact
3.3.1 Constant density ratio r = 0:003
3.3.2 Constant viscosity ratio m = 0:037 and m = 0:07
3.4 Frontier between impact outcomes
3.5 Conclusion
4 Aerodynamic splashing mechanism
4.1 Pre-contact dynamics
4.1.1 Lubrication approximation in the gas layer
4.1.2 Characteristic lengths in the gas layer
4.1.3 Ejection of the small jet
4.1.4 Lift-up of the small jet
4.2 Post-contact dynamics
4.2.1 Central bubble entrapment
4.2.2 Detachment of the thin liquid sheet
4.2.3 Gas entrainment around the contact line
4.3 An aerodynamic splashing mechanism
4.4 Other elects in the impact dynamics
4.4.1 Contact angle elect
4.4.1.1 Impact on a totally-nonwettable surface ( = 180)
4.4.1.2 Impact on a totally-wettable surface ( = 0)
4.4.1.3 Impact on a partially-wettable surface ( = 30; 60; 120; 150)
4.4.2 Gas compressibility elect
4.5 Conclusion
5 Drop impact on a highly-viscous liquid
5.1 Physical description
5.2 Numerical investigations
5.2.1 Wave-like period
5.2.2 Solidifcation period
5.2.3 Comparison with an impact on solid
5.3 Comparison with the experiments
