Efficient local energy dissipation preserving algorithms for the Cahn–Hilliard equation

Highlights

• We derive that Cahn–Hilliard equation possesses a local energy dissipation law (LEDL).

• Based on the observation, three LEDL schemes for the CH equation are derived.

• Three schemes are independent of boundary conditions.

• Such schemes are proven to conserve the LEDL in any local region.

• Our schemes hold the total energy stable laws and the mass laws.

Abstract

In this paper, we show that the Cahn–Hilliard equation possesses a local energy dissipation law, which is independent of boundary conditions and produces much more information of the original problem. To inherit the intrinsic property, we derive three novel local structure-preserving algorithms for the 2D Cahn–Hilliard equation by the concatenating method. In particular, when the nonlinear bulk potential $f(\varphi )$ in the equation is chosen as the Ginzburg–Landau double-well potential, the method discussed by Zhang and Qiao (2012) [50] is a special case of our scheme II. Thanks to the Leibnitz rules and properties of operators, the three schemes are rigorously proven to conserve the discrete local energy dissipation law in any local time–space region. Under periodic boundary conditions, the schemes are proven to possess the discrete mass conservation and total energy dissipation laws. Numerical experiments are conducted to show the performance of the proposed schemes.

Introduction

There has been an increasing emphasis on the free interface problem of multi-phase incompressible flows which is often modeled by the diffuse interface method with a long history [39], [44]. Its main idea is to apply one or more continuous functions to describe the volume fraction in the multi-phase immiscible system and represent the fluid interfaces by a thin and smooth transition layers, see, e.g., [1], [3], [7], [20], [22], [33], [34], [37], [47] and [38], [43] for details. The Cahn–Hilliard equation induced by the variational approach, is a widely used phase-field model.

In this paper, we consider the following Cahn–Hilliard (CH) equation [6]

