A Hele-Shaw-Cahn-Hilliard model for incompressible two-phase flows with different densities

TL;DR

This paper develops a diffuse interface model for incompressible two-phase flows with different densities in a Hele-Shaw cell, deriving it from existing models, proving solution existence, and demonstrating numerical simulations of complex flow phenomena.

Contribution

It introduces a new diffuse interface model for two-phase flows with density differences in Hele-Shaw geometry, extending previous models and analyzing their solutions.

Findings

Derived a simplified diffuse interface model for non-matched densities.

Proved the existence of weak solutions for the model.

Performed numerical simulations showing bubble rise and fingering instabilities.

Abstract

Topology changes in multi-phase fluid flows are difficult to model within a traditional sharp interface theory. Diffuse interface models turn out to be an attractive alternative to model two-phase flows. Based on a Cahn-Hilliard-Navier-Stokes model introduced by Abels, Garcke and Gr\"{u}n (Math. Models Methods Appl. Sci. 2012), which uses a volume averaged velocity, we derive a diffuse interface model in a Hele-Shaw geometry, which in the case of non-matched densities, simplifies an earlier model of Lee, Lowengrub and Goodman (Phys. Fluids 2002). We recover the classical Hele-Shaw model as a sharp interface limit of the diffuse interface model. Furthermore, we show the existence of weak solutions and present several numerical computations including situations with rising bubbles and fingering instabilities.

Equations

$$\mathrm{div}\mathbf{v}$$

$$12{\eta }_1\mathbf{v}=$$

$$12{\eta }_2\mathbf{v}=$$

$${\left[\mathbf{v}\right]}_1^2\cdot \mathbf{\nu }$$

$${\left[p\right]}_1^2$$

$$\mathcal{V}$$

$$\mathrm{div}\mathbf{v}$$

$${\partial }_t({\rho }_{\ast }(\phi )\mathbf{v})+\mathrm{div}({\rho }_{\ast }(\phi )\mathbf{v}\otimes \mathbf{v})-\nabla p$$

$$-\frac{\epsilon }{\mathrm{Ca}}\mathrm{div}(\nabla \phi \otimes \nabla \phi )+\frac{{\overline{\rho }}_2-{\overline{\rho }}_1}{2{\overline{\rho }}_2}\mathrm{div}(m(\phi )\mathbf{v}\otimes \nabla \mu ),$$

$${\partial }_t^{\bullet }\phi$$

$$\mu$$

$$\mathrm{div}\mathbf{v}$$

$$12{\eta }_{\ast }(\phi )\mathbf{v}$$

$${\partial }_t^{\bullet }\phi$$

$$\mu$$

$$\mathrm{div}\mathbf{v}-\frac{\alpha }{\mathrm{Pe}}\Delta \mu$$

$$\rho (c)({\partial }_{t}c+\nabla c\cdot \mathbf{v})-\frac{1}{\mathrm{Pe}}\Delta \mu$$

$$\mathbf{v}+\frac{1}{12\eta (c)}\left(\nabla p+\frac{\mathrm{Ch}}{\mathrm{Ma}}\mathrm{div}(\rho (c)\nabla c\otimes \nabla c)-\rho (c)\hat{\mathbf{g}}\right)$$

$$\mu -f_0^{\prime }(c)+\frac{\mathrm{Ch}}{\rho (c)}\mathrm{div}(\rho (c)\nabla c)-\mathrm{Ma}\alpha p$$

$$\mathbf{v}=u_1{\mathbf{v}}_1+u_2{\mathbf{v}}_2.$$

$$\mathrm{div}\mathbf{v}$$

$${\partial }_t(\rho (\phi )\mathbf{v})+\mathrm{div}(\rho (\phi )\mathbf{v}\otimes \mathbf{v})$$

$$-\sigma \epsilon \mathrm{div}(\nabla \phi \otimes \nabla \phi )+\frac{{\overline{\rho }}_2-{\overline{\rho }}_1}{2}\mathrm{div}(m(\phi )\mathbf{v}\otimes \nabla \mu ),$$

$${\partial }_t^{\bullet }\phi$$

$$\mu$$

$$\mathbf{G}(\phi )=\rho (\phi )g\hat{\mathbf{g}},$$

$$\frac{d}{dt}{\int }_{\Omega }\left(\frac{\rho }{2}{\left|\mathbf{v}\right|}^2+\frac{\sigma }{\epsilon }\Psi (\phi )+\frac{\sigma \epsilon }{2}{\left|\nabla \phi \right|}^2\right)dx+{\int }_{\Omega }\left(2\eta {\left|\mathrm{D}\mathbf{v}\right|}^2+m{\left|\nabla \mu \right|}^2\right)dx={\int }_{\Omega }\mathbf{G}\cdot \mathbf{v}dx,$$

$${\partial }_{\mathbf{\nu }}\phi =\nabla \phi \cdot \mathbf{\nu }=0,{\partial }_{\mathbf{\nu }}\mu =0,\mathbf{v}=0$$

$$x_i$$

$$x_3$$

$$\epsilon ={\epsilon }_{\ast }L,\mu =\frac{\sigma }{L}{\mu }_{\ast },\mathrm{Pe}=\frac{V{L}^2}{\sigma m},$$

$${\partial }_{t_{\ast }}\phi +{\nabla }_{\ast }\phi \cdot {\mathbf{v}}_{\ast }$$

$${\mu }_{\ast }$$

$$0$$

$$v_{j,\ast }$$

$$\phi$$

$${\mu }_{\ast }$$

$${\partial }_3{\phi }_j(x_{1,\ast },x_{2,\ast },0)$$

$${\partial }_3{\phi }_j(x_{1,\ast },x_{2,\ast },1)$$

$${\partial }_3^2{\mu }_0={\partial }_3^2{\phi }_0=0,{\partial }_3^2{\mu }_1={\partial }_3^2{\phi }_1=0.$$

$${\partial }_{t_{\ast }}{\phi }_0+{\partial }_1{\phi }_0{v}_{1,0}+{\partial }_2{\phi }_0{v}_{2,0}$$

$${\mu }_0$$

$${\partial }_{t_{\ast }}{\phi }_0+{\partial }_1{\phi }_0\overline{v_{1,0}}+{\partial }_2{\phi }_0\overline{v_{2,0}}$$

Full text

**Key words. ** Hele–Shaw flows, multi-phase flows, Cahn–Hilliard model, diffuse interfaces, sharp interface limit, isogeometric analysis.

**AMS subject classification. ** 35Q35, 76D27, 76D45, 76T99, 76S05, 35D30.

1 Introduction

Interfaces in fluid flow play an important role in many applications. A mathematical description of such phenomena typically involves highly nonlinear equations due to the unknown interfaces. Two-phase flow in the special case of a Hele–Shaw cell which involves the slow flow of a fluid between two parallel flat plates which are fixed at a small distance apart still contains many ingredients of more complicated systems. Especially interesting instabilities like the Saffman–Taylor fingering instability [48] can occur and this instability has important applications in technology. We refer for example to the analogy of the Hele–Shaw cell to instabilities that appear when one tries to extract residual oil from a porous rock, see [52, 53]. Water is pumped into the porous rocks to direct the oil to the producing wells, but it was observed that a lot of oil remained in the ground when the water appeared at the wells. One explanation of this phenomenon was attributed to the instabilities of the oil-water interface, which allowed the water to flow through the porous rocks without displacing much of the oil.

In a sharp interface description, the Hele–Shaw model is given as follows. The overall domain $\Omega$ is occupied by two fluids, modeled as time-dependent disjoint regions $\Omega_{1}$ and $\Omega_{2}$, which are separated by a time-dependent hypersurface $\Sigma$. Introducing the fluid velocity $\bm{v}$, the viscosities $\eta_{i}$ and densities $\overline{\rho}_{i}$, $i=1,2$, (which can be different in the two phases), the pressure $p$, the gravity vector $\bm{g}=g\hat{\bm{g}}$ with modulus $g$ and unit vector $\hat{\bm{g}}$, the unit normal $\bm{\nu}$ on $\Sigma$ pointing into $\Omega_{2}$, one has to study the following set of equations:

[TABLE]

Here, $\left[\cdot\right]_{1}^{2}$ denotes the jump across the interface, $\sigma$ is the surface tension, $\kappa$ is the mean curvature and $\mathcal{V}$ is the normal velocity. These equations can be derived from more complete models involving the (Navier–)Stokes equations in situation where the flow is slow (small Reynolds number) and is confined between two parallel plates at a small distance apart, see for example [43].

Such a sharp interface description breaks down when the topology of the interface changes. As a remedy, various diffuse interface models have been introduced to describe incompressible two-phase flows. A first model restricted to equal densities was introduced by Hohenberg and Halperin [32], while a first diffuse interface model for two-phase flow allowing for a density contrast was introduced by Lowengrub and Truskinowsky [40]. However, the model in [40] leads to a velocity field which is not divergence-free (solenoidal) although both individual fluids are. Let us remark that Lowengrub and Truskinowsky used a mass averaged velocity field to define their diffuse interface model. More recently, Abels, Garcke and Grün [3] introduced a diffuse interface model with a divergence-free velocity field which also allows for different densities.

We base our derivation of a diffuse interface model for a Hele–Shaw cell on the Cahn–Hilliard–Navier–Stokes model of [3], which in nondimensional form reads as

[TABLE]

