Continuum Approach for Interfaces and Free Surfaces

From Thermal-FluidsPedia

Jump to: navigation, search

There are many engineering problems that involve separated phases, in which the interface between phases is clearly defined but 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.

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 Fig. 5.29. The techniques are as follows:

1. Lagrangian Techniques: grid and fluid move together, interface is directly captured

2. 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

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.

Each of the above methods will be briefly decribed below.


Lagrangian Techniques

In a Lagrangian technique, fluid particles are tracked directly. Therefore, a particle on the interface will also be tracked directly, so the interface will automatically be resolved. The Lagrangian approach is used by Shopov et al. (1990) to model a droplet interacting with a wall. Even though a Lagrangian approach directly captures an interface, an Eulerian approach is often preferable to a Lagrangian approach to solve the bulk flow of a fluid. An Eulerian approach is generally preferable, because fluid flow is usually complex, and fluid particle path lines often intersect and/or disperse, making fluid interaction very difficult to manage with a Lagrangian method. Techniques 2, 3, and 4 all use an Eulerian point of view to solve the governing equations in the bulk fluids with special consideration to track the interface.

Stationary Grid Approach

The second class of numerical techniques for multiphase fluids 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 (1965). 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, (Nichols and Hirt, 1975). Further development of this class of models lead to the Volume of Fluid method (VOF) by Hirt and Nichols (1981); 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, and it 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]


When the volume fraction is between 0 and 1 in a computational cell, that 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


Therefore, (k-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}}}


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 (5.249)

The continuity equation is the same as eq. (3.38), except it 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 }}


The momentum equation is also the same as eq. (3.52), with the same exceptions as the continuity equation and the momentum transfer due to phase change. In the energy equation, the enthalpy is mass averaged instead of volume averaged.

Figure 5.30 (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

{{h}_{eff}}=\frac{1}{{{\rho }_{eff}}}\sum\limits_{k=1}^{\Pi }{{{\varepsilon }_{k}}{{\rho }_{k}}}{{h}_{k}} (5.251)

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}}


The energy equation has the same form as eq. (3.79), except 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 (1980). 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 Fig. 5.30.

As noted before, the first of such methods is the donor-acceptor scheme proposed by Hirt and Nichols (1981). 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 show in Fig. 5.30(b). 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|}


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, and 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 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 (1982), 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. (6), as shown in Fig. 5.30(c). The advection of fluid through each face is calculated in a manner similar to the donor-acceptor scheme. Both the donor-acceptor 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 Fig. 5.30 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. (1992). The curvature is defined as the divergence of the surface normal vector.

{{K}_{k}}=\nabla \cdot {{\mathbf{n}}_{k}}


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)}}}


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)}


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. (2002).

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 other drawback of the VOF method 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.

Phase Interface Fitted Grid

The third method, in which a phase interface fitted grid is employed, can be very useful for multiphase problems in which the interface does not deform greatly or break apart (although it can be handled with complex algorithms). The benefit of these models is that they directly capture an interface; therefore, a phenomenon occurring at a phase interface is not subject to interpolation error, such as the interface location itself. The third class of numerical method used to capture a free surface is the semi-Lagrangian approach. Since a Lagrangian approach is not feasible for most fluid mechanics problems, only the interface is tracked in a Lagrangian manner, while the rest of the mesh can be manipulated to maintain good mesh qualities. Therefore, this type of technique can be deemed “semi-Lagrangian,” and is a front tracking technique. Rice and Faghri (2005a) developed a semi-Lagrangian approach that works effectively on the interface between phases with phase change. They used a transient method in which the present interfacial velocity, , is calculated based on the error in mass flux through each computational face from the previous time step. Therefore the interfacial velocity is lagging by one time-step. The conservation of mass between phase k and phase j is

