A fully-decoupled energy stable scheme for the phase-field model of non-Newtonian two-phase flows

1.   Introduction

In this paper, we focus on the numerical methods for a non-Newtonian fluid dynamical system coupled to the Cahn-Hilliard equation ^[1], which is governed by the following:

where $T > 0$, $\Omega\subset\mathbb{R}^{2}$ is a bounded domain, $f(\phi) = F'(\phi) = \phi^{3}-\phi$, $\textbf{u}$ represents the velocity, $p$ is the pressure, $\mu$ denotes the chemical potential of the binary mixture depending on the phase function $\phi$, and $M > 0$ is the mobility. Here, we consider the extra-stress tensor $\sigma(\textbf{u})$ that obeys the following Carreau's law ^[2]:

where $D(\textbf{u}) = (\nabla\textbf{u}+\nabla \textbf{u}^{T})/2$ represents the strain and the parameters $\nu_{1} > 0$, $\nu_{2} > 0$, and $\kappa > 0$. The fluids are characterized to have a shear-thinning behavior for $1 < r < 2$ and a shear-thickening behavior for $r > 2$.

Throughout the paper, we close the aforementioned system (1.1) with the following initial and boundary conditions:

where $x\in\Omega$ is the spatial variables, and $\textbf{n}$ denotes the unit outward normal to the boundary $\partial\Omega$. Phase-field models have been widely studied by many authors. Bian-Zhang ^[3] established the weak solutions for a degenerate phase-field model via a Galerkin approximation. In view of a physical point, the background of the phase-field model is well studied, see ^[4,5].

For the case $r = 2$ (or $\kappa = 0$) in (1.2), the system (1.1) is reduced to the phase-field models related to the Newtonian-fluid flows. As far as the author knows, for the classical phase field model of two-phase flows the only available energy-stable fully-decoupled scheme was developed in ^[6,7,8], in which a full decoupling structure was achieved by adding a stabilization term to the explicit advection velocity term. However, this decoupling type scheme is only first-order accurate in time and it requires more calculations at each time step. A second-order fully-decoupled scheme was developed in ^[9,10], in which the advection and surface tension terms were discretized through explicit and implicit combination methods, which can inevitably lead to either expensive nonlinear fully-coupled schemes or relatively fast linear fully-coupled schemes. Therefore, the focus of this paper to develop a scheme for the given system that maintains an unconditionally energy stable scheme and a second-order time accuracy, while also achieving full decoupling calculations. The scalar auxiliary variable $V(t)$ is used to develop a stabilized Scalar Auxiliary Variable (SAV) approach, which has been studied by many authors in many equations such as the Allen-Cahn equation, the phase field crystal equation, and the magneto-hydrodynamic equations, see ^[11,12,13]. The reason why $V(t)$ can be used as a stabilizer is that we use the explicit format to deal with $(u, \phi)$, thus producing an error. In order to overcome this disadvantage, the auxiliary variable is introduced, see Lemma 2.2 for the details. Previous studies have considered a fully decoupled scheme with a stabilization term, which was essential for the phase-field model of non-Newtonian two-phase flows ^[14,15]. The numerical approaches to solve the phase-field model of non-Newtonian two-phase flows contain the classical finite difference, and the finite element method, and we will provide the concrete scheme. Inspired by the works of SAV methods ^[16,17] and the invariant energy quadratization (IEQ) methods ^[18,19,20], we develop a novel stabilized-SAV method to solve the coupled system (1.1). The introduction of the new auxiliary variable can provide a sufficient skill to eliminate all the troublesome nonlinear terms, which can be independently solved to achieve the fully-decoupled structure.

The rest of the paper is organized as follows. In the next section, we derive an equivalent system of the phase-field model of non-Newtonian two-phase flows with the SAV approach. In Section 3, we propose an efficient, fully-decoupled and second-order time-accurate scheme to solve the given system. In Section 4, a numerical simulation is performed to confirm the effectiveness of the proposed method.

2.   Equivalent system with the SAV approach

To begin with, in order to study the interaction of the phase-fluid interface problem, the total energy of the dynamic system is defined as a sum of the kinetic energy $E_{kin}(\textbf{u})$ and the interfacial free energy $E_{f}(\phi)$:

Then, the system of phase-field models within the context of non-Newtonian fluids satisfies the following basic energy law. We say $(\phi, \mu, \textbf{u}, p)$ is a solution of system (1.1)–(1.3) if $(\phi, \mu, \textbf{u}, p)$ satisfies (1.1)–(1.3) in a weak sense.

