A coupled VOF/embedded boundary method to model two-phase flows on arbitrary solid surfaces

We present an hybrid VOF/embedded boundary method allowing to model two-phase flows in presence of solids with arbitrary shapes. The method relies on the coupling of existing methods: a geometric Volume of fluid (VOF) method to tackle the two-phase flow and an embedded boundary method to sharply resolve arbitrary solid geometries. Coupling these approaches consistently is not trivial and we present in detail a quad/octree spatial discretization for solving the corresponding partial differential equations. Modelling contact angle dynamics is a complex physical and numerical problem. We present a Navier-slip boundary condition compatible with the present cut cell method, validated through a Taylor-Couette test case. To impose the boundary condition when the fluid-fluid interface intersects a solid surface, a geometrical contact angle approach is developed. Our method is validated for several test cases including the spreading of a droplet on a cylinder, and the equilibrium shape of a droplet on a flat or tilted plane in 2D and 3D. The temporal evolution and convergence of the droplet spreading on a flat plane is also discussed for the moving contact line given the boundary condition (Dirichlet or Navier) used. The ability of our numerical methodology to resolve contact line statics and dynamics for different solid geometries is thus demonstrated.


Introduction
Multiphase flows in complex geometries are present in many environmental contexts and industrial applications, from raindrops impacting on leaves and plants [3,21] to droplet deposition on fog harvesting structures [34,44,27], from multiphase fluid-structure interactions in powerplant [35] to 3D printing.In these situations, three different phases are at least concerned: two fluids, a gas and a liquid in general (any liquid vapor contained in the gas phase, commonly, has a negligible effect on the dynamics); and a solid.This makes these types of problems challenging and difficult to tackle either experimentally or numerically because of the complexity of the 1 arXiv:2402.10185v1[physics.flu-dyn]15 Feb 2024 flow coupled with generally non-flat, heterogeneous or even porous solid structures.This can be even more complex when phase change impacts the dynamics, such as in condensation/vaporization or melting/solidification, for instance in steam generators or ice formation.Numerical modelling tools, combined with theoretical approaches and experimental measurements, have an important role to play for understanding the fundamentals of the fluid dynamical interactions in these applications.Such three-phase numerical simulations are difficult however, even in the case without phase change.Three main difficulties are combined: (i) liquid-gas interface tracking, (ii) interaction between the liquid-gas and a solid that may be geometrically complex and, finally, (iii) liquid-gas and solid interaction including models of contact line dynamics.
(i) When it comes to multiphase flow tracking methods, the challenge is to define sharp interfaces while keeping the accuracy of the interfacial attributes (normal, curvature) and robustness.In past decades, major progress has been achieved and various tracking techniques have been investigated in the literature to handle multiphase flow simulations.
A class of methods for handling fluid-fluid interface on a fixed Eulerian mesh is the fronttracking approach [52].A surface mesh is built on the interface in a Lagrangian way while the Navier-stokes equations are solved on the Eulerian fixed grid.Front-tracking approaches intrinsically maintain the sharpness of the interface making this approach one of the most accurate methods to track interfaces and deal with capillary effects.However, they involve high coding complexity of dynamical re-meshing process and mesh management in two-and three-dimensional flows.Moreover, they do not automatically handle interface rupture or coalescence.
On the other hand, front-capturing methods capture the presence of the two fluid interface using an auxiliary Eulerian variable.Among them, the level-set method [51] is an approach that is very often used.This method uses a signed distance function which is continuous across the interface to locate it.The continuous distribution of phases through the distance function allows accurate estimation of the curvature and normal at the interface.Despite these advantages, the level-set method generally suffers from a lack of mass conservation.The Volume Of Fluid (VOF) method [57] is another front-capturing method for which the interface is located with a "color function" which denotes the volume fraction of fluid in a grid cell.Due in particular to its good mass conservation properties and relative simplicity, the VOF approach is one of the most popular front capturing method even though it requires specific efforts to transport interfaces without introducing numerical diffusion and to compute an accurate curvature at the interface.Geometric schemes such as SLIC [36] and PLIC [57] have been developed to give a sharp interface representation.Various authors [50,15,9] have proposed a combination of both level-set and VOF methods to give a simpler representation of the interface thanks to the level-set method while ensuring natural mass conservation with the VOF approach.More recent developments [11,32,14,24,56] have been focused on unstructured grids for VOF and level-set methods.
(ii) Coupling the liquid-gas system with an arbitrary solid generally requires computational tools to describe the solid boundary in addition to the approach selected to model the multiphase flows.Non-body-conforming methods are a class of methods using a Cartesian grid on which the fluid equations of motion (Navier-Stokes) are solved whereas the solid bodies are tracked using a surface mesh independent from the fixed grid.Here, the Cartesian grid does not conform to the solid body, thus it is necessary to modify the equations in the vicinity of the solid boundary to account for its presence.For these approaches, the main difficulties arise in the numerical cells where the three phases are present and for which no simple numerical reconstruction is available, in contrast with two-phase flows [46].Body conforming methods use this latter principle, fitting the grid to the solid boundary.They might appear appealing at first because they give a sharp description of the solid boundary while fullfilling all the boundary conditions, possibly with a high-order of accuracy.However, the lack of flexibility in handling complex geometries, the computational overhead associated with the remeshing process and the additional equations for the mesh motion make them less attractive.Non-body-conforming methods such as the popular Immersed Boundary Methods (IBMs) provide great advantages over traditional body conforming methods since they allow simple grid generation and discretization of the Navier-Stokes equations while requiring less memory.IBMs were originally proposed by [39] for cardiac flow modelling and are used now for many different applications, mostly for fluid-structure interactions.IBMs can be classified into two groups: • Continuous forcing (CF) approaches [39] use a continuous forcing source term based on a regularized Dirac function in the momentum equation to enforce no-slip boundary condition at the solid frontier.CF methods are very attractive for flows with immersed elastic boundaries (biology, interfacial flows).However, in the rigid limit, forcing terms used in these approaches behave less well with for instance problems of numerical accuracy and stability.They are notoriously non-sharp since the regularized Dirac function kernel requires a compact support on Eulerian grids and therefore spreads the solid boundary over a few cells.
• With discrete forcing approaches (DF), the boundary conditions are accounted for at the level of the discretized momentum equations.Since the forcing procedure is directly linked to the discretization method used, DF approaches are not as straightforward to implement as CF methods.However, the discrete nature of the forcing in DF approaches allow to generate sharp immersed boundaries and do not require any stability constraints in the solid body representation.DF approaches can be classified in two categories.Indirect boundary condition imposition on one side and direct boundary condition imposition on the other side with the ghost-cells finite difference method and the cut cell finite volume approach.
A review of the large variety of IBMs is presented in [33].Due to the sharp representation of the solid boundary in DF approaches, the cut-cell or embedded boundary method [10,23,25,40,48,20], which is second-order accurate and conservative, is used in this work, and a more detailed description will be given in the following sections.
(iii) When coupling solid and multiphase flows, different issues have to be considered at their interface, corresponding to boundary conditions that have to be correctly implemented: when the two fluids are in contact with the solid (wetting or dewetting problems), the coupling problem lies in the correct numerical modelling of the contact angle that the liquid/gas interface adopts at the solid frontier.Microscopic physico-chemical interactions between molecules of the two immiscible fluids and the solid drive the contact line behavior.At the macroscopic scale, one can see this contact line as an apparent contact angle between the liquid/gas interface and the solid substrate.Multiple studies throughout the literature have been carried out to give a modelling of this macroscopic apparent contact angle with either static or dynamical value but generally using simple solid geometries [1,16,28].In the last decade, several attempts have been made for more complex geometries using generally IBMs coupled with either VOF, level-set or front-tracking to model the multiphase flow.A diffuse interface immersed boundary method for the simulation of flows with moving contact lines on curved substrates was proposed by [30], where they considered only 2D problems.The authors of [38] used in their work a VOF/IBM coupling with an implicit DF approach to simulate 3D multiphase flow and solid interaction with contact line dynamics.A VOF coupled to a DF ghost-cells/IBM to enforce contact angle dynamics on 2D arbitrarily solids has been proposed by [37].Another immersed boundary based contact angle method coupled with VOF was developed by [22] to handle complex surfaces of mixed wettabilities with a focus on contact angle models.More recently, [4] proposed a wetting benchmark using an unstructured VOF method including a contact line model, allowing to deal with complex solid geometries.
In this article, we present an hybrid model coupling a geometric volume of fluid (VOF) method representing the liquid-gas interface with a cut-cell embedded boundary method accounting for the solid body.The resulting scheme maintains a sharp interface and solid boundary representation, while ensuring mass conservation.While both methods are nowadays classical and have been extensively validated for two-phase flows on one side and solid on the other side, their coupling remains scarce [23].Different boundary conditions for the velocity fields and numerical schemes are addressed in order to be able to model the contact line dynamics which are crucial in such three-phases flows on complex geometries.
This paper is structured as follows: the second section is dedicated to the description of the governing equations with a detailed description of the embedded boundary method and the coupled VOF/embedded boundary method.The contact angle model developed in this work is also presented in static and dynamical situations.In the third section, validation test cases are proposed to show the ability of our approach to accurately deal with three-phase problems involving a triple line in two-or three-dimensional test cases such as a droplet wetting a cylinder or a flattened droplet on an horizontal or inclined embedded plane.