with suitable boundary and initial conditions. Here, $\varphi$ is an order parameter which represents the difference in the volume fractions, such that $\{\varphi=-1\}$ represents fluid 1 and $\{\varphi=1\}$ represents fluid 2. The function $\rho_{*}(\varphi)=\frac{\overline{\rho}_{2}-\overline{\rho}_{1}}{2\overline{\rho}_{2}}\varphi+\frac{\overline{\rho}_{2}+\overline{\rho}_{1}}{2\overline{\rho}_{2}}$ is the nondimensionalized density of the fluid mixture, $\mathrm{D}\bm{v}=\frac{1}{2}(\nabla\bm{v}+(\nabla\bm{v})^{\top})$ is the symmetric gradient for the volume-averaged velocity $\bm{v}$, $p$ denotes the pressure, $\eta_{*}(\varphi)=\frac{\eta_{2}-\eta_{1}}{2\eta_{2}}\varphi+\frac{\eta_{2}+\eta_{1}}{2\eta_{2}}$ is the nondimensionalized viscosity of the mixture, $\mathrm{Re}$ denotes the Reynolds number, $\mathrm{Ca}$ denotes the capillary number, $\varepsilon>0$ is a (small) parameter related to the thickness of the interfacial regions, $\Psi^{\prime}$ is the derivative of a potential $\Psi$ which has equal minima at $\pm 1$, $\mu$ is the chemical potential, $m(\varphi)$ is a non-negative mobility which, in the case of a constant mobility $m(\varphi)=m$, can be seen as the reciprocal of the Pélect number $\mathrm{Pe}$, $\partial_{t}^{\bullet}\varphi=\partial_{t}\varphi+\nabla\varphi\cdot\bm{v}$ is the material derivative of $\varphi$, and $\bm{G}$ denotes an external body force.

We will show, via a formal asymptotic analysis, that for slow flow in a Hele–Shaw cell geometry the above model leads to a Hele–Shaw–Cahn–Hilliard model

[TABLE]

which inherits a divergence-free velocity field from the Cahn–Hilliard–Navier–Stokes model (1.2). In this paper, we will study the model (1.3) in detail both from an analytical and also from a numerical point of view.

An earlier Hele–Shaw–Cahn–Hilliard model was introduced by Lee, Lowengrub and Goodman [36, 37]. However, they used the Cahn–Hilliard–Navier–Stokes model of Lowengrub and Truskinovsky [40] as a basis and obtained

[TABLE]

where $\mathrm{Pe}$ is the Pélect number, $\mathrm{Ch}$ is the Cahn number, $\mathrm{Ma}$ is the Mach number, $c$ is the mass concentration of fluid 1, so that $\{c=1\}$ represents fluid 1 and $\{c=0\}$ represents fluid 2, $\rho(c)$ is the total density, $\alpha$ is the difference between the reciprocals of the actual mass densities of the fluid, $f_{0}(c)=c^{2}(1-c^{2})$ is a potential with two minima at $c=0$ and $c=1$, $\eta(c)=\eta_{1}c+\eta_{2}(1-c)$ is the interpolation of the two viscosities, and $\hat{\bm{g}}$ is the unit vector of gravity. We refer the reader to Section 2.3 for more details.

It is important to note that the velocity $\bm{v}$ in (1.4) is the mass-averaged velocity, which is in contrast to the volume-averaged velocities in (1.2) and (1.3). One observes that the mass-averaged velocity is not divergence-free and that the pressure $p$ enters the equation for the chemical potential (1.4d). These facts make the analysis and the numerical approximation of this model quite involved. We remark that Lee, Lowengrub and Goodman derived (1.3) from (1.4) in the case where a Boussinesq approximation is valid, i.e., the deviation of $\rho$ from its spatial average needs to be small which basically means that the densities of the two fluids are very close. Our derivation however is valid for any density contrast among the fluids.

We spatially approximate the Hele–Shaw–Cahn–Hilliard equations by means of NURBS-based Isogeometric Analysis [16, 33] as it allows a straightforward construction of the finite dimensional function spaces for high order problems [31, 50]. Indeed, in this paper, we formulate the Hele–Shaw–Cahn–Hilliard model (1.3) in terms of the pressure $p$ and order parameter $\varphi$, thus yielding a fourth order problem in the latter variable. In this respect, our finite dimensional function spaces are built out of globally $C^{1}$-continuous B-spline basis functions of degree $2$ [44]. For the time discretization, we use Backward Differentiation Formulas (BDF) of order $2$ [45] with equal order extrapolation of the unknowns to obtain a semi-implicit formulation of the full discrete problem as e.g. in [22].

Finally, we propose and discuss numerical results for two benchmark problems: the rising bubble and viscous fingering tests [34, 36].

The outline of this paper is as follows: In Section 2 we derive (1.3) from (1.2) by means of a formal asymptotic analysis. In Section 3, we derive the sharp interface limit of (1.3) and prove the existence of weak solutions to (1.3). In Section 4 we present the numerical scheme for (1.3) reformulated in terms of the pressure $p$ and the order parameter $\varphi$, and in Section 5 we present and discuss the numerical results.

2 Derivation of the Hele–Shaw–Cahn–Hilliard model

2.1 A Navier–Stokes–Cahn–Hilliard model for incompressible two-phase flows

We start from the volume-averaged velocity model introduced by Abels, Garcke and Grün in [3]: For fluid $i$, $i=1,2$, let $\rho_{i}$ denote the actual mass density, $\overline{\rho}_{i}$ the density of a pure component, $u_{i}:=\frac{\rho_{i}}{\overline{\rho}_{i}}$ the volume fraction, $\bm{v}_{i}$ the individual velocity, and $\eta_{i}$ the viscosity. The volume-averaged velocity for the fluid mixture is defined as

[TABLE]

We define the order parameter $\varphi$ as the difference in the volume fractions, i.e., $\varphi=u_{2}-u_{1}$, then we obtain the following system of equations:

[TABLE]

Here, $\rho(\varphi)=\frac{\overline{\rho}_{2}-\overline{\rho}_{1}}{2}\varphi+\frac{\overline{\rho}_{2}+\overline{\rho}_{1}}{2}$ is the density of the fluid mixture, $\mathrm{D}\bm{v}=\frac{1}{2}(\nabla\bm{v}+(\nabla\bm{v})^{\top})$ is the symmetric gradient, $p$ denotes the pressure, $\eta(\varphi)=\frac{\eta_{2}-\eta_{1}}{2}\varphi+\frac{\eta_{2}+\eta_{1}}{2}$ is the viscosity of the mixture, $\sigma$ is a constant related to the surface energy density, $\varepsilon>0$ is a (small) parameter related to the thickness of the interfacial regions, $\Psi^{\prime}$ is the derivative of a potential $\Psi$ which has equal minima at $\pm 1$, $\mu$ is the chemical potential, $m(\varphi)$ is a non-negative mobility, $\partial_{t}^{\bullet}\varphi=\partial_{t}\varphi+\nabla\varphi\cdot\bm{v}$ is the material derivative of $\varphi$, and $\bm{G}=\bm{G}(\rho(\varphi))$ denotes a body force which may depend on the density. The example we have in mind refers to the gravitational force and reads:

[TABLE]

where the unit vector $\hat{\bm{g}}$ indicates the direction of gravity and $g$ is the modulus.

The model (2.1) consists of the Navier–Stokes equations coupled with a Cahn–Hilliard system. The capillary forces due to surface tension are modeled by the term $\sigma\varepsilon\,\mathrm{div}\,(\nabla\varphi\otimes\nabla\varphi)$, and the term $\frac{\overline{\rho}_{2}-\overline{\rho}_{1}}{2}\,\mathrm{div}\,(m(\varphi)\bm{v}\otimes\nabla\mu)$ accounts for the effects of non-matched fluid densities. We point out that the simple form for continuity equation (2.1a) is due to the choice of $\bm{v}$ as the volume-averaged velocity, when compared for instance to the approach of Antanovskii [8] and Lowengrub and Truskinovsky [40], where a mass-averaged velocity is used and leads to a more complex expression for the continuity equation.

Furthermore, (2.1) satisfies the energy equality

[TABLE]

when we complement (2.1) with the boundary conditions

[TABLE]

on the boundary $\Gamma$ of the bounded domain $\Omega\subset\mathbb{R}^{d}$, $d=1,2,3$, under consideration. Here $\frac{\rho}{2}\left|\bm{v}\right|^{2}$ denotes the kinetic energy of the fluid mixture, $\frac{1}{\varepsilon}\Psi(\varphi)+\frac{\varepsilon}{2}\left|\nabla\varphi\right|^{2}$ is the Ginzburg–Landau energy density, and its product with $\sigma$ approximates the surface energy density in the limit $\varepsilon\to 0$, cf. [42]. The total energy (consisting of the kinetic energy and the surface energy) is dissipated by viscous stress and diffusion, given by the second integral on the left-hand side. We also obtain energy contributions via the body force $\bm{G}$ in the form of the right-hand side. For the existence of weak solutions to (2.1) we refer to the work of Abels, Depner and Garcke [1, 2].

2.2 Nondimensionalization and the Hele–Shaw approximation

We now follow the procedure outlined in [43, Chapter 4], and consider the Navier–Stokes–Cahn–Hilliard equations (2.1) with the body force $\bm{G}$ given as in (2.2) in a domain $\Omega\subset\mathbb{R}^{3}$ which occupies a region in between two rigid walls, one at $\{x_{3}=0\}$ and one at $\{x_{3}=H\}$, for some $H>0$. To be precise, we assume that $\Omega=\Omega^{\prime}\times(0,H)$ with a domain $\Omega^{\prime}\subset\mathbb{R}^{2}$.

We consider a characteristic length $L$ and a characteristic velocity $V$. We denote by $\delta=\frac{H}{L}\ll 1$ the ratio between the height $H$ and the characteristic length $L$ in the $(x_{1},x_{2})$-directions. We set $T=\frac{V}{L}$ as the characteristic time scale and, due to the geometry of the domain under consideration, we rescale the third component of the spatial variable and the third component of the velocity by $\delta$. That is,

[TABLE]

where the variables with $*$-subscript denote nondimensionalized variables. In the following, we use the notation $\partial_{i}:=\partial_{x_{i,*}}$ for $i=1,2,3$, and $\nabla_{*}=(\partial_{1},\,\partial_{2},\,\partial_{3}\,)^{\top}$. Let us consider a constant mobility $m(\varphi)=m$ and define

[TABLE]

