Extending the isoAdvector Geometric VOF Method to Flows in Porous Media

. We consider the interfacial flow in and around porous structures in coastal and marine engineering. ⋆ During recent years, interfacial flow through porous media has been repeatedly simulated with Computational Fluid Dynamics (CFD) based on algebraic Volume Of Fluid (VOF) methods [1,2]. Here, we present an implementation of a porous medium interfacial flow solver based on the geometric VOF method, isoAdvector [3,4]. In our implementation, the porous medium is treated without resolving the actual pore geometry. Rather, the porous media, pores, and rigid structure are considered a continuum and the effects of porosity on the fluid flow are modelled through source terms in the Navier-Stokes equations, including Darcy-Forchheimer forces, added mass force and accounting for the part of mesh cells that are occupied by the solid material comprising the skeleton of the porous medium. The governing equations are adopted from the formulation by Jensen et al. [1]. For the interface advection using isoAdvector, we also account for the reduced cell volume available for fluid flow and for the increase in the interface front velocity caused by a cell being partially filled with solid material. The solver is implemented in the open source CFD library OpenFOAM ® . It is validated using two case setups: 1) A pure passive advection test case to compare the isolated advection algorithm against a known analytical solution and 2) a porous dam break case by Liu et al. [5] where both numerical and experimental results are available for comparison. We find good agreement with numerical and experimental results. For both cases the interface sharpness, shape conservation as well as volume conservation and boundedness are demonstrated to be very good. The solver is released as open source for the benefit of the coastal and marine CFD community (code repository https://github.com/InterFlowers/porousInterIsoFoam) and as of OpenFOAM-v2112 the new functionality is integrated in the official interIsoFoam solver.


Introduction
Interfacial flow through porous media appears in numerous engineering problems.For instance coastal and marine as well as environmental applications such as waste management facilities.In this work we focus on coastal and marine applications; typical examples from this field are breakwaters, seawalls and other perforated marine structures such as aquaculture cages [7][8][9].Given the complexity of these flows, analytical flow solutions are intractable.The predictive capabilities of experimental model tests, on the other hand, are limited by scaling effects.Numerical simulation provides an alternative engineering tool from which important flow features can be extracted such as pressure, velocity distribution and surface elevation.
Fully resolving the flow in the complex fluid domain inside the porous structure is often too computationally expensive.A more convenient modeling strategy is to treat the porous medium in the volume averaged spirit of the Finite Volume Method (FVM).The first systematic study within an FVM framework was carried out by van Gent [10], who utilized the Volume Averaged Reynolds Averaged Navier Stokes (VARANS) and added explicit source terms to the equations to account for the resistance forces exerted on the fluid by the porous material.Those explicit forces were based on the Darcy-Forchheimer equations.Later Liu et al. [5] conducted experiments in order to validate van Gent's work.More recently, Higuera et al. [2] applied van Gent's approach in their implementation, denoted IHFOAM, of an intrinsic VARANS set of equations along with a two-equation turbulence closure model to take into account the turbulent stresses.For interface advection they used a Multidimensional Universal Limiter with Explicit Solution (MULES) scheme, a description of which is available from Deshpande et al. [11].Similarly, Jensen et al. [1] implemented the porous flow equations in the waves2Foam package [12] also using MULES for interface advection.Their implementation contained a corrected interface advection taking properly into account the reduced cell volume available for fluid flow when it is partially filled with porous material.Turbulence was not modelled and all turbulent effects were assumed incorporated in the body force in the momentum equation.An extensive validation of this implementation was carried out by Jacobsen et al. [13].A recent study performed by Feichtner et al. explores the use of the VARANS approach for simulations of wave interaction with thin perforated structures [14].
In contrast to the previous studies that used algebraic VOF for the interface advection part, in this paper we employ a geometric VOF method called isoAdvector.The algorithm performs very well in terms of advection accuracy, interface sharpness, volume conservation, and boundedness [3].The method is well-documented and tested in multiple areas of application [15][16][17][18].Especially the sharp interface attributes of isoAdvector is a desirable feature when simulating waves interacting with porous structures in coastal and marine engineering.Our newly developed solver, porousInterIsoFoam, is essentially an extension of the existing interfacial flow solver, interIsoFoam, in OpenFOAM enabling it to model flows in and around porous structures.
In the remaining manuscript we briefly describe the theoretical background for the implemented solver.We then present two validation cases using porousInterIsoFoam for 1) pure passive advection of a disc through a porous region and 2) reproducing the porous dam break experiment conducted by Liu et al. [5].In both cases we compare against porousInterFoam, an extension of the interfacial flow solver interFoam that can model flows in porous media.The porousInterFoam solver is a refactored version of the porousWaveFoam solver (part of the waves2Foam package).To ensure fair comparison, the porous-InterFoam and porousInterIsoFoam solver only differ in the interface advection where the former uses MULES and the latter uses isoAdvector.For implementation details, please see the accompanying code repository (https://github.com/InterFlowers/porousInterIsoFoam).

The governing equations
In this section we present the porous interfacial flow equations implemented in our new porousInter-IsoFoam solver.
2.1.The VARANS equations.For the treatment of the momentum equation, we follow the derivation of Jensen et al. [1].By applying a superficial volume average1 to the Reynolds Averaged Navier-Stokes (RANS) equations, the VARANS equations are derived.In these equations, the effect of the porous region on the flow is included based on the Darcy-Forchheimer equation via the linear and non-linear resistance forces and an added mass force proportional to the fluid acceleration.The resistance coefficients are determined according to van Gent [10] and can also account for turbulent effects inside the porous region.
The volume averaged continuity equation is given as where u is the velocity field, the overbar represents ensemble averaging and the angle brackets represent the superficial volume average.The volume-averaged, Reynolds averaged momentum equation can be formulated as Here the added mass coefficient, C m is modelled as C m = 0.34(1 − n)/n, where n is the effective porosity field.This is given by n = V f /V , where V f is the void volume within a cell and V is the total volume of the cell.ρ is the fluid density, ⟨pf ⟩ is the intrinsic volume averaged pressure, g and x are the gravity and position vector, and µ is the dynamic viscosity.The last term, F contains the combined resistance force term exerted by the pore skeleton on the fluid.It is modeled as where a and b are the resistance coefficients determined by van Gent [10].
In the present study we carry out laminar simulations and therefore we do not consider a turbulence closure.Furthermore, the surface tension term is disregarded since the capillary effects are assumed negligible in the cases of interest.This remains a valid approximation for many marine and coastal applications which is the focus of the present implementation.The reader is referred to Jensen et al. [1] where a more exhaustive account of the governing equations is given.

2.2.
Interface advection in porous media.The evolution of the fluid interface can be described by a cell volume integrated form of the continuity equation, where dS is the differential area vector pointing out of the cell volume V while ∂V is the surface enclosing volume V .When employing a VOF method we define the volume fraction of phase-1 fluid within cell C i as where V i is the volume of cell i and H is a three-dimensional Heaviside function equal to one inside the phase-1 region and zero elsewhere.The outcome of this definition is a volume fraction field taking the value α i = 0 in cells filled with phase-2, α i = 1 in cells filled with phase-1, and α i ∈ (0, 1) when an interface is present in the cell.By time integrating the discretized form of Eqn. ( 4) over the interval [t, t + ∆t] one can arrive at where F i the set of faces of cell i.The double integral on the RHS of Eqn. ( 6) is the outward advected phase-1 fluid from cell i through its face f during the time interval [t, t + ∆t].In practice, isoAdvector evaluates this integral, assuming a temporally constant advecting face flux, ϕ f (t), on face f , yielding the following approximation, where S f is the face area and A f (τ ) is the submerged area of face f at time τ .To calculate this area, it is important to know how the interface moves inside a cell during the interval [t, t + ∆t] and hence how it sweeps the face of interest, f .The interface inside a cell is represented by an approximately planar polygonal face with a well-defined center x I i and a unit normal vector n I i .The interface advection velocity U I i is approximated by interpolating the cell centred velocity field to x I i .The aforementioned setup results in a unique description of how the interface travels within a cell and how it sweeps a given face f .Based on that model A f (τ ) is given as a second-degree polynomial in time that can be further integrated as per Eqn.(7).For the sake of clarity, further details of the isoAdvector method are left aside.The reader is referred to Roenby et al. [3] and Scheufler and Roenby [4] where the isoAdvector algorithm is described in detail.
In order to extend isoAdvector to flows in porous media we must make a number of modifications to the existing algorithm.First, the available void volume within cell i is only a fraction n i of the total cell volume.As a result, the volume fraction definition in Eqn. ( 5) should be altered to In this way the volume fraction field retains the property to be in the range 0 to 1 with α i = 1 to refer to a cell with its void volume filled with phase-1 fluid and α i = 0 to a cell with its void volume filled with phase-2 fluid.Subsequently the above volume fraction modification directly affects Eqn.(6), which in order to accommodate the presence of a porous medium becomes The isoAdvector algorithm evaluates the second term on the RHS of Eqn. ( 9) based on geometric operations.Those operations require an estimate for the advection velocity vector U I i at interface centre x I i .As mentioned, we interpolate the cell centred velocity field to x I i to approximate U I i .However, when a porous medium is present, the advection velocity must be adjusted accordingly.This is because the interface in a cell is a Lagrangian surface and thus it will be advected with the superficial velocity, while the cell centred velocity vectors are an intrinsic representation of the velocity field.As a result, the interpolated velocity at the interface center x I i must be divided by the porosity value of the cell.Therefore the adjusted interface advection velocity in cell i will be U I′ i = U I i /n i .
Figure 1.Sketch of the set up for the disc in constant flow through porous region case.
The disc is shown in grey with a dashed outline, while the porous region located in the middle of the figure is depicted in a grey-white pattern.U refers to the constant velocity vector, L is the total length of the domain, H is the total height of the domain, and L p is the length of the porous region in the horizontal direction.

Benchmark cases
Here we present validation of the porousInterIsoFoam solver in two different cases.Initially, a pure passive advection case is simulated, and afterwards, we reproduce the experimental results of Liu et al. [5].These simulations aim to evaluate/illustrate some critical aspects of this new solver, such as interface sharpness, shape preservation, boundedness, and volume conservation.
3.1.Passive Advection of a Disc Through a Porous Region.In this section, the passive advection of a disc through a uniform porous region is simulated.The term passive implies that the flow field is prescribed and unaltered by the presence of the disc throughout the simulation.As a result, this case will allow us to have a closer look at the interface advection algorithm isolated from the complex porous flow phenomena.As depicted in Fig. 1  with porosity n = 0.5.The circular blob with center (x 0 , y 0 ) = (0.5, 0.5) [m] and radius R = 0.25[m] is advected diagonally with a velocity U = (1, 0.5)[m s −1 ] outside the porous region.Inside the porous region, fluid incompressibility dictates a velocity U ′ = U/n, hence a doubling of the interface velocity with the chosen porosity.The result is that the circular blob is stretched as it passes the left boundary into the porous region.The blob regains its circular shape as it leaves the porous region through its right boundary.If we call the left and right limits of the porous domain x L and x R , respectively, then the analytical trajectories for the particles comprising the circular blob can be written We stop the simulation at t = 2.5[s], leaving the disc time to enter, pass through, and exit the porous region.Lastly, the time step selection is adjustable in order to keep the maximum Courant number, maxCo, at 0.25.The Courant number of a cell i with porosity n i is defined as which is identical to the existing definition in OpenFOAM when n i = 1.The porous modification retains the interpretation of the Courant number as the fraction of the total fluid volume in the cell (n i V i ) that is fluxed out of the cell during the time step, ∆t.We note that this porous Courant number definition poses a more strict limitation on the time step size than the non-porous version.In the simulation we also restrict the maximum interfacial Courant number, maxAlphaCo, to 0.25.In accordance with the existing maxAlphaCo implementation in OpenFOAM, the maximum is taken over all cells with 0.01 < α < 0.99.The different stages of the disc during the simulation are illustrated in Fig. 2. As the disc enters the porous region, it elongates in the advection direction, and its advection velocity doubles.When exiting the porous region, the elongated elliptic blob transforms back to its original shape.The interface sharpness is illustrated with the α = 0.01 and α = 0.99 contours.Those contours should be as close as possible to each other with the given mesh resolution.In the porousInterIsoFoam case the contours are observed to be three cells apart during the whole course of the simulation.For porousInterFoam the interface width is approximately 5 cells.The green closed curve in each panel on the figure shows the theoretical interface shape from Eqn. (10).porousInterIsoFoam recovers an almost circular interface shape upon its exit from the porous region whereas porousInterFoam returns a distorted shape.
In Fig. 3a a plot of the volume conservation error is shown.For both solvers the error is observed to be around machine precision with porousInterFoam exhibiting a lower error compared to porousInterIso-Foam.In Fig. 3b the upper bounding error is depicted for both solvers.Here the porousInterIsoFoam error is around machine precision, whereas the porousInterFoam maximum α value settles around 10 −7 below 1.This is due to the smear of the volume fraction field that subsequently leads to a maximum value below one.In Fig. 3c we show the lower bounding error, which stays below machine precision during the whole simulation for both solvers.

Porous Dam Break.
In the second benchmark case we reproduce an experiment of a porous dam break conducted by Liu et al. [5].This particular experiment has been a reference point for many past developments in the field (e.g.[1], [2] and [20]).Here we evaluate the solver behaviour as a whole including the combined interface advection and porous momentum equation.Fig. 4 illustrates the simulation set up.On the left there is a column of water that is released at t = 0[s].Then water flows through the uniform porous block placed in the middle of the domain, until at some point the system finds equilibrium.The full description of the experimental setup can be found in Liu et al. [5].The domain is discretised by a uniform mesh with cell size of ∆x = ∆y ≈ 0.006[m].This simulation is carried out with a maximum interfacial Courant number of 0.5.
The water surface profiles for different time instances are illustrated in Fig. 5.It is clear that porous-InterIsoFoam predicts the water elevation accurately in most of the simulation.There is a deviation in the initial phase, this is presumably due to modelling of the water gate; in the simulation the whole water  column is released instantly while in reality it is a finite process that takes around 0.1[s], as also noted by Liu et al. [5] (similar observations have been made in a different case study by Fekken et al. [21]).Fig. 6 illustrates the interface sharpness capabilities of isoAdvector compared to MULES.For the case of isoAdvector the α = 0.01 and α = 0.99 contours are shown to be separated by three cells which is the best a VOF representation can do.On the other hand, porousInterFoam results in a smeared interface that expands over more cells, while at the same time producing numerical vapor.Both the numerical vapor issue and the smeared interface of MULES has been reported in earlier investigations [22,23].
In Fig. 7a we show the evolution in the volume conservation error over the course of the simulation (i.e.4[s]).Here it is observed that both solvers exhibit excellent volume conservation and retain the error below 10 −13 .Fig. 7b shows that the upper bounding error is kept below 10 −9 for both solvers.As for the lower bounding error Fig. 7c shows that both solvers retain the error below machine precision.

Conclusion
An extension of the interfacial solver interIsoFoam has been derived and implemented.The extended solver, porousInterIsoFoam, is capable of simulating two-phase flows inside and around porous regions.The momentum equation implementation is based on the analysis of Jensen et al. [1], while the geometric VOF method, isoAdvector, is extended to account for a porous medium occupying part of the volume of computational cells.Using isoAdvector for interface advection gives rise to improved shape preservation while retaining excellent volume conservation and boundedness.The solver was validated with a pure advection test case demonstrating the superior shape preservation compared to a similar solver using MULES for interface advection.The porousInterIsoFoam solver was also benchmarked with numerical and experimental data in a porous dam break case where it accurately matched the interface shape measurements while also conserving fluid mass to a high precision with a sharp and well-bounded volume fraction field.The porousInterIsoFoam solver is released as an open source tool capable of facilitating coastal and marine studies of interfacial flows involving porous regions (code repository https://github.com/InterFlowers/porousInterIsoFoam).As of OpenFOAM-v2112 the porousInter-IsoFoam capabilites have been merged into the official interIsoFoam solver (the reader is referred to the v2112 release notes for further details on usage).
the computational set-up consists of a rectangular domain of length L = 5[m] and height H = 3[m].The domain is discretised by a uniform mesh with cell size ∆x = ∆y = 0.05[m].The rectangular domain contains an inner porous region of length L p = 2[m] Initial (b) Entry, t = 1[s] (c) Inside, t = 1.5[s](d) Exit, t = 2[s](e) End, t = 2.5[s]

Figure 2 .
Figure 2. Disc in uniform flow at various phases (Initial, during the entry, inside, during the exit and outside of the porous region).Comparison of the shape preservation and interface sharpness for the two solvers.porousInterIsoFoam results are shown in the top row of panels in grey-scale (α = 0.5 contour is illustrated in black while α = 0.01 and α = 0.99 contours are shown in grey) whereas porousInterFoam results are shown in the bottom row of panels in red-scale (α = 0.5 contour is illustrated in red while α = 0.01 and α = 0.99 contours are shown in light red).In green is illustrated the analytical solution from Eqn. (10).Finally, the limits of the porous region are shown as yellow vertical lines.

Figure 3 .
Figure 3.Comparison of volume conservation error as well as upper and lower bounding of the volume fraction for the two solvers.The two vertical dashed lines mark the time interval where the disc is in contact with the porous region.The vertical dash-doted lines mark the time interval where the disc is fully inside the porous region.

Figure 4 .
Figure 4. Sketch of the set-up for the porous dam break case.The free surface of the initial water distribution is shown in a dashed line, while the porous region located in the middle of the figure is depicted with a grey-white pattern.L = 0.892[m] is the total length of the domain, H c = 0.24[m] and L c = 0.28[m] is the height and width of the initial water column, and L p = 0.29[m] is the length of the porous region in the horizontal direction.Lastly, the height of the short column is h 0 = 0.022[m].

Figure 5 .
Figure 5.Comparison of free surface elevation profiles at 9 different time instances.The experimental data, shown in green dots, are measurements presented in Liu et al. [5] while the numerical values, produced by porousInterFoam, are shown in red line (interFoambased solver using MULES advection scheme).The free surface values acquired by our new implementation are shown in black dashed line.The grey rectangle in the middle of the domain marks the porous region.

Figure 6 .
Figure 6.Comparison of the interface sharpness for the two solvers.Porous dam break at t = 1 [s] .The porousInteIsoFoam results are shown in black palette, α = 0.5 contour is illustrated in black while α = 0.01 and α = 0.99 contours are shown in grey.Results acquired with porousInterFoam are represented by the red palette, α = 0.5 contour is illustrated in red (this line is overlapped by the black α = 0.5 contour) while α = 0.01 and α = 0.99 contours are shown in light red.The porous region is shown in light grey.The rectangle in dashed outline highlights the area where porousInterFoam produces numerical vapor.

Figure 7 .
Figure 7.Comparison of errors in volume conservation and in upper and lower bounding for the two solvers.