Governing equations
We consider two incompressible, Newtonian and non-miscible fluids and a non-deformable solid so that the coupling between the fluids and the solid intervenes only through the boundary conditions.As a first step, we assume that there is no mass and heat transfer at the fluid-fluid interface and that the surface tension between the two fluids σ is constant.The two fluids will be denoted g and l for gas and liquid, with densities and dynamical viscosities ρ g and ρ l , µ g and µ l respectively.
The dynamics of the two fluids are therefore governed by the unsteady incompressible, variable density, Navier-Stokes equations which can be written, in the classical one-fluid formulation [26]: where ū ≡ ū(x, t) is the flow velocity field, P ≡ P (x, t) stands for the pressure, D = ( ∇ū+ ∇ū T )/2 is the rate of deformation tensor, ρ ≡ ρ(x, t) the one fluid density, µ ≡ µ(x, t) the one fluid viscosity, Fσ the surface tension force and ḡ is the acceleration due to the gravity.The surface tension force is modeled using the continuum surface force [8] (CSF) so that Fσ = σκnδ Γ F with σ the surface tension, κ the curvature, n the normal unit vector to the interface and δ Γ F the Dirac distribution function defined by δ Γ = δ(x − xΓ F ) indicating that it is equal to zero everywhere except on the interface i.e x = xΓ F .Within the one fluid formalism, a discontinuous indicator function χ is used to locate the interface: The evolution of the interface is then simply given by the following advection equation: The density ρ and the viscosity µ are fields that also depend on space and time through the relation: Boundary conditions (BC) between the fluids and the solids have now to be added to this set of equations: indeed the one fluid formulation already provides by construction the correct boundary conditions between the two fluids, through the CSF term.At the solid/fluids interfaces, two different types of boundary conditions need to be imposed: on the velocity fields on the one hand, on the indicator function χ on the other hand.As stated above, we consider the solid as a non-deformable body so that its velocity can be simply described by the classical solid body motion Ūs (x, t) = Ū0 + Ω × x.Therefore, since we have viscous fluids, the usual boundary conditions for the fluid velocity fields at the solid interface reads: This boundary condition imposes in fact both the non-penetration of the fluids inside the solid (the continuity of the normal velocity at the interface) and the no-slip boundary condition for the tangential velocity fields because of the viscous stress.This latter condition is often numerically relaxed by introducing a slip at the solid boundary whose origin can be microscopic (mean free path) and/or geometric (effective roughness of the solid surface), leading to the so-called Navier condition: where λ is the slip length, while ūτ is the projection of the velocity tangentially to the solid (with Ūs,τ the tangential component of the body motion).Similarly, for the interface on the solid, we need to impose the correct boundary condition for the contact line leading eventually to a condition on the characteristic function χ.This corresponds to imposing the contact angle between the interface and the solid surface at the contact line.While the contact angle θ s can be known for static configurations, where the Young-Laplace relation holds, it is much more complex to determine when the contact line is moving, a problem which has been much debated in the literature [12,7].For the rest of the present paper, we will simply assume that this dynamical contact angle θ d is given by a model, typically a function of the velocity in the vicinity of the contact line.For validating the numerical scheme we will use the simplest model for the moving contact line consisting in imposing a constant contact angle coupled with a Navier-slip condition for the velocity.

