A Rotational Pressure-Correction Scheme for Incompressible Two-Phase Flows with Open Boundaries

Two-phase outflows refer to situations where the interface formed between two immiscible incompressible fluids passes through open portions of the domain boundary. We present several new forms of open boundary conditions for two-phase outflow simulations within the phase field framework, as well as a rotational pressure correction based algorithm for numerically treating these open boundary conditions. Our algorithm gives rise to linear algebraic systems for the velocity and the pressure that involve only constant and time-independent coefficient matrices after discretization, despite the variable density and variable viscosity of the two-phase mixture. By comparing simulation results with theory and the experimental data, we show that the method produces physically accurate results. We also present numerical experiments to demonstrate the long-term stability of the method in situations where large density contrast, large viscosity contrast, and backflows occur at the two-phase open boundaries.


Introduction
The current work focuses on the motion of a mixture of two immiscible incompressible fluids in a domain that is open on part of its boundary.The domain boundary is open in the sense that the fluids can freely leave or even enter the domain through such boundaries.In particular, we assume that the interface formed between the two fluids will pass through the open portions of the domain boundary.Therefore, the problem will involve truly two-phase outflow or open boundaries.
Two-phase outflows are encountered in many situations: oil plumes in the deep sea, wakes of surface ships, and ocean waves generated by the wind shear are some of the examples.These problems usually involve physically unbounded flow domains.Numerical simulation of such problems will therefore need to truncate the domain to a finite size, and some open/outflow boundary condition (OBC) will be required on the artificial boundary.The presence of fluid interfaces at the open boundary calls for appropriate two-phase open boundary conditions in such problems.
Several challenges are associated with the design of two-phase open boundary conditions.Some of the challenges are common to those encountered with single-phase outflows, for example, the instability associated with strong vortices or backflows at the open boundary (a.k.a.backflow instability) [12].Others are new and unique to two-phase outflows.For example, owing to the presence of fluid interfaces, two-phase outflow problems involve viscosity contrasts, density contrasts, and surface tension at the open boundaries.Large viscosity ratios and large density ratios at the open boundary can cause severe stability difficulties [10].
While outflow/open boundary conditions for single-phase flows have been under intensive investigations for decades, very scarce work exists for two-phase outflows or open boundaries.In [24] the zero-flux (Neumann), convective, and extrapolation boundary conditions originated from single-phase flows have been studied for the two-phase lattice-Boltzmann equation.The zero-flux condition is also employed for the outflow boundary in [1] within the context of a coupled level-set/volume-of-fluids method, and in [32] in the context of a level set method where the outflow boundary involves only a single type of fluid.The outflow condition for two immiscible fluids in a porous medium is discussed in [20].The works of [25,5] have both considered the outflow condition for two-phase compressible flows in one dimension.
In a recent work [10] we have proposed a set of effective two-phase outflow (and also inflow) boundary conditions within the phase field framework.A salient characteristic of these boundary conditions is that they ensure the energy stability of the two-phase system.By looking into the two-phase energy balance relation, we have shown, at the continuum level, that with these boundary conditions the total energy of the two-phase system will not increase over time, even in situations where there are strong vortices or backflows, large viscosity contrast, and large density contrast at the outflow/open boundaries.In [10] we have further developed an algorithm for numerically treating these open boundary conditions based on a velocity-correction type splitting strategy.
Inspired by the two-phase energy balance relation from [10], we propose in the current work several new forms of outflow/open boundary conditions for the two-phase momentum equations within the phase field framework.We also discuss a generalized form of the two-phase open boundary conditions, which includes these new forms and those forms proposed in [10] as particular cases.In addition, we present an algorithm for numerically treating these new open boundary conditions.Different than that of [10], the current algorithm is based on a rotational pressure-correction type strategy for solving the two-phase momentum equations.This algorithm shares one common property with those of [10,13], namely, after discretization it gives rise to linear algebraic systems for the pressure and the velocity that involve only constant and time-independent coefficient matrices, despite the variable density and variable viscosity of the two-phase mixture.Therefore, these coefficient matrices can be pre-computed during pre-processing.This makes the current algorithm computationally very efficient and attractive.
It is commonly observed that, with traditional splitting type schemes, the variable density in the Navier-Stokes equation has resulted in a variable (time-dependent) coefficient matrix for the pressure linear algebraic system after discretization [3,27,15,6,23,28,4,21].This creates a severe computational and performance issue, due to the need for the frequent re-computation of the coefficient matrix and the challenge in efficiently solving the resultant linear algebraic system at large density ratios.Guermond & Salgado [17] have advocated a penalty point of view toward the projection idea, leading to a Poisson type equation for the pressure; see also [18,30].Dong & Shen [13] have proposed a different strategy for coping with the variable density.By a reformulation of the pressure term in the variable-density Navier-Stokes equation, they have developed a scheme which requires the solution of a pressure Poisson equation with constant (time-independent) coefficient matrix; see also [8,9,10,7,11].
The novelties of this paper lie in two aspects: (i) the several new forms of outflow/open boundary conditions for the two-phase momentum equations, and (ii) the pressure-correction based algorithm for numerically treating these two-phase open boundary conditions.On the other hand, we would like to point out that the method for solving the phase field equation employed in the current paper is not new.It was originally developed in [13].
The numerical algorithm presented herein has been implemented using C 0 continuous high-order spectral elements [31,19,34] for spatial discretizations in the current paper.It should however be noted that the algorithm is general and can also be implemented with other spatial discretization techniques.
2 Pressure Correction Scheme for Two-Phase Outflows the constant densities of the two individual fluids, and use µ 1 and µ 2 to denote their constant dynamic viscosities.With the phase field approach, the motion of this mixture can be described by the following system of equations [22,33,13], where x and t are respectively the spatial coordinates and time, u(x, t) is the velocity, p(x, t) is pressure, D(u) = ∇u + ∇u T (the superscript T denotes transpose), f (x, t) denotes some external body force, and ⊗ represents the tensor product.φ(x, t) is the phase field function, −1 φ 1. Regions with φ = 1 denote the first fluid, and the regions with φ = −1 denote the second fluid.The function h(φ) is given by h , where η is the characteristic scale of the interfacial thickness.λ is referred to as the mixing energy density coefficient, and is given by [33] where σ is the surface tension and assumed to be constant in the current paper.The constant γ 1 > 0 is the mobility coefficient associated with the interface.ρ(φ) and µ(φ) are respectively the density and dynamic viscosity of the mixture, given by The function g(x, t) in (1c) is a prescribed source term for the purpose of numerical testing only, and will be set to g(x, t) = 0 in actual simulations.Equation (1c) with g = 0 is the Cahn-Hilliard equation.We assume that the domain boundary consists of three types which are non-overlapping with one another: We refer to ∂Ω i as the inflow boundary, ∂Ω w as the wall boundary, and ∂Ω o as the outflow or open boundary.
On the inflow and the wall boundaries, the velocity u is assumed to be known.In addition, the phase field function φ is also assumed to be known on the inflow boundary.
In the above equations (5a)-(5d), n is the outward-pointing unit vector normal to ∂Ω o , |u| denotes the magnitude of u, µ and ρ are respectively the mixture dynamic viscosity and density given by (3).The function F (φ) is given by , and note that λ 2 ∇φ • ∇φ + F (φ) is the free energy density of the two-phase system [22,33].f b is a prescribed function on ∂Ω o for the purpose of numerical testing only, and will be set to f b = 0 in actual simulations.Θ 0 (n, u) is a smoothed step function whose form is given by [12,10] Θ where U 0 is the characteristic velocity scale, and δ > 0 is a constant that is sufficiently small.δ controls the sharpness of the smoothed step function, and Θ 0 approaches the step function as δ → 0. When δ is sufficiently small, Θ 0 (n, u) essentially assumes the unit value where n • u < 0 and vanishes otherwise.The boundary conditions (5a)-(5d) belong to the following family of boundary conditions for ∂Ω o , where 0 θ 1, β 1 0 and β 2 0 are constant parameters.One can verify that the general form of open boundary condition (7), as δ → 0 and assuming f b = 0, is conducive to the stability of the two-phase energy balance equation given in [10].The boundary conditions (5a)-(5d) are particular cases of (7).For example, (5a) corresponds to (θ, β 1 , β 2 ) = (1/2, 0, 0) and (5d) corresponds to (θ, β 1 , β 2 ) = (0, 1, 0) in (7).This general form also contains the following two open boundary conditions as particular cases, which are proposed in [10] and can be obtained by respectively setting (θ, β 1 , β 2 ) = (1, 0, 0) and (θ, β 1 , β 2 ) = (0, 0, 0) in (7), The above boundary conditions on ∂Ω o are for the momentum equations (1a)-(1b).In addition to them, one also needs to supply appropriate boundary conditions on ∂Ω o for the phase field equation (1c).Note that two independent conditions will be needed on each boundary, due to the fourth spatial order of equation (1c).For the phase field function φ, on the outflow boundary ∂Ω o we will employ the boundary conditions developed in [10] where g a1 and g a2 are prescribed source terms on ∂Ω o for the purpose of numerical testing only, and will be set to g a1 = 0 and g a2 = 0 in actual simulations.The constant D 0 0 is a chosen non-negative constant, and 1 D0 plays the role of a convection velocity at the outflow boundary ∂Ω o .
The boundary conditions for the other types of boundaries (wall and inflow) will be set in accordance with previous works [8,10].We impose a Dirichlet condition for the velocity on the inflow and wall boundaries, where w is the boundary velocity.For the phase field function, we impose the following condition from [10] on the inflow boundary, where φ b denotes the distribution of the phase field function on the inflow boundary, and g b is a prescribed source term for numerical testing only and will be set to g b = 0 in actual simulations.On the wall boundary we employ the contact-angle condition of [8], considering only the effect of the static contact angle, where θ s is the static (equilibrium) contact angle formed between the fluid interface and the wall measured on the side of the first fluid, g c1 and g c2 are two prescribed source terms for the purpose of the numerical testing only and will be set to g c1 = 0 and g c2 = 0 in actual simulations.Finally, we assume that the following initial conditions for the velocity and the phase field function are known where the initial velocity u in and the initial phase field function φ in should be compatible with the above boundary conditions and the governing equations.

