A numerical model of two-phase flow at the micro-scale using the volume-of-fluid method

Article history: Received 3 May 2017 Received in revised form 17 December 2017 Accepted 19 December 2017 Available online 21 December 2017


Introduction
Understanding multiphase flow through porous media is of great importance in a variety of environmental, industrial and engineering applications such as contaminant cleanup [35], fluid separation in fuel cells [48], enhanced oil recovery [29] and carbon dioxide storage in geological porous media [32]. However, accurate modelling and quantification of such flows at the pore level is challenging when interfacial tension effects become dominant. Inaccurate numerical modelling of the interfacial force may upset the balance of forces in the momentum conservation equation and lead to instabilities and unphysical results [54,16,17].
Popular approaches to simulate multiphase pore-scale processes in porous media include pore-network models [13,7,6,53], lattice Boltzmann models [47,40,2], mesh-free Lagrangian particle approaches such as smoothed particle hydrodynamics [36,49], and grid-based methods such as finite difference, finite volume or finite element. The grid-based methods are used in conjunction with various algorithms to model fluid-fluid interfaces, such as front-tracking [18,41], volume-of-fluid [22,42], level-set [39,1] and phase-field [4,3] models. Pore-network models are computationally efficient, but they need to simplify the geometry and physics of the problem. Lattice Boltzmann methods have the advantage of not requiring explicit interface tracking or contact angle models, but numerical instabilities may arise when simulating multiphase flows with large density and viscosity ratios [10]. Lagrangian mesh-free methods are robust, simple to implement and able to handle com-plex interface motions. However, they are not always computationally efficient when dealing with large three-dimensional systems [56].
Grid-based methods are conventionally considered as the preferred approach for most applications, due to their computational efficiency and ability to simulate fluid flow with large density and viscosity ratios [34]. One major challenge in modelling flows with moving interfaces in grid-based methods, however, is the need for an algorithm to model the interface including its location, shape and evolution with time. In front-tracking approaches, such as the height function method [23], the interface is determined explicitly as a sharp boundary and its motion is tracked on the computational domain. The main advantage offered by front-tracking methods lies in their capability to resolve the position of the interface to a high accuracy. However, explicit tracking of the motion of interfaces that undergo large distortions can be a challenging problem. Another class of algorithms represent the interface implicitly by marking the fluids on both sides of the interface, using either virtual particles (marker-and-cell methods) [21,33], or an indicator function such as a volume fraction (volume-of-fluid methods), a level-set function (level-set methods) or an order parameter (phase-field methods). One key advantage of these methods as opposed to front tracking is that they do not require complicated interface tracking algorithms. This is important for modelling two-phase flows through complex geometries such as three-dimensional micro-CT images of porous media in which we usually have large interface motions and interactions. Phase-field models flexibly capture complex interface motions and related physical phenomena by describing the interface as a thin transition region based on thermodynamic arguments. However, resolving the diffusive interface while convecting it in a mass-conservative fashion can be numerically challenging and computationally expensive [31]. Volume-of-fluid methods, unlike level-set methods, are intrinsically mass conservative. In addition, they are more computationally efficient compared to marker-and-cell methods [34]. These advantages, in addition to flexibility for treating complex interface configurations, have made volume-of-fluid a popular and powerful tool for the direct numerical simulation of immiscible two-phase flow [37,5,52,44,14]. In this work, we use a volume-of-fluid approach to capture interfaces.
Different methods have been employed to compute interfacial forces using implicit data from interface capturing approaches such as the volume-of-fluid method on a fixed grid. Continuum models such as the continuous surface force (CSF) method [9] and the continuous surface stress (CSS) method [19] are among the most widely used techniques that treat the singular interfacial force as a body force and concentrate it in the region of the interface using a delta function. The main problem with these models, however, is that they suffer from spurious currents near the interface especially when it comes to simulating flows in which the capillary force is dominant over viscous and gravity effects [28,26]. Another fundamental drawback is that, due to the lack of a sharp boundary between different fluids, these methods may not be able to ensure a zero net capillary force on closed interfaces [50]. The adverse effects of failing to impose this force conservation property have been reported in modelling of free rising bubbles in quiescent liquids and moving droplets in uniform velocity fields [45,42].
These issues have been addressed in many studies. Tryggvason et al. [50] employed a front-tracking method to reconstruct a sharp interface from volume fractions to compute a more accurate interfacial force on fixed grids. Scardovelli and Zaleski [46] used an interface tracking algorithm based on the piecewise linear interface construction method (PLIC) and a Lagrangian advection scheme to balance interfacial forces more accurately and reduce spurious currents to machine accuracy for two-dimensional calculations. Renardy and Renardy [44] developed a parabolic reconstruction of surface tension (PROST) algorithm for three-dimensional cases to model the volume force due to interfacial tension, which effectively reduced spurious currents in comparison to CSF and CSS. Jamet et al. [26] implemented an energy conserving discretization scheme, similar to the phase-field implementation of Jacqmin [25], to achieve a physical equilibrium state without spurious currents. He argued that conserving energy to round-off errors makes it possible to model fluid interfaces accurately, eliminating non-physical currents to machine accuracy, even if momentum is not strictly conserved.
In this paper, we use a face-based reconstruction algorithm to define the fluid-fluid interface as a sharp iso-contour surface to eliminate spurious currents effectively. Moreover, interfacial forces are computed directly on the reconstructed interface elements providing a force conservative interfacial tension model. We use interFoam, a volume-of-fluid based solver, in OpenFOAM [38] as the basis for implementing our numerical method. A modified fully sharpened indicator function is used to present the fluid-fluid interface as a sharp discontinuity [42]. A simple and effective algorithm is devised to reconstruct the interface from volume fractions only at the cell faces for which the gradient of the sharpened indicator function is non-zero. The reconstructed interface is then employed to compute the interfacial force at the centre of cell faces to give a balanced-force algorithm [16]. A primary advantage of this method is its ability to suppress spurious currents effectively and improve the force conservation property of the interfacial tension model even for low mesh resolutions. Herein, the term contour-level surface force (CLSF) method is used to refer to our approach for the treatment of the interfacial force.
In section 2, we present the governing mathematical equations. In section 3, we present the numerical implementation of the governing transport equations, the interface reconstruction algorithm and the interfacial force model in detail. The last section is dedicated to verification studies to ensure the correctness of the proposed numerical approach. The method is validated using several test cases. The process of equilibration of a static droplet in two and three dimensions is simulated to demonstrate the reduction of spurious currents (Section 4.1.1). The accuracy of the numerical method when liquid and solid phases come into contact is validated by modelling a wetting layer resting in the corner of a two-dimensional pore geometry with different corner and equilibrium contact angles (Section 4.1.2). The ability of the model to respect zero net capillary force on closed interfaces is investigated by simulating a micro-scale moving droplet in a uniform velocity field in two and three dimensions (Section 4.2.1). The dynamics of three-phase contact lines is examined by modelling the steady movement of a contact line through two-and three-dimensional micro-channels (Section 4.2.2). Finally, the convergence behaviour of the method with spatial refinement is tested for corner flow of a wetting layer in a square capillary tube (Section 4.3). Here we demonstrate how co-current flow can lead to a greatly enhanced flow conductance when the non-wetting phase has a finite viscosity, as compared to the case where the non-wetting phase is an inviscid gas.

