FULLY DECOUPLED ENERGY-STABLE NUMERICAL SCHEMES FOR TWO-PHASE COUPLED POROUS MEDIA AND FREE FLOW WITH DIFFERENT DENSITIES AND VISCOSITIES

. In this article, we consider a phase field model with different densities and viscosities for the coupled two-phase porous media flow and two-phase free flow, as well as the corresponding numerical simulation. This model consists of three parts: a Cahn–Hilliard–Darcy system with different densities/viscosities describing the porous media flow in matrix, a Cahn–Hilliard–Navier–Stokes system with different densities/viscosities describing the free fluid in conduit, and seven interface conditions coupling the flows in the matrix and the conduit. Based on the separate Cahn–Hilliard equations in the porous media region and the free flow region, a weak formulation is proposed to incorporate the two-phase systems of the two regions and the seven interface conditions between them, and the corresponding energy law is proved for the model. A fully decoupled numerical scheme, including the novel decoupling of the Cahn–Hilliard equations through the four phase interface conditions, is developed to solve this coupled nonlinear phase field model. An energy-law preservation is analyzed for the temporal semi-discretization scheme. Furthermore, a fully discretized Galerkin finite element method is proposed. Six numerical examples are provided to demonstrate the accuracy, discrete energy law, and applicability of the proposed fully decoupled scheme.