Description of the numerical method
The numerical scheme developed here is based on the open source library for partial-differential equations, Basilisk [40,41,42].A detailed description of Basilisk can be found at (http:// basilisk.fr).Consequently, we will only give a brief description of the elements when they are not specific to the coupled VOF/embedded boundary implementations.An effort is devoted to the description of the VOF approach used on one side and the embedded boundary on the other, to fully understand the features of the coupled VOF/embedded boundary method.The section is divided into four main parts: (i) the temporal discretization focusing on numerical solution of the incompressible Navier-Stokes equations, (ii) the VOF approach used to deal with the dynamics of deformable fluid-fluid interface, (iii) the embedded boundary method handling the presence of solids, (iv) the coupled VOF/embedded boundary method with its main features, namely when a contact angle boundary condition is required to ensure the triple point/line definition at the intersection between the fluid-fluid interface and the solid.
In order to discretize the equations, the volume fraction F corresponding to the average value of χ in each computational cell Ω F of characteristic volume V is introduced: Then, the evolution of the interface is given by the following advection equation: The local density ρ and viscosity µ are defined from the local volume fraction F using the standard arithmetic mean:

Temporal discretization
The discretization of the Navier-Stokes equations ( 1) is summarized here.The incompressible Navier-Stokes equations are approximated using a a time-staggered approximate projection method on a Cartesian grid.The advection term is discretized with the explicit and conservative Bell-Colella-Glaz second-order unsplit upwind scheme [6].For the discretization of the viscous diffusion term, a second-order Crank-Nicholson fully-implicit scheme is used.Spatial discretization is achieved using a quadtree (octree in three dimensions) adaptive mesh refinement on collocated grids.The adaptive mesh technique is based on the wavelet decomposition [42,53] of the variables such as the velocity field or the volume fraction and is efficient and critical to the success of some simulations as it allows high resolution only where needed, therefore reducing the cost of calculation.
The following fractional step method is used: • The Navier-stokes equations (1b)-(1c) are first discretized as . ∇ū with ∆t def = t n+1 − t n the time step.The subscript c and f refers to the cell center and cell face respectively.Equation (15) gives a temporary velocity ūn+1 f obtained by an arithmetic averaging from the centered velocity field ( c→f ) to which the acceleration field ān+1 f is added.This acceleration term is simply due to the surface tension forces and the gravity, by considering (ā = σκnδ Γ F + ρḡ).To ensure a consistent discretization, the acceleration term is defined on the faces f like the pressure gradient, in order to minimize the parasitic currents close to the interface, according to the well-balanced surface tension model [41,43].The combined centered acceleration and pressure gradient is also obtained by arithmetic averaging ( f →c ).
• To ensure incompressibility, the temporary velocity (15) is projected so that: with ∇•ū n+1 f = 0.The pressure P n+1 is then obtained using a multigrid iterative method to solve Eq. ( 16).We get the divergence free face velocity ūn+1 f by substituting P n+1 in Eq. ( 16).
• The final centered velocity Eq. ( 17) is obtained using the combined pressure gradient and acceleration interpolated from faces to cells.

Interface advection (VOF)
The geometric VOF PLIC method on an Eulerian fixed grid is used, because of its accuracy to maintain a sharp interface between the fluids while ensuring mass conservation.The integration of the hyperbolic advection equation Eq. ( 8) is done in two steps.First, the explicit interface is locally reconstructed from the volume fraction and then the advection is performed using a split advection method.

Geometric reconstruction of the interface
Consider a discrete representation of the liquid-gas system separated by the interface Γ F , integrated on a regular Cartesian grid, composed of squares (cubes in 3D) cells of size ∆.Each cell is a control volume V, delimited by the closed contour δV and associated with a velocity or pressure variable.The fluid-fluid interface Γ F is represented with a piecewise-linear reconstruction (VOF PLIC), satisfying in each cut cell the following equation for a line (or plane in 3D): where nΓ F is the the local normal to the interface, x the vector position and α F the intercept.The normal nΓ F is estimated using the Mixed-Youngs-Centred (MYC) implementation of [5].By ensuring that the fluid volume contained within a plane of normal nΓ F is equal to F , the intercept α F can be determined in a unique way, using analytical relations [47] already implemented in Basilisk:

VOF advection and geometric flux computation
Once the reconstruction step is completed, the advection equation: is integrated using dimension-splitting.The method proposed by [55] is used to discretize the equation ( 19) which is solved along each dimension successively, with a fractional step method in order to guarantee exact mass conservation for the multi-dimensional advection scheme as long as the velocity field is divergence free.The left-hand-side term is computed using the net fluxes per direction whereas the second term represents the dilatation required to compensate the non-divergence-free flow advection per direction.With an update of the interface from t n− 1 2 to t n+ 1 2 (it is assumed here that the interface lags the velocity/pressure fields by half a time step), we get: where G and H are the face fluxes estimated respectively by F u xf and F u y f .To ensure exact volume conservation, the centered volume fraction field is set to: As shown in Figure 6a, the fluxes G and H are computed using the geometry of the reconstructed interface [36,13].In this illustration, the flux G, corresponding to the horizontal face velocity F u xf and delimited by the dotted line is estimated by the dark red area.This area can be computed using an approach similar to that in the previous section.Once the volume fraction F n+ 1 2 is updated, the local properties ρ n+ 1 2 and µ n+ 1 2 can be evaluated using Eqs.( 9)-( 10).

Computing the surface tension force
In the framework of VOF methods, the continuum surface force model (CSF) [8] appears to be the natural way to compute the surface tension term in the momentum equation.
As mentioned previously (see Sec. 3.1), the surface tension force is computed here using the balanced-force calculation method [41,43].Both pressure gradient and surface tension terms are defined on cell faces to ensure a consistent discretization and eventually reach an exact balance between these two terms.Computing the curvature accurately is more difficult and is done here using the generalized height-function approach of [41].

The embedded boundary or cut cell method
The Cartesian grid embedded boundary method, also known as cut cell method, used in this work has already been implemented in Basilisk and used for moving boundaries [20] and solidification modelling [29].We give here the general lines of this method, originally designed by [25,48].

Solid-fluid interface description
A discrete representation of a solid domain Ω s with boundary Γ s is embedded in a regular Cartesian grid, composed of square (cubic in 3D) cells of length ∆ as shown in Fig. 1b.Each cell is a control volume V, delimited by the closed contour δV consisting of edges f d (faces in 3D).The index d refers to the direction of the cell so that d is the left, bottom, right or upper direction for a 2D cell.The reconstruction of the solid interface Γ S is done using a piecewise linear reconstruction thus satisfying in each cut cell the following equation for a line (a plane in 3D): where nΓ S is the the local normal to the interface, x the position vector and α S the intercept.The subscript F , S refers respectively to the fluid and solid domain.We consider firstly the simplified where D is the dimension of the problem.The complementary solid volume is therefore Similarly a fluid face fraction 0 < f d F < 1 is defined for each face d so that the effective area occupied by the fluid on a cut face is: and the solid area is simply given by With this integral description of the discrete rigid boundary, it is possible to represent arbitrary solid shapes whose size is greater than the cell size ∆.

Embedded fractions in a cut cell
The geometric properties including the normal, the face and the volume fraction of the interface are computed using a signed distance function ϕ(x− xΓ S ) sampled on the vertices of the Cartesian grid cell (see Fig. 1b).This signed distance function represents the solid-fluid interface Γ S and is taken by convention positive "inside" the interface when x < xΓ S and negative outside when x > xΓ S .Figure 1a shows the description of a Cartesian square (cubic in 3D) cell of size unity centered on the origin.We give here some elements of the interface characteristics calculation for a 2D cell.Consider for instance the bottom left and right vertices, respectively ϕ[0, 0] and ϕ[0, 1] to evaluate the bottom line fraction f d = f y F .• For each edge of an interfacial cell (i.e ϕ[0, 0] × ϕ[0, 1] < 0), the bottom line fraction for instance f y F is given by: • The normal nΓ S is found by expressing the circulation along the closed boundary δV of the volume control V: where nd is the outward unit normal of each face f d and nΓ S the non-unit outward normal to the interface Γ S .
• For each edge, if 0 < f d < 1, then α d S = nΓ S • xd .For the bottom edge, the coordinates of xd are for example: The final intercept α S is the average of all face intercepts α S d for all intersections.• Given nΓ S and α S , the interface is fully described and the volume fraction can be calculated using the predefined functions [47] mentioned above so that The generalization to 3D is made quite straightforward by considering the faces of the cell instead of the edges.The whole procedure described above to get the volume fraction C in 2D is now used to obtain the face fraction f d .The rest of the procedure consists in getting the normal, the intercept and finally the volume fraction with some adaptations to account for the third dimension.The complete procedure of the geometric calculation of the interface properties is detailed in http://basilisk.fr/src/fractions.h.From the face fractions, we are also able to define the position of the centroid bΓ S in a cut cell: which is an important quantity for approximating differential operators on interfacial cells.

Discrete operators in a cut cell
Once the embedded fractions and the centroids are computed, it is possible to write a conservative finite volume approximation of discrete operators in the cut cell such as the divergence of a quantity Φ (typically needed to compute the non linear advection term or the viscous term Eq. (1c) in the fluid side for instance): where nd is the outward unit normal vector to the face d of the closed boundary δV within the cell.The challenge is thus to estimate accurately such operator taking into account the correct boundary conditions at the solid/fluid frontier.The quantity ΦΓ S is discretized at the centroid b Γs as well as Φ which is taken at the center of the full or "regular" Cartesian face d rather than at the center of the irregular fluid face, even when the center of the Cartesian cell is outside the fluid domain.Here, it is assumed that Φ is sufficiently smooth to be extended to the cut cell (0 < C < 1), even if its center is included in the solid part of the cell.This leads to a conservative and at least second-order accurate, finite-volume methodology.
The detailed algorithms and demonstration of 2nd order accuracy of the embedded boundary method can be found in [25,40].The interest in using embedded boundary on Cartesian grids over unstructured grid methods or other immersed boundary methods is the generation of simple grids and data structures in 2D or 3D and the straightforward extension to quad/octree meshes.Complete details of how this method is implemented in Basilisk are available in http://basilisk.fr/src/embed.h.
Finally, a crucial issue with cut-cell methods is to implement accurately the correct boundary conditions at the solid boundaries for the velocity fields and the contact angle in the case of multiphase flows.While the non-penetrating condition for the velocity field does not require additional treatment for the mass conservation equation, the conditions for the velocity field needs to be taken into account when estimating the viscous term in the Navier-Stokes equations.We therefore investigate firstly how to impose any boundary condition for the velocity in the next section.

Viscous flux through the discrete rigid boundary in a cut cell
For computing the viscous term in the Navier-Stokes equations, we need, in our formulation, to access the gradients of the velocities in order to evaluate the fluxes associated with the viscous stress tensor.Considering an incompressible flow and constant viscosity, the viscous term of the Navier-Stokes equations Eq. (1c) can be further simplified so that: The right hand term of the relation (33) then gives decoupled scalar Poisson equations for each velocity component.To compute the viscous flux through the discrete rigid boundary Γ S in a cut cell, we then need to write the conservative finite-volume approximation for each velocity component taking the horizontal component of the velocity for instance.Here, dS is the elementary boundary surface and ∇u x •n Γ S is the embedded face gradient that can also be written ∇u x | Γ S .The gradient must be determined at the centroid of each face f d and also at the centroid bΓ S (see 32) of the cut cell interface Γ S (see fig. 1b) .Since the calculation of the flux through the face f d is detailed in [20], we focus here solely on the gradient calculation over the discrete rigid boundary Γ S of centroid bΓ S , which is crucial for the accuracy of the method.We consider both Dirichlet and Navier boundary conditions on the velocity.

Dirichlet boundary condition
If a Dirichlet boundary condition u xΓ or u y Γ (u z Γ in 3D) is required for the velocity on the solid embedded boundary, the following second-order discretization is used to compute the embedded face gradient of the scalar velocity component.Following the methodology presented in [25], we write for the horizontal component: The gradient ∇u x | Γ S or ( ∇u x • nΓ S ) is evaluated at the centroid b Γ S of the cut-cell (see Figure 2: Gradient calculation in the embedded boundary method depending on the normal orientation nΓ S . Fig. 2).The second-order accuracy of the gradient is achieved by computing u x I1 and u x I2 with a third-order quadratic (biquadratic in 3D) interpolation along the direction orthogonal to the direction of the normal nΓs at distances d 1 and d 2 from the centroid b Γs (see http://basilisk.fr/src/embed.h#dirichlet-boundary-condition for more details).

Neumann boundary condition
To impose Neumann boundary condition s Γ on the velocity component u x , the gradient ∇u x | Γ S is thus defined precisely at the interface and we simply have to take for evaluating the viscous flux (34) this value: Both these Dirichlet and Neumann boundary conditions have been validated in previous works [20,29], showing at least second-order convergence for the multigrid Poisson-Helmholtz solver in presence of embedded boundaries.

Navier boundary condition
As for the Dirichlet condition, taking into account the Navier boundary condition Eq. ( 6), is not straightforward since we need to express the embedded gradient on a fixed Cartesian grid while taking into account the local Navier boundary condition.Considering the tangential and normal local coordinates (τ , n) for a 2D system, attached to the solid rigid boundary, we write: ∇ū • n = ∂uτ ∂n ∂un ∂n (37) The transformation between the local system (τ , n) and the fixed Cartesian system ( ēx , ēy ) is simply given by: ∇ū with We can then express the embedded gradient on Γ S for each velocity component: Following the Navier boundary condition in Eq. ( 6), and considering Ūs,τ = 0 (motionless embedded boundary), ∂u τ ∂n and ∂u n ∂n are a priori known at the interface Γ S since: Given the relation u n = 0, we use the second-order gradient calculation Dirichlet method (35) to evaluate ∂u n ∂n .The term ∂u τ ∂n can be approximated using the same principle, But here, the value of the tangential velocity u τ Γ needs to be determined.Following [25] we can write the relation ( 42) where the coefficients a 0 , a 1 and a 2 are directly determined by (42).Combining equations ( 41)-( 43) , we can eliminate ∂uτ ∂n and deduce the following relation to determine Given the value of u τ Γ , the Dirichlet gradient calculation Eq. ( 35) is used again to calculate ∂u τ ∂n .Note that this method degenerates naturally into a Dirichlet boundary condition when λ → 0 and into a free-slip condition when λ → ∞, so that is can be used formally for any value of λ.The generalization to 3D is also straightforward.

Validation of the Navier boundary condition implementation
To validate our method, we solve the (Stokes) Couette flow between two rotating cylinders using the Dirichlet, Navier or free slip boundary condition with embedded solid boundaries for a single phase flow.We set the geometry of the embedded boundary as two cylinders of radii R 2 = 0.5 and R 1 = 0.25 (see fig. 3).Here, the outer cylinder is fixed and the inner cylinder is rotating with an angular velocity ω.For the free-slip boundary condition (outer cylinder free of viscous stresses), the analytical angular velocity u θ is given by: For a Navier boundary condition on the outer cylinder, we have: The simulation is performed in a two-dimensional box of size D = 2R 2 .Adaptive mesh refinement is not activated and we consider only a uniform Cartesian grid with respectively N = 16, 32, 64 and 128 grid points within the box size.Figure 4a illustrates the comparison between the numerical solution for N = 32 and the analytical solution of equation (45) in the free slip case.We observe very good agreement between both numerical and analytical solutions.The same comparison is carried out for the Navier boundary condition considering different values of λ, respectively 0.4, 0.5, 0.8 and 1.The agreement between the numerical and the theoretical solutions is also good whatever the value λ.To further quantify the differences between numerical and analytical solutions, we define the errors:    47), E ∞ Eq. ( 48) and order of convergence.
Considering the free slip case, we can notice in table 1 that the numerical solution converges to the theory at 2nd order for E 1 and E ∞ .For the Navier slip case (table 2), the numerical solution also converges to the theoretical values but is sensitive to the value of λ.As the value of λ is decreased (from 1 to 0.5), the order of convergence slowly drops from 2 to ≈ 1.4 on average for E 1 .The same trend applies for E ∞ in this case.

Hybrid VOF/embedded boundary method
We consider now a three phases system with a liquid-gas interface namely the liquid (Ω L ), the gas (Ω G ) and the solid (Ω S ) (see Fig. 5c illustrating the typical configurations encountered).
To solve such problems, we can a priori simply couple the two methods described in the previous section, the Volume of Fluid one for the liquid-gas, and the embedded boundary method for the solid-fluid system.As stated in (see fig. 5c,5b) the flow must now account for two volume fractions: F to track the liquid-gas interface and C for embedded solid cells.Both volume fractions F and C are defined in the whole domain Ω G ∪ Ω L ∪ Ω S following a priori: whereas  It is important to note that the value of liquid-gas function F in the domain where C = 0 can in fact take any value since it is not used so far in the equations.Moreover, since the embedded boundary method only intervenes in the fluid cells near the cut-cells, the coupling is at work only in this region.Therefore, in the presence of a solid boundary, the solution of the two-phase flow system (Eqs.( 1),( 8)) described in section 2 needs to be adapted only near the solid boundary, since it will not at all be considered in the cells where C = 0. Let us first consider consider the VOF advection, to be performed in the full fluid cells but also in the "mixed" cells (0 < F < 1 and 0 < C < 1, Ω G ∩ Ω L ∩ Ω S ): a simple approximation, which we use here, is to ignore the solid volume fraction.
This ensures robustness of the VOF advection and preserves the conservation of the total fluid mass, if the latter is computed ignoring the volume occupied by the solid in "mixed" cells.
The solid can then be interpreted as "slightly porous in surface" i.e. able to retain some of the fluid, but only within a thin shell (i.e.thinner than the mesh size).So we expect this method to be no more than first-order in convergence.It will be therefore useful to quantify this "local mass absorption" and this is done in the validation section 4 for different test cases.We observe that it can reach 5% in the most severe configurations (smallest contact angles).A rough way to compensate for this mass absorption could be to add or subtract the mean mass absorption from each cell at the fluid interface [38].In that case, one must assume that the mass compensation does not affect the global flow, thus the local mass absorption must be small.A higher-order approach could also be employed by modifying the flux in "mixed" cells as illustrated in Figure 6a to calculate the flux G = F u xf , including the solid volume.Figure 6b shows the horizontal VOF flux calculation in a "mixed" cell.Here, the dotted line delimits the flux going to the right-hand neighbor.Even if the normal face velocity is weighted by the face fraction u xf f d , the total flux estimation includes the solid region.A better estimate of this flux should take into account the solid region area in the calculation.In the example given in Fig. 6b, the exact geometric flux estimation is straightforward as the solid region is a rectangle.In more general cases including 3D configurations, this estimation is less trivial and requires a specific algorithm to account for the possible orientations of general solid geometries.This problem is not addressed here but would be an interesting improvement of the coupled VOF/embedded boundary method.
The other important issue when coupling VOF and embedded boundary is the boundary conditions to impose on the liquid-gas volume fraction F .For the single flow cut cells (where F = 0 or F = 1, with 0 < C < 1), as illustrated on figures 5c and 5b, no special action is required, beside taking into account the correct viscous boundary conditions explained above.On the other hand, the situation is more complicate for mixed cells since the reconstruction of the interface is crucial to estimate accurately the fluid fluxes and the surface forces.This corresponds also to the presence of a triple point/line in the cell as illustrated in Figure 5d.In this configuration a contact angle boundary condition has to be explicitly set because of the non-conforming body description of the fluid and solid region.To do so, we use "ghost cells" taking advantage of the fact that the function F can be arbitrarily set in the full solid cells (where C = 0).The principle of the ghost-cell method is then to impose the values of F in the full solid cells near the mixed cell in order to obtain the correct contact angle for the interface in the mixed cells.

Contact angle boundary condition with the VOF/embedded boundary method
In a three-phase system, the contact angle computation is one of the main challenges since it plays a major role on surface tension effects controlled by the motion of contact lines and therefore the wetting/dewetting of the solid.Our objective is to propose a simple methodology which captures the physics of the wetting by imposing a prescribed contact angle.This angle could either be static θ s or dynamical θ d (t) with an instantaneous dynamical angle determined with an appropriate model [49,28,1].In what follows, the algorithm to enforce this contact angle (noted θ s for simplicity further on, with no loss of generality) is described.The contact angle θ s is thus imposed at the embedded boundary wall by setting the corresponding fluid volume fraction F in Ω S , using the concept of ghost cells [4,38,37].This is done in several steps: • First, according to a prescribed value of θ s , the properly-oriented normal nθ at the liquidgas-solid intersection is determined: For "mixed" cells (0 < C < 1 and 0 < F < 1), nθ is calculated so that the intercept

/ normal a t t h e t r i p l e p o i n t
4 algorithm (see Tab. 3) summarizes the first steps of the contact angle procedure with the calculation of the interface geometric properties nθ and α θ of the cells with three phases.

nθ ) // i n t e r c e p t a t t h e t r i p l e p o i n t c e l l
• We then consider the neighboring ghost cells (C = 0) where the volume fraction F has to be modified to account for the contact angle condition.If the neighbor of a ghost-cell within two cells is a cell with three phases (0 < C < 1 and 0 < F < 1), a coefficient to weight the fluid volume fraction in the "ghost-cells" depending on how well the "mixed" cell is defined, is computed.We assume that a "mixed" cell is better defined when C ≈ F ≈ 0.5.Knowing the values of the intercept α θ and the normal nθ in these cells (see Tab. 4), the prescribed interface (accounting for the prescribed value of θ s ) is extrapolated into the considered "ghost cells" and a new fluid volume fraction (F = F(α θ , nθ )) is computed.The final volume fraction F G in the "ghost cells" is a weighted average of the volume fractions reconstructed from the cells neighboring a cell with three phases.