Mathematical model
The equations of motion for an isothermal, incompressible two-phase flow of Newtonian fluids, can be written using a single-fluid continuum approach as follows:

∂(ρu) ∂t
where u is the velocity vector, p is pressure, g is the gravitational acceleration, f σ v is the interfacial force concentrated on the interface, and ρ and μ are the fluid density and viscosity, respectively, defined as Here α is an indicator function (physically representing the volume fraction of one phase) that has a value of one in the first fluid (subscript 1) and zero in the second fluid (subscript 2). In the volume-of-fluid method, α represents volume fractions of one of the two fluids and varies smoothly between 0 and 1. The interface can be defined as the surface where the value of α is equal to 0.5. The indicator function is evolved using the following advection equation: The interfacial force per unit area, f σ , at any point on the interface, assuming the interfacial tension, σ , is constant, can be obtained as [51] f σ = σ κn, where for an infinitesimal volume element dV enclosing an interface dS bounded by a closed contour C (Fig. 1). Here κ is the local interface curvature, n is the local unit vector normal to the interface, d is the length of an increment along the perimeter of the contour line, C , and the vector m = t × n lies on the interface and is perpendicular to both the normal vector, n, and the unit vector t tangent to the interface along C . Finally, the interfacial force per unit volume, f σ v , can be computed as where δ s is a delta function concentrated on the interface, i.e. is non-zero only at the interface points.
In the single-fluid approach, the fluid-fluid interface is treated as a discontinuity in fluid properties rather than an explicit boundary. Therefore we do not need to prescribe kinematic and dynamic boundary conditions explicitly at the fluid-fluid interface [15]. In our study, apart from bulk phase viscosity, μ, Eq. (3), we do not consider any additional hydrodynamic resistance to flow at the fluid-fluid interface, i.e., zero interfacial shear viscosity [20].
In contact line regions, where the fluid-fluid interface joins a solid boundary, wetting effects must be taken into account. This is achieved by modifying unit vectors normal to the interface at the solid boundary in accordance with a prescribed contact angle, θ (see Section 3.2). The contact angle is an input parameter and, in our numerical tests, it is set equal to the equilibrium contact angle.

Numerical method
The finite volume method, with a collocated arrangement of variables at centres of discrete non-overlapping control volumes [55], is used to numerically solve the governing equations of flow, Eqs. (1), (2) and (4). The flow algorithm in conjunction with a volume-of-fluid based interface tracking method, is taken from the interFoam solver of OpenFOAM [38,27,45]. In this section, first the algorithm employed to solve and couple the transport equations is briefly described for the sake of completeness, see Deshpande et al. [11] and references therein for more details. Then, we present a new numerical algorithm to reconstruct the location of the interface. Finally, the computation of the interfacial force on the reconstructed interface and its transfer to the underlying fixed mesh are discussed.