Two-Phase Momentum Equations: Algorithm and Implementation
The system of equations (1a)-(1c), the boundary conditions (9a)-(12b), and one of the conditions among (5a)-(5d), together with the initial conditions (13) for the velocity and the phase field function, constitute the overall system that need to be solved in numerical simulations.We next consider the numerical solution of this system.Because the phase field equation (1c) is coupled to the momentum equations (1a)-(1b) only through the convection term, it is possible and convenient to treat the momentum equations and the phase field equation individually.Indeed, by treating the convection term in (1c) explicitly, one can de-couple the computation for the phase field function from those for the momentum equations.On can first solve (1c) for the phase field function, and then solve the momentum equations for the pressure and the velocity.
In the following we will first concentrate on the momentum equations (1a)-(1b), together with the associated boundary conditions (5a)-(5d) for ∂Ω o and (10) for ∂Ω i and ∂Ω w .We defer the discussion of the solution to the phase field equation to an Appendix (see the subsequent Section 2.3 for detail).In subsequent discussions of this section we assume that the variables φ and ∇ 2 φ have been computed in appropriate ways and are already available.
To facilitate the following discussions we introduce an auxiliary pressure which will also be loosely called pressure where no confusion arises.Then equation (1a) can be transformed into We further re-write the open boundary conditions (5a)-(5d) or (8a)-(8b) into a unified compact form where The following algorithm is for the equations ( 15) and (1b), together with the boundary conditions (10) on ∂Ω i ∪ ∂Ω w and ( 16) on ∂Ω o .Note that the variables φ and ∇ 2 φ are assumed to be known here, as discussed before.
Let n denote the time step index, and (•) n denote the variable (•) at time step n.We use ũn and u n to denote two slightly different approximations of the velocity at time step n.Define ũ0 = u in , u 0 = u in , φ 0 = φ in .
By enforcing equation ( 15) at time step zero, one can compute the initial pressure P 0 as follows.Let and q ∈ H 1 p0 (Ω) denote the test function.By taking the inner product between ∇q and equation ( 15) and integrating by part, one obtains the weak form about P 0 , where ∂w ∂t 0 is the time derivative at time step zero, which can be numerically computed based on the second-order backward differential formula (BDF2) because the boundary velocity w(x, t) is known on ∂Ω i ∪ ∂Ω w .Ψ represents the projection of ∇ 2 φ 0 into the H 1 (Ω) space, and is given by the following weak form (ϕ denoting the test function), The weak forms (20) and ( 22) can be discretized in space using C 0 spectral elements (or finite elements).We solve equation (20), together with the Dirichlet condition to obtain the initial pressure P 0 , where Ψ is obtained by solving equation (22).Given (ũ n , u n , P n , φ n+1 , ∇ 2 φ n+1 ), where φ n+1 and ∇ 2 φ n+1 are assumed known and result from the algorithm for the phase field equation to be discussed later, we compute ũn+1 , u n+1 and P n+1 , together with an auxiliary field variable ξ n+1 , successively in a de-coupled fashion as follows: For ũn+1 : For ξ n+1 : For P n+1 : The notation employed in the equations (24a)-( 26d) is as follows.Let J (J = 1 or 2) denote the temporal order of the scheme, and χ denote a generic variable.Then in the above equations, χ * ,n+1 is a J-th order explicit approximation of χ n+1 , given by The expression 1 ∆t (γ 0 χ n+1 − χ) denotes an approximation of ∂χ ∂t n+1 by the J-th order backward differentiation formula, where ∆t is the time step size and In equations (25a) and (25b) G n+1 is given by The function E(ρ, n, u) is defined by equation (17).ρ n+1 and µ n+1 are given by equation (3) and by using φ n+1 .In equation (26d) µ min = min(µ 1 , µ 2 ).
In the above equations, ρ m is a chosen constant that must satisfy the condition This condition is critical to the stability of the scheme.The scheme is observed to be unstable if this condition is violated.We will employ ρ m = min(ρ 1 , ρ 2 ) for the numerical simulations in Section 3. ν m is a chosen constant that is sufficiently large, and a reasonable condition is µ 0 in equation ( 24c) is a chosen constant that is sufficiently large.In the presence of open boundaries and when µ 1 = µ 2 , the scheme is observed to be unstable if µ 0 min(µ 1 , µ 2 ).We will use µ 0 max(µ 1 , µ 2 ) in the numerical simulations in Section 3. It is observed that increasing ν m tends to improve the stability.Increasing µ 0 also tends to improve the stability in the presence of open boundaries.Note that the constant µ 0 here should not be confused with the field variable µ 0 = µ(φ 0 ) in equation (21), which represents the distribution of the dynamic viscosity at time step zero.
We would like to make several comments on the above scheme: • The computations for the pressure P n+1 and the velocity ũn+1 are de-coupled in this algorithm, and the velocity u n+1 can be evaluated based on equation (26a) once P n+1 is computed.One can recognize that the overall construction of the scheme resembles a rotational pressure-correction type strategy [16].However, the formulation here contains features that distinguish it from the usual pressure-correction formulations.Most notably, the current scheme involves a discrete equation and associated boundary conditions, (25a)-(25c), about an auxiliary variable ξ n+1 .In addition, the pressure P n+1 from the current scheme resides in the H 1 (Ω) space.In contrast, the pressure from the usual pressure-correction formulations resides in the L 2 (Ω) space (see [16]).
• The variable ξ n+1 is an approximation of the quantity ∇ • ũn+1 .The equation (25a) about ξ n+1 exists only in the discrete sense.It is different than the dynamic equation about ∇ • u at the continuum level.
• The scheme leads to linear algebraic systems involving only constant and time-independent coefficient matrices for the pressure, velocity, and the variable ξ n+1 after discretization.This is due to the reformulations of the pressure term and the viscous term, and the introduction of the constants ρ m and ν m in the scheme.The reformulation of the pressure term for coping with the variable density is proposed by [13].The idea for the treatment of the viscous term for dealing with the variable viscosity can be traced to the early works in the 1970s (e.g.[14]); see also later works in e.g.[2,13].Because only constant and time-independent coefficient matrices are involved, which can be pre-computed during pre-processing, the current scheme is computationally very attractive and efficient.
• In the velocity substep we impose a velocity Neumann-type condition, (24c)-(24d), on the open boundary ∂Ω o .The discrete condition (24c) originates from the open boundary condition (16).But it contains constructions involving the constant µ 0 , which are critical to the stability if open boundaries are present.In the absence of the µ 0 constructions, the computation is unstable when the viscosity ratio of the two fluids becomes large and when the fluid interface passes through the open boundaries.The idea of the µ 0 construction for treating the variable viscosity at the open boundary is first proposed by [10].However, there exists a crucial difference in terms of stability between the current scheme and that of [10].The algorithm of [10] is based on a velocity-correction type strategy, and it is observed that a smaller µ 0 constant tends to improve the stability of that scheme in the presence of open boundaries [10].In contrast, the current scheme is based on a pressure-correction type strategy, and we observe that a larger µ 0 constant tends to improve the stability of the scheme when open boundaries are present.
• In the pressure substep we impose a pressure Dirichlet condition, (26d), on the open boundary ∂Ω o .This discrete condition results essentially from taking the inner product between n and the open boundary condition (16).However, note that it contains an extra term µ min ∇• ũn+1 in the construction.
We employ C 0 continuous spectral elements [31,19,34] for spatial discretizations in the current paper.Let us next consider how to implement the algorithm, represented by (24a)-(26d), using C 0 spectral elements.The formulations presented below with no change also applies to C 0 finite elements.
The main issue with regard to the implementation arises from the terms such as ∇ × ∇ × ũ * ,n+1 , ∇ × ∇ × ũn+1 , and ∇ • G n+1 involved in the algorithm.These terms cannot be directly computed in the discrete space of C 0 elements.Note that the term ∇ 2 φ n+1 itself may also cause difficulty to C 0 elements.However, this term will be computed in a proper fashion using C 0 elements later when discussing how to solve the phase field equation.So here we assume that ∇ 2 φ n+1 is already available in a suitable form.
We will derive weak forms of the algorithm for different flow variables.In the process the terms causing difficulty to C 0 elements will be treated in an appropriate way.
Let ω = ∇ × ũ denote the vorticity.Equation (24a) can be re-written as where G n+1 is given by (29).Let and ϕ ∈ H 1 u0 (Ω) denote the test function.Taking the L 2 inner product between ϕ and equation (32), and integrating by part, we get the weak form about ũn+1 , When deriving the above weak form we have used the equations (24c) and (24d), and the following identity (K denoting a scalar field function) Let ϑ ∈ H 1 p0 (Ω) denote the test function, where H 1 p0 (Ω) is defined in (19).Taking the L 2 inner product between ϑ and equation (25a), and integrating by part, we have where we have used the fact that ϑ ∈ H 1 p0 (Ω), equation (25b), the divergence theorem, and the following identity (K denoting a scalar field function) We note the identity where we have used the fact ϑ ∈ H 1 p0 (Ω), and have repeatedly used the divergence theorem.Then, equation (36) can be transformed into the final weak form about ξ n+1 , Let q ∈ H 1 p0 (Ω) denote the test function.Taking the L 2 inner product between ∇q and equation (26a) and integrating by part, we obtain the weak form about P n+1 , where we have used the divergence theorem, and the equations (26b) and (26c).One can observe that the weak forms ( 34), ( 39) and (40) involve no derivatives of order two or higher, and all the terms can be computed directly with C 0 elements.These weak forms can be discretized in space using C 0 spectral elements in the standard way [19].
Given (ũ n , u n , P n , φ n+1 , ∇ 2 φ n+1 ), our final algorithm for solving the momentum equations therefore consists of the following procedure.We refer to this procedure as AdvanceMomentum hereafter.It produces (ũ n+1 , u n+1 , P n+1 ) as follows: AdvanceMomentum procedure: • Solve equation ( 34), together with the velocity Dirichlet condition (24b) on ∂Ω i ∪ ∂Ω w , for ũn+1 ; • Solve equation ( 39), together with the Dirichlet condition (25c) on ∂Ω o , for ξ n+1 ; • Solve equation ( 40), together with the pressure Dirichlet condition (26d) on ∂Ω o , for P n+1 ; • Evaluate u n+1 based on equation (26a) in the following form: In the above algorithm, when imposing the Dirichlet condition (25c) about ξ n+1 on ∂Ω o and when imposing the pressure Dirichlet condition (26d) on ∂Ω o , it should be noted that with C 0 elements one needs to first project the Dirichlet data computed from these equations into the H 1 (∂Ω o ), and then impose the projected data as the Dirichlet condition.This is because the expressions for the boundary conditions of (25c) and (26d) involve derivatives, which may not be continuous across element boundaries on ∂Ω o for C 0 elements.
One can observe that the AdvanceMomentum algorithm has the following characteristics: (i) The computations for the velocity, the pressure, and the field variable ξ n+1 are all de-coupled; (ii) The computations for the different components of the velocity ũn+1 are de-coupled in (34); (iii) All resultant linear algebraic systems from the algorithm involve only constant and time-independent coefficient matrices, which can be pre-computed.
As discussed in [13], the density ρ n+1 and the dynamic viscosity µ n+1 computed according to equation (3) based on φ n+1 may encounter numerical difficulties when the density ratio between the two fluids becomes very large or conversely very small.This is because the numerically-computed φ may not exactly lie within the range [−1, 1] and may be slightly out of bound at certain spatial points in the domain, because of the interaction between mass conservation and the minimization of the free energy inherent in the Cahn-Hilliard dynamics [13].At large density ratios, the slightly out-of-range values of φ may cause the density or the dynamic viscosity computed from (3) to become negative at certain points, thus causing numerical difficulties.Following [13], when the density ratio becomes large or conversely small (typically beyond 10 2 or below 10 −2 ), we will use the following modified function for computing the mixture density and dynamic viscosity, Table 1: Normalization constants for the flow variables and parameters.