Algorithm 2
1 f o r a l l c e l l s

i n t h e g h o s t c e l l s
3

w e i g h t t o e s t i m a t e t h e volume f r a c t i o n r e c o n s t r u c t i o n
4 f o r a l l n e i g h b o r s w i t h i n 2 c e l l s

then // p o t e n t i a l t r i p l e p h a s e c e l l
6 This contact angle procedure allows to apply the proper boundary condition at the liquid-gassolid intersection by modifying only the values of the fluid volume fraction F inside the "ghost cells".A shown in fig.7b a neighborhood of two "ghost cells"around the triple point is required to obtain an accurate definition of the contact angle θ s .

Contact line dynamics
The motion of the contact line along a no-slip solid surface involves a mathematical paradox.Indeed, the no-slip boundary condition is known to introduce a "stress singularity" at the contact line since the stress generated by the moving contact line diverges with grid refinement.Several authors in the literature have proposed various models to deal with this singularity [1,19,2].Using the VOF method intrinsically removes this contradiction as the volume fraction advection uses the face centered velocity which is not strictly zero, even if a zero velocity is imposed at the solid boundary.Thus, the contact line can move due the so-called "numerical slip" [45].This "numerical slip" is linked to the grid size and tends toward the no-slip limit as the size of the mesh is decreased.A way to control this mesh dependency while relaxing the stress singularity at the contact line is to use slip models.The most common slip model in the literature is the Navier slip law Eq: (6).In what follows, We will use either Dirichlet (35) or Navier Eq: ( 40) boundary conditions coupled with the contact angle calculation (thus taking a fixed θ s ) described above for the contact line dynamics.