Numerical formulation of transport equations
The linearized semi-discrete form of the momentum equation can be written as where A d is a scalar field containing the diagonal elements of the coefficient matrix, corresponding to a system of linear algebraic equations resulting from the discretization of the momentum equation, Eq. (2), the term H(u) accounts for off-diagonal elements of the coefficient matrix and also all source terms apart from pressure gradients and body forces, and F is any body force (in our case, only interfacial forces). A first-order, forward Euler method is used for the temporal discretization, while second-order finite volume schemes, based on Gauss's theorem, are used for the spatial discretizations. The bounded self-filtered central differencing (SFCD) scheme is used to interpolate velocity and the indicator function from cell centres to face centres in the convective terms. Central differencing is used to compute required gradients normal to the cell faces.
By rearranging Eq. (8), we have Imposing the continuity condition, Eq. (1), upon this velocity field, after discretization, leads to the following equation for the pressure: Here, φ f can be considered as an intermediate face-centred flux field that lacks the pressure gradient effect, S f and S f are, respectively, the outward-pointing area vector of face f and its magnitude, f denotes the face-centred field interpolated from the corresponding cell-centred field, ∇ ⊥ f denotes face normal gradient evaluated as a scalar field directly at the face centres, and φ F , f = F f .S f is the interfacial force flux computed directly at the face f , which is the focus of discussion in Once the pressure field is obtained, a corrected face-centred flux field, φ f , is computed using the following equation: Finally, the cell-centred velocity vector for each cell c, u c , can be reconstructed in a conserved way from its face fluxes, φ f , as follows: where S −1 is the inverse of the following second rank symmetric tensor: Here f ∈c is a sum operator over all faces belonging to the cell c. The major problem with advecting the indicator function using Eq. (4), is numerical diffusion, which tends to smear the interface sharpness through time. To alleviate this problem, an additional artificial convection term is added to the advection equation: ∂α ∂t where C α is a compression factor, set to 1 in our simulations, and u r is the relative velocity vector [45,11]. The fluid-fluid interface is advected by solving the system of linear algebra arising from the discrete form of Eq. (15) for each cell c: Here ∂ t represents the first-order Euler method for the temporal discretization, v c is the volume of cell c, φ f is the volumetric flux at face f computed using Eq. (12), and φ r, f is a relative volumetric flux at face f defined as where n f = S f /S f is the unit normal to the face f , and n s is the unit normal to the interface defined as n s = ∇α |∇α| .
Eqs. (10), (12) and (16) are solved iteratively in two loops: (a) an outer-corrector loop to semi-implicitly couple the interface evolution and consequently the interfacial forces with the velocity field, and (b) an inner-PISO loop, to couple pressure and velocity fields, based on the pressure implicit with splitting of operator (PISO) method of Issa [24]. The solution algorithm from time step n to time step n + 1 can be summarized as follows: ii. Solve the pressure equation to update the pressure field (Eq. (10)). iii. Correct the volumetric flux field (Eq. (12)). iv. Reconstruct the cell-centred velocity field (Eq. (13)). e. End of inner-PISO loop.

End of outer-corrector loop.
We use two iterations in both the outer-corrector and the inner-PISO loops, n outerLoop = n P I S O = 2, in our simulations.