{\dot{m}}''={{\rho }_{k}}\left( {{\mathbf{V}}_{k}}-{{\mathbf{V}}_{I}} \right)\cdot \mathbf{n}={{\rho }_{j}}\left( {{\mathbf{V}}_{j}}-{{\mathbf{V}}_{I}} \right)\cdot \mathbf{n}={{{\dot{m}}''}_{actual}}


The interfacial continuity equation represents the actual mass flux at the interface. Instead of implicitly calculating the interfacial velocity of the next time step (n+1), the interfacial velocity from the previous time step (n) is used. By lagging the interfacial velocity, there is an inherent error in mass flux at the interface.

{{m}''}^{n+1}_{error}={{\rho }_{k}}\left( {V}_{k}^{n+1}-{V}_{I}^{n} \right)*{n}-{{{{m}}''}^{n+1}}_{actual}={{\rho }_{k}}\left({V}_{I}^{n+1}-{V}_{I}^{n} \right)*{n}



The normal velocity of phase j at the interface is fixed by continuity. However, the normal velocity of phase k at the interface is not fixed, and fluid is allowed to pass through this interface without regard to the interfacial mass balance. Figure 5.31 demonstrates how the interfaces moves and shows that there is an error in the mass of the system due to the interfacial movement. In order to compensate for the error in mass through each computational face, the new velocity can be found by rearranging eq. (5.258) and integrating over time and area.

\Delta t\int_{{{A}_{I}}}^{{}}{{{\rho }_{k}}\mathbf{V}_{I}^{n+1}\cdot \mathbf{n}\,dA}=\Delta t\int_{{{A}_{I}}}^{{}}{\left( {{{{{m}}''}}^{n+1}}_{error}+{{\rho }_{k}}\mathbf{V}_{I}^{n}\cdot \mathbf{n} \right)dA}



In discretized form, eq. (5.259) becomes

{{\left[ \Delta t\left( {{\rho }_{k}}\mathbf{V}_{I}^{n+1}\cdot \mathbf{n} \right)A \right]}_{i}}={{\left[ \Delta t\left( {{{{{m}}''}}^{n+1}}_{error}+{{\rho }_{k}}\mathbf{V}_{I}^{n}\cdot \mathbf{n} \right)A \right]}_{i}}={{m}_{i}}


The subscript i represents a single face on the interface. Rice and Faghri (2005a) suggest moving the interface at the beginning of each time step based on the error in mass flux and interfacial velocity of the previous time step, so that the mass the interface gains by movement is equal to mi. They specified a direction in which the interface location \left( {{\delta }_{I}} \right) will always be a function of {{\delta }_{I}}=f\left( \eta  \right), and moved the nodes on the interface normal to this direction. Specifying the coordinate system can prevent the interfacial face area of a single cell from becoming zero or negative. Some knowledge of the problem is necessary in order to specify the direction in which the interface can move.

The pressure of phase k at the interface should be defined from the interfacial momentum equation. When a second phase also influences the interface, this interface will move with the interface of phase k. However, the mass flux that passes through this interface should be identically equal to {{{\dot{m}}''}_{actual}}. This side of the interface affects the interfacial movement through the momentum equations. When surface tension effects are important, they are directly applied as an interfacial boundary condition in the interfacial momentum equation. This technique was shown by Rice and Faghri (2005b) to eliminate the parasitic currents caused by surface tension effects that are experienced in the VOF method. Also, this technique directly tracks the interface; therefore, the interface is always directly resolved, no matter what the mesh spacing is. For problems with large fluid distortion, it is very difficult to define a coordinate system of which the interface will always be a function; therefore, this coordinate system will need to be able to move over time. If the interface is largely distorted, advanced remeshing techniques will be needed for the computational domain away from the interface to maintain a good quality of computational cells.

Rice and Faghri’s (2005a, b) technique used a transient solver that allowed a small amount of liquid mass to pass through the interface. During each new time-step, the liquid that passed through the interface during the previous time-step was recaptured in a manner in which mass is directly conserved. Even though this technique is generally not suitable for interfaces undergoing large deformations, such as jet breakup, it has the advantage of directly capturing the interface. Therefore, this technique offers a great advantage in problems involving multiple components, whose interfacial values are not continuous.

