A decoupled finite element method for a modified Cahn-Hilliard-Hele-Shaw system

1.   Introduction

The Cahn-Hilliard-Hele-Shaw system is a very important mathematical model which describes the motion of a viscous incompressible fluid between two closely spaced parallel plates and can be viewed as the simplification of the Cahn-Hilliard-Navier-Stokes system ^[1,2,3]. The model are widely applied in different fields, such as simulations of nonlinear tumor growth and neovascularization ^[4,5,6,7], spinodal decomposition in a Hele-Shaw cell ^[8], and two-phase flow in porous medium ^[9,10], etc.

The Cahn-Hilliard-Hele-Shaw system is a gradient system coupled with fluid motion, which is difficult to solve because of its complex form. For this model, purely explicit methods are limited by strict time step constraints for stability, and completely implicit numerical methods must contend with potentially large systems of nonlinear algebraic equations ^[11]. There have been many effective numerical schemes for the Cahn-Hilliard-Hele-Shaw system. Guo et al. proposed a semi-implicit time integration scheme based on convex splitting technique, and proved the unconditional stability of the fully discrete scheme of the Cahn-Hilliard-Hele-Shaw system ^[12]. S.M. Wise put forward an unconditionally stable finite difference scheme for the Cahn-Hilliard-Hele-Shaw ^[13]. Chen et al. established a finite difference simulation of Gagliardo-Nirenberg-type inequalities to analyze stability and convergence ^[14]. Liu et al. developed a mixed finite element numerical scheme for the Cahn-Hilliard-Hele-Shaw system and proved its unconditional stability ^[15]. Guo carried out a numerical analysis for the Cahn-Hilliard-Hele-Shaw system with variable mobility and logarithmic Flory-Huggins potential ^[16]. The above mentioned works are numerical methods to solve the Cahn-Hilliard-Hele-Shaw system. However, there are few researches on the modified Cahn-Hilliard-Hele-Shaw system.

The modified Cahn-Hilliard equation (also named Cahn-Hilliard-Oono equation) used to suppress phase coarsening in ^[17] is as follows

where $\overline{\phi}_{0}: = \dfrac{1}{|\Omega|}\int_{\Omega}\phi_{0}(x)dx$. More works on the modified Cahn-Hilliard equation can be found in ^{[18,19,20,21]}. For the modified Cahn-Hilliard equation, when $\theta$ = 0, the equation becomes the classical Cahn-Hilliard equation^[22,23,24]. When the modified Cahn-Hilliard is coupled with the Darcy equation, the modified Cahn-Hilliard-Hele-Shaw equation can be obtained. Jia et al. introduced a novel finite element method for the modified Cahn-Hilliard-Hele-Shaw system ^[25], in which the time discretization was based on the convex splitting of the energy functional in the modified Cahn-Hilliard equation. Of course, the above numerical methods are directly solved based on the coupling equation, and the solving process is complicated. To solve this kind of problem, many decoupled methods have been proposed to solve the Cahn-Hilliard-Hele-Shaw system in recent years. Han ^[26] presented a decoupled unconditionally stable numerical scheme for the Cahn-Hilliard-Hele-Shaw system with variable viscosity, in which the operator-splitting strategy and the pressure-stabilization technique were used to completely decouple the nonlinear Cahn-Hilliard equation from pressure. Similar strategies were also adopted in ^[27]. Then, Gao ^[28] studied the fully decoupled numerical scheme of the Cahn-Hilliard-Hele-Shaw model, in which the scalar auxiliary variable method was used to deal with the nonlinear term in the free energy. Similarly, decoupled schemes are also effectively used in other systems and models recently. Zhao et al. ^[29] developed an energy-stable scheme for a binary hydrodynamic phase field model of mixtures of nematic liquid crystals and viscous fluids. A second-order decoupled energy-stable schemes for Cahn-Hilliard-Navier-Stokes equations was suggested in ^[30]. For thermodynamically consistent models, Zhao ^[31] investigated a general numerical framework for designing linear, energy stable, and decoupled numerical algorithms. However, to the best of our knowledge, there are few researches on decoupling methods of the modified Cahn-Hilliard-Hele-Shaw system, it will be the purpose of our paper.

Based on Eqs (1.1)-(1.3), the modified Cahn-Hilliard-Hele-Shaw system with double well potential is given by