Reconstructing the interface
Our motivation in reconstructing the interface is to reduce spurious velocities and calculate the pressure jump across the fluid-fluid interface more accurately within a finite-volume volume-of-fluid framework. The proposed algorithm is facebased, robust and easy to implement. We define the interface as face-based surface elements that consist of connected points where the value of α is equal to 0.5. For the purpose of constructing a sharp interface, we first introduce a new indicator function by clipping α and rescaling it between 0 and 1 as follows [42]: where the clipping factor, C c , controls the sharpness of α c in the vicinity of α = 0.5. In this study, this factor is set to 0.98.
Then, we declare a cell face to be a mixed-face, f mix , when for that face the following condition is satisfied (Fig. 2): where is a small number set to 10 −8 and h is the face normal distance between the centroid of the two cells that share the face f . This technique converts the smooth profile of the indicator function, α, into a sharp discontinuity, α c , without distorting the location of α = 0.5. The sharpness of the discontinuity can be controlled using the clipping factor, C c . As the clipping factor approaches 1, the discontinuity becomes sharper and the profile of α c becomes closer to a stair-stepped approximation. Employing the face normal gradient of α c , instead of α, to mark mixed-faces, i.e. the only faces on which the interfacial force is non-zero, increases the interface locality. Consequently, this improves the balance of the interfacial forces with the pressure gradient which can effectively eliminate spurious currents, as shown for static test cases in Section 4.
The smooth profile of α may get distorted due to the artificial compression of the interface during the advection process, Eq. (16), especially for dynamic test cases. This can adversely affect the accuracy of the interface normal vector computations which require a smooth change in the gradient of the indicator function. Therefore, we define a smoothed indicator function, α s , by interpolating α from cells to faces and back to faces, recursively, for four iterations (i = 0-3): where c denotes the cell-centred field interpolated from the corresponding face-centred field, and C s is the smoothing relaxation factor, which is set to 0.6 in our simulations. This smoothed indicator function is then used to obtain the interface normal vectors. This step may also include an optional prior clipping of α to exclude localized unphysical variations in α close to 0 or 1 in the vicinity of the interface zone, where C c,0 is a constant controlling the length of the clipping that is set to 1.02.
A second order least-squares method is used to estimate the gradient of α s at mesh points, p, which form the mixedfaces: where G −1 p is the inverse of the following symmetric tensor: and ( α s ) i = α s,i − α s,p . (25) Here N f is the number of faces that share the point p, d i is the distance vector between the centroid of face i and the location of point p, ω i is the inverse distance weighting factor calculated as ω i = 1/|d i |, α s,i is the interpolated value of the smoothed indicator function, α s , at face i, and α s,p is the interpolated value of α s at the point p, computed using the following inverse-distance-weighting interpolation method: Using a least-squares method to estimate (∇α s ) p is preferable because it only uses the values of immediate neighbour cells of point p, hence preserving its locality to a greater extent than the wide stencil finite difference approximations at cell centres followed by interpolation. Localized estimations of (∇α s ) p are crucial to compute the interface location, α = 0.5, accurately.
Finally, the interface unit normal for each corner point p is defined as follows: For the points on solid boundaries, fluid-fluid-solid contact points, however, the direction of n p must be modified in accordance with the contact angle, θ , to obtain interface unit normal vectors at the contact line, n p | w (Fig. 3). The interface normals at the contact line, n p | w , are obtained as a linear combination of n p and the unit normal vector to the solid wall, where θ = cos −1 (n w · n p ). (29) Here n w is obtained by interpolating face unit normal vectors of the solid boundary to the mesh points. In the case of stair-like solid boundaries, used in Section 4.1.2, n w is smoothed by recursively interpolating from the mesh points to their adjacent solid boundary faces and back to the mesh points [42], four times in this work. For every mixed-face, an interface element is considered so that a one-to-one correspondence is established between the points of the mixed-face and the interface points. The location of the interface points, where α = 0.5, is calculated from a simple linear extrapolation (see Fig. 2(c)) where x int,p is the position vector of the interface points, and x p is the position of point p on the mixed-face. Once the locations of the interface points are computed, the area vector of each interface element, S int, f , is calculated using the following sum over its points: (31) where N p is the total number of the corner points that form the interface element, and is the centroid of the interface element. The interface element points are ordered in a way that the computed area vector always points towards the same fluid (fluid 1). Employing the least-squares method to estimate the interface normal at corner points alongside the one-to-one correspondence between corner points of interface elements and mixed-faces results in a simple and robust face-based methodology for reconstructing the interface. All that is required is a set of discrete marked faces described by their corner points. We stress that the clipped and smoothed indicator functions, α c and α s , are only used to reconstruct the interface elements, without affecting the original indicator function, α. Moreover, the indicator function, α, is advected in a conservative control-volume framework that guarantees mass conservation to round-off in all test cases; see the test case in Subsection 4.2.2, for an example.

Computing the interfacial force
In this section, we introduce our algorithm for computation of interfacial forces, the contour-level surface force (CLSF) model. Then, we briefly describe the previously-published continuous surface force (CSF) [9] and sharp surface force (SSF) [42] models. In Section 4, we show a comparison of the accuracy of the CLSF, CSF and SSF algorithms.

Contour-level surface force (CLSF) model
In this work, the interface is reconstructed and explicitly represented as discrete elements. The interfacial force for each interface element, f int,σ , can be written as an integral of f σ , Eq. (5), over the surface of the element and converted to a line integral using the Stokes theorem [50]: where C int is the closed boundary along the edges of the interface elements. The interfacial force for each element is computed using this line integral over the edges of the interface element. This offers the unique advantage of not needing to calculate the interface curvature at cell centres and then interpolating them to faces, which improves the balance of forces in the momentum equation [50].
On the discrete level, the interfacial force on an interface element corresponding to face f , f σ ,int, f , is obtained by estimating Eq. (32) using the following sum over the corner points of the element: The computed interfacial forces are then smoothed by looping over corresponding mixed-faces to equally distribute f σ ,int, f between the two grid cells that share the mixed-face, and then interpolate the cell values back to the mixed-face based on an area-weighting method, where and Here, ω i is the area-based interpolation weighting for cell i, N f mix is the number of mixed faces belong to cell i, and | α c | f is the magnitude of the difference of α c between the two cells that share the mixed-face f . This smoothing step is essential to keep the conservation of the net capillary force on closed surfaces within an acceptable limit for our dynamic test cases (see Section 4.2), when applying the interfacial forces onto the underlying fixed grid.
Finally, the interfacial force flux term, φ F , f = F f · S f , Eq. (11), is approximated on the fixed grid faces as follows: where n s is the interface normal computed using Eq. (18). The steps taken in the CLSF method to compute the face-based interfacial force at each time step, can be summarized as follows: 1. Compute a sharp indicator function, α c (Eq. (19)). b. Obtain the interface normals at the vertices using the computed gradients (Eq. (27)). Modify the normal vectors on solid boundaries to include wettability effects (Eq. (28)).
c. Estimate the location of the fluid-fluid interface, where α = 0.5, in the direction of the interface normal vectors (Eq. (30)). d. Form interfacial surface elements using the computed interface location coordinates by holding a one-to-one correspondence between the mixed-faces and their interface elements (Eq. (31)). e. Compute the interfacial force on the interface elements (Eq. (33)). 6. Smooth the computed interfacial forces (Eq. (34)). 7. Use the computed interfacial forces to estimate the interfacial force flux on the corresponding mixed-faces on the fixed grid (Eq. (37)).