where $\mathrm{Pe}$ is the Pélect number. Then, the Cahn–Hilliard part (2.1c)-(2.1d) and the Neumann boundary conditions become

[TABLE]

Since $\bm{v}_{*}$, $\varphi$, and $\mu_{*}$ depend on $\delta$ via the third spatial component, we assume that there exists an asymptotic expansion in $\delta$, i.e.,

[TABLE]

We will substitute these expansions into (2.3) and solve them order by order. On the surfaces $\{x_{3,*}=0\}$ and $\{x_{3,*}=1\}$, as $\nu_{1}=\nu_{2}=0$ we obtain from (2.3c)-(2.3d) for all orders $j=0,1,2\dots$,

[TABLE]

Meanwhile, to orders $\mathcal{O}(\frac{1}{\delta^{2}})$ and $\mathcal{O}(\frac{1}{\delta})$, we obtain from (2.3a)-(2.3b),

[TABLE]

Upon integrating with respect to $x_{3,*}$ and using the conditions (2.4) we have that $\mu_{0}$, $\mu_{1}$, $\varphi_{0}$, and $\varphi_{1}$ are independent of $x_{3,*}$. Then, to order $\mathcal{O}(1)$ we obtain from (2.3a)-(2.3b),

[TABLE]

Integrating the above equations with respect to $x_{3,*}$ from [math] to $1$, and using the condition (2.4) leads to

[TABLE]

where

[TABLE]

denotes the components of the mean velocity $\overline{\bm{v}}=(\overline{v_{1,0}},\overline{v_{2,0}})^{\top}$. In particular, in the limit $\delta\to 0$, we obtain a two-dimensional Cahn–Hilliard system convected by the mean velocity $\overline{\bm{v}}$ and complemented with Neumann boundary conditions on $\partial\Omega^{\prime}$ from (2.3c)-(2.3d).

For the Navier–Stokes part (2.1a)-(2.1b), the continuity equation after the transformation becomes

[TABLE]

From the above computation with the Cahn–Hilliard part, we expect that a scale factor of $\frac{1}{\delta^{2}}$ will appear from the term $\,\mathrm{div}\,(2\eta(\varphi)\mathrm{D}\bm{v})$ in (2.1b). Thus, in order to retain the pressure, the body force and the capillary term in the limit $\delta\to 0$, we set

[TABLE]

and define

[TABLE]

where the capillary number $\mathrm{Ca}$ is the ratio between viscous forces and surface tension, the Reynolds number $\mathrm{Re}$ is the ratio between inertial forces and viscous forces, and the Bond number $\mathrm{Bo}$ is the ratio between gravitational forces and viscous forces. We now nondimensionalize the first component of the momentum equation (2.1b):

[TABLE]

where $\lambda_{\rho}=\frac{\overline{\rho}_{1}}{\overline{\rho}_{2}}$ denotes the density ratio. We point out that, in the case $\overline{\rho}_{2}\geq\overline{\rho}_{1}$, i.e., fluid 2 is the heavier fluid, then the Atwood number $\mathrm{A}:=\frac{\overline{\rho}_{2}-\overline{\rho}_{1}}{\overline{\rho}_{2}+\overline{\rho}_{1}}$ can be expressed as $\mathrm{A}=\frac{1-\lambda_{\rho}}{1+\lambda_{\rho}}$. Similarly, for the second component of the momentum equation (2.1b) we obtain

[TABLE]

Meanwhile, for the third component of the momentum equation (2.1b) we have

[TABLE]

The no-slip boundary condition becomes

[TABLE]

and thus on the surfaces $\{x_{3,*}=0\}$ and $\{x_{3,*}=1\}$ we have

[TABLE]

The procedure to obtain a set of equations from the Navier–Stokes part in the limit $\delta\to 0$ is similar to what we have performed for the Cahn–Hilliard part. In the following, we will only sketch the details. Let $p_{*}=p_{0}+\delta p_{1}+\delta^{2}p_{2}+\text{ h.o.t.}$ denote an asymptotic expansion of the pressure. Due to the fact that $\varphi_{0}$ and $\varphi_{1}$ are independent of $x_{3,*}$, to order $\mathcal{O}(\frac{1}{\delta^{3}})$ we find that (2.10) yields

[TABLE]

and thus $p_{0}$ is independent of $x_{3,*}$. Similarly, thanks to the fact that $\partial_{3}\varphi_{0}=\partial_{3}\varphi_{1}=\partial_{3}\mu_{0}=0$, to order $\mathcal{O}(\frac{1}{\delta^{2}})$ we obtain from (2.8) and (2.9),

[TABLE]

for $i=1,2$. Integrating the above equation with respect to $x_{3,*}$ from [math] to $1$, and using the conditions (2.4) and (2.11) leads to

[TABLE]

for $i=1,2$. Dividing by $\eta_{*}$ and integrating over $s$ from [math] to $1$ leads to the equation for the mean velocity $\overline{\bm{v}}=(\overline{v_{1,0}},\overline{v_{2,0}})$, (recall (2.6)):

[TABLE]

for $i=1,2$. Furthermore, thanks to (2.7) and the condition (2.11), we obtain

[TABLE]

Thus, from the Navier–Stokes part (2.7)-(2.10) we obtain

[TABLE]

where $\overline{\nabla}f=(\partial_{1}f,\partial_{2}f)^{\top}$ denotes the two-dimensional gradient of a scalar function $f$, and $\overline{\,\mathrm{div}\,}\bm{f}=\partial_{1}f_{1}+\partial_{2}f_{2}$ denotes the two-dimensional divergence of a vector function $\bm{f}$. Here, we reuse the notation $\hat{\bm{g}}=(\hat{g}_{1},\hat{g}_{2})^{\top}$.

Dropping the subscripts and combining (2.5) and (2.12) leads to the following nondimensionalized Hele–Shaw–Cahn–Hilliard model:

[TABLE]

where $\,\mathrm{div}\,\cdot$, $\nabla\cdot$, and $\Delta\cdot$ are to be interpreted as the two-dimensional divergence, gradient and Laplace operators, respectively. Using the identity

[TABLE]

and defining the modified pressures

[TABLE]

we obtain two variants of (2.13b):

[TABLE]

In the case where there is no density contrast, i.e., $\overline{\rho}_{1}=\overline{\rho}_{2}$, and the gravitational forces are neglected, the model (2.13) with (2.14b) has been studied by Wang and Zhang in [55] concerning strong well-posedness globally in time for two dimensions and locally in time for three dimensions, and by Wang and Wu in [54] concerning long-time behavior and well-posedness in three dimensions with well-prepared data.

If, in addition, there is no viscosity contrast, i.e., $\eta_{1}=\eta_{2}$, then Feng and Wise established the global existence of weak solutions in two and three dimensions via the convergence of a fully discrete and energy stable implicit finite element scheme in [21]. Uniqueness of weak solutions can be shown if additional regularity assumptions on the solutions are imposed, see [21, Thm. 2.4], and the error analysis of the numerical scheme is performed in [39]. For the convergence analysis of finite difference schemes, we refer the reader to [13, 14, 56].

Meanwhile, Bosia, Conti and Grasselli proved that weak solutions to the Cahn–Hilliard–Brinkman model converge to a weak solution of the Hele–Shaw–Cahn–Hilliard model in [11]. The Cahn–Hilliard–Brinkman model is a related system where an addition term of the form $-\,\mathrm{div}\,(\nu\mathrm{D}\overline{\bm{v}})$ is added to the left-hand side of (2.14b). Here, $\mathrm{D}\overline{\bm{v}}:=\frac{1}{2}(\nabla\overline{\bm{v}}+(\nabla\overline{\bm{v}})^{\top})$ is the rate of deformation tensor and $\nu>0$ is the approximation parameter. Error estimates in terms of $\nu$ between the Cahn–Hilliard–Brinkman model and the Hele–Shaw–Cahn–Hilliard have also been derived in two dimensions. A nonlocal version of the results of [11] has been recently established in [18].

Recently, the asymptotic behavior $\varepsilon\to 0$ of global weak solutions ***We point out that the $L^{2}$ temporal regularity for the time derivative $\partial_{t}\varphi$ (written as $\partial_{t}c^{\varepsilon}$) in [20] may be a typo, cf. Theorem 3.1 below. to the Hele–Shaw–Cahn–Hilliard model (2.13) with (2.14a), and the particular scaling $\frac{1}{\mathrm{Pe}}=\varepsilon^{\alpha}$ for $0\leq\alpha<1$ and $\mathrm{Bo}=0$ has been studied by Fei in [20], which employs the varifold approach of Chen [15]; see also [4, 24, 41] and [5, Appendix A]. In Section 3.3 below we will establish the global in time existence of weak solutions to (2.13) (with the variant (2.14a) and a general body force $\bm{G}(\varphi)$ replacing $\mathrm{Bo}\,\rho(\varphi)\hat{\bm{g}}$) for two and three dimensions.

2.3 Comparison with the Lee–Lowengrub–Goodman model

In this section, we compare the model (2.13) with the model of Lee, Lowengrub and Goodman [36]. In the sequel, we will denote the mass-averaged velocity by $\bm{w}$. Let $c$ denote an order parameter distinguishing the two fluid phases, with $\Omega_{1}:=\{c=1\}$ and $\Omega_{2}=\{c=0\}$. Recalling $\rho_{i}$ and $\bm{v}_{i}$ as the actual mass density and individual velocity of fluid $i$, $i=1,2$, the total density $\rho(c)$ and mass-averaged velocity $\bm{w}$ are defined as

[TABLE]

Let $\eta(c)=\eta_{1}c+\eta_{2}(1-c)$ denote the interpolation of the two viscosities. We introduce the coefficient

[TABLE]

and let $g$ denote the modulus of the gravity vector $\bm{g}=g\hat{\bm{g}}$ with unit vector $\hat{\bm{g}}$. Then, the nondimensionalized Hele–Shaw–Cahn–Hilliard equations of [36, Equ. (2.18)-(2.21)] are

[TABLE]