Validations tests
In this section, we investigate different configurations in order to test and validate the implementation of the hybrid VOF/embedded boundary method.These test cases are designed to challenge the VOF/embedded coupling and particularly the boundary condition at the triple point/line with the geometrical approach detailed in the previous section.
The constant static contact angle θ s is imposed and we start in general with an initial condition where the contact angle is θ i ̸ = θ s so that we investigate the relaxation towards an equilibrium static configuration.Indeed, since the initial contact angle is out of equilibrium, the VOF/embedded dynamics should converge towards a situation where the contact angle is the static one.In most cases, analytical static configurations exist that we will compare to the numerical solution.We will investigate both 2D (plane) and 3D geometries.
We first simply study the "numerical dynamics" that is using no-slip boundary conditions on the solid boundaries and no dynamical contact angle modeling.These conditions are known to mimic a (numerical) Navier slip boundary conditions at the solid boundaries [1,16] and we will investigate a classical dynamical model for the contact line dynamics using a slip length [18].We will therefore consider droplets on various solid geometries: a cylinder, an horizontal plane and an inclined plane with or without gravity.Except when explicitly stated, the density and viscosity ratios for the liquid-gas interface are set to ρ l /ρ g = 1 and µ l /µ g = 1 and the surface tension σ = 1.To ensure the stability of the simulations, the time step is constrained by the oscillation period of the smallest capillary wave given by ∆t = . Moreover, the maximum CFL number is set to 0.5 for all cases.

