Computational methods for free surface flow

The free surface should be a sharp boundary separating the two fluids. There should be no flow through this boundary, i.e.,

${\displaystyle [(v-v_{b})\cdot n]_{fs}=0,}$ or

${\displaystyle {\dot {m}}_{fs}=0}$

where ‘fs’ stands for free surface. This implies that the normal component of the velocity of the fluid at the surface is equal to the normal component of the velocity of the free surface.

The forces acting on the fluid at free surface should be in equilibrium, i.e. the momentum is conserved at the free surface. The normal forces on either side of the free surface are equal and opposite in direction and the forces in tangential direction should be equal in magnitude and direction.

${\displaystyle (n\cdot T)_{l}\cdot n+\sigma K=-(n\cdot T)_{g}\cdot n,}$

${\displaystyle (n\cdot T)_{l}\cdot t-{\frac {\partial \sigma }{\partial t}}=(n\cdot T)_{g}\cdot t,}$

${\displaystyle (n\cdot T)_{l}\cdot s-{\frac {\partial \sigma }{\partial s}}=(n\cdot T)_{g}\cdot s,}$

Here σ is the surface tension, n, t and s are unit vectors in a local orthogonal coordinate system (n,t,s) at the free surface (n is outward normal to the free surface while the other two lie in the tangential plane and are mutually orthogonal). The indices 'l' and 'g' denote liquid and gas, respectively and K is the curvature of the free surface.

${\displaystyle K={\frac {1}{R_{t}}}+{\frac {1}{R_{s}}}}$

with R_t and R_s being radii of curvature along coordinates t and s.

The surface tension σ is force per unit length of a surface element and acts tangential to the free surface.

${\displaystyle f_{\sigma }=\sigma \ dl}$

For an infinitesimally small surface element dS, the tangential components of the surface tension forces cancel out when σ = constant, and the normal component can be expressed as a local force that results in a pressure jump across the surface.

This is a methods which treat the free surface as a sharp interface whose motion is followed. In this method, boundary-fitted grids are used and advanced each time the free surface is moved.
Interface tracking method is useful in situations like calculation of flow around submerged bodies. This is done by linearizing about unperturbed free surface, so a height function is introduced for the free surface elevation relative to its unperturbed state.

${\displaystyle Z=H(x,y,t)}$

This gives the kinematic boundary condition a new form:

∂H/∂t = u_z - u_x∂H/∂x - u_y∂H/∂y

This equation can be integrated and the fluid velocity at free surface can be obtained either by extrapolation from the interior or by using dynamic boundary condition. For the calculation of flow, FV method is widely used. The steps for a fully conservative FV method of this type are:

• momentum equation is solved to obtain velocity ${\displaystyle u}$_{${\displaystyle i}$}^{${\displaystyle *}$} at the current free surface using specified pressure.
• Local mass conservation is enforced in each CV by solving a pressure-correction equation. Mass is conserved both globally and locally, but velocity-correction is produced at free surface giving a non-zero mass flux.
• Position of free surface is corrected to compensate for the non-zero mass flux with the volume flux due to the movement of the each free-surface cell face by enforcing the kinematic boundary conditions.
• Iterate until no further correction is needed, satisfying the continuity and momentum equations.
• Advance to the next time step.[1]

The main problem with the algorithm in this procedure is that there is only one equation for one cell but large number of grid nodes moving. To avoid instability and wave reflection, the method is modified as follows:
From the previous steps, we can calculate the volume of fluid to be flowed in or out of the CV to have mass conservation. To obtain the coordinates of CV vertices at free surface, we have more unknowns and less equations due to single volumetric flow rate for each cell.

Hence the CVs are defined by the cell face centers rather than vertices and vertices are obtained by interpolation. This gives a tridiagonal system for 2D and can be solved using TDMA method. For 3D, the system is block tridiagonal and is best solved by one of the iterative solvers.

In computation of two-fluid flows, in some cases the interface might be too complex to track while keeping the frequency of re-meshing at an acceptable level. Not being able to reduce the frequency of re-meshing in 3D might introduce overwhelming mesh generation and projection costs, making the computations with the interface-tracking technique no longer feasible. In such cases, interface-capturing techniques, which do not normally require costly mesh update steps, could be used with the understanding that the interface will not be represented as accurately as we would have with an interface-tracking technique.[2] Methods which do not define the interface as sharp boundary. A fixed grid extends beyond the free surface over which the computation is performed. To determine the shape of the free surface, the fraction of each cell near the interface is computed that is partially filled.

MAC scheme was proposed by Harlow and Welch in 1965. In this method, a massless particle is introduced at the initial time at the free surface. The motion of this massless particle is followed with the passage of time.

Benefit: This scheme can treat complex phenomena like wave breaking.

Drawback: In three dimensional flow solving the equations governing fluid flow and also following the motion of a large number of markers both simultaneously demands high computational power.

VOF scheme was proposed by Hirt and Nichols in 1981. In this method, fraction of the cell occupied by the liquid phase can be calculated by solving the transport equation.[3] The transport equation is:

∂c/∂t + div(cv) = 0

where c is the fraction of control volume filled. c=1 for completely filled and c = 0 for completely empty control volumes.
So in total, for VOF method, one has to solve three forms of equations, conservation equations for mass, conservation equations for momentum, equation for filled fraction for each control volume.

NOTE: IN INCOMPRESSIBLE FLOWS, ABOVE EQUATION GIVES SAME RESULTS WITH c AND 1 - c MAKING THE ENFORCEMENT OF MASS CONSERVATION A MUST.

Since the higher order schemes are preferred over lower order schemes to prevent artificial mixing of the two fluids, it is important to prevent the overshoots and undershoots due to the condition 0 ≤ c ≤ 1. For such problems, modifications were made to MAC and VOF schemes.

Marker and micro-cell method in which local grid refinement is done according to the following criteria:

only the cells having 0 < c < 1 are refined.

This method is more efficient than MAC scheme because only the cells at the boundary are refined. But in this method, the free surface profile is not sharply defined.

There are some fluid flows which do not come under either of the category, for example, bubbly flows. For the computation of such two-phase flows which do not come under any of the above discussed categories, elements are borrowed from both surface-capturing and surface-tracking methods. Such methods are called hybrid methods. In this method, fluid properties are smeared over a fixed number of grid points normal to the interface. Now, as in interface capturing method, both fluids are treated as single fluid with variable properties. Interface is also tracked as in interface-tracking method to prevent it from smearing by moving the marker particles using the velocity field generated by the flow solver. marker particles are added and removed to maintain the accuracy by keeping the approximate spacing between them equal.