where $f_{0}=c^{2}(1-c^{2})$ has two minima at $c=0$ and $c=1$, and the dimensionless constants $\mathrm{Pe}$, $\mathrm{Ch}$ and $\mathrm{Ma}$ are the Pélect number, the Cahn number and the Mach number, respectively.

Here we point out that the continuity equation (2.13a) and the equation for the chemical potential (2.13d) in the volume-averaged model (2.13) are considerably simpler than their counterparts (2.16a) and (2.16d) in the mass-averaged model (2.16). In particular, the pressure appears explicitly in (2.16d) and compressibility effects may be introduced as the mass-averaged velocity $\bm{w}$ need not be solenoidal. In contrast, these features are not present in (2.13).

3 Analysis of the volume-averaged model

3.1 Sharp interface asymptotics

We now consider the sharp interface asymptotics of the nondimensional model (2.13) (using $\bm{v}$ to denote the averaged velocity $\overline{\bm{v}}$ and $\sigma$ to denote the reciprocal of the capillary number $\mathrm{Ca}$) in the following setting:

Assumption 3.1**.**

• •

We set $\mathrm{Pe}=\frac{1}{\varepsilon}$ and consider a more general function $\bm{G}$ replacing the term $\mathrm{Bo}\,\rho(\varphi)\hat{\bm{g}}$, where $\bm{G}$ depends only on $\varphi$ but not high order derivatives.

• •

We assume that there is a family $(\varphi_{\varepsilon},\bm{v}_{\varepsilon},p_{\varepsilon},\mu_{\varepsilon})_{\varepsilon>0}$ of solutions to (2.13), which are sufficiently smooth. For small $\varepsilon$, the domain $\Omega$ can be divided into two open subdomains $\Omega^{\pm}(\varepsilon)$, separated by an interface $\Sigma(\varepsilon)$, given as the zero-level set of $\varphi_{\varepsilon}$, that does not intersect with $\partial\Omega=\Gamma$.

• •

We assume that $(\varphi_{\varepsilon},\bm{v}_{\varepsilon},p_{\varepsilon},\mu_{\varepsilon})_{\varepsilon>0}$ have an asymptotic expansion in $\varepsilon$ in the bulk regions away from $\Sigma(\varepsilon)$ (the outer expansion), and another expansion in the interfacial region close to $\Sigma(\varepsilon)$ (the inner expansion).

• •

We assume that the zero level sets of $\varphi_{\varepsilon}$ converge to a limiting hypersurface $\Sigma$ moving with normal velocity $\mathcal{V}$ as $\varepsilon\to 0$.

• •

We rescale the potential $\Psi$ such that

[TABLE]

For example, the classical quartic double-well potential $\Psi(s)=\frac{1}{4}(1-s^{2})^{2}$ is rescaled to $\Psi(s)=\frac{3}{2\sqrt{2}}\frac{1}{4}(1-s^{2})^{2}$.

The equations we study are

[TABLE]

The idea of the method is to plug the outer and inner expansions in the model equations and solve them order by order, and in addition we have to define a suitable region where these expansions should match up. For $\alpha=-2,-1,0,1,\dots$, we will use the notation $\eqref{HSCH:div}_{O}^{\alpha}$ and $\eqref{HSCH:div}_{I}^{\alpha}$ to denote the terms resulting from the order $\alpha$ outer and inner expansions of (3.2a), respectively.

3.1.1 Outer expansion

We assume that $(\bm{v}_{\varepsilon},p_{\varepsilon},\varphi_{\varepsilon},\mu_{\varepsilon})$ have the following outer expansions

[TABLE]

To leading order $\eqref{HSCH:pressure}_{O}^{-1}$ we obtain

[TABLE]

and so $p_{-1}$ is constant in the bulk regions. Meanwhile $\eqref{HSCH:mu}_{O}^{-1}$ gives

[TABLE]

The stable solutions to the above equation are the minima of $\Psi$, which yields that $\varphi_{0}=\pm 1$. This allows us to define the bulk fluid domains $\Omega_{1}:=\{\varphi(x)=-1\}$ and $\Omega_{2}:=\{\varphi(x)=1\}$. To leading order we obtain from $\eqref{HSCH:div}_{O}^{0}$

[TABLE]

and to first order we obtain from $\eqref{HSCH:pressure}_{O}^{0}$

[TABLE]

3.1.2 Inner expansions

By assumption, $\Sigma$ is the limiting hypersurface of the zero level sets of $\varphi_{\varepsilon}$. In order to study the limiting behavior close to $\Sigma$ we introduce a new coordinate system, which involves the signed distance function $d(x)$ to $\Sigma$. Setting $z=\frac{d}{\varepsilon}$ as the rescaled distance variable to $\Sigma$, and using the convention that $d(x)<0$ in $\Omega_{1}$, and $d(x)>0$ in $\Omega_{2}$, we see that the gradient $\nabla d$ points from $\Omega_{1}$ to $\Omega_{2}$, and we may use $\nabla d$ on $\Sigma$ to denote the unit normal of $\Sigma$, pointing from $\Omega_{1}$ to $\Omega_{2}$.

Let $\alpha(t,s)$ denote a parametrization of $\Sigma$ with tangential coordinates $s$, and let $\bm{\nu}$ denote the unit normal of $\Sigma$, pointing into $\Omega_{2}$. Then, in a tubular neighborhood of $\Sigma$, for a sufficiently smooth function $f(x)$, we have

[TABLE]

In this new $(t,s,z)$-coordinate system, the following change of variables apply, see [29],

[TABLE]

where $\mathcal{V}$ is the normal velocity of $\Sigma$, $\nabla_{\Sigma}g$ denotes the surface gradient of $g$ on $\Sigma$ and h.o.t. denotes higher order terms with respect to $\varepsilon$. In particular, we have

[TABLE]

where $\kappa=-\,\mathrm{div}\,_{\Sigma}\bm{\nu}$ is the mean curvature of $\Sigma$. If $\bm{v}$ is a vector-valued function with $\bm{V}(t,s,z)=\bm{v}(x)$ for $x$ in a tubular neighborhood of $\Sigma$, then we obtain

[TABLE]

The inner variables of $(\bm{v}_{\varepsilon},p_{\varepsilon},\varphi_{\varepsilon},\mu_{\varepsilon})$ are denoted as $(\bm{V}_{\varepsilon},P_{\varepsilon},\Phi_{\varepsilon},\Xi_{\varepsilon})$ with the inner expansion

[TABLE]

Since the zero level sets of $\varphi_{\varepsilon}$ converge to $\Sigma$, we additionally impose that

[TABLE]

In order to match the inner expansions valid in the interfacial region to the outer expansions of Section 3.1.1 we employ the matching conditions, see [29],

[TABLE]

where $f_{0}^{\pm}(t,x):=\lim_{\delta\searrow 0}f_{0}(t,x\pm\delta\bm{\nu})$ for $x\in\Sigma$. For the pressure, we have

[TABLE]

We will employ the following notation: Let $\delta>0$ and for $x\in\Sigma$ with $x-\delta\bm{\nu}\in\Omega_{1}$ and $x+\delta\bm{\nu}\in\Omega_{2}$, we denote the jump of a quantity $f$ across the interface by

[TABLE]

Then, the expansions of (3.2a), (3.2c) and (3.2d) in terms of the inner variables are

[TABLE]

For the tensor product $\varepsilon\,\mathrm{div}\,(\nabla\varphi\otimes\nabla\varphi)$ we obtain the formula

[TABLE]

so that the expansion of (3.2b) becomes

[TABLE]

3.1.3 Expansions to leading order

To leading order we obtain from $\eqref{HSCH:mu}_{I}^{-1}$

[TABLE]

This is a second order equation in $z$ and together with the conditions $\lim_{z\to\pm\infty}\Phi_{0}(t,s,z)=\pm 1$, and $\Phi_{0}(t,s,0)=0$ we obtain a unique solution $\Phi_{0}(z)$ to (3.14) that is independent of $s$ and $t$, i.e., (3.14) can be viewed as an ordinary differential equation in $z$. For the double-well potential $\Psi(s)=\frac{1}{4}(1-s^{2})^{2}$, the unique solution is given by $\Phi_{0}(z)=\tanh\left(\frac{z}{\sqrt{2}}\right)$. Furthermore, multiplying (3.14) by $\Phi_{0}^{\prime}$, integrating and applying matching conditions (3.6) and (3.7) to $\Phi_{0}$ leads to the so-called equipartition of energy

[TABLE]

By (3.1), we see that

[TABLE]

Then, to leading order $\eqref{HSCH:div}_{I}^{-1}$, we obtain

[TABLE]

which implies that $\bm{V}_{0}\cdot\bm{\nu}$ is independent of $z$. Integrating and applying the matching condition (3.6) to $\bm{V}_{0}$ yields

[TABLE]

Meanwhile, from $\eqref{HSCH:pressure}_{I}^{-2}$ we have

[TABLE]

Taking the scalar product with $\bm{\nu}$ and upon integrating with respect to $z$ leads to

[TABLE]

for some function $\hat{P}$ independent of $z$. Sending $z\to\pm\infty$ and applying the matching condition (3.9) to $P_{-1}$ and (3.7) to $\Phi_{0}$, we see that

[TABLE]

In particular, the constant values of $p_{-1}$ in the bulk phase (see (3.3)) should match. We take $p_{-1}^{\pm}=0$ so that $P_{-1}$ is a function only in $z$ and

[TABLE]

To leading order $\eqref{HSCH:mu}_{I}^{-1}$ gives

[TABLE]

By (3.16), $\bm{V}_{0}\cdot\bm{\nu}$ is independent of $z$, and so upon integrating and apply matching conditions (3.6) to $\Phi_{0}$ and (3.7) to $\Xi_{0}$, we obtain

[TABLE]

This implies that

[TABLE]

3.1.4 Expansions to first order

To first order, we obtain from $\eqref{HSCH:mu}_{I}^{0}$

[TABLE]