Droplet on an embedded cylinder
We investigate in this test case the ability of the coupled VOF/embedded boundary method described above to impose a contact angle on a curved solid geometry.The equilibrium shape of a droplet on a solid circular cylinder of radius R c = 0.5 is investigated here for different static contact angles θ s .As shown in Fig. 8, a droplet of initial radius R 0 = 0.5 wets a cylinder with an initial contact angle θ i = 90 • .In the absence of gravity, the equilibrium shape of the droplet can be obtained analytically.For a two-dimensional droplet, the area corresponding to the equilibrium state is obtained by carving the cylinder out of a circle so that, the initial area S 0 of the droplet is given by: The final radius of the droplet R and the distance between the center of the cylinder and the droplet l c can be computed by finding the roots of the equation with and Figure 8: Schematic representation of the initial configuration of the droplet on a cylinder.
The simulation is performed in a two-dimensional box of size 2D × 2D with D = 2R 0 .It is enclosed by solid walls for all boundaries except the cylinder which is described using the embedded boundary method.The contact angle boundary condition is applied using the coupled VOF/embedded boundary method as described above.The adaptive mesh refinement is not enabled here so that we always have a Cartesian uniform grid with respectively N = 32, 64, 128 and 256 cells and 8, 16, 32 or 64 cells describing the initial radius R 0 of the drop (for instance R = 16∆ for N = 64).Five different contact angles θ s , 30   shows satisfactory qualitative agreement with the analytical solutions for most of the contact angles studied, as shown in Figures 9.More precisely, the agreement is better for hydrophobic contact angles while it is less good for hydrophilic cases (see for instance θ s = 30 • ).Figures 10a-10b shows the analytical and numerical solutions as the grid is refined, for θ s = 30 • and θ s = 120 • .
For θ s = 30 • , the discrepancy between analytical and numerical solutions is clearer.However, it is difficult to state if the numerical solution slowly converges to the analytical one with grid refinement or not.The lack of convergence is related to a numerical pinning of the contact line: depending on the mesh orientation with the embedded boundary, a situation where the VOF flux advection reduces to zero can occur and the contact line cannot relax anymore to the prescribed angle.This pinning effect is clearer for the test case of the droplet on an inclined plane (see sec.4.2).This "pinning effect" can also be interpreted through the dependency of the "numerical Navier slip" condition on the details of the cut-cell configuration: for (nearly) full fluid cells close to the embedded boundary, the numerical Navier slip length is of the order of ∆/2 while for (nearly) full solid cells this numerical slip length tends toward zero (i.e.no-slip).To quantify the accuracy of our method, we also computed the numerical final radius of the droplet shown on figures 11a.Except for the smallest value of θ s = 30 • , the numerical and analytical radius are in good agreement.To quantify it more precisely, we show on figure 11b the relative error on the radius as a function of the mesh size ∆ for θ s = 30 • , 90 • and 150 • , exhibiting roughly the expected first-order convergence.Furthermore, we also compute the relative mass variation during the simulation, plotted on figure 12a for different contact angles θ s for N = 64.We observe that this mass variation converges to an asymptotic value at large times, when the droplet has reached its equilibrium state.
Figure 12b shows this relative mass absorption with respect to the mesh size ∆.We observe again the expected first-order of convergence.The origin of this mass absorption has been previously analyzed and the results are in agreement with the "porosity" of the coupled VOF/embedded boundary method.Indeed, in a cell with three phases, the solid fraction is simply ignored (the solid is assumed to be "porous" at the ∆ scale).

2D sessile droplet on an embedded plane
In this test case, a two-dimensional drop of radius R 0 is released at rest with an initial contact angle θ i = 90 • on an embedded plane with a static contact angle θ s .Two possible orientations of the embedded plane are studied: either 0 • or 45 • .In the first situation, we have an embedded geometry perfectly aligned with the grid, an ideal situation, whereas the latter configuration represents an embedded plane poorly aligned with the Cartesian grid.Initially, the drop is placed in the domain as an half-disk (i.e. the initial contact angle is 90 • , see Fig. 13).From its initial contact angle, the droplet will relax to reach its equilibrium position (the prescribed angle θ s ).
If the effect of gravity is taken into account, the drop is also flattened on a horizontal plane or can move if the inclined plane is considered.The Eötvös number Eo is then used to quantify the relative effect of gravitational forces compared to surface tension forces: Eo = ρ L gR 2 0 /σ.

