Numerical Simulation of Interfaces and Free Surfaces

From Thermal-FluidsPedia

**Related Topics Catalog**
Computational methodologies for forced convection One-Dimensional Steady-State Convection and Diffusion Central Difference Scheme Upwind Scheme Hybrid Scheme Exponential and Power Law Schemes A Generalized Expression of Discretization Schemes Multidimensional Convection and Diffusion Problems Numerical Solution of Flow Field Special Difficulties Staggered grid Pressure Correction Equation The SIMPLE Algorithm Numerical Simulation of Interfaces and Free Surfaces Application of Computational Methods

Many engineering applications involve interfacial phenomena because of the high heat transfer rates that can be achieved. There are many complexities that need to be addressed when modeling an interface, since an interface is generally irregular, involves mass transfer, is three-dimensional, and is not at a fixed location. There are two different approaches when considering an interface: a continuum and a noncontinuum approach. In many problems, the scale of the entire computational domain is much larger than the scale of the interfacial thickness; therefore the interface can be considered as a planar surface and the continuum approach will give highly accurate results. With the development of nanotechnology, the scales of systems are becoming very small. Therefore, neglecting the thickness of the interface is no longer a valid assumption, and a noncontinuum approach is needed. This approach can be achieved through molecular dynamic simulations. A detailed discussion on tracking an interface using the various approaches is given by Faghri and Zhang ^[1]^[2].

Interface Tracking Techniques

There are many engineering problems that involve separated phases in which the interface between phases is clearly defined, but is not at a fixed location. The interface can be assumed to be a planar surface in these situations, and therefore, a continuum approach is valid. In these problems, it is necessary to solve for both the phases as well as the interfacial location. These problems can be solved numerically on an Eulerian or Lagrangian mesh. An Eulerian mesh is stationary and defined prior to the start of a solution. When using an Eulerian mesh, the interface is tracked by solving an additional scalar equation. In the Lagrangian approach, a boundary of the mesh is aligned with the interface, and this boundary moves with the interface.

Interfacial representations for different interface tracking techniques.

When thinking of a multiphase system from a continuum approach, in the bulk region a phase is continuous and is discontinuous at an interface between different phases. In general, the interface is free to deform based on the nature of the flow. Therefore, it is difficult to efficiently capture an interface between phases with just one model. Consequently, there have been strong efforts to track an interface based on several different techniques, each with its own pros and cons. An actual interface, represented by the different numerical techniques, is presented in the figure to the right. The techniques are as follows:
1. Stationary Grid Approach: standard CFD modeling in cells that contain only a single phase; special consideration taken in cells in the vicinity of an interface.
2. Lagrangian Techniques: grid and fluid move together, interface is directly captured.
3. Phase Interface Fitted Grid: interface directly tracked as boundary, rest of grid moves as a function of interfacial movement, “semi-Lagrangian”.
4. Front Tracking Methods: stationary grid is used and modified near the interface, so that the grid is aligned with the interface; combination of 2 and 3.
The above methods are discussed in detail in Faghri and Zhang ^[1]. The stationary grid approach, which is used widely for solving single phase problems, is presented below.

Stationary Grid Approach

The stationary grid approach is based on a stationary grid where the fluid interface is captured directly. The first of such approaches is the marker-and-cell approach originated by Harlow and Welch ^[3]. In this approach, massless particles are introduced into the flow field, and the locations are projected from their interpolated velocities. Cells with a particle are considered to have one phase in them, and cells without a particle do not have that phase in them. An interface is considered where cells with particles are neighbored with cells without particles. For efficiency, this method was extended to only track particles on the surface ^[4]. Further development of this class of models led to the Volume of Fluid method (VOF) by Hirt and Nichols ^[5] in which they used a donor-acceptor method to effectively eliminate numerical diffusion at an interface.

The VOF method is one of the most widely used routines to solve a free surface problem on an Eulerian mesh. In this method, there is one velocity, pressure, and temperature field, which is shared by all of the phases modeled. The interface between the phases is tracked by the volume fraction of phase k, $\varepsilon _{k}$ . The volume fraction equation is the continuity equation of phase k divided by the density of that phase.

