Volume of Fluid: Interface Capturing Method (Multiphase flow)

Reference: Slideshare

Multiphase flow classification

Multiphase flows have phase boundaries that have an inter-phase coupling with mass, momentum, and energy transfer between different phases. They can be classified into 2 categories based on phase interaction, Continuous-Continuous phase interaction & Dispersed-Continuous phase interaction.  

boat | GIF | PrimoGIF
Ocean Wave (Continuous-Continuous phase interaction)
Whiskey Drink GIF - Whiskey Drink Beer GIFs
Bubble motion in Liquid (Continuous-Dispersed phase interaction)

Continuous-Continuous phase interaction flow has a single continuous surface as an interface without any breaks or discontinuities. This continuous-discrete interface separates the two phases. Example: Ocean Wave, where water and air motion has a discrete interface separating both phases.  

Dispersed-Continuous phase interaction flow has distributed discrete interfaces where one phase is dispersed into the continuous phase. Example: Motion of bubble in a glass of Beer, Here the continuous phase is Beer (liquid) and the dispersed phase is bubbles (gas), Cavitation phenomena in the liquid jet is another example.  

2-Phase Flows Modeling

Multiphase Flow in CFD: Basics and Modeling | SimScale Blog
Reference: SimScale

There are many fluid flow applications where we have more than a single phase. One such special application is the interaction of immiscible fluids that have a continuous-continuous phase interaction type of behavior. Dam break, Tank sloshing, Stratified flow, and Liquid jet spray are such examples. In all such applications, we deal with at least 2 phases Liquid-Liquid or Liquid-Gas.

If we want to model such flows using the application of CFD, we have to capture the interface between two phases and track it over time as it advects with the flow under the application of boundary conditions. That’s where the concept of Interface capturing and Interface advection comes into the picture.

Challenges in Interface Capturing & Advection

  1. Enforcing mass, momentum, and energy conservation at the interface.
  2. Interface motion is inherently transient in nature which sometimes requires very small time steps in order to capture the rapidly changing fluid dynamic behavior of the flow.
  3. Modeling discontinuities in properties across the interface, especially large jumps in density and viscosity.
  4. Handling complex topologies and scale separation (When ocean waves approach close to the coastline, waves will start to break into droplets at a level of microscales).
  5. A finer mesh is required in order to capture these droplet breakups with a discrete interface and with the transient motion of these small scales (for example liquid droplets breakup from a continuous liquid jet as shown below) making the problem complicated and computationally expensive to tackle. You can find the few images below showing scale separation phenomena.
  6. The computational numerical method must be robust enough to capture small breaking scale structure and large continuous scale. This is required to be captured from the physics point of view. Enforcing conservation principles (mass, momentum, and energy) makes the system of PDEs stiff, and therefore numerical solver must be robust enough to solve such phase interaction.
Scale separation of ocean wave (NPTEL)
Reference: ResearchGate

Fluid Mechanics aspects with Interface

The approach is very simple, starting with deriving the equation of flow without interface. Later the mathematical representation of the moving/advecting interface is introduced in the equation. Appropriate jump conditions required to couple the equation with the interface is imposed by looking into the aspects of fluid dynamics of the interface.

Finally, we will introduce the “One-Fluid” formulation which is the origin of the Volume of Fluid Method where the interface is introduced as a singular distribution in equations written for the whole flow field.

One-Fluid” formulation has a single fluid with different properties in different regions across the interface. The interface definition is introduced using a singular distribution using a source term in the momentum equation (discussed in the later section of the blog).

Three are the major assumptions while deriving such equations:-

Reference: NPTEL
  • The Continuum hypothesis.
  • The hypothesis of the sharp interface.
  • Intermolecular forces are modeled by introducing the surface tension concept.

The hypothesis of the sharp interface is such that it separates the pair of fluids with different thermodynamic phases, such as solid & liquid or vapor or gas. Properties change across this interface.