Multiplying by $\Phi_{0}^{\prime}$, integrating over $\mathbb{R}$ with respect to $z$ leads to

[TABLE]

Integration by parts, applying the matching conditions (3.6) and (3.7) applied to $\Phi_{0}$, and using that $\Psi^{\prime}(\pm 1)=0$, we see that

[TABLE]

and so the first two terms on the right-hand side of (3.19) are zero. Then, using (3.15) and (3.18), we obtain from (3.19),

[TABLE]

Next, using that $P_{-1}$ and $\Phi_{0}$ depend only on $z$, to first order we obtain from $\eqref{HSCH:pressure}_{I}^{-1}$

[TABLE]

Taking the scalar product with $\bm{\nu}$, integrating and applying the matching condition (3.10) and using (3.15) leads to

[TABLE]

where we used that $\,\mathrm{div}\,_{\Sigma}(\bm{\nu}\otimes\bm{\nu})=-\kappa\bm{\nu}$. Hence, the sharp interface limit of (3.2) is

[TABLE]

Remark 3.1**.**

We point out that the formal asymptotic analysis performed with the degenerate mobility

[TABLE]

will yield the same sharp interface limit (3.21). For more details, we refer to [3, 28].

Remark 3.2**.**

If we use the variant (2.14a) of the velocity equation instead of (3.2b), i.e.,

[TABLE]

then the outer and inner expansions of the pressure $q$ do not require a term scaling with $\frac{1}{\varepsilon}$. That is, we can consider

[TABLE]

as the corresponding outer and inner expansions, respectively. While the analysis for the outer expansions remains unchanged, from the leading order inner expansion $\eqref{HSCH:pressure:alt}_{I}^{-1}$ we obtain after taking the scalar product with $\bm{\nu}$ and integrating, and using (3.19) and (3.20),

[TABLE]

which is the nondimensionalized Young–Laplace law (3.21d) for the modified pressure $q$.

3.2 Sharp interface limit for the mass-averaged model

It turns out that in choosing

[TABLE]

and the rescaling $\mu\mapsto\frac{1}{\varepsilon}\mu$ in the mass-averaged model (2.16), that is,

[TABLE]

will result in a sharp interface limit that coincides with (3.21) when we consider $\bm{G}(\varphi)=\rho(\varphi)\hat{\bm{g}}$. We will briefly sketch the details below.

• •

We consider an outer expansion for the pressure $p=p_{0}+\varepsilon p_{1}+\dots$, that is, $p_{-1}=0$. Then, one obtains to leading order $\eqref{LLG:rescaled:chem}_{O}^{-1}$ that $f^{\prime}_{0}(c_{0})=0$, which yields the solutions $c_{0}=0$ or $1$, and the bulk domains can be defined as $\Omega_{1}=\{c_{0}=1\}$ and $\Omega_{2}=\{c_{0}=0\}$. Then, from $\eqref{LLG:rescaled:div}_{O}^{0}$ and $\eqref{LLG:rescaled:velo}_{O}^{0}$ we obtain

[TABLE]

• •

For the inner expansions, we denote the inner variable of $c$ and $\bm{w}$ by $C$ and $\bm{W}$, respectively, and assume that the $\frac{1}{2}$-level sets of $c_{\varepsilon}$ converges to $\Sigma$, which implies that

[TABLE]

Furthermore, we assume that the inner expansion for the pressure $P_{\varepsilon}$ is given as in (3.4), and we alter the matching conditions (3.9), (3.10) to

[TABLE]

• •

To leading order $\eqref{LLG:rescaled:div}_{I}^{-1}$ we obtain $\left[\bm{w}_{0}\right]_{1}^{2}\cdot\bm{\nu}=0$, and to leading order $\eqref{LLG:rescaled:phase}_{I}^{-1}$ we obtain $\mathcal{V}=\bm{w}_{0}\cdot\bm{\nu}$ whenever $\rho>0$ and $\partial_{z}C_{0}\neq 0$.

• •

To leading order $\eqref{LLG:rescaled:velo}_{I}^{-2}$ we obtain

[TABLE]

Integrating and applying the matching conditions for $P_{-1}$ and $\partial_{z}C_{0}$ yields that $P_{-1}(t,s,z)=-\frac{\sigma}{\Sigma_{r}}(\rho(C_{0})(\partial_{z}C_{0})^{2})(t,s,z)$. Then, substituting this into $\eqref{LLG:rescaled:chem}_{I}^{-1}$ gives

[TABLE]

Together with the conditions $\lim_{z\to\infty}C_{0}(t,s,z)=0$, $\lim_{z\to-\infty}C_{0}(t,s,z)=1$, and $C_{0}(t,s,0)=\frac{1}{2}$ this yields a second order ODE in $z$, which implies that we can choose $C_{0}$ to be a function depending only on $z$, and thus $P_{-1}$ only depends on $z$. Multiplying (3.24) by $C_{0}^{\prime}$, applying the product rule to the second term and using the definition of $\alpha$ leads to $0=(f_{0}(C_{0})-\frac{1}{2}\left|C_{0}^{\prime}\right|^{2})^{\prime}$, and upon integrating yields the equipartition of energy

[TABLE]

• •

Lastly, using the fact that $C_{0}$, $P_{-1}$ are independent of $s$ and $t$, we obtain from $\eqref{LLG:rescaled:velo}_{I}^{-1}$

[TABLE]

Taking the scalar product with $\bm{\nu}$, integrating with respect to $z$ and applying the matching conditions for $C_{0}^{\prime}$, we obtain with the help of the equiparition of energy (3.25) and a change of variables $s=C_{0}^{\prime}(z)$,

[TABLE]

3.3 Global existence of weak solutions

In this section, we investigate the existence of weak solutions to the Hele–Shaw–Cahn–Hilliard model (3.2) with the parameters $\varepsilon=\sigma=1$, and rescaling the viscosity by a factor of $\frac{1}{12}$. For a bounded domain $\Omega\subset\mathbb{R}^{d}$, $d=2,3$, with boundary $\Gamma$ and an arbitrary but fixed terminal time $T>0$, we consider

[TABLE]

Here $a>0,b\geq 0$ are constants, $h$ is a prescribed boundary function. Although (2.13) is derived as a model in two dimensions, we include in our analysis the existence theory for three dimensions, which is applicable to the situation of fluid flow in a porous medium. We also point out that, in the case $b=0$, the pressure $q$ is determined up to a constant, and therefore we prescribe in addition that $\int_{\Omega}q\,dx\,=0$ for the case $b=0$. Before presenting the existence result we introduce the notation and useful preliminaries for this section.

Notation.

We set $H:=L^{2}(\Omega)$, $V:=H^{1}(\Omega)$, $H_{\Gamma}:=L^{2}(\Gamma)$. For a (real) Banach space $X$ its dual is denoted as $X^{\prime}$ and $\langle\cdot,\cdot\rangle_{X}$ denotes the duality pairing between $X$ and $X^{\prime}$. The $L^{2}$-inner product on $\Omega$ and on $\Gamma$ will be denoted by $(\cdot,\cdot)$ and $(\cdot,\cdot)_{\Gamma}$, respectively. For convenience, we use the notation $L^{p}:=L^{p}(\Omega)$ and $W^{k,p}:=W^{k,p}(\Omega)$ for any $p\in[1,\infty]$, $k>0$ to denote the standard Lebesgue spaces and Sobolev spaces equipped with the norms $\|\cdot\|_{L^{p}}$ and $\|\cdot\|_{W^{k,p}}$. In the case $p=2$ we use notation $\|\cdot\|_{H}:=\|\cdot\|_{L^{2}}$, $\|\cdot\|_{H_{\Gamma}}:=\|\cdot\|_{L^{2}(\Gamma)}$, and $\|\cdot\|_{V}:=\|\cdot\|_{H^{1}}$. We denote $\mathbb{R}^{d}$-valued functions and spaces consisting of $\mathbb{R}^{d}$-valued functions in boldface, that is $\bm{H}:=(L^{2}(\Omega))^{d}$ and $\bm{V}:=(H^{1}(\Omega))^{d}$. The mean of an integrable function $f:\Omega\to\mathbb{R}$ is defined as $\overline{f}:=\frac{1}{\left|\Omega\right|}\int_{\Omega}f\,dx\,$, and we denote

[TABLE]

For the velocity, we introduce the space

[TABLE]

i.e., $\bm{H}_{\,\mathrm{div}\,}$ is the closure of the space of all divergence free vector fields in $(C^{\infty}_{0}(\Omega))^{d}$ in the $\bm{L}^{2}$-norm. Integration with respect to the Hausdorff measure on $\Gamma$ will be denoted by $\,d\Gamma\,$.

Useful preliminaries.

We have the Sobolev embedding $V\subset L^{r}$ for any $r\in[1,\infty)$ in two dimensions and $r\in[1,6]$ in three dimensions, and the following compact embeddings in dimension $d$ (see [6, Thm. 6.3] and [23, Thm. 11.2, p. 31])

[TABLE]

for any $q\in[1,\infty)$ in two dimensions and $q\in[1,6)$ in three dimensions. We state the Gagliardo–Nirenberg interpolation inequality in dimension $d$ (see [23, Thm. 10.1, p. 27], [19, Thm. 2.1] and [6, Thm. 5.8]): Let $\Omega$ be a bounded domain with Lipschitz boundary, and $f\in W^{m,r}(\Omega)\cap L^{q}(\Omega)$, $1\leq q,r\leq\infty$. For any integer $j$, $0\leq j<m$, suppose there is $\alpha\in\mathbb{R}$ such that

[TABLE]

If $r\in(1,\infty)$ and $m-j-\frac{d}{r}$ is a nonnegative integer, we in addition assume $\alpha\neq 1$. Under these assumptions, there exists a positive constant $C$ depending only on $\Omega$, $m$, $j$, $q$, $r$, and $\alpha$ such that

[TABLE]

