Hessian recovery based finite element methods for the Cahn-Hilliard Equation
HHESSIAN RECOVERY BASED FINITE ELEMENT METHODS FORTHE CAHN-HILLIARD EQUATION
MINQIANG XU ∗ , HAILONG GUO † , AND
QINGSONG ZOU ‡ Abstract.
In this paper, we propose a novel recovery based finite element method for theCahn-Hilliard equation. One distinguishing feature of the method is that we discretize the fourth-order differential operator in a standard C linear finite elements space. Precisely, we first transformthe fourth-order Cahn-Hilliard equation to its variational formulation in which only first-order andsecond-order derivatives are involved and then we compute the first and second-order derivatives of alinear finite element function by a least-square-fitting recovery procedure. When the underlying meshis uniform of regular pattern, our recovery scheme for the Laplacian operator coincides with the well-known five-point stencil. Another feature of the method is some special treatments on Neumann typeboundary conditions for reducing computational cost. The optimal-order convergence properties andenergy stability are numerically proved through a series of benchmark tests. The proposed methodcan be regarded as a combination of the finite difference scheme and the finite element scheme. AMS subject classifications.
Primary 65N30; Secondary 45N08
Key words.
Hessian recovery, Cahn-Hilliard equation, phase separation, recovery based, su-perconvergence, linear finite element.
1. Introduction.
The phase field model is a powerful tool to characterize in-terfacial problems in which the dynamics of the physical systems are described by agradient flow. The Cahn-Hilliard equation [7] is a famous phase field model intro-duced by Cahn and Hilliard to model the phase separation in binary alloys. Later on,it is widely used to model multiphase flow [4, 11], tumor growth [36, 41], and imageimpainting [5] etc.As a nonlinear parabolic type equation, the analytic solution of the Cahn-Hilliardequation is usually hard to be obtained. Numerical simulation looks like to be the onlyfeasible way to study the physical problems governed by the Cahn-Hilliard equation.During the past several decades, a huge number of numerical methods have beendeveloped in the literature, including finite different methods [8], spectral methods[38], and finite element methods [18, 40, 43]. In this paper, we concentrate on finiteelement methods for the Cahn-Hilliard equation.One of the main difficulties in the numerical solution of the Cahn-Hilliard equa-tion is the discretization of the fourth-order differential operator in a certain finiteelement space. In the literature, finite element methods for fourth-order elliptic equa-tions can be roughly categorized into the following four classes: conforming finiteelement methods [6, 12], nonconforming finite element methods [33], mixed finite ele-ment methods [13], discontinuous Galerkin method [20]. In corresponding, the Cahn-Hilliard equation has been numerically solved by conforming finite element methodsin [15, 19], nonconforming finite element methods in [16, 43], mixed finite elementmethods in [17,22], and discontinuous Galerkin methods [40,42]. All the above meth-ods in the primary form discretize the Cahn-Hilliard equation at least in a quadratic ∗ School of Data and Computer Science, Sun Yat-sen University, Guangzhou 510275, China([email protected]). † School of Mathematics and Statistics, The University of Melbourne, Parkville, VIC 3010, Aus-tralia ([email protected]). ‡ School of Data and Computer Science and Guangdong Province Key Laboratory of Computa-tional Science, Sun Yat-sen University, Guangzhou 510275, China ([email protected])1 a r X i v : . [ m a t h . NA ] N ov nite element space, which means that there are at least six degrees of freedom oneach triangular element. To reduce the complexity, a new class of finite elementmethods, called recovery based finite element methods [10, 25, 28], is proposed tosimulate fourth-order partial differential equations. The key idea of those methodsis to facilitate the simplest element, the continuous linear element, to discretize thefourth-order differential operator. It is known that the second-order derivative of C piecewise linear function is not well-defined and thus usually we can not solve afourth-order differential equation in a C linear finite element space. Such a barrier isalleviated by using the classical gradient recovery operator G h [44, 46] to smooth thediscontinuous piecewise constant function into a continuous piecewise linear functioni [10, 25].In this paper, we will also discretize the Cahn-Hilliard equation only in the sim-plest linear element space. Comparing to the recovery-based FEMs in [10, 25, 28], thedifference here is that we directly recover the Hessian matrix of a linear finite elementfunction instead of recovering its gradient. Note that the Hessian recovery has beenstudied for the purpose of post-processing [1, 23, 37, 39]. In particular, in [23], Guoet al. proposed a new Hessian recovery method and established its complete super-convergence theory on mildly unstructured meshes and ultraconvergence theory onstructured meshes. The Hessian recovery technique in [23] is then applied to solvinga sixth-order PDE in [26]. In this paper, we use a Hessian recovery technique byfirstly recovering a local quadratic polynomial and then taking second order deriva-tives of the recovered polynomial as the second-order derivatives of the linear finiteelement function. We sprucely discover that there is an intrinsic connection betweenthe Hessian recovery method and the finite difference method. In specific, we findthat the Hessian recovery method reproduces the standard five-point finite differencescheme on regular pattern uniform meshes. This means, on the regular pattern uni-form meshes, the new recovery based finite element method is a kind of infusion of thefinite difference method and the finite element method in the sense that we first usethe standard five-point finite difference scheme to discretize the Laplacian operatorand then put it back into the standard linear finite element framework.Different from second-order elliptic equations, the Neumann boundary conditionfor fourth-order partial differential equations is an essential boundary condition inthe sense that the boundary condition should be enforced in their associate solutionspaces. But it looks like impossible to be enforced into a C linear finite element space.In our previous paper [10, 25, 26], it is imposed by the penalty method [14, 45] andthe Lagrange multiplier method [3]. Like it for the second order partial differentialequations, the resulting linear system of the penalty method is usually ill-conditioned.The Lagrange multiplier method shows the potential to overcome such drawback butit introduces additional degrees of freedom. In this paper, we adopt two differentmethods to deal with the Neumann boundary condition. On general unstructuredmeshes, we propose to impose the boundary condition weakly based on a techniquecalled Nitsche’s method [35], which is originally introduced by Nitsche to incorporatethe Dirichlet boundary condition for second order elliptic equations. A Nitsche’svariational formulation for the Cahn-Hilliard equations is presented. It paves theway for implicitly imposing the Neumann boundary conditions in Hessian recoverybased finite element methods. As mentioned in the previous paragraph, the Hessianrecovery method reduces to the standard five-point standard finite difference schemeon uniform meshes. Such key observation enables us to incorporate the Neumannboundary condition into the Hessian recovery operator by the celebrated ghost point ethod in finite difference methods [29].The rest of paper is organized as follows: In Section 2, we present a simpleintroduction to the Cahn-Hilliard equation and review some of its property. In Section3, we first revisit the Hessian recovery method and uncover its relationship with theclassical finite difference method; then, we propose the new recovery based finiteelement method to discretize the spatial variable; the fully discrete formulation isdiscussed through an energy stable time stepping method. The proposed method isnumerically verified and validated using a series of benchmark examples in Section 4.We end with some conclusive remarks in Section 5.
2. The Cahn-Hilliard equation and its variational formulations.
Let Ωbe a bounded polygonal domain with Lipschitz boundary ∂ Ω in R . For a subdomain A of Ω, let P m ( A ) be the space of polynomials of degree less than or equal to m over A and n m be the dimension of P m ( A ) with n m = ( m + 1)( m + 2). We denote by H k ( A ) the Sobolev space with norm (cid:107) · (cid:107) k, A and seminorm | · | k, A .The well-known Cahn-Hilliard equation on a space domain Ω and a certain timeperiod [0 , T ] can be described as below: ∂u∂t = − ε ∆ u + ∆ F (cid:48) ( u ) , in Ω × [0 , T ] ,∂ n u = ∂ n ( − ε ∆ u + F (cid:48) ( u )) = 0 , on ∂ Ω × [0 , T ] ,u ( · ,
0) = u ( · ) , in Ω , (2.1)where n is the unit outer normal vector of ∂ Ω. The unknown function u often indicatesthe concentration of one of the two metal components constituting the alloys, ε is thesize of the interface of two alloys, F is a double well nonconvex function defined as F ( u ) = ( u − . The Cahn-Hilliard equation (2.1) can be viewed as an H − -gradient flow of theGinzburg-Landau free energy functional E ( u ) := (cid:90) Ω (cid:18) ε |∇ u | + F ( u ) (cid:19) dx, (2.2)of which the first part is called the interfacial energy and the second part is called the bulk energy . Thanks to the homogeneous Neumann boundary conditions, it is easy toverify that the mass conservation property ddt (cid:90) Ω udx = 0 , and the energy decay property dE ( u ) dt = −(cid:107) − ε ∆ u + F (cid:48) ( u )) (cid:107) ≤ , ∀ t > , always hold for the solution of the Cahn-Hilliard equation.Let f = F (cid:48) . To implicitly impose the Neumann boundary condition ∂ n u = 0, weintroduce the bilinear form a ( w, v ) = (cid:90) Ω ∆ w ∆ vdz − (cid:90) ∂ Ω ∆ w∂ n vds − (cid:90) ∂ Ω ∂ n w ∆ vds + γ (cid:90) ∂ Ω ∂ n w∂ n vds, ∀ v, w ∈ H (Ω) , (2.3)with γ is a positive stability parameter to be specified in the sequel. It is easy toverify that if u ∈ H (Ω) is the solution of (2.1), then u satisfies (cid:18) ∂u∂t , v (cid:19) + ε a ( u, v ) + ( ∇ f ( u ) , ∇ v ) = 0 , ∀ v ∈ H (Ω) . (2.4) onversely, if u ∈ C (Ω) and u t ∈ C (Ω) which satisfies (2.4), then for all v ∈ H (Ω),we have (cid:18) ∂u∂t + ε ∆ u − ∆ f ( u ) , v (cid:19) − (cid:90) ∂ Ω ∂ n u ∆ vds + γ (cid:90) ∂ Ω ∂ n u∂ n v + (cid:90) ∂ Ω ∂ n ( − ε ∆ u + F (cid:48) ( u )) vds = 0 . (2.5)Choosing v ∈ C ∞ (Ω) in (2.5), we obtain (cid:18) ∂u∂t + ε ∆ u − ∆ f ( u ) , v (cid:19) = 0 . which implies ∂u∂t + ε ∆ u − ∆ f ( u ) = 0. Consequently, (2.5) becomes − (cid:90) ∂ Ω ∂ n u ∆ vds + γ (cid:90) ∂ Ω ∂ n u∂ n v + (cid:90) ∂ Ω ∂ n ( − ε ∆ u + f ( u )) vds = 0 , v ∈ H (Ω) . Since v ∈ H (Ω) is arbitrary, we derive from the above equation that ∂ n u = ∂ n ( − ε ∆ u + f ( u )) = 0It means u is the classical solution of (2.1). From the above reasonings, we obtainthat the solution of (2.4) is a weak solution of (2.1) and we call (2.4) a variationalformulation of (2.1). Remark 2.1.
An alternative variational equation of (2.1) is ( ∂u∂t , v ) + ε a ( u, v ) + ( ∇ f ( u ) , ∇ v ) = 0 , ∀ v ∈ H (Ω) (2.6) where a ( w, v ) = (cid:90) Ω D w : D vdx − (cid:90) ∂ Ω ∂ n w∂ n vds − (cid:90) ∂ Ω ∂ n w∂ n vds + γ (cid:90) ∂ Ω ∂ n w∂ n vds, (2.7) and A : B is the Frobenius norm of × matrices. We observe that here, the bilinearform a ( · , · ) differs from a ( · , · ) by replacing the Laplace operator ∆ with the Hessianmatrix operator D .In both variational formulations, the Neumann boundary conditions ∂ n u = 0 isweakly built into the bilinear formations. The idea is similar to impose the Dirichletboundary condition for the second-oder elliptic equations by Nitsche [35]. We callthose two methods the Nitsche’s method. Remark 2.2.
Originally, the Neumann boundary condition for fourth-order par-tial differential equations is enforced into their solution space. For such purpose, let V = { v ∈ H (Ω) : ∂ n v = 0 on ∂ Ω } . (2.8) The bilinear forms a and a in V reduce to a ( w, v ) = (cid:90) Ω ∆ w ∆ vdx, w, v ∈ V (2.9) and a ( w, v ) = (cid:90) Ω D w : D vdx, w, v ∈ V (2.10) espectively. Correspondingly, the variational formulations become to find u ∈ ( L ([0 , T ]; V ) such that (cid:18) ∂u∂t , v (cid:19) + ε a ( u, v ) + ( ∇ f ( u ) , ∇ v ) = 0 , ∀ v ∈ V, (2.11) or, to find u ∈ ( L ([0 , T ]; V ) such that (cid:18) ∂u∂t , v (cid:19) + ε a ( u, v ) + ( ∇ f ( u ) , ∇ v ) = 0 , ∀ v ∈ V. (2.12)
3. A novel recovery based finite element method.
In this section, we de-sign novel recovery-technique-based finite element methods for Cahn-Hilliard equa-tions. Since our main attention is on novel space discretization techniques, for thetime discretization, we choose a simple one-step energy stable linear scheme proposedin [27,38,43]. Precisely, we make use of a semi-implicit scheme with an extra stabilizedpenalty term added to ensure energy stability. Let the time step size be ∆ t = TN , u ( x ) = u ( x, u n ( x ) ≈ u ( x, n ∆ t ) , n = 1 , , . . . , N , then the stabilized first-ordersemi-implicit method reads as: find u n ∈ S, n = 1 , , . . . N such that for all v ∈ S , (cid:18) u n +1 − u n ∆ t , v (cid:19) + ε a i ( u n +1 , v )+ ( ∇ f ( u n ) , ∇ v ) + κ (cid:0) ∇ ( u n +1 − u n ) , ∇ v (cid:1) = 0 , i = 1 , , , , (3.1)where S = H (Ω) for i = 1 , S = V for i = 3 ,
4. Note that the choice of κ hasa great influence on the stability of (3.1), see [30, 32, 38, 43] for the details. In thispaper, we always take κ = 2. In fact, [30, 32] give you rigorous analysis on the chooseof κ .Next we explain how to discretize (3.1) in a linear finite element space. Let T h be a shape regular triangulation of Ω with mesh size h . The set of all vertices andof all edges of T h are denoted by N h and E h , respectively. We define the standardcontinuous linear finite element space S h on T h by S h := (cid:8) v h ∈ C (Ω) : v h | T ∈ P , ∀ T ∈ T h (cid:9) . (3.2)and denote its nodal basis by { φ z } z ∈N h .To construct our fully discrete schemes on the linear finite element space S h , wefirst introduce a Hessian recovery technique based on least-squares fitting in the firstsubsection. Then we apply the Hessian recovery operator to develop our novel fullydiscrete schemes for the Cahn-Hilliard equation in the second subsection. It isknown that the second order derivative of a function v h ∈ S h equals to 0 in theinterior of each element T ∈ T h and is not well-defined on an edge E ∈ E h . In thefollowing, we propose a least-square type method to calculate the approximate secondorder derivatives of v h . In other words, we will define a Hessian recovery operator H h from S h to S h which maps a function v h ∈ S h to H h v h ∈ S h so that H h v h can beregarded as an approximation of the Hessian matrix of v h in some sense.Since H h v h ∈ S h , to define H h v h , it is sufficient to define the value ( H h v h )( z )for all z ∈ N h . For this purpose, we first construct a local patch associated with z hich is a polygon surrounding the node z . Given a vertex z ∈ N h and a nonnegativeinteger n ∈ N , let the first n layer element patch be L ( z, n ) = { z } , if n = 0 , (cid:83) { T : T ∈ T h , T ∩ L ( z, (cid:54) = φ } , if n = 1 , (cid:83) { T : T ∈ T h , T ∩ L ( z, n −
1) is an edge in E h } , if n ≥ . (3.3)For all z ∈ N h , let n z be the smallest integer such that L ( z, n ) satisfies the rankcondition in the following sense. Definition 3.1.
A surrounding z polygon is said to satisfy the rank condition ifit admits a unique least-squares fitted polynomial p z in (3.4) . We define the local patch associated with z as Ω z = (cid:32)L( z, n z ). Using the verticesin Ω z as sampling points, we fit a quadratic polynomial p z at the vertex z in thefollowing least-squares sense p z = arg min p ∈ P (Ω z ) (cid:88) x ∈∈ Ω z ∩N h | p ( x ) − v h ( x ) | . (3.4)Then, we define the recovered Hessian node value by( H h v h )( z ) = (cid:18) H xxh v h ( z ) H xyh v h ( z ) H yxh v h ( z ) H yyh v h ( z ) (cid:19) = (cid:32) ∂ p z ∂x ( z ) ∂ p z ∂x∂y ( z ) ∂ p z ∂y∂x ( z ) ∂ p z ∂y ( z ) (cid:33) . (3.5)With this definition, we have H h v h = (cid:80) z ∈N h H h v h ( z ) φ z and the symmetric property H xyh = H yxh of the Hessian matrix function H h v h . Moreover, based on H h , we definea discrete Laplacian operator ∆ h : S h → S h as∆ h v h = H xxh v h + H yyh v h . (3.6)Note that in the same way, we can recover the gradient of v h by letting( G h v h )( z ) = ∇ p z ( z ) , ∀ z ∈ N h (3.7)and G h v h = (cid:80) z ∈N h G h v h ( z ) φ z ∈ S h . Note that the Hessian recovery operator H h in (3.5) has been applied to post-process the finite element solution in [23]. According to the numerical results, some-times it might be inefficient or even not convergent as a post-processing operator.Since it involves a relatively smaller number of neighbourhood vertices in its stencil,here we choose it as our pre-processing operator. In fact, H h can be regarded as aspecial finite difference operator of the second order on uniform meshes and of thefirst order on general unstructured meshes. To elucidate this basic idea, we consider aspecial case that T h is a regular pattern uniform triangular mesh. In this case, for aninterior node z = z i , the local patch Ω z i is defined as the polygon z · · · z , see Figure3.1a and thus the sampling points in Ω z i include z i , z i , · · · , z i with z i = z i . Usingthese seven sampling points, we fit a quadratic polynomial p z i in the least-squaressense and take second-order differentiation, which produces( H xxh u )( z i ) = 1 h ( u − u + u ) , ( H xyh u )( z i ) = 12 h (2 u − u + u − u − u + u − u ) , ( H xxh u )( z i ) = 1 h ( u − u + u ); here the values u j is defined by u j = u ( z i j ). By the definition of the discreteLaplacian operator (3.6), we have(∆ h u )( z i ) = 1 h ( u + u − u + u + u ) . (3.8)The formula (3.8) implies the discrete Laplacian operator on regular pattern uniformmeshes is the well-known five-point-finite-difference stencil of the Laplace operator,as illustrated in Figure 3.1b. By the standard approximation theory, we have (cid:107) ∆ u − ∆ h u I (cid:107) , Ω ≤ Ch (cid:107) u (cid:107) , Ω . (3.9)This exciting discovery implies that the Hessian recovery operator can be regarded asan extension of the classic second-order difference operator on regular meshes to thedifference operator on non-uniform meshes and thus it can be used to design discreteschemes for higher-order differential equations. z i z i z i z i z i z i z i (a) − − − − (b) Fig. 3.1: Illustration of Hessian recovery on uniform mesh: (a) local patch ; (b)Discrete Laplace operator.
Remark 3.1.
For a boundary vertex z ∈ ∂ Ω , there are other approaches toconstruct the local patch Ω z , see [24] for the details. To present our fully discrete schemes, we firstintroduce the discrete bilinear form a i,h ( · , · ) on the linear finite element space S h . Forall w h , v h ∈ S h , we define a ,h ( w h , v h ) = (cid:90) Ω ∆ h w h ∆ h v h dz − (cid:90) ∂ Ω ∆ h w h ( G h v h · n ) ds − (cid:90) ∂ Ω ( G h w h · n )∆ h v h ds + γ (cid:90) ∂ Ω ( G h w h · n )( G h v h · n ) ds, (3.10)and a ,h ( w h , v h ) = (cid:90) Ω H h w h H h v h vdz − (cid:90) ∂ Ω ( n T H h w h n )( G h v h · n ) ds − (cid:90) ∂ Ω ( G h w h · n )( n T H h v h n ) ds + γ (cid:90) ∂ Ω ( G h w h · n )( G h v h · n ) ds, (3.11) here γ = Ch with C a sufficiently large positive constant.The fully discrete Hessian recovery based finite element method for the Cahn-Hillliard equation (2.1) reads as : find { u nh } n ≥ ∈ V h such that for all v h ∈ S h , (cid:18) u n +1 h − u nh ∆ t , v h (cid:19) + ε a i,h ( u n +1 h , v h ) + κ ( ∇ u n +1 h − ∇ u nh , ∇ v h ) + ( ∇ f ( u nh ) , ∇ v h ) = 0 . (3.12)Note that both the schemes in (3.12) work on general unstructured meshes. More-over, our numerical experiments indicate that there is no essential difference betweenthese two schemes. We also observe that the stiffness matrices corresponding to bothschemes (3.12) are symmetric and positive definite, so both schemes are stable anduniquely solvable.Next, we present a simple fully discrete scheme derived from the variational for-mulation (2.11) on uniform meshes. We define the finite element space S h, ⊂ S h as S h, = { v h ∈ S h : ∇ h v h ( z ) · n = 0 , ∀ z ∈ N h ∩ ∂ Ω } , (3.13)where ∇ h is a discrete gradient operator so that S h, be a discrete analogous of V .Note that in the continuous linear finite element space S h , ∇ v h · n is not well definedon a vertex of T h , so at each boundary vertex, we use a central finite difference schemeinstead of ∇ v h ( z ) to define ∇ h v h ( z ).The key part in construction of a simpler scheme on the uniform meshes isbased on the fact the discrete Laplacian operator ∆ h reduces to the five-point-finite-difference stencil at an interior vertex, as illustrated in (3.8). Now, we suppose ∆ h is the discrete Laplacian operator on the finite element space S h, . Different fromthe general treatment of the recovery on the boundary as introduced in the previ-ous subsection, we borrow the idea of ghost point method from the finite differencemethod [29]. In specific, at every boundary vertex z i , we introduce one or more ghostpoints. Then, we still have the fact that the the discrete Laplacian operator ∆ h isjust the five-point finite difference stencil at the boundary vertex z i but it involves thevalue of the finite element function at the ghost points. To eliminate it, we combinethe discrete boundary condition ∇ h v h ( z i ) · n = 0 which also involves the same ghostpoints.To illustrate idea, we consider a typical boundary vertex z i , as illustrated inFigure 3.2a. In that case, z i is a boundary vertex with three neighbour mesh vertices z i , z i , z i . To apply the five-point finite difference scheme at z i , we introduce a ghostpoint z i , as the red dot point in Figure 3.2a, and the discrete Laplacian ∆ v h ( z i ) is∆ v h ( z i ) = 1 h ( v h ( z i ) + v h ( z i ) − v h ( z i ) + v h ( z i ) + v h ( z i )) , (3.14)which involves the ghost point finite element function value v h ( z i ). By the definitionof the finite element space S h, , at the boundary vertex z i , we also ∇ h v h ( z i ) · n = 12 h ( v h ( z i ) − v h ( z i )) = 0 . (3.15)Using (3.14) and (3.15) to eliminate v h ( z i ), we obtain∆ v h ( z i ) = 1 h ( v h ( z i ) + 2 v h ( z i ) − v h ( z i ) + v h ( z i ))) , (3.16) hich only depends on the value at the vertices in N h . In other word, we haveembedding the discrete Neumann boundary condition into the discrete Laplacianoperator ∆ h . Similarly, we can explicitly construct the discrete Laplacian operator∆ h at a corner boundary vertex, which may need two ghost points as plotted inFigure 3.2b. In this case, the computation of ∆ h does not need to use an implicitleast-squares fitting process. z i z i z i z i z i (a) z i z i z i z i z i (b) Fig. 3.2: Illustration of ghost point method: (a). One ghost point; (b) Two Ghostpoints.
Remark 3.2.
The key processing is to build the discrete Neumann boundarycondition into the discrete Laplacian operator ∆ h . Such process is only possible for theuniform meshes. For general unstructured meshes, we can define a similar discretefinite element space as S h, = { v h ∈ S h : G h v h ( z i ) · n = 0 , ∀ z i ∈ N h ∩ ∂ Ω } . Butthe discrete boundary condition G h v h ( z i ) · n = 0 can not be embedded into discreteLaplacian operator ∆ h . We have to use other methods like the penalty method [14,45]and the Lagrange multiplier method [3] to impose the discrete Neumann boundarycondition G h v h ( z i ) · n = 0 , ∀ z i ∈ N h ∩ ∂ Ω . However, these two methods perform badlyfor the Cahn-Hilliard equation. Then the discrete bilinear a ,h ( · , · ) on S h, as a ,h ( w h , v h ) = (cid:90) Ω ∆ h w h ∆ h v h dx, ∀ v h , w h ∈ S h, . (3.17)A simple fully discrete for (2.1) on uniform meshes is to find { u nh } n ≥ ∈ S h, suchthat for all v h ∈ S h, , (cid:18) u n +1 h − u nh ∆ t , v h (cid:19) + ε a ,h ( u n +1 h , v h ) + κ ( ∇ u n +1 h − ∇ u nh , ∇ v h ) + ( ∇ f ( u nh ) , ∇ v h ) = 0 . (3.18)Since the bilinear form a ,h ( · , · ) does not involve the computation of the gradientrecovery operator G h either, the scheme (3.18) is very computationally efficient andaccurate. Moreover, the scheme (3.18) can be regarded as a mixture of the finitedifference method and the finite element method since we first use the finite differenceoperator ∆ h to recover the second order derivatives of a linear finite element function nd then bring them back to the framework of the finite element method. Namely, thescheme (3.18) sheds some light on using the finite difference operators to constructsimple finite element methods for higher-order partial differential equations.It may worth mentioning that the gradient recovered method proposed in [25]for fourth-order problems use the minimum number of degrees of freedom among thefinite element spaces, see for details. However, our numerical experiments show thata direct application of the gradient recovered method to the Cahn-Hilliard equationleads to an unstable scheme. Moreover, compared with the gradient recovered method,the present Hessian recovered method uses the same number of total degree butits stiffness matrix is more sparse than the one derived from the gradient-recoveredmethod.
4. Numerical Experiments.
In this section, we present several numerical ex-amples to demonstrate the properties of our proposed methods.Except for the last numerical example, the domain Ω of the problems in thissection is chosen as the unit square [0 , . In our experiments, we will adopt twodifferent types of meshes: the uniform and unstructured meshes. Our uniform meshesare generated by first dividing Ω into m congruent subsquares and then splittingeach subsquare into two right-angled triangles, see Figure 4.1a. Our unstructuredmeshes are generated by the first partition of the domain with the Delaunay meshgenerator EasyMesh [34] to obtain the first level mesh and then uniformly refines eachtriangle in the first level mesh several times, see Figure 4.1b. On a uniform mesh,we will use the fully discrete scheme (3.18), while on an unstructured mesh, we willuse the scheme (3.12) with i = 1 , κ = 2 , C = 1 to solve numerically the Cahn-Hilliardequations. (a) (b) Fig. 4.1: (a) A uniform mesh; (b) An unstructured mesh.Throughout this section,we define discrete interface energy and discrete buckenergy at time t n respectively as : E n = (cid:15) (cid:88) τ ∈T h (cid:90) τ | G h u nh | d x d y, E n = 14 (cid:88) τ ∈T h (cid:90) τ (( u nh ) − d x d y. Moreover, we denote different kinds of numerical errors by e = (cid:107) u − u h (cid:107) , Ω , e = (cid:107)∇ u − ∇ u h (cid:107) , Ω ,e ,r = (cid:107)∇ u − G h u h (cid:107) , Ω , e = (cid:107) D u − H h u h (cid:107) , Ω , nd use r = log( e h /e h )log(2) to indicate the convergence rates. We consider the non-homogeneous Cahn-Hilliardequation (cid:26) ∂u∂t = − ε ∆ u + ∆( u − u ) + g, in Ω × [0 , T ] ,∂ n u = ∂ n ∆ u = 0 , on ∂ Ω × [0 , T ] , (4.1)with the parameter (cid:15) = 0 .
1. The initial solution u and g are chosen such that theexact solution is u ( x, y, t ) = e − t cos( πx ) cos( πy ) . To compute the convergence rates with respect to the space meshsize h , we fix thetime step size ∆ t = 10 − and study the convergence order of the numerical solutionat T = 0 . i = 1 , C = 1 , κ = 2) on unstructured meshes. The numerical results by (3.18)and (3.12) are depicted in Table 4.1 and Table 4.2, respectively. From these twotables, we observe that for both schemes, the L -norm errors converge with order2 while the H -seminorm errors converge with order 1 which are both optimal fora linear finite element method. We also observe that the recovered H -seminormerror is superconvergent of O ( h ) while the recovered H -seminorm errors convergesoptimally with order 1.To test the convergence rate of the scheme (3.18) with respect to the time dis-cretization, we fix the spacial mesh size h = 1 / T = 0 .
01 with different time step ∆ t are shown in Table 4.3. The numericalresults evidently indicate that the scheme (3.18) is of first order in time, which isconsistent with the first-order semi-implicit scheme.Table 4.1: Spatial errors and convergence rates by scheme (3.18) for Example 1 h e r e r e ,r r e r /
16 1 . × − . × − . × − . × − /
32 4 . × − . × − . . × − . × − /
64 9 . × − . × − . . × − . × − /
128 2 . × − . × − . . × − . × − /
256 6 . × − . × − . . × − . × − Table 4.2:
Spatial errors and convergence rates by scheme (3.12) for Example 1 dof e r e r e ,r r e r
513 1 . × − . × − . × − . × − . × − . × − . . × − . . × − . × − . × − . . × − . . × − . × − . × − . . × − . . × − Example 2:
We consider the Cahn-Hilliard equation (2.1) with the parameter (cid:15) = 0 . u = cos( πx ) cos( πy ).We use the simple scheme (3.18) to compute the numerical results. As the exactsolution of Example 2 is unknown, we use the computable quantity u h − u h to replacethe “true error” e = u − u h in our real computations. As in the previous example, able 4.3: Temporal errors and convergence rate by scheme (3.18) for Example 1 ∆ t − − / − / − / re . × − . × − . × − . × − we fix ∆ t = 10 − and T = 0 . h = 1 /
128 and T = 0 .
01. The numerical results with differenttime step ∆ t are presented in Table 4.5. We observe that the scheme (3.18) has afirst-order accuracy in time discretization.Table 4.4: Spatial errors and convergence rates by scheme (3.18) for Example 2 h e r e r e ,r r e r /
16 1 . × − . × − . × − . × − /
32 3 . × − . . × − . . × − . × − /
64 8 . × − . . × − . . × − . × − /
128 2 . × − . . × − . . × − . × − Table 4.5:
Temporal errors and convergence rate by scheme (3.18) for Example 1 ∆ t − − / − / − / re . × − . × − . × − . × − In this subsection, we numerically solve theCahn-Hilliard equation to show the spinodal decomposition: a process or phenomenonto rapid unmix a mixture of liquids or solids from one thermodynamic phase, to formtwo coexisting phases. In the following examples, we apply the scheme (3.18) on theuniform triangular mesh with the space stepsize h = 1 /
128 and time stepsize ∆ t =10 − . Since the numerical results computed by the scheme (3.12) on unstructuredmeshes are similar to those by (3.18), they will be not reported here. Example 3:
We consider the Cahn-Hilliard equation suggested in [43] where theparameter (cid:15) = 0 .
02 and the initial value is given by u = 10 − sin πx h sin πy h , ( x, y ) ∈ (0 , h ) × (0 , h ) . We depict the phases at six different times in Fig. 4.2. Note that typical phasetransition phenomena can be clearly observed from these pictures. Moreover, we findthat under a small perturbation( u is small near the origin), the spinodal decompo-sition occurs and then coarsens, and after a period of evolution, the two coexistingphases become stable. Note that, compared to the subsequent motion, the initial sep-aration occurs over a very small time scale. Moreover, the evolution of the energies,including bulk energy and interfacial energy, is shown in Fig. 4.3a, the developmentof the mass is displayed in Fig. 4.3b, the maximum-norm of the numerical solutionis illustrated in Fig. 4.3c. Apparently, the presented method almost preserves the a) t=0 (b) t=0.01 (c) t=0.1(d) t=0.5 (e) t=1 (f) t=10 Fig. 4.2:
Example 4, spinoidal decomposition at six fixed time. (a) (b) (c)
Fig. 4.3:
Example 4 (a):
Energies evolution , (b):
Mass evolution , (c):
Evolution of themaximum-norm of the solution. properties of energy dissipation and mass conservation, while the numerical solutionitself is uniformly bounded.
Example 4:
We consider the Cahn-Hilliard equation suggested in [2] where theinitial date u is a random value field which is uniformly distributed between − (cid:15) is set to be 0 .
02. We depict the phase evolution of Example 4 inFig.4.4. The process of phase evolution is similar to that in Example 3. That is, thespinodal decomposition takes place very early, and after a brief period of evolution,the separation becomes very slow. Fig. 4.4 is also in good agreement with the onepresented in [2] by using the C virtual element method. The discrete energies andmass are shown in Fig. 4.5a and Fig. 4.5b. The maximum norm of the approximatesolution is displayed in Fig. 4.5c. These numerical results reveal that our numericalscheme is energy stable. a) t=0 (b) t=0.01 (c) t=0.1(d) t=0.5 (e) t=1 (f) t=10 Fig. 4.4:
Example 4, spinoidal decomposition at six fixed time. (a) (b) (c)
Fig. 4.5:
Example 4 (a):
Energies evolution , (b):
Mass evolution , (c):
Evolution of themaximum-norm of the solution.
In this subsection, we focus on tracking the evo-lution of initial data’s interfaces, including a cross-shaped, an elliptic-shaped and twocircles-shaped interfaces between phases. In all examples, we use the uniform meshwith mesh size h = 1 /
128 and time step size ∆ t = 5 × − . Example 5:
We consider the Cahn-Hilliard equation suggested in [9] with (cid:15) = 0 .
01 and the initial value u ( x, y ) = . , if 5 | ( y − . − ( x − . | + | ( x − . − ( y − . | < , . , if 5 | ( x − . − ( y − . | + | ( y − . − ( x − . | < , − . , otherwise . From Fig. 4.6, we observe that the cross-shaped interface evolves toward a steadycircular interface. From Fig. 4.7, we observe that the mass is well preserved, theenergy is dissipative and the maximum norm of the approximate solution is controlled. a) t=0 (b) t=0.005 (c) t=0.01(d) t=0.05 (e) t=0.1 (f) t=1 Fig. 4.6:
Example 5
Evolution of a cross-shaped interface at six temporal frames.
Comparing Fig. 4.6 with the numerical results given in [9] computed by the mixedFEM, they are in good agreement. (a) (b) (c)
Fig. 4.7:
Example 5 (a)
Energies evolution (b)
Mass evolution , (c)
Evolution of the maximumnorm of the approximate solution.
Example 6:
We consider the Cahn-Hilliard equation with a piecewise constantinitial data u whose jump set has a shape of an ellipse: u ( x, y ) = (cid:40) . , if 81( x − . + 9( y − . < , − . , otherwise . and the parameter is set to be (cid:15) = 0 . a) t=0 (b) t=0.003 (c) t=0.05(d) t=0.1 (e) t=0.3 (f) t=1 Fig. 4.8:
Example 6
Evolution of a cross-shaped interface at six temporal frames. (a) (b) (c)
Fig. 4.9:
Example 6 (a)
Energies evolution (b)
Mass evolution , (c)
Evolution of the maximumnorm of the approximate solution.
Example 7:
We consider the Cahn-Hilliard equation with (cid:15) = 0 .
025 on thedomain [ − , and the initial value u ( x ) = tanh( 1 √ (cid:15) min { (cid:112) ( x + 0 . + y − . , (cid:112) ( x − . + y − . } ) . We generate six snapshots at six fixed time points in Fig. 4.10. This graphclearly indicates that the two circle interfaces gradually evolve into one circle, whichis consistent with the maximum-norm results obtained in [21]. Numerical resultsdepicting the mass, energies and solution’s evolution are shown in Fig. 4.11.
5. Conclusion.
We designed a C linear finite element method to solve theCahn-Hilliard equations. This method has a minimum total degree of freedoms and isvery simple in implementation. A series of numerical examples indicate that the new a) t=0 (b) t=0.001 (c) t=0.005(d) t=0.01 (e) t=0.05 (f) t=0.1 Fig. 4.10:
Example 7
Evolution of a cross-shaped interface at six temporal frames. (a) (b) (c)
Fig. 4.11:
Example 7 (a)
Energies evolution (b)
Mass evolution , (c)
Evolution of the maxi-mum norm of the approximate solution. method is stable, efficient and is able to capture some important physical features suchas energy decay and mass conservation during the phase evolution process governedby the Cahn-Hilliard equations. Meanwhile, the numerical results reveal that ournovel method has the optimal convergence orders.Ongoing research topics include a theoretical analysis of the proposed methodand an extension of the presented method to 3D Cahn-Hilliard equations and/orother high order differential equations. For 3D Cahn-Hilliard equation, we adoptsimilar stabilization method [31] for time discretization.
REFERENCES171]
A. Agouzal and Y. Vassilevski , On a discrete Hessian recovery for P finite elements , J.Numer. Math., 10 (2002), pp. 1–12.[2] P. F. Antonietti, L. Beir˜ao da Veiga, S. Scacchi, and M. Verani , A C virtual elementmethod for the Cahn-Hilliard equation with polygonal meshes , SIAM J. Numer. Anal., 54(2016), pp. 34–56.[3] I. Babuˇska , The finite element method with Lagrangian multipliers , Numer. Math., 20(1972/73), pp. 179–192.[4]
V. E. Badalassi, H. D. Ceniceros, and S. Banerjee , Computation of multiphase systemswith phase field models , J. Comput. Phys., 190 (2003), pp. 371–397.[5]
Andrea L. Bertozzi, Selim Esedolu, and Alan Gillette , Inpainting of binary images usingthe Cahn-Hilliard equation , IEEE Trans. Image Process., 16 (2007), pp. 285–291.[6]
S. C. Brenner and L. R. Scott , The mathematical theory of finite element methods , vol. 15of Texts in Applied Mathematics, Springer, New York, third ed., 2008.[7]
John W Cahn and John E Hilliard , Free energy of a nonuniform system. i. interfacial freeenergy , The Journal of Chemical Physics, 28 (1958), p. 258.[8]
H. D. Ceniceros and A. M. Roma , A nonstiff, adaptive mesh refinement-based method forthe Cahn-Hilliard equation , J. Comput. Phys., 225 (2007), pp. 1849–1862.[9]
F. Chave, D. A. Di Pietro, F. Marche, and F. Pigeonneau , A hybrid high-order method forthe Cahn-Hilliard problem in mixed form , SIAM J. Numer. Anal., 54 (2016), pp. 1873–1898.[10]
H. Chen, H. Guo, Z. Zhang, and Q. Zou , A C linear finite element method for two fourth-order eigenvalue problems , IMA J. Numer. Anal., 37 (2017), pp. 2120–2138.[11] Y. Chen and J. Shen , Efficient, adaptive energy stable schemes for the incompressible Cahn-Hilliard Navier-Stokes phase-field models , J. Comput. Phys., 308 (2016), pp. 40–56.[12]
P. G. Ciarlet , The finite element method for elliptic problems , vol. 40 of Classics in AppliedMathematics, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA,2002. Reprint of the 1978 original [North-Holland, Amsterdam; MR0520174 (58
P. G. Ciarlet and P.-A. Raviart , A mixed finite element method for the biharmonic equation ,(1974), pp. 125–145. Publication No. 33.[14]
R. Courant , Variational methods for the solution of problems of equilibrium and vibrations ,Bull. Amer. Math. Soc., 49 (1943), pp. 1–23.[15]
Q. Du and R. A. Nicolaides , Numerical analysis of a continuum model of phase transition ,SIAM J. Numer. Anal., 28 (1991), pp. 1310–1322.[16]
C. M. Elliott and D. A. French , A nonconforming finite-element method for the two-dimensional Cahn-Hilliard equation , SIAM J. Numer. Anal., 26 (1989), pp. 884–903.[17]
C. M. Elliott, D. A. French, and F. A. Milner , A second order splitting method for theCahn-Hilliard equation , Numer. Math., 54 (1989), pp. 575–590.[18]
C. M. Elliott and S. Larsson , Error estimates with smooth and nonsmooth data for a finiteelement method for the Cahn-Hilliard equation , Math. Comp., 58 (1992), pp. 603–630,S33–S36.[19]
C. M. Elliott and S. Zheng , On the Cahn-Hilliard equation , Arch. Rational Mech. Anal., 96(1986), pp. 339–357.[20]
G. Engel, K. Garikipati, T. J. R. Hughes, M. G. Larson, L. Mazzei, and R. L. Taylor , Continuous/discontinuous finite element approximations of fourth-order elliptic problemsin structural and continuum mechanics with applications to thin beams and plates, andstrain gradient elasticity , Comput. Methods Appl. Mech. Engrg., 191 (2002), pp. 3669–3750.[21]
X. Feng, Y. Li, and Y. Xing , Analysis of mixed interior penalty discontinuous Galerkinmethods for the Cahn-Hilliard equation and the Hele-Shaw flow , SIAM J. Numer. Anal.,54 (2016), pp. 825–847.[22]
X. Feng and A. Prohl , Error analysis of a mixed finite element method for the Cahn-Hilliardequation , Numer. Math., 99 (2004), pp. 47–84.[23]
H. Guo, Z. Zhang, and R. Zhao , Hessian recovery for finite element methods , Math. Comp.,86 (2017), pp. 1671–1692.[24]
H. Guo, Z. Zhang, R. Zhao, and Q. Zou , Polynomial preserving recovery on boundary , J.Comput. Appl. Math., 307 (2016), pp. 119–133.[25]
H. Guo, Z. Zhang, and Q. Zou , A C linear finite element method for biharmonic problems ,J. Sci. Comput., 74 (2018), pp. 1397–1422.[26] , A C linear finite element method for sixth order elliptic equations , 2018.arXiv:1804.03793v2.[27] Y. He, Y. Liu, and T. Tang , On large time-stepping methods for the Cahn-Hilliard equation ,Appl. Numer. Math., 57 (2007), pp. 616–628.[28]
B. P. Lamichhane , A finite element method for a biharmonic equation based on gradient ecovery operators , BIT, 54 (2014), pp. 469–484.[29] R. J. LeVeque , Finite difference methods for ordinary and partial differential equations , Soci-ety for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2007. Steady-stateand time-dependent problems.[30]
D. Li and Z. Qiao , On second order semi-implicit Fourier spectral methods for 2D Cahn-Hilliard equations , J. Sci. Comput., 70 (2017), pp. 301–341.[31] ,
On the stabilization size of semi-implicit Fourier-spectral methods for 3D Cahn-Hilliardequations , Commun. Math. Sci., 15 (2017), pp. 1489–1506.[32]
D. Li, Z. Qiao, and T. Tang , Characterizing the stabilization size for semi-implicit Fourier-spectral method to phase field equations , SIAM J. Numer. Anal., 54 (2016), pp. 1653–1681.[33]
L. S. D MORLEY , The triangular equilibrium element in the solution of plate bending prob-lems , Aerosp. Quart., 19 (1968), pp. 149–169.[34]
B. Niceno , EasyMesh: A two-dimensional quality mesh generator . , 2001.[35] J. Nitsche , ¨uber ein Variationsprinzip zur L¨osung von Dirichlet-Problemen bei Verwendungvon Teilr¨aumen, die keinen Randbedingungen unterworfen sind , Abh. Math. Sem. Univ.Hamburg, 36 (1971), pp. 9–15. Collection of articles dedicated to Lothar Collatz on hissixtieth birthday.[36] J. T. Oden, A. Hawkins, and S. Prudhomme , General diffuse-interface theories and an ap-proach to predictive tumor growth modeling , Math. Models Methods Appl. Sci., 20 (2010),pp. 477–517.[37]
M. Picasso, F. Alauzet, H. Borouchaki, and P.-L. George , A numerical study of someHessian recovery techniques on isotropic and anisotropic meshes , SIAM J. Sci. Comput.,33 (2011), pp. 1058–1076.[38]
J. Shen and X. Yang , Numerical approximations of Allen-Cahn and Cahn-Hilliard equations ,Discrete Contin. Dyn. Syst., 28 (2010), pp. 1669–1691.[39]
M.-G. Vallet, C.-M. Manole, J. Dompierre, S. Dufour, and F. Guibault , Numericalcomparison of some Hessian recovery techniques , Internat. J. Numer. Methods Engrg., 72(2007), pp. 987–1007.[40]
G. N. Wells, E. Kuhl, and K. Garikipati , A discontinuous Galerkin method for the Cahn-Hilliard equation , J. Comput. Phys., 218 (2006), pp. 860–877.[41]
S. M. Wise, J. S. Lowengrub, H. B. Frieboes, and V. Cristini , Three-dimensional multi-species nonlinear tumor growth—I: Model and numerical method , J. Theoret. Biol., 253(2008), pp. 524–543.[42]
Y. Xia, Y. Xu, and C.-W. Shu , Local discontinuous Galerkin methods for the Cahn-Hilliardtype equations , J. Comput. Phys., 227 (2007), pp. 472–491.[43]
S. Zhang and M. Wang , A nonconforming finite element method for the Cahn-Hilliard equa-tion , J. Comput. Phys., 229 (2010), pp. 7361–7372.[44]
Z. Zhang and A. Naga , A new finite element gradient recovery method: superconvergenceproperty , SIAM J. Sci. Comput., 26 (2005), pp. 1192–1213 (electronic).[45]
O. C. Zienkiewicz, R. L. Taylor, and J. Z. Zhu , The finite element method: its basis andfundamentals , Elsevier/Butterworth Heinemann, Amsterdam, seventh ed., 2013.[46]