$\frac{1}{\rho _{k}}\left[ \frac{\partial }{\partial t}\left( \varepsilon _{k}\rho _{k} \right)+\nabla \cdot \left( \varepsilon _{k}\rho _{k}\mathbf{V} \right)=\sum\limits_{j=1\left( j\ne k \right)}^{\Pi }{{\dot{m}}'''_{jk}} \right]$	(1)

where $Π$ is the number of phases. When the volume fraction is between 0 and 1 in a computational cell, the cell is considered to be an interfacial cell. For a two-phase system, the phases are k and j. When the volume fraction is 1, that cell is occupied by only phase k, and when it is 0, that cell is occupied by only phase j. The sum of the volume fractions is unity,

$\sum\limits_{k=1}^{\Pi }{\varepsilon _{k}}=1$	(2)

Therefore, ( $Π - 1$ ) volume fraction equations need to be solved. The fluid properties, such as density, viscosity, and thermal conductivity, are calculated by their volume weighted average.

$\Phi _{eff}=\sum\limits_{k=1}^{\Pi }{\varepsilon _{k}\Phi _{k}}$	(3)

The overall continuity equation in conjunction with the VOF method is:

$\frac{\partial }{\partial t}\left( \rho _{eff} \right)+\nabla \cdot \left( \rho _{eff}\mathbf{V} \right)=0$	(4)

The continuity equation uses the effective density and the velocity field is shared by both phases; therefore, subscript k is dropped. The momentum equation is

$\frac{\partial }{\partial t}\left( \rho _{eff}\mathbf{V} \right)+\nabla \cdot \left( \rho _{eff}\mathbf{VV} \right)=\nabla \cdot {\tau }'_{eff}+\rho _{eff}\mathbf{X}+\sum\limits_{k=1}^{\Pi }{\sum\limits_{j=1\left( j\ne k \right)}^{\Pi }{\left\langle {\dot{m}}'''\mathbf{V}_{jk} \right\rangle }}$	(5)

In the energy equation, the enthalpy is mass averaged instead of volume averaged.

$h_{eff}=\frac{1}{\rho _{eff}}\sum\limits_{k=1}^{\Pi }{\varepsilon _{k}\rho _{k}}h_{k}$	(6)

Therefore the energy equation, neglecting pressure effects and viscous dissipation, is:

$\frac{\partial }{\partial t}\left( \rho _{eff}h_{eff} \right)+\nabla \cdot \left( \rho _{eff}\mathbf{V}h_{eff} \right)=-\nabla \cdot \mathbf{{q}''}_{eff}+\sum\limits_{k=1}^{\Pi }{\sum\limits_{j=1\left( j\ne k \right)}^{\Pi }{{\dot{m}}'''_{jk}h_{k}}}+\dot{{q}'''}_{eff}$	(7)

(a) An actual interface between two phases, (b) an interfacial representation using the Donor-Acceptor scheme, and (c) a piecewise linear reconstruction scheme with the VOF method.

In the energy equation effective properties are incorporated and a latent heat term due to phase change is added.

The continuity, momentum, and energy equations can be solved by standard solution procedures, such as the SIMPLE class of algorithms in Patankar's works ^[6]^[7]. However, if the volume fraction equation is solved using a standard implicit or explicit scheme, the interface will quickly lose resolution due to numerical diffusion. The numerical diffusion will lead to inaccurate solutions and/or a solution that will not converge. Therefore, special consideration must be taken on interfacial cells to construct the interface so that the advection transport of fluid is representative of the physical problem. An actual interface between two phases, as well as the corresponding interface represented by two special interpolation schemes, is presented in the figure to the right.

As noted before, the first of such methods is the donor-acceptor scheme proposed by Hirt and Nichols ^[5]. If a cell is an interfacial cell, $0<\varepsilon _{k}<1$ , the fluid will be rearranged in the cell to be on one face, as shown in (b) of the figure to the right. The face on which the fluid will be rearranged will depend on the normal direction of the interface. The normal direction of the interface, with respect to phase k, can be calculated by the gradient of the volume fraction:

$\mathbf{n}_{k}=-\frac{\nabla \varepsilon _{k}}{\left\| \nabla \varepsilon _{k} \right\|}$	(8)

The component in which the normal is the greatest will occur where the fluid is perpendicular to the interface while the interface is either horizontal or vertical. Once this is done, one cell is designated as a donor, while its neighbor is designated as an acceptor. The amount of fluid leaving the donor cell is exactly equal to the amount of fluid entering the acceptor through each computational face. Also, the maximum amount of fluid that can leave a cell is equal to either the amount of fluid in that cell or the amount that would make another cell fill with fluid.

A more refined interface interpolation scheme was developed by Youngs ^[8], in which the interface was approximated as piecewise linear in each cell. The normal of the reconstructed interface in each cell is the same as the normal calculated in eq. (8), as shown in (c) of the above figure. The advection of fluid through each face is calculated in a manner similar to the donor-acceptor scheme. Both the donor-acceptor scheme and piecewise linear interpolation can only be run in transient simulations, and the time step must be kept small enough so that the fluid near the interface will only advance by one cell at a time. Also, the resolution of the interface is limited to the grid spacing of the computational mesh. Any surface waves that are smaller than the spacing of the mesh will be smoothed out, as seen in the above picture on the left side of the interface. Surface tension effects are important in many free surface problems. The surface tension effects can be applied as a body force in the momentum equations. This method is called the continuum surface force (CSF) model and was proposed by Brackbill et al. ^[9]. The curvature is defined as the divergence of the surface normal vector.

$K_{k}=\nabla \cdot \mathbf{n}_{k}$	(9)

The volume force due to surface tension, $F σ$ , is

$F_{\sigma }=\sum\limits_{k=1}^{\Pi }{\sum\limits_{j=1}^{k-1}{-\sigma _{jk}\frac{\varepsilon _{j}\rho _{j}K_{k}\nabla \varepsilon _{j}+\varepsilon _{k}\rho _{k}K_{j}\nabla \varepsilon _{k}}{\frac{1}{2}\left( \rho _{j}+\rho _{k} \right)}}}$	(10)

If only two phases are present, the body force reduces to

$F_{\sigma }=\sigma _{jk}\frac{\rho _{eff}K_{k}\nabla \varepsilon _{j}}{\frac{1}{2}\left( \rho _{k}+\rho _{j} \right)}$	(11)

It can be seen that the body force is proportional to the cell density. It is important to note that when surface tension forces are large compared to other flow characteristics, numerical inaccuracies of the surface tension forces in the CSF model create artificial currents called parasitic currents. Much work has been done to eliminate parasitic currents, such as the second gradient method proposed by Jamet et al. ^[10].

Despite the adverse effects of parasitic currents, the VOF method has been widely used and can give reasonably accurate results for a wide range of applications. It is robust in its handling of free surface problems with large interface distortion and can easily handle problems in which the free surface breaks apart, such as droplet formation. One drawback of the VOF method, however, is that the interface resolution is limited to the grid spacing. Therefore, a refined mesh is needed anywhere the interface is going to travel. This refinement can lead to many computational cells in regions of the mesh resided in by the interface for a short period of time, which will increase the total computational time of the solution. Therefore, advanced remeshing algorithms are needed for problems of this type.

References

↑ ^1.0 ^1.1 Faghri, A., and Zhang, Y., 2006, Transport Phenomena in Multiphase Systems, Elsevier, Burlington, MA.
↑ Faghri, A., Zhang, Y., and Howell, J. R., 2010, Advanced Heat and Mass Transfer, Global Digital Press, Columbia, MO.
↑ Harlow, F.H., and Welch, J.E., 1965, “Numerical Calculation of Time-Dependent Viscous Incompressible Flow,” Physics of Fluids, Vol. 8, pp. 2182-2189.
↑ Nichols, B.D., and Hirt, C.W., 1975, “Methods for Calculating Multidimensional, Transient Free Surface Flows Past Bodies,” Proc. of the First International Conf. On Num. Ship Hydrodynamics, Gaithersburg, MD, Oct. 20-23.
↑ ^5.0 ^5.1 Hirt, C.W. and Nichols, B.D., 1981, “Volume of Fluid (VOF) Method for the Dynamics of Free Boundaries,” Journal of Computational Physics, Vol. 39, pp. 201-225.
↑ Patankar, S.V., 1980, Numerical Heat Transfer and Fluid Flow, Hemisphere, Washington, DC.
↑ Patankar, S.V., 1991, Computation of Conduction and Duct Flow Heat Transfer, Innovative Research.
↑ Youngs, D.L., 1982, “Time-Dependent Multimaterial Flow with large Fluid Distortion,” numerical Methods for Fluid Dynamics, K.W. Morton and M.J. Baines, eds. Academic Press, pp. 273-285.
↑ Brackbill, J.U., Kothe, D.B., and Zemach, C., 1992, “A Continuum Method for Modeling Surface Tension,” Journal of Computational Physics, Vol. 100, pp. 335-354.
↑ Jamet, D., Torres, D. and Brackbill, J.U., 2002, “On the Theory and Computation of Surface Tension: The Elimination of Parasitic Currents through Energy Conservation in the Second-Gradient Method,” Journal of Computational Physics, Vol. 182, pp. 262-276.

Retrieved from "http://www.thermalfluidscentral.org/encyclopedia/index.php/Numerical_Simulation_of_Interfaces_and_Free_Surfaces"

Numerical Simulation of Interfaces and Free Surfaces

From Thermal-FluidsPedia

Interface Tracking Techniques

Stationary Grid Approach

References

Views

Personal tools

Navigation

Search