We recall the Poincaré inequalities (see for instance [51, Equ. (1.35), (1.37a) and (1.37c)]): There exist positive constants $C_{p}$ depending only on $\Omega$ such that, for all $f\in V$,

[TABLE]

For fixed $b>0$ and a given function $\varphi$, we introduce the operators $\mathcal{N}_{b,\varphi}:V\to V^{\prime}$ and $\mathcal{N}_{0,\varphi}:V\cap L^{2}_{0}\to V^{\prime}_{0}$ by

[TABLE]

Under a boundedness assumption on $\eta$ (see $(\mathrm{A2})$ below), the Lax–Milgram theorem and the Poincaré inequality (3.29) yield that the inverse operator $\mathcal{N}_{b,\varphi}^{-1}$ is well-defined and stable under perturbations. I.e., for any $g\in V^{\prime}$, there exists a unique $u\in V$ such that

[TABLE]

for some positive constant $C$ not depending on $g$ and $u$. Furthermore, given $g_{1},g_{2}\in V^{\prime}$ and the corresponding unique solution $u_{1},u_{2}\in V$ it holds that

[TABLE]

Similarly, using the Poincaré inequality (3.28) with zero mean, the inverse operator $\mathcal{N}_{0,\varphi}^{-1}:V^{\prime}_{0}\to V\cap L^{2}_{0}$ is also well-defined and stable under perturbations.

Assumption 3.2**.**

1. $(\mathrm{A1})$

We assume that $\Omega\subset\mathbb{R}^{d}$, $d=2,3$, is a bounded domain with $C^{3}$-boundary $\Gamma$. 2. $(\mathrm{A2})$

We assume that $\eta\in C^{0}(\mathbb{R})$, $\bm{G}\in C^{0}(\mathbb{R};\mathbb{R}^{d})$ and

[TABLE]

for some positive constants $\eta_{0}$, $\eta_{1}$, $G_{0}$ and $G_{1}$. 3. $(\mathrm{A3})$

We assume that $h\in L^{2}(0,T;H_{\Gamma})$ and $\varphi_{0}\in V$. 4. $(\mathrm{A4})$

The potential $\Psi\in C^{2}(\mathbb{R})$ is nonnegative and satisfies

[TABLE]

for positive constants $c_{0}$, $c_{1}$, $c_{2}$, $c_{3}$, $c_{4}$, $c_{5}$, and exponent $r\geq 2$ for two dimensions and $r\in[2,6)$ for three dimensions.

Theorem 3.1** (Existence of weak solutions).**

Under Assumption 3.2, for

[TABLE]

there exists a quadruple of functions $(\varphi,\mu,q,\bm{v})$ with

[TABLE]

such that $\varphi(0)=\varphi_{0}$ and

[TABLE]

for a.e. $t\in(0,T)$, and for all $\phi\in V$ and $\bm{\zeta}\in\bm{H}$.

Note that by the compact embedding

[TABLE]

the initial value $\varphi(0)$ makes sense as a function in $H$ and thus the initial condition $\varphi_{0}$ is attained. Furthermore, the boundary condition (3.26f) can be attained by choosing $\bm{\zeta}=\eta(\varphi)^{-1}\nabla\phi$ in (3.31a), leading to

[TABLE]

We further point out that the temporal regularity for $\partial_{t}\varphi$ and the pressure $q$ have been similarly observed in the work of [11, 26, 35].

Proof.

The proof is based on a Galerkin approximation. We consider the set of eigenfunctions of the Neumann-Laplacian $\{w_{i}\}_{i\in\mathbb{N}}$ which forms an orthonormal basis of $H$. In [26, §3] it has been shown that $\{w_{i}\}_{i\in\mathbb{N}}$ is also a basis of $H^{2}_{N}$. Let $W_{k}:=\mathrm{span}\{w_{1},\dots,w_{k}\}$ denote the finite dimensional subspace spanned by the first $k$ eigenfunctions, and let $\Pi_{k}$ denote the orthogonal projection into $W_{k}$. We consider a Galerkin ansatz $(\varphi_{k},\mu_{k},q_{k},\bm{v}_{k})_{k\in\mathbb{N}}$ which satisfy $\varphi_{k}=\sum_{i=1}^{k}\alpha_{ik}(t)w_{i}\in W_{k}$,

[TABLE]

where $\eta_{k}:=\eta(\varphi_{k})$, and $q_{k}$ satisfies an elliptic problem whose weak formulation reads as

[TABLE]

Let us define the linear functionals $F_{b,\varphi_{k}},F_{0,\varphi_{k}}\in V^{\prime}$ by

[TABLE]

for all $\zeta\in V$, where $\mu_{k}$ is as defined in (3.33). Then, we may express $q_{k}$ as

[TABLE]

where the operators $\mathcal{N}_{b,\varphi_{k}}$ and $\mathcal{N}_{0,\varphi_{k}}$ are defined in (3.30). Taking the inner product of (3.32) with $w_{j}$, $j=1,\dots,k$, and substituting (3.33), (3.34) and (3.36) leads to a system of nonlinear ODEs for the coefficients $\{\alpha_{ik}(t)\}_{1\leq i\leq k}$. The right-hand side depends continuously on the coefficients $\{\alpha_{ik}(t)\}_{1\leq i\leq k}$. Applying the theory of ordinary differential equations yields the existence of $t_{k}\in(0,T]$ such that the resulting ODE system has a solution $\bm{\alpha}_{k}=(\alpha_{ik})_{1\leq i\leq k}\in C^{0}([0,t_{k});\mathbb{R}^{k})$ that is absolutely continuous. We may define $\mu_{k}$ by the equation (3.33), then $q_{k}$ is defined by (3.36) and $\bm{v}_{k}$ is defined by (3.34).

We now derive a priori estimates for the Galerkin ansatz $(\varphi_{k},\mu_{k},q_{k},\bm{v}_{k})$. In the following, the constant $C>0$ may vary line to line, but it is independent of $k$. For convenience, we denote

[TABLE]

Note that by the assumption $\varphi_{0}\in V$, the growth assumptions $(\mathrm{A4})$ on $\Psi$ and the Sobolev embedding $V\subset L^{p}$ for $p\geq 2$ in two dimensions and $p\in[1,6]$ for three dimensions, there exists a constant $C$ such that

[TABLE]

First estimate.

Substituting $\zeta=q_{k}$ in (3.35), and taking the inner product of (3.32) with $\mu_{k}$, the inner product of (3.33) with $\partial_{t}\varphi_{k}$, and the inner product of (3.34) with $\bm{v}_{k}$, summing and integrating from [math] to $s\in(0,T]$ leads to

[TABLE]

By the growth conditions for $\bm{G}$ in $(\mathrm{A2})$ we see that

[TABLE]

Using the lower bound for $\Psi$ in $(\mathrm{A4})$, we have

[TABLE]

and thus, by the lower bound on the viscosity $\eta$, we obtain from (3.37)

[TABLE]

Applying Young’s inequality to the last term on the right-hand side of (3.38), and then applying Gronwall’s inequality (see [27, Lem. 3.1]) leads to

[TABLE]

For the case $b=0$, we obtain (3.38) without the terms in the $H_{\Gamma}$-norm. Furthermore, the a priori estimate (3.39) guarantees that we can extend the Galerkin ansatz to the whole of $[0,T]$, and thus $t_{k}=T$ for all $k\in\mathbb{N}$.

Second estimate.

Integrating (3.33) and using $(\mathrm{A4})$ leads to

[TABLE]

By (3.39) we have that the mean $\overline{\mu_{k}}$ is bounded uniformly in $L^{\infty}(0,T)$, and thus by the Poincaré inequality (3.28) and the boundedness of $\|\nabla\mu_{k}\|_{L^{2}(0,T;\bm{H})}$, we have

[TABLE]

Third estimate.

We may view (3.33) as an elliptic equation for $\varphi_{k}$:

[TABLE]

Then, the argument in [26, §4.2] yields that $\{\varphi_{k}\}_{k\in\mathbb{N}}$ is bounded uniformly in $L^{2}(0,T;H^{3})$. We will omit the details and refer the reader to [26].

Fourth estimate.

Substituting $\zeta=q_{k}$ in (3.35) leads to

[TABLE]

Case (i) $b>0$.

For this case, Young’s inequality gives $b\|h\|_{H_{\Gamma}}\|q_{k}\|_{H_{\Gamma}}\leq\frac{ba}{2}\|q_{k}\|_{H_{\Gamma}}^{2}+\frac{b}{2a}\|h\|_{H_{\Gamma}}^{2}$, which leads to

[TABLE]

Case (ii) $b=0$.

For this case, we obtain

[TABLE]

For both cases, we obtain an a priori estimate of the form

[TABLE]

By $(\mathrm{A2})$ and $(\mathrm{A3})$, we have that $\bm{G}(\varphi_{k})\in L^{\infty}(0,T;\bm{H})$ and $h\in L^{2}(0,T;H_{\Gamma})$. Thus, we expect that the temporal regularity of $\nabla q_{k}$ will be no greater than the temporal regularity of the product $\mu_{k}\nabla\varphi_{k}$. Let $s\in[1,\infty)$ for two dimensions, then by the Gagliardo–Nirenberg inequality (3.27), we see that

[TABLE]

By Hölder’s inequality and Sobolev embedding, we obtain for two dimensions,

[TABLE]

and so $\mu_{k}\nabla\varphi_{k}\in L^{r}(0,T;\bm{H})$ for $\frac{4}{3}\leq r<2$ in two dimensions. For three dimensions, we obtain analogously

[TABLE]

and so $\mu_{k}\nabla\varphi_{k}\in L^{\frac{8}{5}}(0,T;\bm{H})$. Thus, from (3.43), and using the Poincaré inequality (3.29) for the case $b>0$ or the condition $\overline{q_{k}}=0$ and the Poincaré inequality (3.28) for the case $b=0$, we obtain that

[TABLE]

Fifth estimate.

Using (3.44), in three dimensions, for an arbitrary test function $\zeta\in L^{\frac{8}{3}}(0,T;V)$ we have