$$\{\begin{array}{c}\frac{\partial \varphi }{\partial t}=\mathrm{\Delta }\mu ,\hfill \\ \mu =f^{\prime }(\varphi )-\kappa \mathrm{\Delta }\varphi ,\hfill \end{array}$$

where ϕ is an unknown real-valued phase function, which stands for the relative concentration of one phase, μ is the chemical potential which is the variation of the free energy, κ is a parameter measuring the strength of the conformational entropy. Here the free energy is chosen as

$$E[\varphi ]=\underset{\mathrm{\Omega }}{\int }(\frac{\kappa }{2}|\mathrm{\nabla }\varphi {|}^2+f(\varphi ))d\mathbf{x},$$

where Ω is the domain that the fluid occupies, and $f(\varphi )$ is the bulk energy density. For the choice of nonlinear potential $f(\varphi )$, Ginzburg–Landau double-well type potential,

$$f(\varphi )=\frac{\gamma }{4}{(1-{\varphi }^2)}^2,\varphi \in [-1,1],$$

and Logarithmic Flory–Huggins potential [41],

$$f(\varphi )=\gamma (\varphi -{\varphi }^2)+(1-\varphi )\ln (1-\varphi )+\varphi \ln \varphi ,\varphi \in [0,1],$$

are two widely used nonlinear bulk potentials in the literature, where γ measures the strength of the repulsive potential. As is well known, the CH equation, with suitable boundary conditions, possesses the mass conservation law (MCL)

$$I(t):=\underset{\mathrm{\Omega }}{\int }\varphi (\mathbf{x},t)d\mathbf{x}\equiv I(0),$$

and the energy dissipation law (EDL)

$$\frac{dE}{dt}=-\underset{\mathrm{\Omega }}{\int }|\mathrm{\nabla }\mu {|}^{2}d\mathbf{x}.$$

To date, extensive mathematical studies have been carried out for the CH equation, see, e.g., [7], [8], [9], [10], [19], [21], [31], [32], [37], [40], [42], [47] for details. For the analytical front, the steady state solution, the existence and the uniqueness of the system, we refer to [13], [24] and the references therein. For the numerical front, various numerical methods have been considered in the literature. Elliott et al. [12] discussed the physical background and its numerical solution by the Galerkin finite element method. Later, Feng and Prohl [17] developed the mixed finite element method which was convergent with quasi-optimal order in time and optimal order in space. Sun [40] proposed the finite difference method by the method of reduction of order, which is conditionally stable and not mass conserving. Choo et al. [5] gave a conservative nonlinear finite difference scheme, which was proven to preserve the total mass, yet the energy-based stability was not discussed. Subsequently, Furihata [16] gave a conservative finite difference scheme, which was proven to be unconditionally stable in the sense of energy decay. Chen and Shen [8] developed a semi-implicit Fourier spectral scheme, in which the linear part was treated implicitly and the nonlinear part was calculated explicitly. Later, Zhu et al. [49] proposed the semi-implicit Fourier spectral method for the CH equation with a variable mobility. Shen and Yang [42] obtained energy stability method without any a priori assumption on the numerical solution. Zhang and Qiao [50] discussed a finite difference scheme and proposed an adaptive time-stepping technique to quickly solve the 2D CH equation. He et al. [31] proposed a large time-stepping methods. Guo et al. [27] used the convex splitting skills, whose idea was proposed more than two decades earlier by Elliott et al. and Eyre in Ref. [14], [15]. Li et al. [35], [36] developed the semi-implicit Fourier spectral method, which was proved unconditional energy stability for modified energy functionals by the new stabilization techniques. Recently, the new technique called “Invariant Energy Quadratization” (IEQ) which was successfully applied to different phase-field models by authors [30], [46], [48], was extended to handle the CH equation in Ref. [25], [26]. The efficient iterative solution of linear system was investigated for the coupled Navier–Stokes Cahn–Hilliard system, in which the block-triangular preconditioners were introduced in Ref. [2]. Gomez et al. gave a good reference book for stable time integration of phase-field equations in Ref. [23]. We notice that most existing methods and their stabilities which are discussed in the literature are based on the global energy dissipation property which is defined on the global space region and depends on the suitable boundary conditions. Here in this paper, we explore to construct novel schemes basing on our observation that CH equation has a local energy dissipation law which is valid in any time–space region and independent of the boundary conditions. It is known that boundary condition of the system is very important in many practical problems.

Proposition 1

The system (1.1) possesses the local energy dissipation law (LEDL)

$${\partial }_t(\frac{\kappa }{2}|\mathrm{\nabla }\varphi {|}^2+f(\varphi ))-\mathrm{\nabla }\cdot (\kappa {\varphi }_t\mathrm{\nabla }\varphi +\mu \mathrm{\nabla }\mu )+|\mathrm{\nabla }\mu {|}^2=0,$$

which can be derived by multiplying μ and ${\varphi }_t$ to the first line and the second line of (1.1) respectively, and adding the two resluting formulas together. Compared with (1.4), the LEDL (1.5) has an additional flux term which can supply much more information of the original system. Moreover, the LEDL is independent of boundary conditions. Therefore, it is valid in any time space domain for the CH equation with any boundary condition. Under suitable boundary conditions, the LEDL (1.5) implies the total EDL (1.4). This observation motivates us to construct numerical schemes preserving the LEDL (1.5) for the CH equation.

Actually, how to construct numerical schemes preserving certain local invariant quantities for the continuous dynamical systems is brought into sharp focus in recent years. The numerical methods preserving some structural properties of the system are sometime referred as geometric integrators or structure-preserving numerical methods. So far, on construction of the local structure-preserving schemes for some classical PDEs, many results have been achieved [28], [45] which further show that numerical methods with local property supply richer information of the original system than that of schemes with global property.

In this paper, we aim to construct three novel numerical schemes to preserve discrete analogues of the local energy dissipation law (1.5) for the 2D Cahn–Hilliard equation which are achieved through the concatenation of basic algorithms suggested by Wang et al. [45], including the implicit midpoint method, the leap-frog method and the discrete variational derivative method proposed by Furihata et al. [11], [18] for the nonlinear term, on the associated first-order subsystems. The corresponding discrete local energy dissipation laws are rigorously proved by the Leibnitz rules, with no dependence of the boundary conditions. While imposed of the periodic boundary condition or the homogeneous boundary condition, the proposed schemes thereby possess the global EDL as well as MCL.

The remainder of this paper is arranged as follows. Some operator definitions and theirs properties are given in Section 2. In Section 3, we construct three novel numerical schemes preserving the LEDL for solving the CH equation. When the equation is imposed on appropriate boundary conditions, the schemes will conserve the global MCL and EDL. Numerical experiments are presented in Section 4 and some concluding remarks are given in the final section.