Continuous surface force (CSF) model
In the CSF model [9], the interfacial force is estimated as where δ s is a delta function concentrated on the interface, n s is the unit normal interface, and κ = −∇ · n s is the interfacial curvature.
On the discrete level, to have a balanced-forced flow algorithm similar to Francois et al. [16], Eq. (38) is used to approximate the interfacial force flux term, φ F , f , on the faces of the fixed-grid, as follows: where interfacial unit normals, n s , are computed at cell centres using Eq (18), and δ sf is approximated as

Sharp surface force (SSF) model
In the SSF model, the same approach is employed as the CSF model, Eq. (38), however, the delta function in Eq. (39), δ sf , is approximated using the normal gradient of the sharp indicator function, α c , Eq. (19), rather than α [42]: This implies that the SSF method employs a more compact delta function in comparison to the CSF method, while the interfacial normal vectors, n s , and consequently the interfacial curvature, κ, are estimated with a wide stencil as in the CSF method. We show that even though this approach may help to reduce spurious currents in some static test cases, it has the same limitations as CSF: large errors in pressure jump and consequently an unbalanced (non-zero) interfacial force on closed surfaces [45].

Model validation
In this section we use several two-and three-dimensional test cases to validate our numerical method. We first consider examples where the two phases are in static equilibrium (Section 4.1). Then, the accuracy of the method is examined using dynamic test cases (Section 4.2). Finally, to explore applications in porous media, corner flow through a three-dimensional non-circular capillary tube is investigated (Section 4.3).
In all simulations presented herein, for liquid-liquid systems we use two liquids, both with viscosities of 10 −3 Pa·s and densities of 1000 kg/m 3 and for gas-liquid systems a gas phase with a viscosity of 10 −5 Pa·s and a density of 1 kg/m 3 .
The interfacial tension is constant and taken to be σ = 0.03 N/m. We compare the results of the CLSF model to analytical solutions and also the CSF and SSF models.

Static test cases
A common test to validate two-phase flow models with interfacial tension, is simulating fluid-fluid systems at static equilibrium conditions where the velocity field must remain zero [41,44]. The discrete balance of the pressure gradient with the interfacial force can be examined by measuring the value of the maximum velocity in the computational domain. In addition, the accuracy of the interfacial force model can be quantified by measuring the error in pressure jump across the interface.

Static droplet
First, we simulate a two-dimensional circular droplet in the absence of gravity as either (a) a gas-liquid system (gas is surrounded by liquid) or (b) a liquid-liquid system. A static droplet with a radius of R = 10 μm is centred in a rectangular domain having side lengths of L = 40 μm. The computational domain is discretized using a uniform Cartesian grid with a resolution of R/δx = 8, where δx is the side length of uniform grid cells. Since the interfacial force is treated explicitly, small time steps, δt, are taken to satisfy the time step constraint δt ≤ where ρ = ρ 1 +ρ 2 2 [9].
where μ and ρ are, respectively, viscosity and density of the liquid, and D is the diameter of the droplet.
It can be observed that Ca reduces to machine precision, 10 −15 , in both systems using the CLSF method. For the SSF method, the velocity field also dies out to machine precision but with a slower convergence rate. On the other hand, for the CSF method, the maximum dimensionless error in the velocity field remains of order 10 −3 for the gas-liquid and the liquid-liquid systems, which obviously is unfavourable when using the numerical method to model multiphase flow at low capillary numbers.  Table 1. One can see that the structure and location of spurious vortices for the CSF method are different from those of the SSF and CLSF methods. For the SSF and CLSF methods, the spurious currents spread over a smaller area of the domain, due to employing a more compact delta function approximation than the CSF method, and are insignificant away from the interface to such an extent that they are effectively not present. This clearly shows the importance of implementing the interfacial force with a sharp force model to eliminate spurious currents efficiently.
The error in velocity is measured using the following L ∞ error norm: while the relative error of the pressure jump is estimated as The pressure jump across the droplet, according to the Young-Laplace equation, is P exact = σ /R = 3000 Pa.  Table 1 Effect of mesh resolution on the error in maximum velocity, L ∞ (u), Eq. (44), and pressure jump across the interface, E( P ), Eq. (45), for the liquid-liquid two-dimensional static droplet at a dimensionless time t d = 0.25, Eq. (43), using the CSF, SSF and CLSF methods. In Table 1, we investigate the results for the error in maximum velocity, L ∞ (u), and pressure jump, E( P ), Eqs. (44) and (45) respectively, as a function of mesh resolution for the two-dimensional static droplet (liquid-liquid system) at t d = 0.25. These results show that the CLSF method can predict the pressure jump accurately, with about 1% error, even at low resolutions, R/δx = 4 and R/δx = 8, while for the CSF and SSF methods the error is about 16% for these cases. This accuracy at low resolution is essential for modelling large problems and especially for three-dimensional calculations where the resolution may be low due to constraints imposed by high simulation costs. Further refinements to R/δx = 16 and R/δx = 32 lead to further reduction of the error in pressure jump for all methods. However, the magnitude of the error is still high, above 10%, for the CSF method. For the SSF method the error in pressure jump decreases dramatically from about 16% to 5% as the interface is resolved to R/δx = 16, whereas for the CSF method the reduction in the error is not that significant. It can also be observed that mesh refinement does not decrease the L ∞ norm of the spurious currents significantly in the CSF method, while for the SSF and CLSF methods spurious currents reduce to very small values, O (10 −11 ).
Next, we consider a three-dimensional droplet of radius R = 10 μm positioned in the centre of a uniform cubic mesh of resolution R/δx = 8. The system is taken to be the liquid-liquid system and the theoretical pressure jump across the interface is P exact = 2σ /R = 6000 Pa. Fig. 6 shows a visualization of the final shape of reconstructed fluid-fluid interface used to compute interfacial forces. Note that each interface element corresponds uniquely to one mixed-face on the fixed grid. Fig. 7 depicts the evolution of the dimensionless maximum velocity for the three-dimensional test case using the CSF, SSF and CLSF methods. The spurious current capillary number is reduced efficiently to a very small value, O (10 −13 ), with the SSF and CLSF methods. On the other hand, similar to the two-dimensional case, Ca persists at the order of 10 −3 , for the CSF method. The pressure profile for a cross-section passing through the centre of the droplet for the three methods are shown in Fig. 8. Improvements in both the accuracy and profile of the computed pressure with the CLSF method are evident. It yields a sharp pressure profile with a relative pressure jump error of 2%, while the CSF method leads to a smooth pressure distribution with a high relative pressure jump error of approximately 12%. For the SSF method, the pressure profile is not quite sharp, some transitional grid cells can be observed in comparison to the CLSF method, and the pressure jump error is 8%.
These results suggest that the CLSF method can predict the sharp pressure jump across the interface accurately, even for low resolution grids. In addition, it can virtually eliminate spurious currents in a stable manner for both the two-and three-dimensional static droplets.    to obtain an analytical equation for the wetting layer area, A l . It is obtained by subtracting the grey area from the shaded area, Eq. (48). a = R l sin η, b = R l cos η/ tan η, h = R l cos η and η = θ + γ . Corner angle, 2γ , is the same for all corners. θ is the equilibrium contact angle and L is the characteristic length (fixed and equal to 20 μm) for the geometry.

