A coupled volume-of-fluid and level set method in interDyMFoam solver

Ahmad Baniabedalruhman1

1Yarmouk University, Irbid, Jordan

Vibroengineering PROCEDIA, Vol. 30, 2020, p. 210-213. https://doi.org/10.21595/vp.2020.21342
Received 19 February 2020; accepted 27 February 2020; published 2 April 2020

Copyright © 2020 Ahmad Baniabedalruhman. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Creative Commons License

InterDyMFoam is a solver for calculating immiscible two-phase flow in an open source OpenFOAM software. The solver uses the volume-of-fluid method (VOF) for capturing the fluid-fluid interface. The solver is modified to work with the coupled VOF and level set (LS) method. In the proposed method, the VOF is used to capture the interface since it is a mass conserving and the LS is used to calculate the curvature and physical quantities near the interface. A bubble rising in column liquids is used to compare the solvers. The proposed method gives an accurate calculations for the physical properties in two phase flow problems and handles the topological change of the interface.

Keywords: openFoam, interDyMFoam, volume of fluid, level set, bubble rising.

1. Introduction

The study of two-phase flow is very important because their various industrial applications. Their applications includes spray, lubrication and food stuff processing. Capturing the interface in incompressible two-phase flow is a challenging problem in computational fluid dynamics. The most popular methods for capturing the interface are the level set and volume-of-fluid. In the level set method [1], the signed distance function φ is used to describe the interface, where the value of φ is negative in the continuous phase, positive in the disperse phase and zero at the interface. Moreover, the level set function is smooth, giving an accurate curvature k. The VOF [2] uses the volume fraction function α instead of the level set function. The volume fraction α takes the value of one in the dispersed phase and zero in the continuous phase. In addition, α is between zero and one in the cells that contain the interface. The VOF has a challenge in calculating the curvature k because the volume fraction function is discontinuous. On the other hand, the VOF is a mass conserving which is not the case in the level set method. In order to get an accurate calculation of the two-phase flow, a coupled VOF and LS method is used. In other words, the mass conservation in the VOF is combined with calculating k and other physical properties using the LS method. The interDyMFoam solver in OpenFOAM for 3D is modified to work with 2D planar and axisymmetric flows [3, 4]. The interDyMFoam for 2D is achieved to use the coupled VOF and LS method instead of using the VOF method only.

2. Main result

The solver interDyMFoam solves a two-phase flow problem for velocity u, pressure p and interface location. The continuity, momentum, and volume fraction equations are respectively given by:

u = 0 ,
( ρ u ) t + ρ u u = - p + η γ ˙ + F σ + ρ g ,
α t + ( α u ) = 0 ,

where ρ is the density, η is the viscosity, γ˙=u+(u)T is the rate-of-strain tensor and g is the gravity. Moreover, Fσ=σknδ represents the surface force tension, where σ is the interfacial tension, n is the unit normal vector on the interface and δ is the Dirac delta function.

In the VOF, the density and viscosity are respectively calculated as ρα=αρ1+(1-α)ρ2 and ηα=αη1+(1-α)η2 where the subscript represents the fluid phase. In the coupled VOF and LS method, only Eq. (3) is solved to find the volume fraction function. In addition, the LS function is initialized from the VOF function α as, φ0=1.5Δx(α-0.5), where Δx is the mesh cell size. Note that the initial value is positive in the disperse phase and negative in the continuous phase. The next step is to solve the re-initialization equation to assign the LS function to a signed distance function:

φ τ = s i g n ( φ 0 ) ( 1 - | φ | ) ,

where τ is the artificial time. The LS function at the beginning of the iterations is φ(x,0)=φ0(x).

The unit normal and the curvature are respectively computed using the LS function as, n=φ/|φ| and k=n. Therefore, the surface tension force is computed using the LS method as, Fσ=σk(φ)nδ(φ), where δ(φ) is defined as:

δ φ = 0 ,                                                                                     φ > ϵ , 1 2 ϵ 1 + cos π φ ϵ ,     φ ϵ .  

Also, the density and viscosity are calculated using the LS function as, ρ(φ)=ρ2+(ρ1-ρ2)H(φ) and η(φ)=η2+(η1-η2)H(φ), where H is the Heaviside function and defined as:

H φ = 0 ,                                                                                                         φ < - ϵ , 1 2 1 + φ ϵ + 1 ϵ sin π φ ϵ ,     φ ϵ ,     1 ,                                                                                                         φ > ϵ .    