Overall Method for Two-Phase Flow Simulations
Let us now consider the numerical solution of the phase field equation (1c), together with the boundary conditions (11a) and (11b) for ∂Ω i , (12a) and (12b) for ∂Ω w , and (9a) and (9b) for ∂Ω o .In a previous work [13], we have developed an algorithm for the phase field equation (1c).This algorithm computes the phase field function φ n+1 and ∇ 2 φ n+1 (both in H 1 (Ω) space) by solving two Helmholtz type equations in a successive but un-coupled fashion.We will employ this algorithm for the phase field equation in the current work.For the sake of completeness, we provide a summary of this algorithm for solving the phase field equation together with the boundary conditions in the Appendix of this paper, and it is referred to as the AdvancePhase procedure (see the Appendix).
Our overall method for simulating incompressible two-phase flows is a combination of the algorithm presented in Section 2.2 for the momentum equations and the algorithm in the Appendix for the phase field equation.Specifically, given (ũ n , u n , P n , φ n ), the overall discrete formulation of the method consists of equations (54a)-(54h) (in the Appendix), (24a)-(24d), (25a)-(25c), and (26a)-(26d).With C 0 spectralelement spatial discretizations, we go through the developments discussed in Section 2.2 and in the Appendix to obtain the weak forms for the field variables.The final solution procedure is composed of the following steps: • Compute φ n+1 and ∇ 2 φ n+1 based on the AdvancePhase procedure discussed in the Appendix.
It can be observed that this method has the following characteristics: (1) The computations for all the flow variables and auxiliary variables are completely de-coupled; (2) All the resultant linear algebraic systems after discretization involve only constant and time-independent coefficient matrices, which can be pre-computed; (3) Within each time step, the method involves only the solution of individual Helmholtz-type (including Poisson) equations; (4) The method is suitable for large density ratios and large viscosity ratios, which will be demonstrated using numerical simulations in Section 3.