Lemma 2.1. Let $(\phi, \mu, \textbf{u}, p)$ be the solution to the system (1.1)–(1.3). Then, the following energy law holds:

Moreover, suppose that the initial values $\textbf{u}_{0}$ and $\phi_{0}$ satisfy $E_{tot}(\textbf{u}_{0}, \phi_{0}) < \infty$; then, we have the following energy dissipation law:

Proof. From (1.1)–(1.3), making use of integration by parts and using the fact that $\nabla\cdot\textbf{u} = 0$, then we have the following:

where $\delta E_{tot}/\delta \phi$ denotes the variational derivative. Therefore, the proof of the desired energy law (2.2) is completed. Besides, the estimate (2.3) can be easily obtained by integrating (2.2) from $0$ to $t$.

Remark 2.1. When deriving (2.2), it is worth noting that three nonlinear terms do not present any difficulties because they are all canceled out. More precisely, we have the following two identities:

where the first one is due to integration by parts, the second one is due to the divergence-free condition of $\textbf{u}$. The two identities mean that the nonlinear terms do not contribute either to the total free energy or the energy diffusion rate. In the next section, when developing a decoupling type scheme, we will take advantage of this often-overlooked feature of "zero-energy-contribution".

Now, we develop a stabilized-SAV approach to design a linear scheme as follows. A scalar auxiliary variable is defined by

with an energy function

where $\beta$ is an appropriate stabilized parameter and $C_{0} > 0$ is a constant to guarantee that the radicand is positive. It is easy to check that

where we denote $H(\phi) = \frac{\delta E}{\delta \phi}$. Therefore, we can rewrite the energy of (2.1) as the following formulation:

Based on (2.4), (2.5), and (2.6), due to $\nabla\cdot\textbf{u} = 0$, then the system (1.1) can be reformulated as follows:

with the following initial conditions

The boundary conditions in (1.3) still hold due to the new variable $V$, which is independent of the spatial variable.

Remark 2.2. We emphasize that the newly transformed SAV system (2.1) is exactly equivalent to the original system (2.3), since (2.6) can be easily obtained by integrating $(2.10)_{5}$ with respect to time. Thus, the new energy law (2.9) for the transformed system is the same as the energy law (2.1) for the original system. We state it as the following lemma.

Lemma 2.2. Let $(\phi, \mu, \textbf{u}, p)$ be the solution to the system (2.10). Then, the solution satisfies the following energy dissipation law:

namely,

Proof. Taking the $L^{2}$ inner product of $(2.10)_{1}$ with $\mu$, $(2.10)_{2}$ with $-\phi_{t}$, $(2.10)_{3}$ with $\textbf{u}$, $(2.10)_{4}$ with $p$, multiplying $(2.10)_{5}$ by $2V(t)$, and adding the resulting equations together, we directly obtain (2.11) and (2.12).

The next proposed scheme is required to formally obey the new energy dissipation law (2.9) in the discrete sense.

3.   Numerical scheme

In this section, we provide a numerical scheme. The method used here has been widely studied, see ^{[21,22,23,24]} for the details. We assume a uniform partition of the time interval $[0, T]$ as $t_{k} = k \Delta t, k = 0, 1, \cdots, K$ with a time step $\Delta t = T/K$. For a smooth function $\varphi$, the approximation of $\varphi$ at time $t_{k}$ is denoted as $\varphi^{k} = \varphi(t_{k})$. Moreover, for convenience, the following notations will be used in the sequence:

Now, it is ready to build up a numerical scheme to discretize the new system (2.8) by using the second-order Crank-Nicolson (CN) scheme. It reads as follows.

Step 1: Find $(\phi^{k+1}, \mu^{k+1}, \widetilde{\textbf{u}}^{k+1})$ such that

Step 2: Find $(\textbf{u}^{k+1}, p^{k+1})$ such that

Step 3: Update $V^{k+1}$ as follows

The boundary conditions of the schemes (3.1)–(3.6) are given by the following:

Now, the energy stability of the developed schemes (3.1)–(3.6) is shown as follows.

Theorem 3.1. The scheme (3.1)–(3.6) is unconditionally energy stable, that is, the following discrete energy law holds:

where $\mathcal{E}^{k+1}$ is the discrete version of the total energy given by the following:

Proof. We take the $L^{2}$-inner product of (3.1) and (3.2) with $\Delta t\mu^{k+1/2}$ and $-(\phi^{k+1}-\phi^{k})$; respectively, by exploiting integration of parts and the boundary conditions (3.7) and adding the resulting equations together, we obtain the following:

