Difference between revisions of "Bubble growth detachment influence BC domain size by Michael Alletto"

From OpenFOAM Wiki
Jump to navigation Jump to search
 
(4 intermediate revisions by the same user not shown)
Line 203: Line 203:
 
The inlet value of α was set to 0 since air is released through the pipe at the bottom. For the pressure prgh the fixedFluxPressure boundary condition was used for the inlet. This boundary condition sets the pressure gradient such that the flux on the boundary is that specified by the velocity boundary condition (see https://www.openfoam.com/documentation/guides/latest/api/classFoam_1_1fixedFluxPressureFvPatchScalarField.html\#details).  
 
The inlet value of α was set to 0 since air is released through the pipe at the bottom. For the pressure prgh the fixedFluxPressure boundary condition was used for the inlet. This boundary condition sets the pressure gradient such that the flux on the boundary is that specified by the velocity boundary condition (see https://www.openfoam.com/documentation/guides/latest/api/classFoam_1_1fixedFluxPressureFvPatchScalarField.html\#details).  
  
Regarding the outlet a zero gradient boundary condition was used for α and a pressureInletOutletVelocity for the velocity (see https://www.openfoam.com/documentation/guides/latest/doc/guide-bcs-outlet-pressure-inlet-outlet.html). This boundary conditions set the velocity to zero gradient if the flow is pointing outside the domain and computes it according to the flux φ if it is pointing inside the domain. Since we do not know the value of the velocity in case the flow is pointing inside the domain, we prefer this boundary to the inletOutlet boundary condition where one has to explicitly specify the inlet value.  
+
Regarding the outlet a zero gradient boundary condition was used for α and a pressureInletOutletVelocity for the velocity (see https://www.openfoam.com/documentation/guides/latest/doc/guide-bcs-outlet-pressure-inlet-outlet.html). This boundary conditions set the velocity to zero gradient if the flow is pointing outside the domain and computes it according to the flux φ if it is pointing inside the domain. Additionally also a inletOutlet boundary condition is analyzed as outflow condition for the velocity. The inlet value of the velocity is set to a zero vector (0 0 0).  
For the pressure prgh tree types of boundary conditions are testes: prghPressure, prghTotalPressure and fixedValue.  
+
For the pressure prgh three types of boundary conditions are testes: prghPressure, prghTotalPressure and fixedValue.  
 
The prghPressure condition provides a static pressure condition for prgh. It means that the static pressure is specified and prgh is calculated accordingly. For details see https://www.openfoam.com/documentation/guides/latest/api/classFoam_1_1prghPressureFvPatchScalarField.html\#details.
 
The prghPressure condition provides a static pressure condition for prgh. It means that the static pressure is specified and prgh is calculated accordingly. For details see https://www.openfoam.com/documentation/guides/latest/api/classFoam_1_1prghPressureFvPatchScalarField.html\#details.
 
The prghTotalPressure conditions is similar to the prghPressure. Also here one specifies the static pressure a the patch. The prgh value for this condition is however reduced by the dynamic pressure 0.5 &rho; u<sup>2</sup> compared to the prghPressure condition.  
 
The prghTotalPressure conditions is similar to the prghPressure. Also here one specifies the static pressure a the patch. The prgh value for this condition is however reduced by the dynamic pressure 0.5 &rho; u<sup>2</sup> compared to the prghPressure condition.  
Line 298: Line 298:
  
  
[[File:Table_bubble_growth_detachment_bc_ds.png|450px|center|cases simulated]]
+
[[File:TableNew_bubble_growth_detachment_bc_ds.png|700px|center|cases simulated]]
  
  
Line 306: Line 306:
  
  
In the next figure the temporal evolution of the center of gravity in z-direction of the air is shown. It is evident that for the two cases with a no slip side boundary where the bubble leaves the domain (labeled as case B and case C), the center of gravity in z-direction suddenly decreases when the bubble is advected through the outlet. For case A where a fixedValue is used for prgh at the outlet, the center of gravity in z-direction of the air phase remains constant when the air bubble reaches the top of the domain. Regarding the detachment time of the air bubble similar values are obtained for this three cases. The time when the bubble detaches from the bottom surface is clearly visible by the steep gradient of the center of gravity of the air phase.  
+
In the next figure the temporal evolution of the center of gravity in z-direction of the air is shown. It is evident that for the tree cases with a no slip side boundary where the bubble leaves the domain (labeled as case B, case C and case F), the center of gravity in z-direction suddenly decreases when the bubble is advected through the outlet. For case A where a fixedValue is used for prgh at the outlet, the center of gravity in z-direction of the air phase remains constant when the air bubble reaches the top of the domain. Regarding the detachment time of the air bubble similar values are obtained for this four cases. The time when the bubble detaches from the bottom surface is clearly visible by the steep gradient of the center of gravity of the air phase.  
For the two cases where we varied the pressure condition at the side boundary (case D and case E), using a total pressure condition leads to a similar solution with the cases where a no slip wall was used. For a fixedValue condition for prgh at the side wall we can observe an increased detachment frequency after the air bubble detaches the first time from the bottom surface. The predicted detachment time is too small for all simulations. We will see in the next tutorials if we can improve the results and if possible analyze the reason of the discrepancies between simulation and experiment. A first cautions guess of the reason of the discrepancies between simulation and experiment is that the curvature computation using the gradient of the volume fraction &alpha; is only zero-order accurate (see [7]). For this reason it is the preponderating error in the simulations where the surface tension place a crucial role effecting the flow behavior and the curvature is computed by the gradient of the volume fraction.
+
For the two cases where we varied the pressure condition at the side boundary (case D and case E), using a total pressure condition leads to a similar solution with the cases where a no slip wall was used. For a fixedValue condition for prgh at the side wall we can observe an increased detachment frequency after the air bubble detaches the first time from the bottom surface. The predicted detachment time is too small for all simulations. The best agreement with the experiments could be achieved when an inletOutlet boundary condition is applied to the velocity at the outlet where the latest detachment is obtained. We will see in the next tutorials if we can improve the results and if possible analyze the reason of the discrepancies between simulation and experiment. A first cautions guess of the reason of the discrepancies between simulation and experiment is that the curvature computation using the gradient of the volume fraction &alpha; is only zero-order accurate (see [7]). For this reason it is the preponderating error in the simulations where the surface tension place a crucial role effecting the flow behavior and the curvature is computed by the gradient of the volume fraction.
  
  
 
[[File:cOGa_varyBC_growth_detachment_bc_ds_alletto.png|450px|center|evolution of the center of gravity of the air phase]]
 
[[File:cOGa_varyBC_growth_detachment_bc_ds_alletto.png|450px|center|evolution of the center of gravity of the air phase]]
 
  
 
===Varying the domain size===
 
===Varying the domain size===

Latest revision as of 14:22, 29 April 2021

Go back to Collection by authors.


Bubble growth detachment influence BC domain

You can download the case files [1] here.

Introduction

In a resent OpenFOAM version (see https://www.openfoam.com/releases/openfoam-v1706/numerics.php) a new interphase advection scheme was released. It is a geometric volume of fluid method. Geometric volume of fluid method means that a geometric model for the interphase position and orientation within a cell is established and advected within a time step. One can imagine that the interphase between two fluids within a cell is described as a plane which has a given position of the plane center and orientation of the normal vector. Two different method to reconstruct this plane are currently available within interIsoFoam: isoAlpha and plicRDF. isoAlpha calculates the surface dividing the two fluids based on the α iso-value which cuts the cell according to the volume occupied by the two fluids. For details see [6,7]. plicRDF uses a piece wise linear radial distance function to compute the orientation and position of the interphase within a cell. The gradient of the distance function gives the normal of the interphase. By moving the interphase position within the cell the volume fraction of the two fluids divided by the interface can be adjusted to the α value. For details see [6]. This reconstruction method was released with OF2006 (see https://www.openfoam.com/releases/openfoam-v2006/solver-and-physics.php).

In a resent paper [2] performed an extensive comparison and validation of the two reconstruction methods of interIsoFoam. They compared interIsoFoam with interFoam and another volume of fluid code (Basilisk). The main findings of the paper were that: (i) plicRDF causes the smallest parasitic currents in a static bubble, (ii) regarding a rising 2D bubble plicRDF delivers slightly better results than MULES with a bit smaller computational times and (iii) plicRDF was the only OpenFOAM version capable to predict the spiraling movement on a coarse grid of a 3D rising bubble. Note that the test cases of the paper can be downloaded here: https://wiki.openfoam.com/Collection_by_authors\#Lionel_Gamet.

Inspired by the novelty of the plicRDF interface reconstruction method and the promising results of [2], the purpose of this tutorial is to analyze the solver interIsoFoam with the plicRDF reconstruction method a bit more in detail. For this purpose the experiment described in [1] is taken as reference. [1] measured the shape of a axial symmetric air bubble exiting a circular orifice. The orifice was submerged by a liquid. This kind of test case is often used as preliminary test case for solvers which are meant to be used later to simulate boiling (see [3,4]).

In this first tutorial of a series we will analyze the influence of the boundary conditions and domain size on the flow and detachment time of a rising bubble. The next tutorials will have a look at the influence of the grid size, time step size, number of inner and outer iterations, settings of the parameter influencing the interphase advection scheme and finally also the way the curvature of the interphase is computed.

Setup of the experiment

In the next figure the experimental setup used in [1] is sketched. The experimental facility comprises a horizontal surface from which the bubble growth takes place. This surface is submerged in a 20 mm water column within a square 50 mm glass tank. The growth orifice has a diameter of 1.6 mm, is manufactured from stainless steel and is inserted into the submerged surface. With a high speed camera [1] measured the evolution and detachment of the gas bubble released from the orifice. The measured detachment time was 0.684 s. This is our main quantity of interest against we compare our simulations.

setup of the experiment


Having analyzed the experimental setup, we see that we have clear defined physical boundaries: Air enters the domain by a circular pipe. From the description of the experiment it is known that the velocity profile was parabolic like in a fully developed laminar pipe flow. The water was placed in a quadratic glass tank: From this we known that the setup was confined by walls. The height of the liquid column is also known. From the description of the experiment we also know that the rising bubble was rotational symmetric. Unfortunately the tank wall are not rotational symmetric (they are quadratic actually). For this reason the problem is not completely rotational symmetric.

Knowing the physical boundaries of the experiment we want to reproduce, we have to decided the boundary of our numerical simulation. Being this a tutorial, we want to keep the computational afford small. For this reason, we want to exploit the rotational symmetric bubble shape and decide to use a rotational symmetric grid. With this approach we can perform a 2D simulation which reduces the computational time considerably.

Furthermore we decide that the height of the domain should be equal to the liquid column since at the upper boundary the boundary condition is known. The last decision which has to be taken is the lateral extent of the computational domain. If we decide to cut the computational domain before the lateral wall, we need to carefully choose the type of boundary condition which mimics best the real condition at this artificial boundary.


Setup

The next figure shows the setup chosen for our simulation. It consists in a wedge type computational domain (it is how a rotational symmetric problem is solved in OpenFOAM). The z-coordinates points against the direction of gravity and the x-coordinate points in radial direction. At the lateral boundary we set an empty boundary condition. At the axis we set an axis boundary condition. At the inflow we will set a parabolic velocity profile to be in line with the experiment. The bottom wall boundary is also fixed. Since the top boundary conditions was not trivial to set (at least for me) we will vary it in order to see its influence on the simulation results. The boundary conditions at the outer radius is a artificial boundary conditions: Since the experimental setup was not rotational symmetric (the squared tank in which the bubble is released) we artificially restrict the possible solution by using an axis symmetric configuration. For this reason we will also vary the type of this boundary condition and its axial position to see its influence on the solution.


setup of simulation

Mesh generation

The next listing shows the blockMeshDict. The main parts are taken from here: https://openfoamwiki.net/index.php/Main_ContribExamples/AxiSymmetric. Note that the line mergeType points should be added in the file. Otherwise zero area faces are created at the axis which cause errors in the checkMesh utility. The extension of the domain in z and x-direction can be specified by the entries Lx and Lz. The position of the points delimiting the domain boundary are calculated via the eval functionality of OpenFOAM. The number of points in x- and z-direction a specified by the variables nx and nz. The points are distributed with equal spacing.


/*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: v2012 | | \\ / A nd | Website: www.openfoam.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile {

   version     2.0;
   format      ascii;
   class       dictionary;
   object      blockMeshDict;

} // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

scale 0.001;

mergeType points; // Wedge geometry - Merge points instead of topology

Lx 10.0; Lz 20.0;

nx 50; ny 1;

nz 100;

Px #eval{ $Lx * cos(pi()*2.5/180) }; Py #eval{ $Lx * sin(pi()*2.5/180) }; MPy #eval{ -$Lx * sin(pi()*2.5/180) };


vertices ( (0 0 0) ($Px $Py 0) ($Px $Py $Lz) (0 0 $Lz) ($Px $MPy 0) ($Px $MPy $Lz) );


blocks (

   hex (0 4 1 0 3 5 2 3) ($nx  1 $nz) simpleGrading (1 1 1)

);

edges ( );

boundary (

    front
    { 
          type wedge;
          faces  
          (
              (0 1 2 3)
          );
     }
    back
    { 
          type wedge;
          faces  
          (
              (0 3 5 4)
          );
     }
    tankWall
    { 
          type wall;
          faces  
          (
              (1 4 5 2)
          );
     }
    wallBottom
    { 
          type patch;
          faces  
          (
              (0 4 1 0)
          );
     }
    outlet
    { 
          type patch;
          faces  
          (
              (3 5 2 3)
          );
     }
    axis
    { 
          type empty;
          faces  
          (
              (0 3 3 0)
          );
     }

);

mergePatchPairs ( );

// ************************************************************************* //


Simulation setup

The computational setup was inspired by the setup used by [1,4] which simulated the same case as examined in this tutorial.

For the inlet we prescribed a parabolic velocity profile which is in agreement with the specified mass flow (see [4]). In order to achieve a parabolic inlet profile, a codedFixedValue boundary condition was used. The next listing shows the applied inlet boundary condition for the velocity.

inlet

   {
       name    parabolicInlet;
       type    codedFixedValue;
       code
       #{
           // Uncomment for testing on non-windows systems [fragile]
           // IOobject::scopeSeparator = '_';
           const scalar Q(4.17e-08);
           const scalar maxR(gMax(this->patch().Cf().component(0)));
           scalar vMax;
           vMax = 2.*Q/(M_PI*pow(maxR,2));
           scalarField vz(vMax*(1. -  pow(this->patch().Cf().component(0),2)/pow(maxR,2)));


           vectorField v(this->patch().size(), vector(0, 0, 0));
           v.replace(vector::Z, vz);
           operator==(v);
       #};
       value           $internalField;


   }

The inlet value of α was set to 0 since air is released through the pipe at the bottom. For the pressure prgh the fixedFluxPressure boundary condition was used for the inlet. This boundary condition sets the pressure gradient such that the flux on the boundary is that specified by the velocity boundary condition (see https://www.openfoam.com/documentation/guides/latest/api/classFoam_1_1fixedFluxPressureFvPatchScalarField.html\#details).

Regarding the outlet a zero gradient boundary condition was used for α and a pressureInletOutletVelocity for the velocity (see https://www.openfoam.com/documentation/guides/latest/doc/guide-bcs-outlet-pressure-inlet-outlet.html). This boundary conditions set the velocity to zero gradient if the flow is pointing outside the domain and computes it according to the flux φ if it is pointing inside the domain. Additionally also a inletOutlet boundary condition is analyzed as outflow condition for the velocity. The inlet value of the velocity is set to a zero vector (0 0 0). For the pressure prgh three types of boundary conditions are testes: prghPressure, prghTotalPressure and fixedValue. The prghPressure condition provides a static pressure condition for prgh. It means that the static pressure is specified and prgh is calculated accordingly. For details see https://www.openfoam.com/documentation/guides/latest/api/classFoam_1_1prghPressureFvPatchScalarField.html\#details. The prghTotalPressure conditions is similar to the prghPressure. Also here one specifies the static pressure a the patch. The prgh value for this condition is however reduced by the dynamic pressure 0.5 ρ u2 compared to the prghPressure condition. For details see https://www.openfoam.com/documentation/guides/latest/api/classFoam_1_1prghTotalPressureFvPatchScalarField.html\#details

For one of the side boundary condition analyzed we used a no slip condition for the velocity and a fixedFluxPressure for prgh. Even a no slip condition is a quite crude assumption of an artificial boundary in general, this decision was inspired by the simulations of [1,4]. The velocity induced in the domain by the rising bubble should decreases with the distance to the bubble. If the side boundary of the domain is far enough from the axis, the induced velocity should decrease nearly to zero and a no slip condition is an approximation which is not in big contradiction with the expected flow picture at this boundary. Since applying a no slip boundary condition at an artificial boundary is a condition which can produce a solution which is quite far away from the expected solution, we further analyze a fixed value condition for the pressure at the side boundary and also a total pressure boundary condition. The velocity for this two pressure condition was set to pressureInletOutletVelocity in order to allow the velocity to adjust to the pressure.

Regarding the bottom wall a no slip boundary conditon was used for the velocity, a fixedFluxPressure for prgh and we specified a static contact angle of 20° for α. The same value for the static contact angle was used by [1,4].

The next listing shows the coded function object which computed the geometric center in z-direction of the air phase and writes this value in a text file. The main purpose of this function object is to detect the detachment time of the air bubble: When the bubble detaches from the bottom we will see a strong gradient of the z-position of the geometric center of the air bubble.

   centerOfGravityAir
   {
       libs        ("libutilityFunctionObjects.so");
       type coded;
       // Name of on-the-fly generated functionObject
       name writeCenterOfGravityAir;
       writeControl  timeStep;
       
       executeControl timeStep;
       
       codeOptions
       #{
           -I$(LIB_SRC)/finiteVolume/lnInclude \
           -I$(LIB_SRC)/OpenFOAM/lnInclude
       #};


       codeInclude
       #{
           #include "volFieldsFwd.H"
           #include "OFstream.H"
       #};
       
       codeData
       #{
             autoPtr<OFstream> outputFilePtr;
       #};
           
       codeRead
       #{
          
            outputFilePtr.reset(new OFstream("centerOfGravityAir.dat"));
          
       #};
       
       codeExecute
       #{
       
            const volScalarField& alpha1 = mesh().lookupObject<volScalarField>("alpha.water");
            
            scalar alphaTot (0.0);
            scalar cOG (0.0);
            
            
            forAll(alpha1,cellI)
            {
                alphaTot += 1.0 - alpha1[cellI];
                cOG += (1.0 - alpha1[cellI])*alpha1.mesh().C()[cellI].z();
            }
            
            reduce(alphaTot, sumOp<scalar>());
            reduce(cOG, sumOp<scalar>());
            
            cOG /= (alphaTot + VSMALL);
           
            if (Pstream::myProcNo() == 0)
            {
                outputFilePtr() << mesh().time().timeName() << " " << cOG << endl;
            }
            
            Info << "cog = "  << mesh().time().timeName() << " " << cOG << endl;
            
       #};
       
       codeWrite
       #{
       #};
   }

Results

Varying the upper and outer axial boundary condition

In this section the results obtained by changing the pressure boundary condition at the top and the boundary condtions at the outer axial boundary are shown. This study was performed since finding the correct conditions was not trivial (at least for me).

In order to provide an overview of the simulations performed, the next table shows the combination of boundary conditions used at the top boundary and at the outer axial boundary. The denomination of the cases is also displayed.


cases simulated


In the next figure snapshots of the α values are shown for different time steps. For this simulations a no-slip wall at the outer radial boundary was used. The left column shows the results obtained by applying a prghPressure boundary condition at the top and the right column shows the results when applying a fixed value at the top for prgh. Note that for a prghTotalPressure similar results to a prghPressure boundary condition is obtained. It is evident that the air bubble exists the pipe located at the bottom surface and slowly stretches with time. When the surface tension is not strong enough to hold the bubble it brakes and moves towards the exit. This qualitative picture is the same for all boundary condition applied at the top. The main difference between the simulations where different boundary conditions for prgh at the top are applied, is the latest stage of the simulation: for a fixedValue the bubble does not leave the domain (see bottom most pictures) and for a prghPressure the bubble leaves the domain. Thanks to Henning Scheufler who suggested to change the type of boundary condition in order to allow the bubble to leave the domain (see https://develop.openfoam.com/Development/openfoam/-/issues/2016).

alpha different time steps


In the next figure the temporal evolution of the center of gravity in z-direction of the air is shown. It is evident that for the tree cases with a no slip side boundary where the bubble leaves the domain (labeled as case B, case C and case F), the center of gravity in z-direction suddenly decreases when the bubble is advected through the outlet. For case A where a fixedValue is used for prgh at the outlet, the center of gravity in z-direction of the air phase remains constant when the air bubble reaches the top of the domain. Regarding the detachment time of the air bubble similar values are obtained for this four cases. The time when the bubble detaches from the bottom surface is clearly visible by the steep gradient of the center of gravity of the air phase. For the two cases where we varied the pressure condition at the side boundary (case D and case E), using a total pressure condition leads to a similar solution with the cases where a no slip wall was used. For a fixedValue condition for prgh at the side wall we can observe an increased detachment frequency after the air bubble detaches the first time from the bottom surface. The predicted detachment time is too small for all simulations. The best agreement with the experiments could be achieved when an inletOutlet boundary condition is applied to the velocity at the outlet where the latest detachment is obtained. We will see in the next tutorials if we can improve the results and if possible analyze the reason of the discrepancies between simulation and experiment. A first cautions guess of the reason of the discrepancies between simulation and experiment is that the curvature computation using the gradient of the volume fraction α is only zero-order accurate (see [7]). For this reason it is the preponderating error in the simulations where the surface tension place a crucial role effecting the flow behavior and the curvature is computed by the gradient of the volume fraction.


evolution of the center of gravity of the air phase

Varying the domain size

The next figure shows the influence of the domain size on the detachment time of the bubble. In the left figure a no slip wall for the velocity and a fixedFluxPressure condition for prg was applied at the outer radius of the domain. In the right figure a pressureInletOutletVelocity for the velocity and a prghTotalPressure condition for prg was applied at the outer radius of the domain. Note that the same resolution was applied for all simulations.

Regarding the simulations with a no-slip wall at the outer axial boundary we see that the detachment time decreases slightly with the size of the domain. The simulation where the outer axial boundary was placed at 20 mm and 25 mm from the center show nearly the same detachment time.

Regarding the simulation with a totalPressure condition at the outer axial boundary is evident that no obvious trend of the detachment time on the domain size can be observed. The detachment time is similar for all configurations simulated. This indicate that the size of the smallest domain is sufficient to have a solution which is fairly independent of the domain size.

evolution of the center of gravity of the air phase


Lessons learned

  • for a new case to solve a careful evaluation of the influence of the domain size and also of the applied boundary condition should be perform
  • applying a fixedValue for prgh at the outlet leads the air phase to remain inside the domain
  • applying a prghPressure condition ( which is equivalent to apply a fixed value for the static pressure) allows the bubble to leave the domain
  • the combination totalPressure for prgh at the top boundary and pressureInletOutletVelocity for U at the outer axial boundary gives similar results as the combination fixedFluxPressure for prgh at the top boundary and no-slip for U at the outer axial boundary
  • for a lateral domain size of roughly 10 bubble radii, the solution is roughly independent of the position of the boundary

References

[1] A Albadawi, DB Donoghue, AJ Robinson, DB Murray, and YMC De- lauré. On the analysis of bubble growth and detachment at low capillary and bond numbers using volume of fluid and level set methods. Chemical Engineering Science, 90:77–91, 2013.


[2] Lionel Gamet, Marco Scala, Johan Roenby, Henning Scheufler, and Jean-Lou Pierson. Validation of volume-of-fluid openfoam R isoad- vector solvers using single bubble benchmarks. Computers & Fluids, 213:104722, 2020.


[3] Anastasios Georgoulas, Manolia Andredaki, and Marco Marengo. An enhanced vof method coupled with heat transfer and phase change to characterise bubble detachment in saturated pool boiling. Energies, 10(3):272, 2017.


[4] Anastasios Georgoulas, P Koukouvinis, Manolis Gavaises, and Marco Marengo. Numerical investigation of quasi-static bubble growth and de- tachment from submerged orifices in isothermal liquid pools: The effect of varying fluid properties and gravity levels. International Journal of Multiphase Flow, 74:59–78, 2015.


[5] Johan Roenby, Henrik Bredmose, and Hrvoje Jasak. A computational method for sharp interface advection. Royal Society open science, 3(11):160405, 2016.

[6] Henning Scheufler and Johan Roenby. Accurate and efficient surface reconstruction from volume fraction data on general meshes. Journal of computational physics, 383:1–23, 2019.


[7] Henning Scheufler and Johan Roenby. Twophaseflow: An openfoam based framework for development of two phase flow solvers. arXiv preprint arXiv:2103.00870, 2021.