On a novel fully-decoupled, linear and second-order accurate numerical scheme for the Cahn–Hilliard–Darcy system of two-phase Hele–Shaw flow

Abstract

We construct a novel fully-decoupled and second-order accurate time marching numerical scheme with unconditional energy stability for the Cahn–Hilliard–Darcy phase-field model of the two-phase Hele–Shaw flow, in which, the key idea to realize the full decoupling structure is to use the so-called “zero-energy-contribution” function and design a special ordinary differential equation to deal with the nonlinear coupling terms between the flow field and the phase-field variable. Compared with the existing decoupling type schemes, the scheme developed here is more effective, efficient and easy to implement. At each time step, one only needs to solve a few fully-decoupled linear equations only with constant coefficients. We also strictly prove the solvability and unconditional energy stability of the scheme, and implement various numerical simulations in 2D and 3D to show the efficiency and stability of the proposed scheme numerically.

Introduction

In this paper, we consider the numerical approximations of the Cahn–Hilliard–Darcy phase-field model of the two-phase Hele–Shaw flow system (CHD, for short). The term Hele–Shaw (or Hele–Shaw cells) is commonly used to describe the motions of viscous fluids which is limited by two parallel plates with a very small gap. For this type of flow, the fluid motion conforms to the mechanical principles in porous media, and the boundary conditions are defined by the pressure. Starting with the pioneering work in [1], the Cahn–Hilliard equation describing the two-phase flows is coupled with the Darcy equations to simulate the motion of the two-phase fluid flow confined in two flat plates, in which both fluids are assumed to be inertialess. Similar models had also been used to describe tumor growth, see [2], [3]. Although the CHD model has received extensive attention, as a complex nonlinear system, to develop easy-to-implement, accurate, and unconditional energy stable schemes for solving it is still very challenging. The main difficulty lies on how to discretize the following three terms, (i) the coupling of the velocity and phase-field function through the nonlinear advection and surface tension; (ii) the stiffness term related to the interfacial width due to the nonlinear double-well potential; and (iii) the coupling between the velocity and pressure through the divergence-free condition in the Darcy equation.

The difficulty referred in (ii) means that for a nonlinear term that causes stiffness, simple implicit or explicit type discretization may lead to strict stability constraints on the time step (cf. [4], [5]), resulting in practical inefficiency. We recall that many successful efforts have been made in the direction of designing appropriate discretization for nonlinear terms to obtain unconditionally energy stable schemes. Some well-known methods include, but are not limited to: the convex splitting [6], [7], [8], [9], [10], [11], the linear stabilization [5], [12], [13], the Invariant Energy Quadratization (IEQ) [14], [15], [16], [17] and its various version of the Scalar Auxiliary Variable (SAV) [18], [19], [20], nonlinear implicit derivative [21], nonlinear quadrature [22], [23], [24] methods, etc. Therefore, it can be considered that the difficulty (ii) has been solved well. The difficulty mentioned in (iii) can also be easily solved by the direct method [25], [26], or the projection-type method [27], [28], [29], [30], [31], [32], [33], [34]. Therefore, we can also consider that it has been solved well.

Regarding the difficulty referred to in (i), we recall that there is also a successful scheme developed in [11], [28] that can achieve the full decoupling structure, second-order time accuracy, and unconditional energy stability. Its key idea is to use the convex-splitting method to deal with the double-well potential and the way of implicit–explicit combination to deal with the advection and surface tension. For this scheme, we temporarily ignore its relatively expensive computational cost due to the nonlinearity and only discuss the decoupling method developed therein. It is worth noting that the implicit–explicit combination used in [11], [28] to deal with advection and surface tension is a very common method that has been widely adopted to deal with the same terms that appear in the Cahn–Hilliard–Navier–Stokes system (referred to as CH–NS) for two-phase flow, cf. [10], [11], [13], [14], [34], [35], [36], [37]. However, for the CH–NS system, this method can only generate a fully-coupled scheme. Hence, a natural question arises: when dealing with the same terms in the same way, what causes the completely different result, that is, the fully-decoupled scheme for the CHD model and the fully-coupled scheme for the CH–NS model? The answer is contained in the special format of the Darcy equation, where the pressure gradient and velocity satisfy an explicit linear relationship. Therefore, the fluid velocity can be given explicitly by the pressure, so as to achieve a fully-decoupled scheme. However, in the Navier–Stokes equation, there is no similar explicit relationship between the velocity and pressure gradient at all. Therefore the way of the implicit and explicit combination of advection and surface tension will inevitably lead to a fully-coupled scheme.

Therefore, in this paper, we aim to construct a novel full decoupling scheme for the CHD model. This scheme does not require an explicit relationship between the velocity and pressure, and can also achieve the fully-decoupling and second-order accuracy. We expect the scheme to be linear and unconditionally energy stable as well. To achieve the full decoupling structure, we take advantage of a well-known but often overlooked special property, the so-called “zero-energy-contribution” feature [38], [39], [40], [41], [42] satisfied by the two coupling terms advection and surface tension. When deducing the energy law, if a proper test function is selected to perform the inner product, these two coupling terms will cancel each other out. In other words, the contribution of these two terms to the energy is zero. We apply this property and develop a novel way to deal with these terms, by adding a specially designed ordinary differential equation (ODE) to the system that contains the inner product of these nonlinear terms. This ODE is trivial at the continuous level because all the terms in it are zero. But after discretization, it can help eliminate all the troublesome nonlinear terms to obtain unconditional energy stability. Besides, the nonlocal variable can be used to decompose all equations into two sub-equations which can be solved independently. As a result, the fully-decoupled feature is achieved.

After combining this novel decoupling method with the projection type method for the Darcy fluid equation, and the SAV method that is used to linearize the nonlinear double-well potential, we arrive at an effective numerical scheme with all the desired properties (linear, fully-decoupling, second-order accurate in time, and unconditionally energy stable). Although the developed scheme is not the first and unique fully-decoupled scheme for solving the CHD model, compared with the scheme in [11], [28], it has two advantages. First, the scheme in [28] uses the convex-splitting method, and its nonlinear nature requires more time expense (cf. Fig. 4.7 (b)), while the scheme developed here is linear and very effective in practice because we only need to solve a few fully-decoupled elliptic equations with constant coefficients at each time step. Second, the scheme developed in [28] utilizes the special property of the Darcy equation, so it is difficult to generalize to other type of fluid coupling system, while the new decoupling method can be used to construct efficient numerical schemes for many coupled nonlinear models, not only the Darcy equation considered in this paper but also the Navier–Stokes equations [14], [43], magnetic field [44], electric field [45], [46], heat equations [17], [47], etc.

The rest of the article is organized as follows. We briefly describe the CHD model and derive its energy structure in Section 2. In Section 3, a fully decoupled numerical scheme is constructed, and then its implementation is described in detail. The solvability and unconditional energy stability are proved rigorously as well. In Section 4, we implement the developed scheme and test its accuracy/stability by simulating plenty of 2D and 3D examples. Some concluding remarks are given in Section 5.