Introduction
Multi-phase flow appears in a wide range of applications, from fluid dynamics and fluid-surfactant dynamics to solid physics and biomedical sciences.Multi-phase flow in karstic geometry, as a typical multi-phase flow in coupled free flow region and porous medium region, plays a crucial role in the research of contaminant transport in groundwater system [59,68,83,86].Related work can be extended to study other important problems with this type of coupled flows, such as oil recovery in petroleum engineering [42,88,96], cardiovascular blood flow [37], and proton exchange membrane fuel cell [84].Therefore, it is worthwhile to develop efficient numerical schemes for studying these coupled multi-phase flow problems.
In the past two decades, scientists have published research results about single-phase models of the coupled free flow and porous media flow [7,17,21,28,53,60,77], which are composed of Stokes/Navier-Stokes equation, Darcy's law, and three interface conditions coupling these two parts.There have been abundant numerical attempts for these models, for instance, domain decomposition methods [10,18,19,25,31,33,48], discontinuous Galerkin method [41,54,65,76], two-grid/multi-grid methods [3,14,70], and many others.However, there have been much less works about the important multi-phase models of the coupled free flow and porous media flow which are much more challenging than the single phase models in this area.
Phase-field has become a popular modeling technique in simulating multi-phase flow [11,12,64,66,67,73,89,95] as a viable tool to describe the interfacial motion between two immiscible fluids.By this technique, the interface is described by a smooth layer with finite width so that the singularities across the interface is removed [34,56,80].One of its attractive features is that it automatically captures the interface transition in numerical simulation based on the energy law [57,62].
In recent years, the phase field modeling for two-phase flow across the interface between free flow region and porous medium region has started to attract researchers' attention [5,26,27,30,39,49,50,82].Among the models proposed in these works, we will consider the modeling ideas in [27,39,49,50].Han et al. [49] deduced a family of phase field models for the two-phase flows in karstic geometry based on a fundamental diffuse interface model with constant density/viscosity, which preserves some physical property, e.g., energy law.Then, for these models, Chen et al. [27] and Han et al. [50] systematically studied the existence and uniqueness of the weak solution, and proposed and analyzed some unconditionally stable numerical algorithms for the involved Cahn-Hilliard-Stokes-Darcy system with constant density/viscosity.Gao et al. [39] developed a partially decoupled numerical method for the Cahn-Hilliard-Navier-Stokes-Darcy (CHNSD) system with constant density/viscosity.Gao et al. [40] developed another type of partially decoupled numerical method in a simplified framework for CHNSD model with different densities/viscosities based on the idea of a single Cahn-Hilliard equation on the whole domain (including both the porous media region and the free flow region).Hence it does not consider the four phase field interface conditions but only the three flow field interface conditions in the design of the partially decoupling scheme.Due to the large scale computation of the realistic applications of the CHNSD system with different densities/viscosities, a fully decoupled scheme is in great needs to further improve the computational efficiency.Therefore, in this article, we will develop a fully decoupled scheme, including the decoupling of the Cahn-Hilliard equations in the porous media region and the free flow region through the four phase field interface conditions, the decoupling of the flow fields through the three flow interface conditions, the decoupling between the phase field and flow field by operator splitting, the decoupling between the velocity and pressure by rotational pressure projection, and the invariant energy quadratization (IEQ) method to handle the nonlinear potential term.
The CHNSD model with different density/viscosity considered in this article consists of a coupled Cahn-Hilliard-Navier-Stokes system with different density/viscosity describing the dynamics of free flow fluid in conduit, a coupled Cahn-Hilliard-Darcy (or Cahn-Hilliard-Hele-Shaw) system with different density/viscosity describing the evolution of porous media flow in matrix, and seven interface conditions describing the transition of two-phase flow across the interface between the conduit region and the matrix region.Three interface conditions are for the flow field and the other four interface conditions are for the phase field.For this sophisticated CHNSD model with different densities and viscosities, here are a few major numerical challenges: (i) the coupling of free flow region and porous media region demands high computing costs; (ii) different densities and viscosity cause the strong coupling and large derivative across the interface between two fluids; (iii) the relative slow motion in matrix requires long time numerical simulation; (iv) the coupling between the phase function and the flow (in both regions) mandates complicated computation; (v) the inherent stiffness in the phase field model associated with the thin transition layer calls for a carefully treatment to capture the evolution of interface between binary immersible fluids.Even though methods dealing with one of the challenges mentioned above may exist, a straightforward combination of the existing methods cannot address these challenges posed by such a complicated coupled system.It is our intention to develop efficient and stable numerical schemes for simulating this CHNSD model with different density/viscosity.
Techniques can be found in the literature for efficiently dealing with the stiffness in diffuse interface models, for example, the convex splitting method [6,43,44,94], the stabilized semi-implicit approach [15,35,90], the invariant energy quadratization (IEQ) method [24,58,72,91,93,97,98], the scalar auxiliary variable (SAV) approach [74,81,92,99], the exponential time differencing method [63,87], etc.The idea of the IEQ is inspired by a Lagrange multiplier method [47].An essential effect of the IEQ method is to semi-implicitly treat the nonlinear term for phase field type models by making free energy quadratic, which yields a linear system at each time step.The IEQ method requires the nonlinear potential to be bounded from below; nevertheless, most nonlinear potentials usually satisfy this constraints.
In this article, we develop energy-law preserving time-stepping schemes for solving the CHNSD system with different density/viscosity, by using the following ideas.To decouple the computation for fluid equations in free flow region and in porous media region, we adopt a novel weak formulation that weakly enforces interface conditions imposed at the boundary between these two regions.The ideas of the IEQ approach and operator slitting technique will be employed to derive a partially decoupled method.The pressure correction method and stabilization strategy are exploited to design a completely decoupled energy-stable numerical scheme.Following these ideas, a fully decoupled numerical method is developed by cooperating the grad-div technique.The proposed methods solve the two Cahn-Hilliard systems in their own subdomains, respectively, instead of directly solving one Cahn-Hilliard equation on the whole domain [27,40].The unconditional stability is proved for the semi-discretization schemes of the two decoupled schemes.It's worth to be pointed out that the decoupling technique for fluid and phase field variables is not trivial because of the different densities considered in the momentum equation for binary fluids.To the best of our knowledge, this is the first work to develop the fully decoupled numerical method for CHNSD system with different densities and viscosities.
The rest of the article is as follows.In Section 2, we present a Cahn-Hilliard-Navier-Stokes-Darcy (CHNSD) model for two-phase flows associated with different densities and viscosities in free fluid region together with porous media.We also prove the energy law of this model.In Section 3, we develop two unconditionally stable decoupled numerical methods by weakly enforcing the domain interface boundary conditions, and discuss their energy dissipation law.In Section 4, we present a fully discrete Galerkin finite element method for this CHNSD model.In Section 5, we present several tests to demonstrate the accuracy and energy law of numerical simulations generated by the proposed schemes, and distinctive features of two-phase flows in the free flow region and the porous media region.