Front Tracking Methods

The fourth class is a combination of the second and third class of multiphase and interface modeling techniques, and is called the front tracking technique. One form was developed by Glimm et al. (2001). This method uses a fixed grid that is modified only near the interface. The interface is modeled as a separate moving grid. The interface can move, and it separates the stationary grid cells into two cells near the interface, representing each fluid. Therefore, each fluid is treated separately. Another front tracking technique developed by Tryggvason et al. (2001) also uses a separate grid to track the interface; however, the phases are considered together, and a single set of governing equations is solved for the whole field. Therefore, the stationary cells that lie below the interface have a shared density and velocity field. Also, since the interfacial area is free to expand or contract, the subgrid that tracks the interface must have the ability to add or subtract nodes as needed. Further consideration of this method is discussed in the following section.

The main concept of the front tracking method is to identify a fluid with a Heaviside function, H. This function is one at a location where a particular fluid exists and zero elsewhere. The interface is located where the H function changes from zero to one, which is the location where the delta function, δ, is nonzero. H can be expressed as a function of the δ function.

H\left( \mathbf{x},t \right)=\int_{\Delta V}{\delta \left( \mathbf{x}-{{\mathbf{x}}_{I}} \right)dV}


The integration is performed over an entire volume, V, with δ being nonzero only on the interface between the phases, AI, at location xI. The calculation of the fluid properties, Φ, is straightforward with the Heaviside function.

\Phi \left( \mathbf{x},t \right)={{\Phi }_{k}}H\left( \mathbf{x},t \right)+{{\Phi }_{j}}\left( 1-H\left( \mathbf{x},t \right) \right)


Note, that when these values are integrated over a finite volume, the properties used in the calculation are volume averaged. The calculation of the interface is done in a Lagrangian fashion.

\frac{\text{d}{{\mathbf{x}}_{I}}}{\ \text{d}t}\cdot \mathbf{n}={{\mathbf{V}}_{I}}\cdot \mathbf{n}


To write the equivalent expression using a kinematic equation, the gradient of H must be defined

\nabla H=\int_{{{A}_{I}}}^{{}}{\mathbf{n}\delta \left( \mathbf{x}-{{\mathbf{x}}_{I}} \right)\text{d}s}


where s is a surface that is integrated over AI. Therefore, the kinematic equation using H to represent the interface is:

\frac{\partial H}{\partial t}=-{{\mathbf{V}}_{I}}\cdot \nabla H=-\int_{{{A}_{I}}}^{{}}{{{\mathbf{V}}_{I}}\cdot \mathbf{n}\delta \left( \mathbf{x}-{{\mathbf{x}}_{I}} \right)\text{d}s}


To calculate the interfacial velocity, the overall continuity equation can be considered.

\frac{\partial }{\partial t}\left( \rho  \right)=-\nabla \cdot \left( \rho \mathbf{V} \right)


The evaluation of fluid properties and flow values uses the H function as shown in eq. (17). The overall continuity equation can be written as:

\frac{\partial }{\partial t}\left( {{\rho }_{k}}H+{{\rho }_{j}}\left( 1-H \right) \right)+\nabla \cdot \left( {{\rho }_{k}}{{\mathbf{V}}_{k}}H+{{\rho }_{j}}{{\mathbf{V}}_{j}}\left( 1-H \right) \right)=0


Using the chain rule, the continuity equation is broken into three regions: the bulk regions for phase k and phase j, and the interface. The interfacial continuity equation is:

\left( {{\rho }_{k}}-{{\rho }_{j}} \right)\frac{\partial }{\partial t}\left( H \right)=-\left( {{\rho }_{k}}{{\mathbf{V}}_{k}}+{{\rho }_{j}}{{\mathbf{V}}_{j}} \right)\cdot \nabla \left( H \right)