Fluid Dynamics aspects with Interface

Across the interface of 2-phases, we can potentially have momentum and energy exchange. From the numerical point of view, we need to understand the fluid dynamic interaction of 2 phases across the interface such that mass conservation, momentum conservation, and any other appropriate constraint are properly imposed on the interface.  

Let’s consider a thin control volume dV with boundary dS across a portion of interface S. The thickness of the control volume is tending to zero such that no accumulation of property takes place inside it. 

 U1 & U2 are the velocity of the fluid in regions Fluid 1 and Fluid 2 respectively.  V is the normal velocity of the interface itself.

So, the jump condition at the interface is :-

The application of mass conservation on the interface leads to a jump condition across the interface which is basically “the jump in velocity across the interface is 0“, in a simple way both phases have to have the same velocity at the interface! This condition must be satisfied when modeling the 2-phase flow numerically.  

Momentum conservation at the interface brings the purview of surface tension effect acting on the interface due to both phase’s interaction. Generally, Surface tension is modeled as a source term in the momentum equation that is discussed later in the One-Fluid Formulation section below. Here we will see the physical and mathematical description of surface tension force.

The interface is not a thermodynamically optimum region because the molecules on both sides of the interface, prefer to be at Fluid1 or Fluid2 by the virtue of surface tension. From the mechanical point of view, Surface tension would have the effect of stretching the interface (by pulling the line segments that construct the interface from the numerical point of view). The work done during this stretching will increase the Interfacial Free Energy.

Let’s apply momentum conservation across the CV.

The above expression gives the mathematical description of the surface stress tensor acting on the interface relating to surface tension coefficient and interface curvature. We can split the expression into normal and tangential stress conditions at the interface.

Jump condition, Normal stress, and Tangential stress condition must be properly imposed on the interface while modeling the fluid interface between 2 phases numerically. 

Now, we know what conditions must be satisfied at the interface in the numerical simulation in order to capture the fluid dynamic physics at the interface accurately. 

Marker Function

If we are capturing a sharp interface then there would be a jump in the fluid properties value. Marker function is a concept where we define the presence of one of the phases using the Heaviside step function H(x,y) and try to capture sharp interface.  

referenece: openFOAM

Instead of identifying the interface by explicitly specifying the location of it at every point on an interface using a height function (link), a Marker Function is defined in the whole domain where we are marking the different phases using some kind of a functional value of the Heaviside step function.  

In a mathematical sense, it is called a sharp interface because there is a step jump in the value of H as we move from Fluid1 to Fluid2.  

One-Fluid Formulation

In the sharp interface approach, the density, viscosity, and any other fluid property experience a jump across the interface. With the “One-Fluid Formulation” approach, it is possible to write one set of governing equations for the whole computational domain occupied by various phases.   

Here, various phases are treated as one-fluid with variable fluid properties that abruptly change across the phase boundary.

Now, In order to take account of the extra force (surface tension) acting at the phase interface, It is necessary to add singular terms in the equation as counterparts of the jump condition. So, basically, we do not explicitly enforce the jump condition, but we introduce a source term in the momentum equation that gives a similar effect to the presence of interface location. In a spatial sense, this source term can become active in the region close to the interface by using the Dirac delta function (which has a spike-like mathematical behavior).

The derivation of the “one-fluid” formulation equation is exactly the same as the basic governing equations, except that we have to add the surface tension as a body force to the momentum equation.

Extra source term has surface tension force coupled with Dirac delta function which controls that the value of surface tension will be zero in the bulk fluid and at interface only it will get a finite value.

The mathematical description of the extra source term that models the surface tension force at the interface allows us to solve the governing equation.

Advecting the Fluid Interface

When the governing equations are solved on a fixed grid, using the One-Fluid Formulation approach for the whole flow field, the different phases must be identified in some way. This is done by using a “Marker Function” that takes the different values in the different regions of the domain.