where $\Omega\in \mathbf{R}^{d} (d = 2, 3)$. $\phi$ is the concentration field, $\mathbf{u}$ is the advective velocity, $\varepsilon > 0$ is the constant to measure the thickness of the transition layer between the two phases, and $\mu$ is the chemical potential. $f(\phi)$ is the derivative of the double well potential $F(\phi)$, $\xi$ is an auxiliary variable. $p$ and $\gamma$ represent the pressure and the dimensionless surface tension parameter, respectively. $\mathbf{n}$ is the unit outer normal of the boundary $\partial\Omega$. when $\theta$ = 0, the equation becomes the classical Cahn-Hilliard-Hele-Shaws equation. With regard to the double well potential corresponding to $f(\phi)$ in Eq (1.2), the following $\check{F}(\phi)$ can be taken^[32,33,34]

Correspondingly, the derivatives of $\check{F}(\phi)$ can be split as follows

$F(\phi)$ and $f(\phi)$ are replaced by $\check{F}(\phi)$ and its derivative $\check{f}(\phi)$, which are still recorded as $F$ and $f$ for simplicity. Typically, the free energy functional of a modified Cahn-Hilliard-Hele-Shaw system with double well potential is given by

In this paper, a decoupled finite element scheme for a modified Cahn-Hilliard-Hele-Shaw system with double well potential is proposed. The temporal discretization is based on the convex splitting of the energy functional in the modified Cahn-Hilliard equation, and the spacial discretization is carried out by the mixed finite element method. The computation of the velocity $\mathbf{u}$ is separated from the computation of the pressure $p$ by using an operator-splitting strategy. A Possion equation is solved to update the pressure at each time step. We prove that the the proposed scheme is unconditionally stable in energy, and the error analyses are obtained. Finally, the numerical results verify the theoretical analysis. The rest of this article is structured as follows.The finite element discrete scheme of the Cahn-Hilliard-Hele-Shaw system combing with the convex splitting is given in Section 2; The theoretical preparations and stability of the proposed numerical scheme are proved in Section 3; The error analyses of the proposed scheme are addressed in Section 4; Some numerical examples are given to verify the previous theory in Section 5, and the conclusion is given in Section 6.

2.   Mathematical model

2.1.   Semi-discrete scheme

Let $L^{2}(\Omega)$ is a space of square integrable function and $H^{k}(\Omega), H_0^{k}(\Omega)$ denote the usual Sobolev spaces. $L^2{(\Omega)}$ inner product and its norm are denoted by $(u, v) = \int_{\Omega}u(x)v(x)dx$, $\|\phi\| = \|\phi\|_{L^{2}(\Omega)} = \sqrt{(\phi, \phi)}$. The weak formulation of the modified Cahn-Hilliard-Hele-Shaw system with double well potential can be written as

where,

Let $N$ be a positive integer and $0 = t_{0} < t_{1} < \cdot\cdot\cdot < t_{N} = T$ be a uniform partition of $[0, T]$, where $t_{i} = i\tau$, $i = 0, 1, \cdot\cdot\cdot, N-1$, $\tau = \dfrac{T}{N}$.

The semi-discrete scheme of the modified Cahn-Hilliard-Hele-Shaw system with double well potential is as follows. For $n\geq0$, find $\{\phi^{n+1}, \mu^{n+1}, \xi^{n+1}, p^{n+1}\}$ such that

where the velocity is given by

Combing with the idea of the literatures ^[26,35], the computation of the modified Cahn-Hilliard equations (2.2)-(2.4) are decoupled from Eq (2.5) after substituting $\mathbf{u}^{n+1}$ into Eq (2.2), since the pressure is explicit in Eq (2.6). The velocity $\mathbf{u}^{n+1}$ is regarded as an intermediate velocity by using the incremental projection method similar to the Navier-Stokes equation. The real velocity $\tilde{\mathbf{u}}^{n+1}$ is obtained from the intermediate velocity and satisfies

Then Eq (2.6) and Eq (2.7) are added together to obtain the original Eq (1.7). If the divergence operator is applied to both side of Eq (2.7), the real velocity $\tilde{\mathbf{u}}^{n+1}$ will vanished. We have

2.2.   Fully discrete scheme

Let $T_{h} = \{K\}$ be a regular partition of the domain $\Omega$ that is divided into triangles with the size $h = \max\limits_{0\leq i\leq N}h_{i}$. $S_{h}$ is a piecewise polynomial space, which is defined as

where $P_{k}(x, y)$ is a polynomial of degree at most r.

Let us denote