Spreading droplet without gravity
First, the simulation is performed for the case where gravity is zero (Eo = 0).In the absence of gravity, the drop is hemispherical and its shape is controlled by surface tension only.As the total drop volume V is constant, it is possible to geometrically calculate at equilibrium the radius of the circle R f , the radius of spreading r f and the height of the drop h f .
The simulation is still performed in a two dimensional box of size 2D × 2D with D = 2R 0 .
Figure 13: Schematic representation of the initial and equilibrium shapes of a droplet on a flat surface with static contact angle θ s and Eo = 0.
All the boundaries are no slip walls and the contact angle boundary condition is applied as described above.The diameter D is chosen as the characteristic length and set to 1 in order to have W e ≈ 1 and Re = 1/µ L .The adaptive mesh refinement is not activated here and the numerical domain is a Cartesian regular grid with the number of cells N = 64, 128, 256 and 512.Only 32 cells describe the initial diameter D of the drop for N = 64 (and thus R 0 = N ∆/4 accordingly).We study the case where the contact angle θ s is varied between 15 • and 165 • .The comparison between the numerical and the analytical equilibrium droplet shapes shows very satisfactory qualitative matching (Fig. 14a) for the horizontal embedded configuration.In the inclined plane situation, the numerical shape of the equilibrium droplet does not match perfectly the analytical one specially in the most hydrophobic or hydrophilic situation.This issue is clearer when comparing the analytical droplet height h f and contact radius r f to the numerical results.Figure 15a presents the numerical droplet height h f and contact radius fr f for the equilibrium shape compared with the analytical results of Eqs.(56b)-(56c), respectively as a function of θ s .For the horizontal embedded configuration; all the considered contact angles show a good   It is less obvious for the inclined plane where both the numerical height h f /R 0 and radius r f /R 0 diverge from the analytical data for extreme angles (15 • , 30 • and 15 • , 165 • ).This problem, already mentioned in the previous test case (see sec.4.1), is linked to the poor grid alignment with the embedded boundary.With a 45 • inclined plane, the VOF interface might cross the embedded boundary and the Cartesian grid in such a way that the VOF flux from one cell to another is close to zero and thus the interface gets pinned.It is to note that the pinning of the contact line is perfectly symmetric between the hydrophobic and the hydrophilic case in that specific configuration (45 • ).
We also define a quantitative error by the mean difference of the volume fraction of the equilibrium drop shape and the volume fraction computed with the corresponding analytical shape: Figure 16a and 16b show the average error E F defined in equation ( 57) for the horizontal and the inclined plane.Surprisingly, a non-monotonic convergence of E F is observed with the grid refinement whatever the angle considered and the embedded configuration.In both situations, the maximum error is obtained for extreme angles (15 • and 165 • ) with a higher level of error (by about a factor 10) for the inclined plane though, mostly for the reasons set out above.The mass absorption of the method is also examined and reported in Fig. 17a and Fig. 17b for the horizontal and 45 • inclined embedded plane and for different constant angles θ s .It appears that for the extreme angles namely 15 • and 165 • the maximum mass absorption is observed.As in the previous test case (see sec.4.1), this is well quantified and linked to the "porosity" of the coupled VOF/embedded boundary method.The spurious currents evolution for the spreading droplet on the flat solid surface is also quantified through the capillary number defined as: Figure 18 shows the evolution of the capillary number (58) as a function of time for three different contact angle θ s = 45 • , θ s = 75 • , θ s = 120 • , two grid resolutions N = 32 and N = 64.For all combination of parameters, the maximum velocity eventually converges to zero within machine precision.

Influence of the gravity on a droplet resting on an embedded horizontal plane
The influence of gravity (Eo > 0) on the droplet shape is now presented.While gravity tends to spread the drop, the surface tension keeps the drop as spherical as possible while the contact angle at the boundary must be satisfied.For small Eötvös number (Eo << 1), the drop shape is controlled by surface tension and the maximum height of the drop h 0 is given by equation (56b).At larger Eötvös number (Eo >> 1), gravity influences the drop shape which starts to form a puddle.The maximum height h ∞ of the drop is then proportional to the capillary length scale and given by the following expression [16,4,38]:

3D sessile droplet on an embedded horizontal plane
We investigate now the same spreading of a spherical droplet on an horizontal plane but in three dimensions.A drop is initialized as a half-sphere of radius R 0 = 0.25 with θ i = 90 • .The analytical equilibrium radius of the droplet is given by: The contact angle is varied between 30 • and 150 • .The drop oscillates and eventually relaxes to its equilibrium position θ s .Here, the adaptive mesh refinement is enabled so that the maximum level of refinement corresponds to R = 16∆.The comparison between the numerical and the analytical equilibrium droplet shapes shows satisfactory qualitative matching (Fig. 20b), even for such a coarse grid.This demonstrates that the algorithm developed for the contact angle boundary condition is also compatible with adaptive mesh refinement.

Dynamics of the droplet on the embedded horizontal plane
So far we have studied the asymptotic regime of a configuration to validate that our simulations converge numerically to the equilibrium state.No specific treatment was made for the dynamics and we simply relied on numerical relaxation.A physically-correct description of the problem must involved two coupled numerical ingredients: the computation of the viscous stress in the mixed cell and the moving contact line.In particular, we know that without any correction, the contact line dynamics should not converge since in the continuum limit, the no-slip boundary condition does not allow the contact line to move [17].
In the following test case, the dynamics of the spreading droplet is therefore studied.The objective is to first emphasize the grid dependency when using a no-slip boundary condition on an embedded solid thus characterizing the effect of the "numerical slip" and then introduce a Navier slip boundary condition Eq. ( 6) in order to (partially) regularize the singularity at the contact line.The contact line behavior as a function of the boundary condition has been debated in many studies for the past decades but remains a numerical challenge especially in the context of hybrid method such as the VOF/embedded boundary in this work.Note that for the Navier boundary condition, grid convergence can be reached if the full hydrodynamic problem is solved in the slip length region [1,28,19].Here, we consider again the spreading of a semicircular droplet on an embedded solid boundary.We investigate the dynamics of the contact angle when the droplet relaxes to its equilibrium shape given by θ s .The physical parameters are identical to those in section 4.2.The contact angle is fixed at θ s = 60 • and the mesh size is varied from N = 32 to N = 256.Gravity is ignored in this simulation.Either a no-slip (λ = 0) or a Navier slip (λ ̸ = 0) condition is imposed along the embedded solid using the procedure described before Eqs.(35) and (40).In this simulation, λ is chosen as an adjustable parameter.Figure 4.2 shows the evolution of the droplet radius r with the dimensionless time for different values of λ and a fixed grid N = 32.For all values of λ, a convergence of the droplet radius r is finally reached.However, as the resolution of λ is increased, we observe a quick convergence of the radius.At least 4 cells are required to describe the slip length λ and a convergence is reached when λ = 16∆ for N = 32.This value of λ = 16∆ 32 is therefore adopted for the rest of the simulations and for all the grids.Considering now the grid convergence, we compare both the no-slip case and the Navier slip case.In both cases, the radius converges to the terminal radius r f .As expected, the radius r does not converge with mesh refinement fig.22a in the no-slip case.The contact line velocity decreases as the mesh resolution is increased.It is a consequence of the numerical slip length which is larger for coarser meshes: as the slip decreases with finer meshes, so does the contact line velocity.On the other hand, the use of the Navier slip boundary condition changes the spreading dynamics of the droplet.By introducing a Navier slip length larger than the grid cell, the results tends to be independent from mesh refinement.Moreover a faster drop spreading is observed in that situation.The Navier boundary condition coupled with a constant contact angle is then sufficient to achieve grid convergence.