As the fluid moves and interface shape & location change, the marker function must be updated and the interface must be reconstructed from the marker function.  

A fluid element has the value of marker function associated with it in the course of its movement which is mathematically represented by the material derivative of the marker function is zero. The individual fluid volume carries along with them the marker function value as they move through the flow field. In other words, the fluid particle no matter where it moves in the flow field does not lose the value of the marker function that it is carrying along with them. That’s what is called Material Derivative.   

The equation shown above governs the advection (motion) of the fluid interface.

Method of using marker function on the interface are called front capturing method where marker function is used to demarcate the 2 phases.  

Looking at the above image, we can see the exact location of the sharp interface which is where the Heaviside function drops sharply from 1 to 0. The red color region is occupied by volume fraction C which is an integrated effect of the Heaviside function. Heaviside function fills the cell using volume fraction by maintaining the area conservation in cell j.

Area under the H curve in jth cell = Area under the C curve in jth cell

Advecting the Volume Fraction Equation

Let’s assume flow to be incompressible and flow is being advected at a constant velocity u in the positive x direction. We can write advection equation for Volume fraction C in conservative form.

First of all, We have to understand how does interface advect in the flow field then only we will be able to discretize the volume fraction advection equation in order to capture advection of the interface.

For simplicity, We will consider the 1D grid as already shown above.  

At time t, the sharp interface is sitting at the jth cell because that’s where Heaviside function H is standing. Therefore cells j-2 & j-1 are completely filled with the volume fraction C. Cell j is partially filled with volume fraction because the Heaviside function hasn’t crossed the cell j yet. As time progressed, the Heaviside function move toward the right (positive-x) direction because flow velocity is u>0 for a 1D grid. 

Once the interface/ Heaviside function reaches “j+1/2” face of the jth cell. The volume fraction C will fill the jth cell completely.

Keeping the above physics in mind we can discretize the volume fraction advection equation in 1D.  

It boils down to appropriately choosing the value of volume fraction at the cell j+1/2 & j-1/2 faces in order to capture the advection of the interface correctly.

If we are going to use a first-order upwind scheme to discretize the above equation, It will introduce numerical diffusion and it will artificially diffuse the sharp interface which is physically not correct.  

Bubble rise in a liquid column is one of the applications of mathematical modeling of interface advection. The bubble has a volume fraction is 1 where the region is red, then there is a sharp jump to the blue region which has a volume fraction of 0. That region is occupied by water. Bubble with its interface is being advected in the liquid column.      

Now if we are trying to model the sharp interface of the bubble & we found that our discretization scheme is artificially diffusion the interface, what it means is that artificially the bubble will get elongated whereas physically it should not be.  

That’s where the Numerical scheme plays an important role in predicting the phase boundary/interface correctly as they advect with the flow.

Volume of Fluid (VOF) Method

The volume of Fluid is one of the famous Interface tracking methods. It consists of reconstruction schemes/techniques of the interface from the volume fraction data stored in different finite volume cells.

Source: NPTEL

The motivation of the VOF method comes from one fluid formulation approach where we solve a single set of the Continuity, Momentum, and Energy equation in the whole computational domain occupied by various phases. We are solving an addition equation, the Volume fraction advection equation that governs the interface advection and computes the volume fraction of different phases in each cell of the computational domain. The volume fraction data let us define different phase properties in different regions.   

Reference: NPTEL

The continuity equation looks similar to the standard continuity equation, but the definition of density is coupled with volume fraction as shown in the above image. Volume fraction C = 1 is the region of fluid 1 and C = 0 is the region of fluid 2. In computational cells with volume fraction C = 1, the density will take the fluid 1 density value and in Computational cells with volume fraction C = 0, the density will take the fluid 2 density value. In cells on which interface is present, the value of C will lie between 0<C<1 so at the interface it will take some averaged value of density. The same analogy can be applied to Viscosity.   