Static contact line
A two-dimensional domain with a star-shaped solid boundary is considered in this subsection. The domain is partially filled with a non-wetting phase while a wetting phase with the same fluid properties resides in the corners (Fig. 9). The grid cells which lie in a circle of radius 40 μm, centred in the computational domain, are initialized with the non-wetting phase while the rest (corner cells) contain the wetting phase. The aim is to test the ability of our model to simulate static immiscible fluid-fluid interfaces in contact with staircase-like solid walls, on a uniform mesh, with different corner and equilibrium contact angles, 2γ and θ , respectively. The accuracy and stability of the method are quantified in terms of the infinite error norm in velocity field and the normalized error in the pressure jump across the interface, where P exact can be calculated from where R l is the radius of the interface curvature in the corners. The wetting phase cross-sectional area, A l , as a function of R l (Fig. 7(c)) is: where η = θ + γ . The layer area, A l , is estimated numerically and plugged in Eq. (48) to obtain R l , which is in turn used in Eq. (47) to obtain P exact . Table 2 lists the simulation results for the CSF, SSF and CLSF methods at t d = 2.5, for different mesh resolutions, L/δx, where L is the characteristic length of the geometry, see Fig. 9, which is taken to be 20 μm. One can observe that the errors resulted from the CSF and SSF methods do not follow a consistent behaviour. The accuracy of the results for these two methods is unacceptably high, and apparently dependent on the values of the corner and equilibrium contact angles, which dictate the final configuration of the contact line with respect to the solid wall. On the other hand, the CLSF method reduces the spurious currents to machine accuracy and predicts the pressure jump across the interface to better than 8% in all cases. These results show that the CLSF method is able to model a static contact line accurately in a liquid-liquid-solid system even with a stair-case-like discretization of solid boundaries.

Dynamic test cases
To complete our numerical investigations, we study two sets of dynamic test cases including (a) a moving droplet in a uniform velocity field (in two and three dimensions), and (b) two-and three-dimensional moving contact lines in microchannels. The moving droplet test case is selected to test the ability of our method to conserve zero total force on a Table 2 Errors in maximum velocity, L ∞ (u), Eq. (44), and normalized error in pressure jump across the interface, E n ( P ), Eq. (46), for different mesh resolutions (L/δx), corner angles (2γ ), and equilibrium contact angles (θ ). closed interface as the indicator function is advected through the domain. This property is crucial, especially for capillarydominated flows, where even a small error in the computation of interfacial forces can lead to a non-physical solution. The moving meniscus test case is considered to check the accuracy and stability of the new method to model the dynamics of immiscible fluid displacement in the presence of a moving contact line.