[TABLE]

This implies that $\{\Pi_{k}(\bm{v}_{k}\cdot\nabla\varphi_{k})\}_{k\in\mathbb{N}}$ is bounded in $L^{\frac{8}{5}}(0,T;V^{\prime})$. Then, from (3.32) we have that $\{\partial_{t}\varphi_{k}\}_{k\in\mathbb{N}}$ is bounded in $L^{\frac{8}{5}}(0,T;V^{\prime})$. For two dimensions, we have for any $s\in[1,\infty)$,

[TABLE]

and so $\{\Pi_{k}(\bm{v}_{k}\cdot\nabla\varphi_{k})\}_{k\in\mathbb{N}}$ and $\{\partial_{t}\varphi_{k}\}_{k\in\mathbb{N}}$ are bounded in $L^{p}(0,T;V^{\prime})$ for $\frac{4}{3}\leq p<2$.

Compactness.

The above a priori estimates and the application of [49, §8, Corollary 4] yield the existence of a relabelled subsequence $(\bm{v}_{k},q_{k},\varphi_{k},\mu_{k})_{k\in\mathbb{N}}$ such that

[TABLE]

for some function $\xi\in L^{p}(0,T;V^{\prime})$ and

[TABLE]

To deduce that $(\varphi,\mu,q,\bm{v})$ is a weak solution of (2.13) that satisfies (3.31), we argue as follows: Fix $j\in\mathbb{N}$ and $\delta\in C^{\infty}_{c}(0,T)$, multiplying (3.32), (3.33) with $\delta w_{j}$, and integrating in time leads to

[TABLE]

where we used $(\Pi_{k}(\Psi^{\prime}(\varphi_{k})),w_{j})=(\Psi^{\prime}(\varphi_{k}),w_{j})$. On one hand we see that

[TABLE]

On the other hand, the strong convergence of $\nabla\varphi_{k}$ to $\nabla\varphi$ in $L^{2}(0,T;W^{1,s})$ and the fact that $w_{j}\in H^{2}$ shows that

[TABLE]

and so $\delta w_{j}\nabla\varphi_{k}\to\delta w_{j}\nabla\varphi$ strongly in $L^{2}(0,T;\bm{H})$. Together with the weak convergence of $\bm{v}_{k}$ in $L^{2}(0,T;\bm{H})$, we obtain

[TABLE]

Equating (3.46) and (3.47) leads to

[TABLE]

Passing to the limit $k\to\infty$ in (3.45), using the above weak/weak* convergences yields

[TABLE]

We refer to [25, §3.1.2] for the details on how to pass to the limit in the term with $\Psi^{\prime}$. Meanwhile, substituting $\zeta=w_{j}$ in (3.35), then multiplying with $\delta$ and integrating over time, we obtain

[TABLE]

Due to the a.e. convergence of $\varphi_{k}$ to $\varphi$ in $Q$, and the continuity of $\eta$ and $\bm{G}$, we have that $\eta(\varphi_{k})^{-1}\to\eta(\varphi)^{-1}$ and $\bm{G}(\varphi_{k})\to\bm{G}(\varphi)$ a.e. in $Q$. Furthermore, by the boundedness of $\eta$, applying Lebesgue’s dominated convergence theorem yields

[TABLE]

Meanwhile, from the strong convergence $\varphi_{k}\to\varphi$ in $L^{2}(0,T;H)$ we find that

[TABLE]

Then, using the growth assumption $(\mathrm{A2})$ for $\bm{G}$ and the generalized Lebesgue dominated convergence theorem ([47, Thm. 1.9, p. 89], [7, Thm. 3.25, p. 60]), it holds that

[TABLE]

By the Gagliardo–Nirenberg inequality (3.27) we find that

[TABLE]

Thus, from the boundedness of $\{\varphi_{k}\}_{k\in\mathbb{N}}$ in $L^{\infty}(0,T;V)\cap L^{2}(0,T;H^{3})$, we see that $\{\nabla\varphi_{k}\}_{k\in\mathbb{N}}$ is bounded in $L^{3}(0,T;\bm{L}^{3})$. Furthermore, using the strong convergence of $\varphi_{k}$ to $\varphi$ in $L^{2}(0,T;W^{2,s})$ for $s\in[1,6)$, and (3.51) for $m=6$, we obtain

[TABLE]

This implies that

[TABLE]

Then, combining (3.51), (3.52), (3.53) and the weak convergences for $q_{k}$ and $\mu_{k}$, after passing to the limit $k\to\infty$ in (3.50) we obtain

[TABLE]

Next, multiplying (3.34) by $\delta\eta(\varphi_{k})(w_{j_{1}},\dots,w_{j_{d}})^{\top}=:\delta\eta(\varphi_{k})\bm{\zeta}_{j}$ for $1\leq j_{1},\dots,j_{d}\leq k$, passing to the limit $k\to\infty$ yields

[TABLE]

Since (3.48), (3.49), (3.54) and (3.55) hold for arbitrary $\delta\in C^{\infty}_{c}(0,T)$, we infer that $(\varphi,\mu,q,\bm{v})$ satisfies (3.31) with $\phi=w_{j}$ and $\bm{\zeta}=(w_{j_{1}},\dots,w_{j_{d}})^{\top}$. Using that $\{w_{j}\}_{j\in\mathbb{N}}$ is basis of $H^{2}_{N}$ and $H^{2}_{N}$ is dense in $V$, we deduce that (3.31) holds for arbitrary $\phi\in V$ and $\bm{\zeta}\in\bm{H}$. This concludes the proof. ∎

4 Numerical approximation

We briefly describe the numerical approximation of the Hele–Shaw–Cahn–Hilliard problem (2.13) with the variant (2.14a). In particular, by recalling that $\mu=\mu(\varphi)=\frac{1}{\varepsilon}\Psi^{\prime}(\varphi)-\varepsilon\Delta\varphi$ from (2.13d) and $\bm{v}(q,\varphi)=-\frac{1}{12\eta(\varphi)}\left(\nabla q-\mathrm{Bo}\,\rho(\varphi)\hat{\bm{g}}-\frac{1}{\mathrm{Ca}}\mu(\varphi)\nabla\varphi\right)$ from (2.14a), we reformulate the dimensionless problem (2.13) in terms of the (modified) pressure $q$ and order parameter $\varphi$ and endow it with suitable initial and boundary conditions as:

[TABLE]

where $\overline{\Gamma}_{N}\cup\overline{\Gamma}_{D}=\Gamma\equiv\partial\Omega$, $\mathring{\Gamma}_{N}\cap\mathring{\Gamma}_{D}=\emptyset$, and $f_{N}$ is a suitable function. We remark that Problem (4.1) is time-dependent, nonlinear, and it involves a fourth order differential operator in (4.1b). Then, we rewrite for convenience the dimensionless density and viscosity as $\rho(\varphi)=\Theta_{1}\,\varphi+\Theta_{2}$ and $\eta(\varphi)=\Lambda_{1}\,\varphi+\Lambda_{2}$, respectively, where $\Theta_{1}:=\frac{1}{2}\left(1-\frac{\overline{\rho}_{1}}{\overline{\rho}_{2}}\right)$, $\Theta_{2}:=\frac{1}{2}\left(1+\frac{\overline{\rho}_{1}}{\overline{\rho}_{2}}\right)$, $\Lambda_{1}:=\frac{1}{2}\left(1-\frac{\eta_{1}}{\eta_{2}}\right)$, and $\Lambda_{2}:=\frac{1}{2}\left(1+\frac{\eta_{1}}{\eta_{2}}\right)$. We recall that for $\varphi=-1$, we obtain the pure phase labeled “$1$”, while $\varphi=1$ refers instead to the pure phase “$2$”.

Let us now introduce the function spaces $\mathcal{V}=\left\{w\in H^{1}(\Omega)\ :\ w=0\ \text{ on }\Gamma_{D}\right\}$ and $\mathcal{H}:=H^{2}_{N}=\left\{w\in H^{2}(\Omega)\ :\partial_{\bm{\nu}}w=0\ \text{ on }\Gamma\right\}$, then, by suitably using integration by parts, the weak formulation of (4.1) reads as follows: Find, for all $t\in(0,T)$, $q\in\mathcal{V}$ and $\varphi\in\mathcal{H}$, such that

[TABLE]

hold for all $\psi\in\mathcal{V}$ and $\vartheta\in\mathcal{H}$ with $\varphi(t=0)=\varphi_{0}$ in $\Omega$.

4.1 Spatial approximation

For the spatial approximation of (4.2) we use NURBS-based Isogeometric Analysis (IGA) [16, 33]. Indeed, in (4.2) we look for a solution $\varphi\in\mathcal{H}\subset H^{2}(\Omega)$ for all $t\in(0,T)$, i.e., we need $H^{2}(\Omega)$-conformal finite dimensional test and trial function spaces, say $\mathcal{H}_{h}$, which are comprised of globally $C^{1}$-continuous basis functions. This requirement can be straightforwardly fulfilled by using B-splines (or NURBS) basis functions [44] of degree $p\geq 2$; we refer the interested reader to [9, 17, 31, 38, 50] for an overview of high order PDEs – including phase field problems – solved by means of NURBS-based IGA.

We introduce the bivariate B-splines basis $\left\{N_{A}(\mathrm{\mathbf{x}})\right\}_{A=1}^{n_{bf}}$ and we write the approximate pressure and order parameter as

[TABLE]

respectively, with the control variables $\left\{q_{A}(t)\right\}_{A=1}^{n_{bf}}$ and $\left\{\varphi_{A}(t)\right\}_{A=1}^{n_{bf}}$ being time-dependent. By introducing the B-splines space $\mathcal{N}_{h}=\text{span}\left\{N_{A},\ \ A=1,\ldots,n_{bf}\right\}$, we define the finite dimensional spaces $\mathcal{V}_{h}:=\mathcal{V}\bigcap\mathcal{N}_{h}$ and $\mathcal{H}_{h}:=\mathcal{H}\bigcap\mathcal{N}_{h}$. Then, the semi-discrete formulation of (4.2) reads as follows: Find, for all $t\in(0,T)$, $q_{h}\in\mathcal{V}_{h}$, $\varphi_{h}\in\mathcal{H}_{h}$, such that