Cahn-Hilliard-Navier-Stokes-Darcy model
In this section, we briefly present the Cahn-Hilliard-Navier-Stokes-Darcy (CHNSD) model with different densities and viscosities based on the existing Cahn-Hilliard-Stokes-Darcy model [49,50] and CHNSD model [39] for matched density in karstic geometry.After deriving a weak formulation in cooperate with the interface conditions for this coupled model, we show that the model obeys a dissipative energy law at the PDE level.

Model system
Consider the temporal domain [0,  ] ( > 0).Let Ω = Ω  ⋃︀ Ω  ⊂ R  , ( = 2, 3) be a bounded domain, where Ω  is the free flow region and Ω  is the porous media region.These two subdomains have Lipschtiz continuous boundaries Ω  and Ω  , respectively, which form the interface Γ := Ω  ∩ Ω  .Denote the unit normal vector from porous media Ω  to free fluid region Ω  by   and its opposite vectors by   .A geometry in two dimensions is depicted in Figure 1.
We consider the double well potential for the Helmholtz free energy  () = 1 4 ( 2 −1) 2 , and let  () =  ′ (), where  is the thickness of transition layer between two immiscible fluids.Let   denote the phase function with respect to Ω  ( = , ) taking distinct values ±1 to indicate two different fluid phases and varying smoothly across the diffuse interfacial region.The mobility variable is represented by   ( = , ) as a function of order parameter .Variable   ( = , ) is the chemical potential related to phase variable   .The parameter  is the elastic relaxation time of mixing interface.
In this paper, we consider the two-phase fluids with different densities  and viscosities  between fluids as follows We denote the external gravitational forces by  arising from density difference corresponding to the gravity vector  = , where the constant  represents the gravity acceleration and  describes the unit vector directing upward.
We assume that the two-phase flow in the porous media region Ω  is described by the Cahn-Hilliard-Darcy (CHD) system as follows: ) where  ∈ [0,  ],   denotes the advective velocity,   denotes the hydraulic head.
We assume that the two-phase flow with different densities and viscosities in the free fluid region Ω  is described by the Cahn-Hilliard-Navier-Stokes (CHNS) system as follows: ) where  ∈ [0,  ],   and   represent the fluid velocity and the pressure in conduit, respectively.T(  ,   ) = 2D(  )−  I prescribes the stress tensor with the deformation tensor D(  ) = (∇  +∇    )/2 and the identity matrix I.
In order to accomplish the exchange of fluids across interface Γ between the conduit and the matrix, we consider the following interface conditions: ) ) ) )  where   ( = 1, • • • ,  − 1) are the unit tangential vectors along the interface and ∏︀ is the permeability of the porous media, ∏︀ = K(  ) [14,22,49].Equations (2.11)-(2.14)describe the continuity of the phase variable, chemical potential function as well as their derivative with respect to normal direction [49,50].Equation (2.15) stands for the continuity of velocity in normal direction and guarantees the mass conservation.Equation (2.16) is known as Lions interface condition describing the balance of forces [28,49,51] arising from the momentum equation written in divergence form [20,32]. Equation (2.17) is the well-known Beavers-Joseph-Saffman-Jones (BJSJ) interface condition [7,16].
We assume that the coupled CHNSD system satisfies the following boundary conditions ) where Γ  = Ω  ∖Γ and Γ  = Ω  ∖Γ, as well as initial conditions   (0, , ) =  0  (, ),  = , ,   (0, , ) =  0  (, ). (2.20) Without loss of generality, K is assumed to be a bounded, symmetric and uniformly positive definite matrix throughout the paper.For the sake of simplicity, we ignore the body force term  as presented in (2.2) and (2.7) in the rest part of this paper except otherwise specified.Following the idea in [45,79], denoting  = √ , we deduce by using the mass conservation equation Remark 2.1.We follow the idea presented in [38,78] for (2.21) and (2.22).There is another alternative thermodynamically consistent model for impressible two-phase flow with different densities [1,80], which obeys the conservation property where  = 2−1 2   ∇  .In this case, one may modify the interface condition (2.16) to obtain energy dissipation law of the two-phase flow since interface conditions are critical for the mathematical modeling of multi-domain two-phase model.In this paper, we focus on the efficient decoupled numerical schemes of CHNSD model associated with (2.21) for different densities.
In order to explicitly handle with the nonlinear term in double well potential, we introduce new variables and then  (  ) = 1      .Therefore, the Cahn-Hilliard equation in matrix and conduit can be transformed into the following equivalent system:
Then the weak formulation (2.30)-(2.37)and the CHNSD model (2.2)-(2.17)are equivalent.In particular, the solution of the weak formulation (2.30)-(2.37)satisfies the strong interface conditions (2.11), (2.12), and the weak interface conditions as follows, for  = , , ) ) Proof.On one hand, one can easily derive the weak formulation (2.30)-(2.37)from the CHNSD system (2.2)-(2.17),by utilizing the Green's formula and the interface/boundary conditions.On the other hand, one can obtain (2.2)-(2.10)after applying integration by parts to the weak formulation while choosing the test functions In the following, we focus on the equivalence on the interface.First, applying integration by parts and using (2.2) and boundary condition for   in (2.19), one can easily conduct (2.39) from (2.30).