Drop impact on a fiber
An additional validation study for a 3D drop impact on a circular fiber is presented here.The parameters used in the simulations are inspired by the experiment of [31] and the related numerical simulation of [54] where they studied a liquid droplet impacting a thin solid fiber in a gaseous environment.In both works, the static contact angle was chosen to be 10 chosen as the characteristic length and velocity of the problem.The circular fiber has a diameter Df so that the ratio Df /D ≈ 0.485, which is slightly larger than the reference fiber diameter (Df /D ≈ 0.416 see [54]).We consider two cases where the impact velocity is varied so that we get the non-dimensional parameters, Re = ρ l U D µ l , W e = ρ l U 2 D σ , F r = U 2 gD , given in table.5. We chose here to work with a static contact angle θ s = 15 • .The domain size is 4D 3 discretized with  an uniform mesh of 128 3 , resulting in Df /∆ ≈ 14.7.
Figures 23 and 24 show snapshots of impact dynamics on a fiber.For the first case, the  drop is captured by the fiber whereas the drop detaches from the fiber in the second case as in the experimental and numerical references.Our 3D simulations show good agreement with both the experimental and numerical observations of the aforementioned references.We do not perform here a more quantitative validation since it would require a better characterization of the experimental conditions, which is beyond the scope of this work.

Conclusion and discussion
In this work, a hybrid VOF/embedded boundary for the modelling of two-phase flows interacting with arbitrary solid geometries is presented.A conservative second-order embedded boundary method is used here to simulate complex solid shapes while preserving the accuracy at the fluid-solid interface for Dirichlet/Neumann like boundary conditions.A specific effort has been devoted to propose a second-order accurate Navier/free slip boundary condition compatible with the embedded boundary method.Thanks to the geometric Volume Of Fluid method (VOF), the fluid-fluid interface is tracked sharply in a conservative way.
To account for the contact angle boundary condition at the intersection between the solid and the fluid-fluid interface (triple point/line), an algorithm using fluid "ghost cells" within the solid region is proposed, allowing to reproduce contact line dynamics in the vicinity of the triple point/line.
Several test cases have been investigated to validate our method, namely the spreading of a droplet in various configurations (cylinder, horizontal or tilted plane, with or without gravity).In the majority of cases, a very satisfactory agreement with reference solutions is observed.The 3D spreading droplet case using Adaptive Mesh Refinement (AMR) also shows satisfactory results overlapping with the analytical solution.
For some configurations (tilted plane, cylinder) and for some values of the contact angle, a pinning of the contact line has been identified.The alignment of the Cartesian grid with the solid leads in certain situations to a zero VOF flux advection resulting in a pinning of the contact line.This is an inherent problem of the method but one could imagine different solutions to overcome pinning situations, for instance imposing a slip length at the solid to enforce the motion of the contact line.
The dynamics of the spreading droplet have also been studied using either Dirichlet or Navier slip boundary conditions on the embedded boundary.Even if a time convergence has been observed in both configurations, a significant difference on the grid convergence is emphasized, as shown in previous studies in the literature.We showed that grid convergence can be obtained only using a Navier slip condition as long as the numerical slip length is well resolved by the grid.
A perspective of this study would be to test different contact angle models in the literature since the proposed algorithm is also compatible with a contact angle varying in space and time.A limitation, common also in other methods, is the case of "shallow" angles (θ s < 10 • or θ s < 170 • ) which would require special treatment.
The implementation of the numerical methods and test cases are freely accessible as part of the Basilisk open-source library (see http://basilisk.fr/sandbox/tavares/test_cases/ ).

Figure 1 :
Figure 1: Stencils used in Basilisk and their respective indexing on a 2D grid: (a) Vertex Φ and face velocity ūf stencils, (b) Cartesian grid cut by the discrete rigid boundary to form irregular control volumes

Figure 3 :
Figure 3: Schematic representation of the Taylor-Couette flow with Navier, free slip or Dirichlet boundary condition imposed at the outer fixed cylinder.

Figure 4 :
Figure 4: Comparison between the numerical and analytical angular velocity u θ as a function of the radius r, N = 32, (a) Free-slip and no-slip cases, (b) Navier slip case with different values of λ = 0.4, 0.5, 0.8 and 1. See also http://basilisk.fr/sandbox/tavares/test_cases/slip.c

Figure 5 :
Figure 5: Rigid boundary (Ω S ) embedded in a two-phase flow (Ω L ), (Ω G ): (a) Continuous solid, liquid and gas interfaces characteristics, (b) Volume fraction F and C in the discrete solid liquid-gas system (without intersection), (c) Continuous solid, liquid and gas interfaces with intersection, (d) discrete solid, liquid and gas system with intersection.

Figure 6 :
Figure 6: Geometric estimation of the flux G: (a) in a two-phase cell, (b) in a mixed (or threephase) cell.

Figure 7 : 1 1 f o r a l l c e l l s 2 i
Figure 7: Contact angle algorithm description: (a) The normal nθ is reconstructed following Eq.(51) (b) The volume fraction in the ghost cells F G is calculated (Tab.4), and the boundary condition at the contact line for the prescribed value of θ s is eventually imposed.α θ = F(F, nθ ) can be determined following the procedure described in section 3.2.1.The Algorithm 1 1 f o r a l l c e l l s 2 i f ( 0 < C < 1 and 0 < F < 1 ) then // p ot e n t i a l t r i p l e p h a s e c e l l

5 end i f 6 end f o r Table 3 :
Interface geometric characteristics calculation at the triple point.

Figure 10 :
Figure 10: Convergence of the equilibrium droplet shape with respect to the grid refinement for θ s = 30 • and 120 •

Figure 11 :Figure 12 :
Figure 11: (a) Dimensionless droplet radius R/R 0 evolution for the equilibrium shape of a droplet on a cylinder for N = 64, (b) Relative error on the final droplet radius for θ s = 30 • , 90 • , 150 • .

Figure 15 :
Figure 15: Dimensionless droplet height h f /R 0 and radius r f /R 0 evolution for the equilibrium shape of a droplet on a flat surface without gravity Eo = 0.

Figure 16 :
Figure 16: Average error on the volume fraction E F with different grids refinement, (a) Horizontal plane, (b) 45 • Inclined plane.

Figure 17 :Figure 18 :
Figure 17: Relative mass absorption for the sessile droplet, for different contact angles θ s for N = 64 for the two different orientations of the substrate: (a) Horizontal plane; (b) 45 • Inclined plane.

Figure 21 :
Figure 21: Evolution of the radius r as a function of the dimensionless time τ for different values of λ, N = 32.

Figure 22 :
Figure 22: Grid convergence using the no slip boundary condition λ = 0 (left) and the Navier slip boundary condition λ = 16∆ (right): Evolution of the droplet radius as a function of the dimensionless time τ (a) no slip, (b) Navier slip, Capillary number evolution (c) no slip, (d) Navier slip, Velocity gradient dū/dn at τ = 0.5 (e) no slip, (f) Navier slip.

Figure 23 :
Figure 23: A drop falling at Re = 175 and θ s = 15 • is captured by the fiber.

Figure 24 :
Figure 24: A drop falling at Re = 830 and θ s = 15 • detaches from the fiber.

Table 2 :
Navier slip case with λ = 0.5 and 1: Error E 1 , E ∞ and order of convergence.

Table 4 :
Fluid volume fraction calculation inside the ghost cells.
• , 60 • , 90 • , 120 • and 150 • have been studied.The comparison between the numerical and the analytical equilibrium droplet shapes

Table 5 :
• and two different impact velocities were studied.Here, the drop has a diameter of D and an impact velocity of U , ρ l /ρ g µ l /µ g Non dimensional parameters for the 3D drop impacting a circular fiber.