[TABLE]

hold for all $\psi_{h}\in\mathcal{V}_{h}$ and $\vartheta_{h}\in\mathcal{H}_{h}$ with $\varphi_{h}(t=0)=\varphi_{0,h}$ in $\Omega$, where $\varphi_{0,h}$ is the $L^{2}(\Omega)$ projection of the initial condition $\varphi_{0}$ onto the space $\mathcal{N}_{h}$.

4.2 Time discretization

The time discretization of (4.3) is based on Backward Differentiation Formulas (BDF) [30, 45] with equal order temporal extrapolations based on Newton–Gregory backward polynomials [12, 46]. Using this semi-implicit formulation yields a fully discrete problem which can be solved in a computationally efficient and accurate manner; see for example [22] and [10] for the use of the BDF scheme together with NURBS-based IGA spatial approximations of the PDEs.

We partition the time interval $[0,T]$ into $N_{t}$ subintervals of equal size $\Delta t=\frac{T}{N_{t}}$ yielding the discrete time instances $t_{n}=n\,\Delta t$ for $n=0,\ldots,N_{t}$. Then, we denote with $q_{h}^{n}$ and $\varphi_{h}^{n}$ the approximations of the pressure $q_{h}$ and order parameter $\varphi_{h}$ at the time $t_{n}$. The approximation of $\partial_{t}\varphi_{h}$ in (4.3) by a $\sigma$-order BDF scheme is

[TABLE]

For example for $\sigma=1$, we have $\alpha_{\sigma}=1$ and $\varphi_{h}^{n,\text{BDF}\sigma}=\varphi_{h}^{n}$ for $n\geq 0$; instead, for $\sigma=2$, $\alpha_{\sigma}=\frac{3}{2}$ and $\varphi_{h}^{n,\text{BDF}\sigma}=2\varphi_{h}^{n}-\frac{1}{2}\varphi_{h}^{n-1}$ for $n\geq 1$. Then, replacing the derivative $\partial_{t}\varphi_{h}$ in (4.3b) with the $\sigma$-order BDF approximation, while the other time dependent terms are evaluated at the time instance $t_{n+1}$ (i.e., terms involving $q_{h}^{n+1}$ and $\varphi_{h}^{n+1}$), yields a nonlinear fully discrete problem at each time instance (for example, for $\sigma=1$ we have the backward Euler scheme).

In order to obtain a semi-implicit fully discrete problem, the nonlinear terms depending on $q_{h}^{n+1}$ and $\varphi_{h}^{n+1}$ are replaced by extrapolations of order $\sigma$ by means of the Newton–Gregory backward polynomials, say $q_{h}^{n+1,\sigma}$ and $\varphi_{h}^{n+1,\sigma}$, respectively. For example, for $\sigma=1$, these are $q_{h}^{n+1,\sigma}=q_{h}^{n}$ and $\varphi_{h}^{n+1,\sigma}=\varphi_{h}^{n}$ for $n\geq 0$; while, for $\sigma=2$, we have $q_{h}^{n+1,\sigma}=2\,q_{h}^{n}-q_{h}^{n-1}$ and $\varphi_{h}^{n+1,\sigma}=2\,\varphi_{h}^{n}-\varphi_{h}^{n-1}$ for $n\geq 1$. For a BDF scheme of order $\sigma$, the semi-implicit formulation of the fully discrete problem reads as follows: Find, for all $n\geq\sigma-1$, $q_{h}^{n+1}\in\mathcal{V}_{h}$, $\varphi_{h}^{n+1}\in\mathcal{H}_{h}$, such that

[TABLE]

hold for all $\psi_{h}\in\mathcal{V}_{h}$ and $\vartheta_{h}\in\mathcal{H}_{h}$ with $\varphi_{h}^{0}=\varphi_{0,h}$ in $\Omega$.

5 Numerical results

We present some numerical results for the Hele–Shaw–Cahn–Hilliard model for incompressible flows. Specifically, we solve two benchmark problems: the rising bubble test, for which a less dense fluid rises into a more dense one in presence of a gravitational field as e.g. in [34, 36], and the viscous fingering test, for which a less viscous fluid is injected into a more viscous one [48].

For both the tests, we use the (dimensionless) computational domain $\Omega=(0,0.5)\times(0,1)$. For the spatial approximation, we consider NURBS-based IGA with globally $C^{1}$-continuous B-splines basis functions of degree $p=2$ – as described in Section 4.1 – with $32,768$ equally-sized mesh elements, yielding the dimensionless mesh size $h=\frac{1}{256}$ and a total of $n_{bf}=33,540$ B-splines basis functions. For the time discretization, we use the semi-implicit formulation (4.4) with the BDF scheme of order $\sigma=2$; the time step size $\Delta t$ and $T$ are specified later for the two tests.

5.1 Test $1$: rising bubble

For this test, we set $\Gamma_{D}=\left\{(\mathrm{x},\mathrm{y})\in\Gamma\ :\ \mathrm{y}=0\right\}$ and $\Gamma_{N}\equiv\Gamma\backslash\Gamma_{D}$ by referring to the boundary conditions in (4.1) with $f_{N}=0$. Then, we set $\rho_{2}=1$, $\eta_{1}=\eta_{2}=0.1$, $\sigma=10^{-5}$, $g=9.80665$, the characteristic length $L=1$, $\delta=0.1$, $V=0.1$, and $\varepsilon=2\,h=\frac{1}{128}$. In this manner, we have $\mathrm{Ca}=10^{6}$ and $\mathrm{Bo}=0.980665$. For the time discretization, we set $\Delta t=2.5\cdot 10^{-5}$ for $T=4.0\cdot 10^{-1}$. As initial condition for the order parameter, we choose $\varphi_{0}(\mathrm{x},\mathrm{y})=\frac{1}{2}\zeta_{+}(\mathrm{x},\mathrm{y})\zeta_{-}(\mathrm{x},\mathrm{y})-1$, with

[TABLE]

which yields the set-up in Fig. 1 (top-left), where the blue color is associated to $\varphi=-1$ – the pure phase labeled “1” corresponding to the “heavy” fluid – while the red color is associated to $\varphi=+1$ – the pure phase “2” corresponding to the “light” fluid.

We start by considering the case $\rho_{1}=5$, for which $\Theta_{1}=-2$ and $\Theta_{2}=3$. We report in Figs. 1, and  2 the time evolution of the order parameter, which highlights the formation and rising of the bubble of light fluid, including topological changes. Correspondingly, we report in Figs. 3, and  4 the evolution of the computed velocity field $\bm{v}(q_{h}^{n+1},\varphi_{h}^{n+1})$; as we can observe, relatively high magnitudes of the velocity occur at pinch-off and when the curvature of the interface is significant.

We also consider the case where the density of the heavier fluid is larger, say $\rho_{1}=20$ (for which $\Theta_{1}=-\frac{19}{2}$ and $\Theta_{2}=\frac{21}{2}$) yielding the result highlighted in Fig. 5 with the velocity field.

5.2 Test $2$: viscous fingering

We set $\Gamma_{N}\equiv\Gamma$ with $\Gamma_{D}=\emptyset$ and, in order to enforce the injection of the fluid into the domain, $f_{N}=-V$ on $\Gamma_{N,b}=\left\{(\mathrm{x},\mathrm{y})\in\Gamma\ :\ \mathrm{y}=0\right\}$, $f_{N}=V$ on $\Gamma_{N,t}=\left\{(\mathrm{x},\mathrm{y})\in\Gamma\ :\ \mathrm{y}=1\right\}$, and $f_{N}=0$ on $\Gamma_{N}\backslash\left(\Gamma_{N,b}\bigcup\Gamma_{N,t}\right)$, for some injection velocity $V>0$. In this case, in order to obtain a well-posed problem, we prescribe the values of the control coefficients of the pressure field (approximated by IGA) $q_{32,319}=q_{32,320}=0$ for all $n\geq 1$. Then, we choose $\rho_{1}=\rho_{2}=1$, $\eta_{2}=1$, $g=0$, $L=1$, $\delta=0.1$, $V=50$, and $\varepsilon=2\,h=\frac{1}{128}$, for which $\mathrm{Bo}=0$. For the time discretization, we set $\Delta t=2.5\cdot 10^{-6}$ for $T=10^{-3}$. The initial condition is

[TABLE]

which yields the set-up in Fig. 6 (top-left); we recall that the blue color is associated to $\varphi=-1$ – the phase “1” indicating the more viscous fluid – while the red color is associated to $\varphi=+1$ – the phase “2” indicating the less viscous fluid.

We set $\eta_{1}=50$ and $\sigma=10^{-5}$, for which $\mathrm{Ca}=5.0\cdot 10^{8}$. The time evolution of the computed order parameter and velocity are reported in Figs. 6 and 7, respectively, which highlight the insurgency of the viscous fingering phenomenon.

In Fig. 8 we compare the order parameters at different time instances obtained for the values of the viscosity $\eta_{1}=10$, $20$, and $50$ for $\sigma=10^{-5}$, thus yielding $\mathrm{Ca}=10^{8}$, $2.0\cdot 10^{8}$, and $5.0\cdot 10^{8}$, respectively. We remark that the more viscous the fluid “1”, the longer the fingers.

Finally, in Fig. 9 we show the order parameter at different time instances obtained for the viscosity $\eta_{1}=50$ and values of the surface tension $\sigma=10^{4}$, $10$, and $10^{-5}$ thus yielding $\mathrm{Ca}=0.5$, $5.0\cdot 10^{2}$, and $5.0\cdot 10^{8}$, respectively. We observe that, the smaller the surface tension, the longer the fingers.