The corresponding fully discrete scheme have the following expression, find $\{\phi^{n+1}_{h}, \mu^{n+1}_{h}, \xi^{n+1}_{h}, p^{n+1}_{h}\}\in S_{h}\times S_{h}\times\hat{S}_{h}\times\hat{S}_{h}$, such that

where the velocity is given by

3.   The analysis of stability

Definition 3.1. ^[36] The Ritz projection operator $R_{h}(\Omega)$: $\phi\in H^{1}(\Omega)\rightarrow S_{h}$ satisfies

and have the following estimates,

Definition 3.2. ^[36] Define the operator $T_{h}:\hat{H}^{-1}\rightarrow \hat{H}^{1}$ through the following variational problems, given $\zeta\in\hat{H}^{-1}$, find $T_{h}(\zeta)\in\hat{H}^{1}$ such that

Lemma 3.1. ^[12,15] Let $\zeta, \varphi\in\hat{H}^{-1}$ and set

where $(\cdot, \cdot)_{-1, h}$ defines an inner product on the $\hat{H}^{-1}$ and its corresponding $H^{-1}$ norm is written as

Consequently, for $\forall\chi\in\hat{H}^{1}, \zeta\in\hat{H}^{-1}$,

Furthermore, the following Poincar$\acute{e}$ inequalities holds,

Definition 3.3. ^[12,15] Define $\mathbf{W}: = \{\mathbf{u}\in\mathbf{L}^{2}(\Omega)|(\mathbf{u}, \nabla q), \forall q\in H^{1}(\Omega)\}$. The projection operator $\mathcal{P}:\mathbf{w}\in \mathbf{L}^{2}(\Omega)\rightarrow \mathbf{W}$ is defined as

where $p\in\dot{H}^{1}: = \{\phi\in H^{1}(\Omega)|(\phi, 1) = 0\}$ is the unique solution to

Lemma 3.2. ^[12,15] Projection operator $\mathcal{P}$ is linear and satisfies the following properties

and

Definition 3.4. ^[12,15] Define $\mathbf{W}_{h}: = \{\mathbf{u}_{h}\in\mathbf{L}^{2}(\Omega)|(\mathbf{u}_{h}, \nabla q_{h}) = 0, \forall q_{h}\in S_{h}\}$. The projection operator $\mathcal{P}_{h}:\mathbf{w}\in\mathbf{L}^{2}(\Omega)\rightarrow \mathbf{W}_{h}$ is defined as

where $p_{h}\in\hat{S_{h}}$ is the unique solution to

Lemma 3.3. ^[14,15] Projection operator $\mathcal{P}_{h}$ is linear and satisfies the following properties

and

Lemma 3.4. ^[12,15] Suppose that $\mathbf{w}\in\mathbf{H}^{q}(\Omega)$ with the compatible boundary condition $\mathbf{w}\cdot\mathbf{n} = 0$ on $\partial\Omega$ and $q\in\mathbf{H}^{q+1}(\Omega)$, then

Theorem 3.1. Let $\{\phi^{n+1}_{h}, \mu^{n+1}_{h}, p^{n+1}_{h}, \xi^{n+1}_{h}\}$ be the unique solution of Eqs (2.9-2.12). Define

Then for any $h, \tau, \varepsilon > 0, n\geq0$, scheme (2.9)-(2.12) satisfies the following property,

Proof. Taking $v_{h} = \tau\mu^{n+1}_{h}$ in Eq (2.9), one has

In Eq (2.10), $f_{1}(\phi^{n+1}_{h}) = 2\phi^{n+1}_{h}$, $f_{2}(\phi^{n}_{h}) = (\phi^{n}_{h})^{3}-3\phi^{n}_{h}$. For $f_{2}(\phi^{n}_{h})$, through Taylor expansion

where $\eta$ is a number between $\phi^{n}_{h}$ and $\phi^{n+1}_{h}$, we have

Then, choosing $w_{h} = -(\phi^{n+1}_{h}-\phi^{n}_{h})$ and using the fact that $(a, a-b) = \dfrac{1}{2}[a^{2}-b^{2}+(a-b)^{2}]$ give

Replacing $\psi_{h}$ by $-T_{h}(\phi^{n+1}_{h}-\phi^{n}_{h})$ in Eq (2.11). By Eq (3.1) in definition 3.1, Eq (3.5) in lemma 3.1 and $(a, a-b) = \dfrac{1}{2}[a^{2}-b^{2}+(a-b)^{2}]$, one obtains