Representative Numerical Tests
In this section we demonstrate the accuracy of the method presented in Section 2 and its capability for coping with two-phase outflows and open boundaries using several two-phase flow problems.These test problems are in two dimensions.They involve two-phase open boundaries, and large contrasts in densities and dynamic viscosities of the two fluids.Simulation results will be compared with the experimental data and with the exact physical solutions from theory to demonstrate that the method developed herein produces physically accurate results.
We first briefly mention the normalization of the governing equations and physical parameters, which has been discussed at length in previous works [8,10].Let L denote the characteristic length scale and U 0 denote the characteristic velocity scale.In Table 1 we have listed the normalization constants for different physical variables and parameters.For instance, the non-dimensional mixing energy density coefficient is given by λ ρ1U 2 0 L 2 based on this table.When the flow variables and parameters are normalized as given by the table, the forms of the governing equations and the boundary conditions will remain unchanged upon normalization.In the following discussions all the flow variables and physical parameters are given in non-dimensional forms unless otherwise noted, with the understanding that they have all been properly normalized.

Convergence Rates
The goal of this section is to study the convergence behavior of the method from Section 2, and to demonstrate its spatial and temporal convergence rates using a contrived analytic solution to the two-phase governing equations.
The setup of the problem is as follows.Figure 1(a) shows the rectangular domain ABCD for this problem, 0 x 2 and −1 y 1.We consider the following analytic expressions for the flow variables where (u, v) are the x and y velocity components, and A, B, a, W , a 1 , b 1 and W 1 are prescribed constants to be specified below.It is evident that the u and v expressions satisfy the equation (1b).The external force f (x, t) in (15) and the source term g(x, t) in (1c) are chosen such that the analytic expressions in (43) satisfy the equations ( 15) and (1c).
For the boundary conditions, on the sides AD, AB and BC we impose the Dirichlet condition (10) for the velocity with the boundary velocity w chosen according to the analytic expressions of (43), and we impose the contact-angle conditions (12a)-(12b) for the phase field function, in which θ s = 90 0 and g c1 and g c2 are chosen such that the φ expression in (43) satisfies the equations (12a) and (12b).On the side CD we impose the open boundary condition (16), in which f b is chosen such that the analytic expressions in (43) satisfy ( 16), and we impose the conditions (9a)-(9b) for the phase field function, in which D 0 = 0 and g a1 and g a2 are chosen such that the φ expression in (43) satisfies the equations (9a) and (9b).For the initial conditions (13) we choose u in and φ in according to the analytic expressions in (43) by setting t = 0.
We partition the domain along the x direction using two quadrilateral spectral elements of the same size as shown in Figure 1(a).The system of governing equations ( 15), (1b) and (1c) is integrated over time with the algorithm presented in Section 2 from t = 0 to t = t f (t f to be specified below).Then we compute and monitor the errors of the simulation results at t = t f against the analytic solution given in (43).The parameters for this problem are listed in Table 2.
In the first group of tests we fix the final integration time at t f = 0.1 and the time step size at ∆t = 0.001 (100 time steps).Then we vary the element order systematically between 2 and 20. Figure 1(b) shows the L 2 errors of the velocity, pressure and the phase field function at t = t f as a function of the element order.The results correspond to the open boundary condition (5b) on the side of the domain CD.It can be observed that the numerical errors decrease exponentially as the element order increases (when below order 10).As the element order increases beyond 12, the error curves level off due to the saturation by the temporal truncation error.
In the second group of tests we fix the final integration time at t f = 0.1 and the element order at a large value 18, and then vary the time step size systematically between ∆t = 1.953125 × 10 −5 and ∆t = 0.01.In Figure 1(c) we plot the L 2 errors of the flow variables as a function of ∆t in logarithmic scales.A slope of 2 has been observed in the error curves when the time step size becomes small.The results of these tests demonstrate that the method developed in Section 2 has a spatial exponential convergence rate and a temporal second-order convergence rate.