A dissipative energy law
For the solution to the coupled CHNSD system, we consider its energy in the following form: which will not grow as suggested by the theorem below.
Theorem 2.3.The smooth solution (  ,   ,   ,   ) of the coupled CHNSD system with different densities (2.2)-(2.20)satisfies the following energy dissipation law where  is defined by Proof.We firstly consider the conduit part.By letting 37), applying the interface conditions (2.11) and (2.12), and summing up the resultants together, we have Using integration by parts we can see Then, for the matrix part, choosing  =   in (2.30) and  =   in (2.31), adding the resultants together and applying (2.38), we obtain We add (2.53) and (2.54) to obtain

Energy-stable time-stepping methods
In this section, we present and analyze a partially decoupled time-stepping scheme and a fully decoupled scheme for solving the CHNSD system (2.30)-(2.37).First, for a preparation and by following the ideas in [69], we have the following lemma for estimating the interface terms.Lemma 3.1.There exists a constant  such that, for where then Also, in discussions from now on, we let the time step size be Δ =  +1 −   =   ,   be the th discrete time level, and  =  Δ .

A partially decoupled scheme based on the pressure correction method
Following the works in [39], we decouple CHNSD system into two subsystem by the modification of the interface conditions.Then we use a stabilization technique to decouple the CHNS system and the CHD system into two subsystems, respectively.Consequently, we obtain four sub-systems.Finally, we employ the pressure projection correction method [38,78] in dealing with Navier-Stokes equations to derive a linearized system with energy-preserving property.These essential components help us obtain a decoupled, linearized, and energy stable scheme for the temporal discretization as follows.
Remark 3.2.In the scheme above, we imposed the following interface condition for the intermediate velocity where ℰ   is defined as +1  is given by Proof.We first consider the conduit part.By (2.51), we have we obtain 1 2 where the Green's formula is used for term In order to deal with the ‖   ⋆ ‖ 2 term, we rewrite (3.10) as and take the inner product of (3.23) with Δ  ⋆ .Then, by the identity (3.21), we obtain 1 2 Letting  = Δ +1  in (3.8), and applying (2.11) and (3.10), we get    15), and then taking the summation, applying (2.12), and (3.21), we obtain We take the inner product of (3.6) with  ℎ = ΔK −1  +1  , let  = Δ +1  in (3.6), and add the resultants together to obtain we obtain By the triangle inequality, we have Adding (3.33)-(3.35),we derive If we impose  ≥ 2 C which only depends on the geometry of Ω  and Ω  , one can obtain (3.17).Thus, we complete the proof of Theorem 3.5.