Next, we take inner product of Eq (2.13) with $\dfrac{\tau}{\gamma}\mathbf{u}^{n+1}_{h}$ to get

Now, taking $q_{h} = \dfrac{\tau}{\gamma}p^{n}_{h}$ and using the fact that $(a-b, 2b) = a^{2}-b^{2}-(a-b)^{2}$ in Eq (2.12), we arrived at

To deal with the $\dfrac{\tau}{2\gamma}\|\nabla p^{n+1}_{h}-\nabla p^{n}_{h}\|^{2}$ in Eq (3.24), replacing $q_{h}$ with $(p^{n+1}_{h}-p^{n}_{h})$ in Eq (2.12) and using Cauchy-Schwarz inequalities, the following estimation can be obtained

Combining Eqs (3.23)-(3.25), it can be written as

Summing Eqs (3.20)-(3.26), one concludes that

Since $f_{2}'(\phi) = 3(\phi^{2}-1)\leq0$, $\phi\in[-1, 1]$, there is $\dfrac{f_{2}'(\eta)}{2}\|\phi^{n+1}_{h}-\phi^{n}_{h}\|^{2}\leq0$ by Taylor expansion. Therefore,

The proof is completed.

Corollary 3.1. Suppose that $\Xi({\phi^{0}_{h}})\leq C_{0}$, there is a constant $C > 0$ independent of $\tau$ and $h$, such that the following estimates hold for any $\tau, h > 0$,

Proof. Summing the Eq (3.19) from $i = 0\; \rm{to}\; N$, we get

The proof is completed.

4.   Error estimates

In this section, we assume that the weak solution $\{\phi, \mu, \xi, p\}$ satisfies the following regularity

For the convenience of subsequent analysis, we introduce some notations,

Lemma 4.1. ^[36] Suppose the $\{\phi, \mu, \xi, p\}$ is the solution to Eq (2.1), the following estimate holds

Theorem 4.1. Suppose the solutions of the initial problem Eq (2.1) and the fully discrete scheme Eqs (2.9)-(2.12) are $\{\phi, \mu, \xi, p\}$ and $\{\phi^{n+1}_{h}, \mu^{n+1}_{h}, \xi^{n+1}_{h}, p^{n+1}_{h}\}$, respectively. Then for any $h\;, \tau\; > \; 0$, the following estimate holds

Proof. Subtracting Eqs (2.9)-(2.12) from Eq (2.1) at $t = n+1$, one has

We choose $v_{h} = \hat{e}^{n+1}_{\mu}$ in Eq (4.3), $w_{h} = -\delta_{\tau}\hat{e}^{n+1}_{\phi}$ in Eq (4.4), $\psi_{h} = -T_{h}(\delta_{\tau}\hat{e}^{n+1}_{\phi})$ in Eq (4.5), $q_{h} = \varepsilon\hat{e}^{n+1}_{p}$ in Eq (4.6) and sum them to get

where we denote

Next, we estimate $M_{i}$. According to the poincar$\acute{e}$ inequality, the Cauchy-Schwarz inequality, the Young inequality and lemma 4.1, one obtains

Using Eq (3.7) in lemma 3.1, the Young inequality and Eq (3.3) in definition 3.1, we have

Similarly, according to lemma 3.1, the Schwarz inequality, the Young inequality, and Eq (3.3) in definition 3.1, we can estimate $M_{3}$ as follows,

As for $M_{4}$, there is $f_{1}(\phi^{n+1})-f_{1}(\phi^{n+1}_{h}) = 2(\phi^{n+1}-\phi^{n+1}_{h})$ for $f_{1}(\phi) = 2\phi$, and $f_{2}(\phi^{n+1})-f_{2}(\phi^{n}_{h})\leq C(\phi^{n+1}-\phi^{n}_{h})$ for $f_{2}(\phi) = \phi^{3}-3\phi$. Then, according to lemma 3.1, the Young inequality, definition 3.1 and Taylor extension $\|\nabla\tau\delta_{\tau}\phi(t)\|^{2}\leq C\tau^{2}$, the following inequality is established

To deal with $M_{5}$, we denote $b(\phi, \mathbf{u}, v): = (\phi\mathbf{u}, \nabla v)$. Referring to the method in ^[15], $M_{5}$ can be analyzed as

where $D: = \|\phi^{n}_{h}\|^{4}_{L^{\infty}}+1\leq C$. Therefore,