Capillary Wave
The goal of this section is to demonstrate the physical accuracy of our method using a two-phase capillary wave problem, whose exact physical solution is known from the literature [26].The problem involves two fluid phases, density contrast, viscosity contrast, gravity and the surface tension effects.We have considered this problem in a previous work [13].It should be noted that the algorithm tested here is different from that of [13].
Here is the setting of the problem.We consider two immiscible incompressible fluids in an infinite domain.The lighter fluid occupies the top half of the domain, and the heavier fluid occupies the bottom half.The gravity is in the vertical direction and points downward.Without loss of generality we assume that the first fluid is lighter than the second one (ρ 1 ρ 2 ).At t = 0, the interface formed between the two fluids is perturbed by a small-amplitude sinusoidal wave from its equilibrium horizontal position, and starts to oscillate.The goal of this problem is to study the behavior of the interface over time.
Prosperetti [26] reported an exact standing-wave (but time-dependent) solution to this problem under the following condition: The two fluids may have different densities and dynamic viscosities, but their kinematic viscosities must match.The relation of the capillary-wave amplitude over time has been provided.We will simulate this problem under the same condition, and compare with the exact physical solution from [26].
Specifically, we consider a computational domain as depicted in Figure 2 (non-dimensionalized), 0 x 1 and −1 y 1.The un-perturbed equilibrium position of the fluid interface coincides with the x-axis.We assume that the initial perturbation profile of the interface is given by where λ w = 1 is the wave length of the perturbation profile, and H 0 = 0.01 is the initial amplitude of the capillary wave.Note that the capillary wave-length λ w is chosen to be the same as the domain dimension in the x direction, and that the initial capillary amplitude H 0 is small compared to the domain dimension in the y direction.We employ the algorithm developed in Section 2 to solve the governing equations ( 15) and (1b)-(1c), where the external body force in ( 15) is set to f = ρg r and g r is the gravitational acceleration.boundary conditions, in the horizontal direction we assume that it is periodic at x = 0 and x = 1.At the bottom of the domain (y = −1), we assume a solid wall in the simulations, and impose the Dirichlet condition (10) with w = 0 for the velocity, and impose the boundary conditions (12a)-(12b) with g c1 = g c2 = 0 and θ s = 90 0 .On the top side (y = 1) we assume that the domain is open, and impose the open boundary condition ( 16) with f b = 0 for the momentum equation, and impose the open boundary conditions (9a)-(9b) with g a1 = g a2 = 0 and D 0 = 0 for the phase field function.We employ the following initial velocity and phase field function in the simulations We discretize the domain using 240 quadrilateral elements, with 10 elements in the x direction and 24 elements in the y direction.The elements are uniform along the x direction, and are non-uniform along the y direction, clustering around the region −0.012 y 0.012.We have used an element order 14 for all the elements.The non-dimensional time step size is fixed at ∆t = 2.5 × 10 −5 in the simulations.
We choose the physical parameters for this problem in accordance with those in [13].A summary of the values for the physical and numerical parameters in this problem is provided in Table 3.Note that while ρ 2 and µ 2 are varied in different cases in the simulations, the relation µ2 ρ2 = µ1 ρ1 is maintained according to the condition of the exact physical solution by [26].
Let us compare the simulation results with the exact physical solution given by [26]. Figure 3 shows the time histories of the capillary amplitude H(t) from the simulation and from the exact solution [26] at several density ratios.Figures 3(a)-(d) respectively correspond to the density ratios ρ2 ρ1 = 2, 50, 200, and 1000.These results are obtained using the open boundary condition (5b) at the upper domain boundary.It can be observed that the fluid interface fluctuates about its equilibrium position with the amplitude attenuated over time.The oscillation frequency decreases with increasing density ratios between the two fluids.One can further observe that the time-history curves from the simulations almost exactly overlap with those from the physical solution given by [26] for all density ratios.The insets of Figure 3(b) and Figure 3(c) are the blow-up views of the curves, which show that the difference between the simulation and the exact physical solution is small.These results indicate that our method presented in Section 2 has produced physically accurate results for the capillary wave problem.