A fully decoupled scheme based on grad-div stabilized method
The time-stepping scheme (3.4)-(3.15)presented in the previous subsection decouples the Cahn-Hilliard-Navier-Stokes system in the free fluid region and the Cahn-Hilliard-Darcy system in the porous medium region.But the velocity and pressure in conduit are still coupled, while artificial boundary conditions is imposed to directly decouple the pressure and velocity resulting in the loss of accuracy [46].In this subsection, we revise the partially decoupled scheme for a better efficiency.Specifically, we modify the standard pressure correction scheme in rotational version and introduce a grad-div stabilization [8,9,29,36,61] to completely decouple the velocity and pressure for Navier-Stokes equation while retaining the unconditionally stability stated as a discrete energy law.
In the first three steps of this new scheme, we follow the first three steps of the algorithm above to solve (3.4), (3.5), (3.7), (3.8) Step 5. Find  +1 ∈   , such that where  = 1 2 min{ 1 ,  2 }.Step 6. Find  +1 ∈   , such that where ν = min{ where ℰ    is defined as and Proof.We focus on discussion of the decoupling method for Navier-Stokes equation in conduit part.By taking the test function  = Δ +1  in (3.37), combining with (3.20), and the identity (3.21), we obtain 1 2 Taking  = νΔ  in (3.39) and taking the  2 inner product of (3.40) by  = Δ∇ •  +1  at time   , then adding these two equations, using the identity (3.21), we have

Fully discrete numerical schemes
In this section, we present a fully discrete scheme by further discretizing the underlying spaces in the fully decouple time stepping scheme described by (3.4)-(3.10)and (3.37)- (3.42).Specifically, we let ℑ ℎ be a quasiuniform triangulation of domain Ω, on which, we introduce finite element spaces  ℎ ,  ℎ and  ℎ with  = , .
Here we further assume the finite element spaces  ℎ and  ℎ satisfy an inf-sup condition for the divergence operator in the following form: There exists a constant  > 0 independent of ℎ such that the LBB condition inf Y. GAO ET AL.
Step 8. Find  +1 ℎ ∈  ℎ and  +1 ℎ ∈  ℎ such that Remark 4.1.Comparing the above totally decoupled numerical algorithm with the numerical scheme in [40], the major differences include that (1) our method solves two separate Cahn-Hilliard equations in free fluid flow and porous medium regions, respectively, while, in [40], the two Cahn-Hilliard equations are viewed as a single Cahn-Hilliard on the whole domain.Thus, the treatment of interface conditions associated with phase function are also different from those in [40]; (2) the strategy in dealing with the stiffness of nonlinear term  () is different from those in [40], namely, the IEQ method is utilized to linearize the nonlinear system in this paper while a stabilized semi-implicit strategy is exploited in [40]; (3) the rotational pressure projection and grad-div stabilization are utilized to break the coupling of velocity and pressure for Navier-Stokes equation, while the artificial compression method is employed to decouple the coupling in [40].Therefore, the proposed scheme in this paper is fully decoupled and more efficient, while the one in [40] is partially decoupled.

Numerical example
In this section, several numerical examples are presented to illustrate the properties of proposed CHNSD model and developed numerical methods, such as the order of convergence, discrete energy decay law, and the evolution of a droplet under the influence of discontinuous permeability field, curve interface and buoyancy forces.
In all examples, we employ the  2 - 1 Taylor-Hood elements for the Navier-Stokes equation and quadratic elements for the Darcy equation.We also utilize the quadratic elements for Cahn-Hilliard equation in both subdomains, i.e., free flow region and porous medium, respectively.where A uniform triangular mesh and a uniform time partition are used in this simulation.Numerical errors are listed in Table 1 for the approximations to velocity and phase variable and pressure in fluid flow generated by decoupled linearized numerical scheme at final time  = 5 which is consistent with the expected optimal convergence of this scheme.
Example 2 (Convergence for varying permeability).In order to validate the effects of proposed numerical method for CHNSD model with varying permeability parameters, we take T(  ,   ) = ∇  −   I and K = I.The exact solution is given by [71] )︂ cos(), Table 1.Numerical errors at  = 5 and convergence rates of fully decoupled numerical method with Δ = 0.01ℎ.    2 to 4, one can see that the error norms with respect to variables   ,   ,   ,   , and   have the satisfied order of convergence when the hydraulic conductivity K varies.
Example 3 (Shape relaxation and energy dissipation).In this example, we simulate the evolution of a initial cross shape phase function in domain 1], in order to Table 4. Numerical relative errors at  = 5 and convergence rates of fully decoupled numerical method with  = 0.001 and Δ = 0.005ℎ.In the rest of the simulations, we shall use a spatial mesh adaptive strategy from [52] to improve the accuracy for capturing the interface layer without increasing computational costs excessively.Specifically, based on the observations in [56], adaptive mesh strategy uses at least four grid elements across the transition interface between two-phase fluids.In our computations for the rest of the examples, the uniform fundamental mesh has a mesh size of ℎ = 1 32 on which this mesh adaptive strategy related to the value of |∇| is applied around the interface.
The initial shape and adaptive mesh of phase function are shown in a cross shape is depicted Figure 2. Figure 3 shows the dynamical morphotype of the simulated phase function from the initial cross shape relaxing to a circular shape under the effect of surface tension associated with the density ratio  1 :  2 = 1 : 5 and uniform time partition Δ = 0.005, which is in a good agreement with the benchmark study in [2,55].The discrete mass and energy of numerical results are plotted in Figure 4.It is clear that the discrete mass is conserved and the discrete energy-decay property is observed.We would like to point out that an asymmetrical behavior of the morphotype on the two sides of the sharp interface between the porous media and the free flow region.This phenomenon shows that the relaxation in the free flow region is faster than that of the porous media region, due to the faster flow speed of the free flow.
The adaptive mesh and phase variable on enlarged adaptive grid are illustrated in Figure 5.As shown in this figure, the evolution of interface can be accurately captured under the adaptive strategy utilized in this paper.And the adaptive meshes are well concentrated near the interface region between the binary phase fluids at each time step.
As expected from the literature for this type of benchmark study [99], Figure 6a displays the discrete energy and evolutions of phase variables under Δ = 0.005.This figure clearly shows that the four individual square droplets firstly evolve to four circle droplets, then start the coalescence, and finally merge into a single circle, which remains unchanged over time.The discrete energy of numerical algorithms is plotted in Figure 6b associated with different time step sizes.From Figure 6, the discrete energy-decay property is indeed satisfied and the discrete energy achieves the minimum value.This observation numerically verifies theoretical results of the unconditional stability of discrete energy.In fact, all the cases presented in Figure 6b relax the four separate square shapes to the circular shape.Since the procedure is similar to Figure 6a, we omit the detailed pictures here.
Figures 8-10 show the snapshots of interfaces between binary fluids from left to right across the interface  = 1.We clearly see that the mobile moves toward the center region with high hydraulic conductivity, for example at times  = 0.5 and  = 0.6 for Case I depicted in Figure 8, and deforms into a flatter shape at the front propagating in porous media.The mobile can seek pathways toward high permeability across the interface  = 1 as expected due to the effect of hydraulic conductivity in porous media region for all three cases.After the deformed droplet smoothly goes across the interface  = 1, the droplet is further deformed in the porous media and then reached the steady shape.
Figure 11 shows the streamline and the contours of norm of velocity with the brighter color indicating higher speed of the flow.We observe that the velocity is larger inside the trapezoidal region of case I, inside the two horizontal cracks of case II, and outside the two vertical rocks of case III, respectively.The fluid flows toward the high permeability areas for all cases.The expected flow patterns are observed for different permeability values.(5.5) Figure 12 shows the dynamics of the droplet falling and turning into an ellipse in time under the effect of gravity force and surface tension.Here the green line indicates the interface between two subdomains.The droplet has a more pronounced deformation once it crosses the interface  = 1 entering free flow region.We also consider the curve interface between two subregions.Figure 13 describes the process of the droplet falling across the curve interface.Comparing Figures 12 and 13, the interface morphology shows a significant effect on the droplet deformation, which implies the importance of interface between free flow region and porous medium to the evolution of droplet.
Example 7 (Buoyancy-driven flow).In this experiment, we further validate the applicability of proposed numerical scheme by simulating a rising bubble in a heavier medium with respect to different density variations.The computational domain is chosen as Ω = [0, 1] × [0, 2] associated with the interface boundary Γ = [0, 1] × {1} between the porous media region Ω  = [0, 1] × [1,2]   Figure 14 shows the evolution of bubble from conduit to matrix across the interface for the density ratio  1 :  2 = 1 : 50, where  = 1 (red color) represents the lighter fluid and  = −1 (blue color) represents the heavier fluid.We can clearly observe the deformation of the bubble under the effect of surface tension, the sharp interface between the two subdomains, and the porous media.We can see that the numerical simulations are consistent with the results in [27].
Furthermore, evolutions of two coaxial bubbles [13,85] are considered with the following initial condition where ( 1 ,  1 ) = (0.5, 0.25) and ( 2 ,  2 ) = (0.5, 0.6).The other parameters are same with the one bubble rising case.The dynamics of rising bubbles are recorded in Figure 15 which demonstrates the coalescence and ascension of two bubbles.It can be observed that the two merging droplets eventually evolute into the stable shape.Similar numerical behaviour was reported in [23].The reasonable results further validate the numerical schemes.