According to the definition 3.3, definition 3.4, lemma 3.4, Taylor expansion $\|\nabla\tau\delta_{\tau}p^{n+1}\|^{2}\leq C\tau^{2}$, the Cauchy-Schwarz inequality and the Young inequality, the following error estimation formulation holds

Combining Eqs (4.7)-(4.15) gives

For $\|\delta_{\tau}\hat{e}^{n+1}_{\phi}\|^{2}_{-1, h}$, taking $\alpha v_{h} = T_{h}(\delta_{\tau}\hat{e}^{n+1}_{\phi})$ in Eq (4.3) and using a similar idea as $M_{5}$, the following inequality can be obtained,

Therefore, it follows that

Then, combining Eq (4.16) with Eq (4.18) and multiplying by $2\tau$, one has

Finally, we take the appropriate $\alpha(0 < \alpha\leq \dfrac{1}{2})$ and add the above estimates from $i=0$ to $n$. When $0 < \tau\leq\dfrac{\alpha\varepsilon^{2}}{8}$, according to the discrete Gronwall inequality, one concludes that

The proof is completed.

5.   Numerical analysis

In this part, some numerical examples are used to verify the correctness and validity of the theoretical analysis. Next, let us take the initial conditions $\phi_{0} = 0.24*cos(2\pi x)cos(2\pi y)+0.4*cos(\pi x)cos(3\pi y)$, and the domain of the calculation is $[0, 1]\times[0, 1]$.

5.1.   The spatial convergence order

For Tables 1 and 2, the parameters are chosen as follows, $\tau$ = 0.01, $T = 0.1$, $\varepsilon$ = 0.14 and mesh steps $h = \frac{1}{16}, \frac{1}{32}, \frac{1}{64}, \frac{1}{128}$. The spatial convergence orders of relative error $\|\hat{e}_{\phi}\|_{H^{1}}$ are close to 1, which is consistent with the convergence order obtained from theoretical analysis. Moreover, different $\theta$ and $\gamma$ have little effect on the corresponding convergence order.

5.2.   The temporal convergence order

For Tables 3 and 4, the parameters are chosen as follows, $\varepsilon = 0.01$, $T = 0.1$, $h = \tau =$0.0625, 0.03125, 0.015625. The temporal convergence orders of relative error $\|\hat{e}_{\phi}\|_{H^{1}}$ are close to 1, which is consistent with the convergence order obtained from theoretical analysis.

5.3.   Energy dissipation

Let us test the energy dissipation of our proposed scheme. The energy functional Eq (1.13) of the modified Cahn-Hilliard-Hele-Shaw system Eqs (1.4)-(1.10) can be discreteized as

Correspondingly, the modified energy of the fully discrete scheme Eqs (2.9)-(2.13) is defined as

For the test, the parameters are chosen as follows: $T = 5$, $\tau = 0.001$, $h = \dfrac{1}{64}$, $\varepsilon = 0.4$, $\gamma = 0.5$. In Figure 1, we can see that the energy functional is non-increasing for $\theta$ = 0, 0.1, 1.

5.4.   Spinodal Decomposition

In this part, we present the phase separation dynamics that is called spinodal decomposition in the modified Cahn-Hilliard-Hele-Shaw system. In the simulation, the computational domain is chosen as $[0, 1]\times[0, 1]$, the parameters are chosen as follows: $\varepsilon = 0.05$, $\gamma = 0.45$, $\tau = 0.0001$. Then, let us take the initial condition

where $rand()\in [0, 1]$. The process of coarsening is shown in the following figures. From figures 2-19, we can see that the contours of $\phi$ are gradually coarsened over time. However, the profiles obtained by different $\theta$ are similar at the same time $T$. From left to right, the coarsening processes of $\theta = 50, 200$ are not obvious compared with the coarsening processes of $\theta = 0$. We know the bigger $\theta$ can suppress the coarsening process.

6.   Summary

In this paper, a decoupled scheme of the modified Cahn-Hilliard-Hele-Shaw system is studied. In our scheme, the velocity and pressure are decoupled, and a Possion equation is solved to update the pressure at each time step. Unconditional stability of the scheme in energy is proved. The convergence analysis are addressed in the frame of finite element method. Furthermore, the theoretical part is verified by several numerical examples. The results show that the numerical examples are consistent with the results of the theoretical part.

Acknowledgments

The work is supported by the the Provincial Natural Science Foundation of Shanxi (No. 201901D111123) and Key Research and Development (R & D) Projects of Shanxi Province (No. 201903D121038).

Conflict of interest

The authors declare no conflicts of interest in this paper.