Bouncing Water Drop on Superhydrophobic Surface
The goal of this section is to further evaluate and demonstrate the accuracy of the method developed here by comparing simulation results with the experimental measurement.The test problem considered in this section involves large density ratio, large viscosity ratio, and superhydrophobic walls (i.e.contact angle 150 0 ).A similar problem but under a different condition has been considered in a previous work [8].We consider a rectangular domain (see Figure 4 , where L = 4mm.The domain is periodic in the horizontal direction at x = ± L 2 .The top and bottom of the domain are two superhydrophobic solid walls.If the air-water interface intersects the walls, the contact angle is assumed to be 170 0 .The domain is initially filled with air.A water drop, initially circular with a radius R 0 = L 4 , is suspended in the air.The center of the water drop is initially located at a height H 0 above the bottom wall, that is, (x 0 , y 0 ) = (0, H 0 ), where (x 0 , y 0 ) is the coordinate of the center of mass of the water drop.The gravity is assumed to be in the −y direction.At t = 0, the water drop is released, and falls through the air, impacting and bouncing off the bottom wall.The objective of this problem is to simulate and study the behavior of the water drop.
The physical properties of the air, water and the air-water interface employed in this problem are listed in Table 4.The air and the water are respectively assigned as the first and the second fluids in the simulations.We use L as the characteristic length scale, and choose the characteristic velocity scale U 0 = √ g r0 L, where g r0 = 1m/s 2 .The problem is then non-dimensionalized according to Table 1.

Time
Center of Mass (Normalized y coordinate) To simulate the problem we discretize the domain using 150 equal-sized quadrilateral elements, with 10 and 15 elements in the x and y directions respectively.We use an element order 14 for all elements in the simulations.The algorithm presented in Section 2 is employed for marching in time, with a non-dimensional time step size ∆t = 2.5 × 10 −5 .In the horizontal direction we employ periodic boundary conditions for all flow variables.At the top and the bottom walls, we impose the Dirichlet condition (10) with w = 0 for the velocity, and impose the contact-angle boundary conditions (12a)-(12b) with g c1 = g c2 = 0 and θ s = 10 0 for the phase field function.Note that θ s in (12b) is the angle measured on the side of the first fluid, that is, the air for the current configuration.We employ the following initial velocity and phase field function distributions where X 0 = (x 0 , y 0 ) is the initial coordinate of the center of mass of the water drop.The values for the physical and numerical parameters of this problem are summarized in Table 5.
Let us first look into the dynamics of this air-water two-phase system.Figure 4 shows a temporal sequence of snapshots of the air-water interface.The initial height of the water drop is H 0 = 3.2mm above the bottom wall.The interface is visualized by the contour levels φ = 0 at different time instants.Upon release, the water drop falls through the air (Figures 4(a)-(b)), and impacts the bottom wall (Figure 4(c)).One can observe a notable deformation of the water drop upon impact of the wall.Subsequently, the water drop bounces off the bottom wall (Figure 4(d)) and rises through the air, reaching a maximum height (Figure 4(e)).Then the drop falls through the air again and impacts the bottom wall a second time (Figures 4(f)-(h)).This process repeats several times, and the water drop eventually settles down on the bottom wall.
We have monitored the motion of the center of mass of the water drop for different values of the initial drop height, ranging from H 0 = 1.6mm to 4mm.The drop center of mass is defined by where Ω w (t) is the domain occupied by the water drop at time t and demarcated by the contour level φ = 0.

© + x
Figure 6: Comparison of restitution coefficient as a function of impact velocity between current simulations and the experiment [29].H 0 is the initial height of the water drop.a number of times in all these cases.One can also discern the oscillation in the drop shape in later time, before it completely settles down on the wall.For larger values of the initial drop height, we notice that occasionally the water drop can reach a maximum height after a bounce that is quite close to that before the bounce; see for instance the second and third peaks in the curve for H 0 = 3.2mm of Figure 5.This is likely due to the fact that a larger initial drop height tends to cause a more pronounced deformation of the water drop upon impact and a more pronounced oscillation in the drop shape after the bounce-off into the air (see Figures 4(c)-(g)).The elastic energy associated with the drop deformation can be converted to the kinetic energy associated with the motion of the drop center of mass in subsequent impact and lift-off, thus resulting in a maximum height close to that before the bounce.
We have computed the restitution coefficient based on the time histories of the center of mass of the water drop.We follow [29] and define the restitution coefficient C res by where H and H ′ respectively denote the maximum heights of the water drop above the bottom wall before and after the bounce.We also follow [29] and estimate the impact velocity of the water drop V imp by where g r is the gravitational acceleration.
In Figure 6 we plot the restitution coefficient C res as a function of the impact velocity V imp from the current simulations.For comparison, we have also shown the restitution coefficient data from the experiment of [29].The restitution coefficients corresponding to different initial drop heights H 0 from the simulations have been included in this figure.The drop size in the simulations is a little larger than that in the experiment of [29].The bulk of the restitution coefficients from current simulations agree well with the experimentally determined values.On the other hand, some differences can be observed, especially for the data points corresponding to the first few bounces with larger initial drop-height values.We observe that for such cases the restitution coefficients from the simulation tend to be a little smaller than the bulk of the experimental data.This is likely due to the larger drop deformation upon impact and stronger drop-shape oscillation after the bounce-off, associated with a larger initial drop height and a larger impact velocity.The elastic energy associated with the drop deformation may reduce the maximum height the drop can reach after the bounce, and thus results in a smaller restitution coefficient.The outlying data point, with a large restitution coefficient, from the simulation case with an initial drop height H 0 = 3.2mm corresponds to the second and the third peaks in the time-history curve in Figure 5, which has been discussed in a previous paragraph.
The above comparison indicates that the simulation results obtained using our method overall are in good agreement with the experimental measurement.

Air Jet in Water with Two-Phase Open Boundaries
The goal of this section is to demonstrate the effectiveness of the open boundary conditions and the numerical algorithm from Section 2 for two-phase outflow problems.The test problem considered in this section involves open boundaries where the two fluids may leave or enter the domain, large density contrast, and large viscosity contrast.The fluid interface passes through the open domain boundary in this problem.
We consider the long-time behavior of an air-water two-phase flow, in which a train of air bubbles continually forms at a wall inside the water and then moves out of the domain due to buoyancy.This flow problem has been considered in a previous work [10].It should be noted that the open boundary conditions and the numerical algorithm being tested here are different.
Specifically, we consider the flow domain shown in Figure 7 , where L = 3cm.The bottom of the domain is a solid wall, while the other three sides (top, left and right) are all open, where the fluid can freely leave or enter the domain.The domain is initially filled with water, and the gravity is along the vertical direction pointing downward.The bottom wall has an orifice in its center, with a diameter d = 6mm.A stream of air is continuously injected into the domain through the orifice.The air velocity has a parabolic profile at the orifice, with a centerline value U 0 = 17.3cm/s.The bottom wall has a neutral wettability, that is, if the air-water interface intersects the wall the contact angle would be 90 0 .Our objective is to simulate and study the long-time behavior of this system.The two-phase open boundaries coupled with the large density ratio between air and water make this problem very challenging to simulate.
The physical parameters concerning the air, water and the air-water interface have been provided in Table 4.We treat the air and the water as the first and the second fluids, respectively.L and U 0 are employed respectively as the characteristic length and velocity scales.Normalization of the problem then proceeds according to Table 1.
The flow domain is discretized using 600 quadrilateral spectral elements, with 20 and 30 elements in the x and y directions respectively.An element order 12 has been used for all elements in the simulations.At the bottom wall, we impose the velocity Dirichlet condition (10)   (12a)-(12b) with g c1 = g c2 = 0 and θ s = 90 0 .At the air inlet we impose the velocity Dirichlet condition (10), in which w has zero horizontal component and its vertical component takes a parabolic profile with a centerline value U 0 ; for the phase field function, we impose the boundary conditions (11a)-(11b), in which g b = 0 and where R = d 2 = 3mm is the radius of the orifice, and H(x, a) is the heaviside step function taking unit value if x a and vanishing otherwise.On the top, left and right sides of the domain, we impose the open boundary condition (16) with f b = 0 for the momentum equation; for the phase field function we impose the boundary conditions (9a)-(9b) with g a1 = g a2 = 0.For the initial conditions, we have used an instantaneous snapshot of the velocity field and the phase field function from the simulation of [10].Because long-time simulations have been performed, the initial velocity and phase field distributions have no effect on the long-time behavior of the system.
We apply an external pressure gradient in the y direction (− ∆P L ) to balance the weight of water in the simulations, i.e.
where ρ w is the water density and g r is the magnitude of the gravitational acceleration.Table 6 summarizes the physical and numerical parameter values in the simulations for this problem.The D 0 in the open boundary condition (9b) for the phase field function is determined based on a preliminary simulation with D 0 = 0. Preliminary simulations indicate that the air bubbles have a non-dimensional convection velocity about 2.0 ∼ 3.0 at the upper domain boundary.Because 1 D0 plays the role of a convection velocity, we therefore use an outflow dynamic mobility 1 D0U0 ≈ 2.5 in the simulations.A non-dimensional time step size 1.5 × 10 −6 has been employed for the current problem.
Let us first demonstrate the long-term stability of the computation.We have performed long-time simulations of this problem using different open boundary conditions.Figure 8 shows a window of the time histories of the average vertical velocity magnitudes V avg (t), where v is the y velocity component and V Ω = Ω dΩ is the volume of the domain.Results in Figures 8(a)-(d) are obtained using the open boundary conditions (5a)-(5d), respectively.One can make two observations.First, the average velocity magnitude V avg fluctuates over time about some constant mean level and its time history signal exhibits a quasi-periodic nature.This indicates that the flow is at a statistically stationary The dynamics of this air-water flow is illustrated by Figure 9, in which we show a temporal sequence of snapshots of the air-water interface in a time-window between t = 16.9397 and t = 17.2022.The fluid interface is visualized using the contour level φ(x, t) = 0 in the plots.These results are obtained with the open boundary condition (5b), corresponding to the time history in Figure 8(b).These plots demonstrate the process of free air bubbles generated at the wall rising through water and crossing the upper domain boundary to migrate out of the domain.Figures 9(a)-(e) show the leading air bubble passing through the upper open boundary of the domain.They demonstrate that the boundary condition and the numerical algorithm we developed in Section 2 can effectively allow the fluid interface to pass through the open/outflow boundary in a smooth fashion.Simultaneously, one can observe that the trailing free bubble rises through the water, and that a new air bubble is forming at the bottom wall (Figures 9(b)-(h)).Subsequently, the air bubble at the wall breaks free, and the above process will repeat itself.We further illustrate the flow dynamics using instantaneous velocity distributions.Figure 10 is a temporal sequence of snapshots of the velocity fields at identical time instants as those of the interfacial plots of Figure 9.One can observe that a significant flow field is induced in the regions occupied by the air bubbles, and that a particularly strong velocity field exists inside the free air bubble as it initially breaks free from the wall; see the region of the trailing free bubble in Figures 10(a)-(b).On the other hand, the velocity field in the water region is in general quite weak.As the air bubble rises through the water, a pair of vortices forms in the water region trailing the air bubble; see the region behind the second air bubble in Figures 10(e)-(h).These vortices can induce a backflow on portions of the outflow/open boundary after the air bubble passes through (Figure 10(h)).
The above results illustrate one state of the flow.We observe that this air-water flow can exhibit another state, in which the flow characteristics are somewhat different than those seen above.In Figure 11 we show another window in the time history of the average magnitude of the vertical velocity, obtained also with the open boundary condition (5b).The flow evidently is at a statistically stationary state.Contrasting this figure with Figure 8(b), which is computed using the same boundary conditions, we can observe that the velocity-history curves have qualitatively different characteristics in these figures.
This different flow state is further illustrated by the temporal sequence of snapshots of the air-water interface shown in Figure 12, which covers a time window between t ≈ 23 and t ≈ 23.5 in the history plot of Figure 11.These results correspond to the open boundary condition (5b).The plots clearly show the breakaway the air bubble from the wall (Figures 12(a)-(c)) and the bubble motion across the domain and the upper open boundary (Figures 12(d)-(k)).The crucial difference, when compared with Figure 9, lies in the following.When multiple free bubbles are present in the domain, the interaction between the leadingbubble wake and the trailing bubble appears to have caused the trailing bubble to accelerate and nearly catch up with the leading one; see Figures 12(e)-(j).This has also induced significant deformations in the trailing bubble (Figures 12(j)-(l)), and caused it to subsequently break up (Figures 12(m)-(o)).As the free  bubbles (and their daughter bubbles) quickly move out of the domain, one can observe that another bubble is forming, but still attached to the wall (Figure 12(o)).Consequently, the flow domain will be depleted of free bubbles for a period of time beyond the time instant corresponding to Figure 12(o), until the air bubble attached to the wall breaks free.This scenario is more similar to the one discussed in [10], but is quite different from that shown by Figure 9. From Figures 12(i 16) is critical to the stability for this problem.We observe that the computation using an open boundary condition without this term is unstable for this problem, that is, due to the backflows induced by the vortices at the outflow boundary.It is observed that increasing ν m in the algorithm tends to improve the stability, and that a larger µ 0 in (24c) for the numerical treatment of the open boundary condition also improves the stability for the current pressure-correction based scheme.This observation concerning µ 0 seems different from the trend observed in [10], which is for a velocity-correction based algorithm.

Concluding Remarks
We have presented several new open boundary conditions for two-phase outflows, and a rotational pressurecorrection based algorithm for solving the two-phase momentum equations in conjunction with the proposed open boundary conditions.These techniques are then combined with a solver for the phase-field equation to form an efficient and effective method for simulating incompressible two-phase flows involving open/outflow boundaries.
The two-phase open boundary conditions presented here are inspired by the two-phase energy balance discussed in the previous work [10].The current work has provided a generalization and several new forms for the open boundary condition beyond those developed in [10] for the momentum equations.
The algorithm presented herein for the two-phase momentum equations is based on a rotational pressure correction-type strategy for de-coupling the velocity/pressure computations.More importantly, the current algorithm results in velocity and the pressure linear algebraic systems with constant and time-independent coefficient matrices after discretization, despite the variable nature of the mixture density and mixture viscosity.Therefore, these coefficient matrices can be pre-computed during pre-processing.In a previous work [13] we have developed a velocity correction-based algorithm for the variable-density Navier-Stokes equations that possesses similar properties (leading to constant coefficient matrices for pressure/velocity linear systems); see also subsequent applications and further developments based on that algorithm in [8,10,9].The algorithm developed herein in a sense can be considered as the pressure-correction counterpart to the scheme of [13].The implementation of the algorithm presented herein is suitable for C 0 spectral elements, and with no change it also applies to conventional finite elements.It should be noted that the rotational pressure correction formulation embodied in the current algorithm has a difference to the usual pressure correction formulations (see e.g.[16]), in that apart from the velocity/pressure we have introduced a discrete equation and the corresponding boundary conditions for another field variable ξ n+1 in the algorithmic formulation.
The numerical treatments for the open boundary conditions proposed herein involve imposing a discrete Neumann type condition on the outflow boundary at the velocity substep, and two discrete Dirichlet type conditions on the outflow boundary at the substeps for ξ n+1 and pressure respectively.The discrete velocity-Neumann and the pressure-Dirichlet conditions on the outflow boundary stem largely from the continuous open boundary condition.But they contain modifications and additional terms that are essential to the stability of the algorithm.
To demonstrate the physical accuracy of the method developed herein, we have considered the capillary wave problem and compared quantitatively the numerical solution with the two-phase exact physical solution by [26] for a range of density ratios (up to 1000).The comparisons show that our method has produced physically accurate results.We have also considered the bounce of a water droplet on a superhydrophobic surface, and compared the restitution coefficients from the simulations and the experimental measurement of [29].The simulation results are in good agreement with the experimental data. We

Figure 1 :
Figure 1: Spatial/temporal convergence rates: (a) Mesh and boundary conditions; (b) Numerical errors versus element order showing spatial exponential convergence (with fixed ∆t = 0.001); (c) Numerical errors versus ∆t showing temporal second-order convergence rate (element order fixed at 18).On the face CD the open boundary condition (5b) is used.

Figure 2 :
Figure 2: Configuration for the capillary wave problem.

Figure 5 :
Figure 5: Time histories of the water-drop center of mass (y coordinate) corresponding to several different initial drop heights.
we show the time histories of the y component (normalized) of the drop center of mass for several values of the initial drop height H 0 .It can be discerned that the water drop bounces off the bottom wall

Figure 7 :
Figure 7: Configuration of the air jet in water problem.
)-(k) we can again observe that our method allows the air bubble and the air-water interface to cross the open domain boundary in a smooth fashion.The air jet in water problem is a stringent test to the open boundary conditions.The presence of twophase open boundary, combined with the large density ratio between air and water, makes this problem extremely challenging to simulate.The results of this section show that the two-phase open boundary conditions and the numerical algorithm developed in the current work are effective for two-phase outflows with large density and viscosity contrasts at the outflow boundaries.The E(ρ, n, u) term in the open boundary condition ( have further simulated the air jet in water problem to test the effectiveness of the open boundary conditions and the numerical algorithm for two-phase problems involving outflow or open boundaries.This problem involves large density ratio, large viscosity ratio, and backflows or vortices at the two-phase open boundary.The results demonstrate the long-time stability of the method presented herein.It is also shown that our boundary conditions allow the fluid interface to pass through the open domain boundary in a smooth and seamless fashion.

Table 2 :
Parameter values for convergence tests.

Table 3 :
Parameter values for the capillary wave problem.

Table 4 :
Physical properties of air and water.

Table 5 :
Physical and numerical parameter values for the bouncing water drop problem.
with w = 0 and the boundary conditions

Table 6 :
Physical and numerical parameter values for the air jet in water problem.