By taking the $L^{2}$-inner product of (3.3) and (3.4) with $\Delta t\widetilde{\textbf{u}}^{k+\frac{1}{2}}$, respectively, adding the desired results and using (3.5), one arrives at the following:

By taking the inner product of (3.4) with $\frac{\Delta t^{2}}{4}\nabla(p^{k+1}+p^{k})$ in the $L^{2}$ space and using (3.5) again, we have

We multiply (3.6) with $2\delta tV^{k+1/2}$ to obtain the following:

Summing the above equations (3.9)–(3.12) together leads to the proof of the desired result (3.8). The proof is completely finished.

4.   Numerical result

In this section, we aim to investigate the performance of our proposed method for the non-Newtonian two-phase flow. The solution of the systems (1.1)–(1.3) shall be a sought subject to the following initial conditions:

where $\tanh(x) = (e^{x}-e^{-x})/(e^{x}+e^{-x})$ and $\varepsilon = 0.02$. The boundary conditions are $\phi|_{\partial \Omega} = \textbf{u}_{\partial \Omega} = 0$.

We set the computational domain $\Omega = [0, 2]^{2}$ and use the finite element method to discretize the space variables (using $P_2$ for $\phi$ and $P_2-P_1$ for $(\textbf{u}, p)$); the parameters are given as $M = 1, \nu_{1} = 0.03, \nu_{2} = 0.01$, and $\kappa = 1$. Figure 1 shows the $L^{2}$ errors of the phase-field variable, the fluid velocity, the pressure between the numerical solution, and the "exact solution" (which is computed by a very small time step $\delta t = 1/2^{14}$) at time $t = 0.8$ with a different $r$. We can see that the obtained temporal convergence rate is $\mathcal{O}(\Delta t^{2})$.

In the next simulation, we take the temporal step as $\Delta t = 0.01$ in the computations. In Figure 2(a), the evolution of the total energy curves with different parameters $r$ is plotted, from which it is easy to see that the energy is decreasing in time, which is in a good agreement with the theoretical results (see Theorem 3.1). Figure 2(b) exhibits the pressure profile at a cross section $y = 1.0$ along the $x$-axis at time $t = 0.20$. Moreover, one can clearly see that the parameters $r$ can affect the profile of the pressure, which gives rise to different patterns of the phase field. Figure 3 displays the snapshots of the phase field variable $\phi$ and the flow velocity vector $\textbf{u}$ at the same time nodes, where the different patterns are observed under different values of $r$, that is, the parameter $r$ in non-Newtonian fluids can affect the spatial shapes of a substance during a surface phase transition. Subsequently, this gives us a means to design and perform optimal control experiments for the geometric shapes of materials in engineering practice.

5.   Conclusions

In this paper, a novel fully-decoupled, linear and second-order time accurate scheme to solve the phase-field model of non-Newtonian two-phase flows was introduced. The developed scheme was based on a stabilized SAV approach.

For the classical phase-field model of two-phase flows, the only available energy-stable fully-decoupled scheme was developed and a fully-decoupled structure was achieved by adding a stabilization term to the explicit advection velocity term. However, this decoupling type scheme was only first-order accurate in time and it required more calculations at each time step. For the second-order fully-decoupled scheme, the advection and surface tension terms were discretized through explicit and implicit combination method, which can inevitably lead to either expensive nonlinear fully-coupled schemes or relatively fast linear fully-coupled schemes. Therefore, the focus of this paper was developing a scheme for the given system that maintained a unconditionally energy stable scheme and a second-order time accuracy, while also achieving fully-decoupled calculations. We strictly proved the unconditional energy stability of the scheme and conducted numerical simulations to show the accuracy and stability of the proposed scheme. Moreover, we observed that the parameter $r$ in non-Newtonian fluids could affect spatial patterns during phase transitions, which directly enabled us to design and perform optimal control experiments in engineering processes.

Author contributions

Wei Li: Conceptualization, Investigation, Writing-original draft; Guangying Lv: Conceptualization, Validation, Writing-original draft, Writing-review and editing. All authors have read and approved the final version of the manuscript for publication.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work was supported by the NSF of China grants 12171247, 11771123, Jiangsu Provincial Double-Innovation Doctor Program JSSCBS20210466 and Qing Lan Project, and Postgraduate Research and Practice Innovation Program of Jiangsu Province (No. KYCX21 0932).

Conflict of interest

The authors declare that they have no conflict of interest.