Because we have expressions for density and viscosity coupled with the volume fraction, we don’t have to solve two continuity equations for two different phases. We will manage with this single continuity equation but the fluid properties guiding this single continuity equation will change as per the volume fraction field. 

For the momentum equation, again we have to solve a single momentum equation where the fluid properties are the function of volume fraction. On the RHS, we have an additional source term to account for Surface Tension force at the interface between two phases. 

Recalling the mathematical description of the surface tension force that was already discussed in the “Fluid Dynamics aspects with Interface” section. The surface tension force in the RHS of momentum equation is coupled with the delta Dirac function. Delta Dirac function ensures the value of surface tension force is zero in the bulk of the fluid and has a finite value only at the interface.  

We can see surface tension force is coupled with interface unit normal vector which can be calculated by computing volume fraction gradient. We are not going into the volume fraction gradient calculation in this article, there are many gradient schemes available in the literature. What is important to note here is, the accurate calculation of volume fraction C because that decides the accuracy of the direction of the unit normal, If that is not calculated properly then it will induce numerical errors in the calculation of the surface tension force term. 

As time progresses, the Interface advects and its location change. The volume fraction advection equation ensures this advection. Let’s see the discretization of the volume fraction advection equation. For simplicity, we will only consider the 1D domain.  

The below discretization scheme ensures advection of the sharp interface. For the 1D domain, If u > 0 then the flux F through the left face “j-1/2” of the jth cell is always u and the flux through the right face “j+1/2” of the jth cell is calculated using the below a condition which basically checks the exact location of interface and computes the flux from the right face.   

Again, The above discretization scheme ensures advection of the sharp interface without any artificial diffusion.

It is important to note that, In the 1D case, the orientation of the interface is kept vertical meaning the volume fraction is filling the cell from left to right. Now for the 2D case, we can orient the interface in a horizontal manner and the interface can be filled from top and bottom.  

It is easy to do such calculation for the 1D case but it gets more complicated when we do it for 2D & 3D cases.  

When the interface orientation is kept in a straight verticle line (like in the 1D case) then this method is known as SLIC Method (Simple Line Interface Calculation), It is one of the Interface reconstruction methods.  

Interface Reconstruction Techniques

Interface reconstruction techniques are methods of reconstructing the interface based on the interface advection. As time goes on, the Interface is going to advect, it may deform and change its shape, it may even rapture. Interface advection is governed by the volume fraction equation which is solved using the Volume of Fluid method. During the iterative process of numerical simulation, after solving the volume fraction advection equation, we will obtain the new value of volume fraction in each and every cell of the domain. We need to reconstruct the interface such that it governs the newly obtained volume fraction field.  

As we reconstruct the interface, we have to recompute the interface unit normals so that we can calculate the surface tension force that is required to close the momentum equation.  

We are going to discuss 3 Interface Reconstruction Techniques: SLIC Method, Hirt & Nichol’s Method, PLIC (Piecewise Linear interface Calculation).

SLIC Method

We have already seen above how the SLIC method works in 1D, When it is applied for 2D & 3D case, It is done using time splitting. In 2D, we take care of advection along x & y separately. For x-direction, the interface is aligned vertically based on the volume fraction field. Then, the interface is advected by doing time integration of fluxes in the volume fraction advection equation. Once that is done, the same process is repeated for the y-direction by aligning the interface horizontally. Scheme used for time integration of the fluxes is already shown in the previous section.  In 3D, the interface is additionally advected along z-direction. 

Above, we have a 2D structured grid on which we can see the original interface. This interface will deform & advect to a new location as time goes on. We are solving the volume fraction advection equation that governs this physics and we are reconstructing the interface using the SLIC method. As we can see, there are dotted horizontal and verticle lines in each cell where the interface is located. This is how the interface is reconstructed using SLIC Method. After every numerical time-step, the Interface is advected to a new location and we will obtain a new volume fraction field by solving the volume fraction advection equation. We will have to reconstruct the new location and shape of the interface as per the new volume fraction field using the SLIC method after every time-step.  