The above equations are solved using the finite volume method (FVM). In FVM the computational domain is divided into a finite number of control volumes. Then, the equations are integrated over a control volume to get integral form of the equations. The integral equations are discretized using different methods to get the linear system of algebraic equations. In OpenFOAM, the dynamic mesh refinement works only for hexahedral cells. The boundary conditions for the LS equation are zeroGradient Which means that the derivative normal to the boundary is zero. For the other equations, there are several boundary conditions as zeroGradient and fixed value. For more details see [5].

3. Numerical test

The bubble rising in a column liquid is studied to compare the two solvers. The computational domain is a channel of length 1 m and hight 2 m. The bubble is a circle of radius 0.25 m and center (0.5, 0.5). The bubbles at time t= 1 and t= 2 s are shown in Figs. 1-2, where the bubbles on the left result from the interDyMFoam solver while the bubbles on the right result from the coupled interDyMFoam solver. The interface is smeared in the interDyMFoam solver. However, the coupled interDyMFoam gives a sharper interface.

The velocity and pressure from coupled interDyMFoam for this 2D problem are compared along the horizontal line y= 0 at time t= 1 s. Fig. 4 shows the velocity graphs for the same case with three different meshes: coarse, standard, and fine mesh using the coupled interDyMFoam solver. The solutions are almost identical outside the bubble and behave similarly inside the bubble. Also, the graph show a zero velocity outside the bubble which is expected since the top and bottom walls are fixed. The velocity decreased before the bubble, increased inside the bubble, and then decreased after the bubble. Note that the graphs from the standard and fine meshes are almost identical which gives the convergence. Similarly, the pressure graph, Fig. 3, shows the same convergence. They have the same behavior overall, as the jump in pressure is decreasing by increasing the mesh resolution using coupled interDyMFoam solver. However, the solutions from coupled interDyMFoam solver are accurate even on the coarse mesh.

Fig. 1. The bubble at time t= 1 s using a) interDyMFoam and b) coupled interDyMFoam

 The bubble at time t= 1 s using a) interDyMFoam and b) coupled interDyMFoam


 The bubble at time t= 1 s using a) interDyMFoam and b) coupled interDyMFoam


Fig. 2. The bubble at time t= 2 s using a) interDyMFoam and b) coupled interDyMFoam

 The bubble at time t= 2 s using a) interDyMFoam and b) coupled interDyMFoam


 The bubble at time t= 2 s using a) interDyMFoam and b) coupled interDyMFoam


Fig. 3. Pressure on three different meshes using coupled interDyMFoam solver

 Pressure on three different meshes using coupled interDyMFoam solver

Fig. 4. Velocity on three different meshes using coupled interDyMFoam solver

 Velocity on three different meshes using coupled interDyMFoam solver

4. Conclusions

The interDyMFoam solver is modified to allow of using VOF and LS methods. The VOF is a mass conserving and the LS function is continuous to calculate the curvature. The modified solver is tested on the bubble rising in a liquid column. Moreover, the two solvers are compared on the same case test. The modified solver gives a sharper interface while the original solver smears the interface. Also, the convergence is obtained by comparing the pressure and velocity on three different meshes.


Author thankful to the Deanship of Scientific Research and Graduate Studies, College of Science at Yarmouk University for funding to present this work.


  1. Osher S., Fedkiw R. Level Set Methods and Dynamic Implicit Surfaces. Springer-Verlag, New York, 2003. [Publisher]
  2. Deshpande S. S., Anumolu L., Trujillo M. F. Evaluating the performance of the two-phase flow solver interfoam. Computational Science and Discovery, Vol. 5, 2012, p. 014016. [Publisher]
  3. Baniabedalruhman A. Dynamic Meshing Around Fluid-Fluid Interfaces with Applications to Droplet Tracking in Contraction Geometries. Michigan Technological University, 2015. [CrossRef]
  4. Feigl Kathleen, Baniabedalruhman Ahmad, Tanner Franz X., Windhab Erich J. Numerical simulations of the breakup of emulsion droplets inside a spraying nozzle. Physics of Fluids, Vol. 28, 2016, p. 123103. [Publisher]
  5. Jasak H. Error Analysis and Estimation in the Finite Volume Method with Applications to Fluid Flows. Ph.D. Thesis, Imperial College, University of London, 1996. [CrossRef]
  6. Jadidi Mehdi, Tembely Moussa, Moghtadernejad Sara, Dolatabadi Ali Coupled level set and volume of fluid method in OpenFoam with application to compressible two-Phase flow. 22nd Annual Conference of the CFD Society of Canada, Toronto, Canada. [CrossRef]
  7. Shu Bitan, Dammel Frank, Stephan Peter Implementation of the level set method into OpenFOAM for capturing the free interface in incompressible fluid flows. OpenFOAM International Conference, Old Windsor, 2007. [CrossRef]