Moving droplet in a uniform velocity field
We model the steady movement of a liquid droplet of radius R = 10 μm in a uniform velocity field and with zero gravity. Theoretically, the droplet should move with the same velocity as the underlying uniform velocity field, and the pressure jump across the interface is P = σ /R = 3000 Pa and P = 2σ /R = 6000 Pa for two and three dimensions, respectively. First, we investigate the effect of the velocity field, U flow , on the movement of the droplet in a two-dimensional domain with a uniform grid of resolution R/δx = 16. Fig. 11 depicts the predicted dimensionless droplet velocity, Ca droplet = U droplet μ/σ , for different dimensionless flow velocities, Ca flow = U flow μ/σ , as a function of a dimensionless time, t d : here U = U flow and D is the diameter of the droplet. It is clear that the CLSF method is able to predict the movement of the droplet accurately for capillary numbers as low as 10 −5 . Some relatively small fluctuations in the droplet velocity can be observed for Ca flow = 10 −4 and Ca flow = 10 −5 . Ca flow = 10 −2 is the lowest capillary number for which the CSF and SSF methods can predict the movement of the droplet accurately. For capillary numbers lower than 10 −2 , the error in the computation of interfacial forces in the CSF and SSF methods causes the violation of the zero net force constraint on the closed interface. This error manifests itself as large spurious currents which leads to the deviation of the droplet velocity from the theoretical value. Fig. 12 shows the initial configuration of the droplet and velocity vector plots and their computed values for Ca flow = 10 −4 at t d = 6 using different methods. One can observe that for the CSF and SSF methods the droplet is not able to move due to high spurious currents caused by error accumulation in front of the droplet. The related indicator function and pressure field are shown in Fig. 13 for the CLSF method. The pressure jump across the interface of the droplet accurately matches the analytical value of 3000 Pa. Next, we test the same moving droplet in three dimensions with a uniform grid of resolution R/δx = 16 for the CSF, SSF and CLSF methods. Fig. 14 shows the dimensionless average velocity of the droplet, Ca droplet , as a function of dimensionless time t d for three different velocities using the CSF, SSF and CLSF methods. Similar to the two-dimensional simulations, the CLSF formulation can model the movement of the droplet for two orders of magnitude lower capillary number (Ca flow = 10 −4 ) compared to that of the CSF and SSF methods (Ca flow = 10 −2 ). Fig. 15 shows snapshots of the velocity vector plot, the indicator function profile and computed pressure field through a centre plane at dimensionless time t d = 3. In the CLSF method the predicted pressure jump across the interface is highly accurate.
Finally, we compute the relative error in the average velocity of the two-and three-dimensional moving droplets, Ū droplet . Fig. 16 depicts the error in Ū droplet over a period time of t d = 2, as a function of mesh resolution, for different capillary numbers, The solution exhibits a positive convergence for relatively high capillary numbers, Ca flow = 10 −2 and Ca flow = 10 −3 , for both two-and three-dimensional test cases. However, the order of convergence decreases as the capillary number decreases and we need finer meshes to obtain a solution within an acceptable error range. For the two-and three-dimensional moving

Steady contact-line displacement through a micro-channel
One of the most important characteristics of immiscible two-phase flow through porous media is the presence of threephase contact lines where the three phases (fluid-fluid-solid) are in contact. Therefore, testing the ability of our numerical method to model contact line dynamics is crucial.
Here, a two-dimensional micro-channel of width D = 20 μm, as shown in Fig. 17, is considered to investigate the steady movement of a contact line. A no-slip boundary condition is applied at the solid walls while the velocity at the inlet is set to a uniform fixed value in the x-direction. We apply a zero-gradient boundary condition for the pressure on all the boundaries and for the velocity at the outlet. The solid boundaries are assumed to be perfectly flat and homogeneous having a fixed equilibrium contact angle of θ = 30 • . First, we investigate the contact line dynamics when the wetting phase displaces the non-wetting phase for several steady injection velocities, U inj , at the inlet, using a uniform grid of resolution D/δx = 32. A sample visualization of the velocity vector field and the shape of the interface (α = 0.5) at dimensionless time t d = 0.5, for the case when wetting fluid is injected at an injection capillary number of Ca inj = U inj μ/σ = 10 −3 , is shown in Fig. 17. The dimensionless time is defined based on Eq. (49) where U = U inj and D is the width of the channel. Fig. 18(a) shows the dimensionless average velocity of the interface, Ca interface , as a function of dimensionless time t d for a range of injection capillary numbers, Ca inj . The interface translates accurately and its average displacement velocity is the same as the injection velocity. This implies that the numerical method is mass conservative and the effect of spurious currents on the overall displacement of the interface is insignificant. However, in our simulations, we observe that the contact line does not translate perfectly smoothly and experiences a periodic stick-slip motion, which introduces some small velocity fluctuations. The error resulted from these numerical spikes is quantified by using the following L 2 error norm: and depicted as a function of dimensionless time t d in Fig. 18(b). The results suggest that the L 2 error norm increases as the capillary number decreases; however, the numerical errors do not introduce significant errors in the overall shape and movement of the interface. Moreover, Fig. 19 presents a mesh sensitivity analysis for the case with Ca inj = 10 −3 . It can be observed that the amplitude of the numerical spikes in the L 2 error norm decreases as the grid is refined.
Although the assumption of perfectly smooth solid walls can be useful to simplify the study of contact line dynamics, it should not be forgotten that solid surfaces in porous rocks are usually rough. Therefore, the roughness of the solid surface can lead to physical stick-slip behaviour at the micro-scale level. As shown in Fig. 20, the solid boundaries of the micro-channel are roughened using a triangular wave pattern characterized by roughness amplitude A = 225 nm and period P = 2 μm. The flow domain is discretized using the mesh generator from OpenFOAM [38] that uses a uniform initial mesh  and snaps the boundary points to the solid wall, here the triangular wave roughness, so that all the solid boundary faces align with the solid wall. We also ensure that each edge of boundary faces connect to only one internal face. Finally, we model the displacement of a meniscus through a three-dimensional square micro-channel, as shown in Fig. 22. We set the equilibrium contact angle to θ = 60 • and keep the other simulation parameters the same as in the two-   These results show that the CLSF method can be applied to model capillary-dominant two-phase flow in the presence of solid boundaries where the dynamics of the contact line plays an important role.