Conclusions
In this paper, a fully decoupled finite element scheme with discrete energy law is developed for two-phase flows with different densities and viscosities in the free fluid region and porous media region, which are coupled by suitable interface conditions.An energy law is derived for the novel weak formulation.The CHNS equation and CHD equation on the two different subdomains are decomposed base on the seven interface conditions.IEQ technique is utilized to deal with the stiffness of the phase function associated with the interfacial width in the phase function.The application of these ideas result in a CHNS equation and a CHD equation that can be resolved separately at each time step.We further decouple the CHNS equation into CH equation and Navier-Stokes equation in conduit and decouple the CHD equation into CH equation and Darcy equation in the matrix.Furthermore, the modified pressure correction and grad-div stabilization technique are utilized to deal with the difficulties arising from the decoupling of velocity and pressure in Navier-Stokes equation.At the end, the proposed fully-decoupled scheme only needs to solve a sequence of linear equations at each discrete time level.The full discretization is constructed in the framework of the finite element method.The discrete energy law is analyzed for the proposed method.The efficiency and unconditional stability of the proposed schemes are verified by numerical simulations.We have also investigated the dynamic behavior of challenging model scenarios to validate the mathematical model and numerical method. )

Figure 1 .
Figure 1.A sketch of the porous median domain Ω  , fluid domain Ω  , and the interface Γ.

Figure 4 .
Figure 4.The evolution of discrete energy and mass conservation of cross shape.

Figure 6 .
Figure 6.The evolution of discrete energy of merging drops.(a) Dynamic of phase variable.(b) Energy evolution.