Applying eqs. (18) and (19) to eq. (21) yields:

\int_{{{A}_{I}}}^{{}}{\left( {{\rho }_{k}}-{{\rho }_{j}} \right){{\mathbf{V}}_{I}}\cdot \mathbf{n}\delta \left( \mathbf{x}-{{\mathbf{x}}_{I}} \right)\text{d}s}=\int_{{{A}_{I}}}^{{}}{\left( {{\rho }_{k}}{{\mathbf{V}}_{k}}+{{\rho }_{j}}{{\mathbf{V}}_{j}} \right)\cdot \mathbf{n}\delta \left( \mathbf{x}-{{\mathbf{x}}_{I}} \right)\text{d}s}


The interface velocity can be calculated using eq. (21). The momentum equation is:

\frac{\partial \rho \mathbf{V}}{\partial t}+\nabla \cdot \left( \rho \mathbf{VV} \right)=\nabla \cdot \tau +\rho \mathbf{g}+\int_{A{{}_{I}}}^{{}}{\sigma \kappa \mathbf{n}\delta \left( \mathbf{x}-{{\mathbf{x}}_{s}} \right)\text{d}s}


The last term represents the surface tension effects where the twice the mean curvature is κ. Also, the energy equation can be written as:

\frac{\partial }{\partial t}\left( \rho {{c}_{p}}T \right)+\nabla \cdot \left( \rho \mathbf{V}{{c}_{p}}T \right)=-\nabla \cdot \mathbf{{q}''}-\int_{{{A}_{I}}}^{{}}{{\dot{m}}''\left[ {{h}_{kj}}+\left( {{c}_{p,k}}-{{c}_{p,j}} \right)\delta \left( \mathbf{x}-{{\mathbf{x}}_{s}} \right) \right]}\text{d}s


Note that pressure effects and viscous dissipation are neglected. The last term represents the latent heat release in a phase change process.

The interface is made of a separate mesh that goes through the stationary mesh. It generally does not lie on top of the points in which the conservation laws are applied. Therefore, the interfacial values are calculated by grid interpolation, and mass and momentum of the advection front is not necessarily conserved. The error in mass and momentum conservation at the interface can be reduced to an acceptable level by refining the mesh.


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.

Glimm, J., Li, X.L., Liu, Y., and Zhao, N., 2001, “Conservative front tracking and level set algorithms,” Proceedings of National Acedemic of Science: Applied Mathemitics, Vol. 98, pp. 14198-14201.

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.

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.

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.

Patankar, S.V., 1980, Numerical Heat Transfer and Fluid Flow, Hemisphere, Washington, DC.

Rice, J., and Faghri, A., 2005a, “A New Computational Method to Track a Liquid/Vapor Interface with Mass Transfer, Demonstrated on the Concentration Driven Evaporation from a Capillary Tube, and the Marangoni Effect,” Proceeding of the 2005 ASME International Mechanical Engineering Congress and Exposition, Nov. 5-11, Orlando, FL.

Rice, J., and Faghri, A., 2005b, “A New Computational Method for Free Surface Problems,” Proceeding of the 2005 ASME Summer Heat Transfer Conference, July 17-22, San Francisco, CA.

Shopov, P.J., Minev, P.D., Bazhekov, I.B., and Zapryanov, Z.D., 1990, “Interaction of a deformable bubble with a rigid wall at moderate Reynolds numbers,” Journal of Fluid Mechanics, Vol. 219, pp. 241-271.

Tryggvason, G., Bunner, B., Esmaeeli, A., Juric, D., Al-Rawahi, N., Tauber, W., Han, J., Nas, S., and Jan, Y., 2001, “A Front-Tracking Method for the Computations of Multiphase Flow,” Journal of Computational Physics, Vol. 169, pp. 707-759.

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.