Application to porous media
We now study two-phase flow problems typically encountered in subsurface rock geometries at the pore scale, namely flow of a wetting phase along the corners of a rock surface [8]. The dynamics of the corner flow of the wetting phase, especially at low capillary numbers, plays an important role in multiphase displacement process and also in the trapping and remobilisation of the non-wetting phase [30,43,12]. Therefore we examine the capability of our method to model the imbibition of a wetting phase in the corners of non-circular capillary tubes.
A 60 μm square capillary tube with a length of L = 300 μm is considered (Fig. 24) to study the spatial convergence of the dimensionless flow resistance factor of the wetting layer, β [43], where μ and u x are the viscosity and axial average flow velocity of the wetting phase, respectively, ∂ p/∂ x is the pressure gradient in the wetting phase along the axial direction, and R l is the radius of the interface curvature in the corner, calculated based on Eq. (48), with γ = π/4 and θ = π/6 as discussed in Section 4.1.2.
The initial liquid-liquid interface is a circular arc of radius r 0 = 30 μm. The flow domain is discretized by a uniform Cartesian mesh of different resolutions. The no-slip boundary condition is applied at the solid walls. At the inlet and outlet boundaries, zero-gradient boundary conditions are applied for pressure and indicator function, α. Both fluids, with the same properties, viscosities of 10 −3 Pa·s and densities of 1000 kg/m 3 , are injected at the inlet boundary in the x-direction, using a cyclic boundary condition; the velocity profile of the outlet boundary is replicated on the inlet boundary, with an initially uniform velocity of U inj = 0.0003 m/s. In addition, no surface shear viscosity between the fluid phases and no gravity are considered. The interfacial tension is σ = 0.03 N/m and the equilibrium contact angle is θ = 30 • . In the results that follow, the simulations are conducted for the whole computational domain. However, the flow resistance is calculated and analyzed for only one corner of the tube. Fig. 25 depicts the evolution of β as a function of time for five mesh resolutions. Through refining the mesh, the computed β converges to a value of about 52. Fig. 25(b) depicts the error in β, at time t = 10 ms, defined as: where β D/δx=54 is the evaluated dimensionless flow resistance on a mesh of resolution D/δx = 54. Here we observe a positive convergence behaviour for β.
The computed dimensionless flow resistance in this test case, β = 52, is significantly lower than the value reported in Ransohoff and Radke [43] for the same geometry, β = 290.7. The main reason is that the dimensionless flow resistance in Ransohoff and Radke [43] is calculated for a gas-liquid system in which gas is assumed to be an inviscid non-wetting phase. Therefore, the stress exerted by the non-wetting phase on the wetting phase is ignored. In our example, the co-current flow of the non-wetting phase boosts the flow of the wetting phase, leading to a much lower flow resistance.
We simulate wetting layer flow on a mesh of resolution D/δx = 18 for a case where the non-wetting phase is a gas with viscosity μ = 10 −5 Pa·s and density ρ = 1 kg/m 3 . Fig. 26 compares the computed dimensionless resistance factor for the gas-occupied and liquid-occupied pores as a function of time. Clearly, through reducing the stress exerted by the non-wetting phase on the wetting phase in the gas-occupied system, β is increased to approximately 250, which is closer to β = 290.7 found by Ransohoff and Radke [43].
These results confirm that the numerical method can be used to model corner flow in non-circular capillary tubes to quantify their flow resistivity coefficients. However, the dynamics of corner flow can also be affected by fluid properties, here the viscosity of the non-wetting phase. A more extensive analysis of the effect of fluid and geometric properties on the corner flow will be the subject of further investigation.

Conclusions
A simple and efficient numerical method for modelling two-phase flow in porous media at the micro-scale has been devised. We compute interfacial forces using a sharp surface force model that can efficiently eliminate the problem of spurious currents. The explicit representation of the interface used to compute interfacial forces leads to an accurate estimation of capillary pressure even for low grid resolutions. This can be very important to model multiphase flows directly on micro-CT images of porous media for which simulations with a very fine grid resolution are prohibitively expensive using currently available computational resources. In our test cases we have shown, for instance, how the consideration of co-current flow of two viscous fluids leads to a much lower flow resistance in layers than computed previously assuming that the centre of the pore space was gas-filled. This opens up the opportunity for future work where high-resolution pore-scale simulations are used to quantify capillary pressure and flow conductance for input into larger scale flow models.