Hirt & Nichol’s Method

In this method, we only advect the interface in one direction either horizontal or vertical, and do the flux calculation. That one direction is decided based on the normal vector of the interface. If the interface normal vector is more aligned to the horizontal direction then the interface is reconstructed along the x horizontal axis otherwise the interface is reconstructed along the y vertical axis.  

We can see in the above image only 1 front (dotted lines) is used in order to reconstruct the interface based on the original interface orientation towards the x or y-axis.  

PLIC Method

Piece-wise Linear Interface Calculation method is more accurate than the SLIC & Hirt Nichol method. Here, the interface is approximated by a straight line segment in each cell, but the line can be oriented arbitrarily with respect to the coordinate axis. The orientation of the line is determined by the interface unit normal vector. Once the interface in each cell is constructed, the fluxes from one cell to another are computed by geometrical consideration.  

We can see the dotted lines are slanted, trying to capture the exact shape of interface by using a piece-wise linear construction.  

Geometric Interface Reconstruction using PLIC Method 

We can see the interface advection equation written above in terms of the Heaviside function. When we apply it in an area-weighted sense for the volume fraction function C & then discretize it, we will see the equation works out as shown below:-  

We can integrate above equation over a time deltat, then we will obtain the equation as shown below:-

Above is the discretized equation, if we introduce it in all cells of the domain & sum them up over all the cells with appropriate BC, what happens is, that internal fluxes would cancel out in pairs. Flux out of the face ‘i+1/2,j’ from the jth cell is flux going into the cell ‘j+1’ through that face, that’s how they cancel each other.  

Finally, we will be able to show that when we sum it over all the cells, we are able to maintain the conservation of the total area (in 2d) of the volume fraction.  

Conservation of volume fraction

Consevation of volume fraction states that when the interface advects and deforms, the volume fraction field will change over time, but the total volume of each phase will remain constant. This is an extremely important constraint that must be satisfied.

Whatever numerical scheme that we use to approximate these fluxes, it must not ever cause any unphysical overshoots i.e. C >1 or undershoots i.e. C<0. This will lead to numerical error and therefore non-satisfaction with the conservation principles and it can trigger numerical instabilities as well.  

Below, we can see a 2D grid example where the interface is reconstructed using the PLIC technique based on the volume fraction distribution, let’s say on the Nth time-step.

m is the outward unit normal vector of the interface that will locally change as interface shape is continuously changing.  

It is known that the outward pointing normal is equal to the negative gradient of the volume fraction field. Now we can define the inteface in the manner shown below:-

Now in order to move the interface as per the new volume fraction field at N+1th time-step, the Geometrical approach is used where the interface line is moved along the normal direction through the adjustment of parameter alpha until the area under the interface at N+1th time-step for that particular cell equals to h2Cijn+1.  

We are trying to say that, we have obtained the value of Cijn+1 by solving our discretized equation and now we are trying to figure out how to move the interface such that it fills exactly that much volume fraction. For example, let’s say in cell i,j we had volume fraction C = 0.60 at n time level and now we have obtained the volume fraction C = 0.7 at n+1 time level. That means at time-level n, 60% volume of the cell was filled with phase 1. Now at time-level n+1, 70% volume of the cell has gotten filled due to interface advection. Now we have to move the interface location such that it fills the cell exactly to 70%. That is what interface reconstruction is all about. We are saying that the interface can be moved by adjusting the parameter alpha.

Geometrically, how it is done can be shown with a diagram!

Let’s focus on a single cell and interface with a certain volume fraction filled with the shaded region.

The above equation can be written in the form as shown below:-

That’s how geometrically we reconstruct the interface!

End Note!

Interface Reconstruction becomes a complex and tedious task when it is required to be done in a 3D case. There are many advanced Interface Reconstruction techniques available in the literature.  

Leave a Reply