On linear schemes for a Cahn–Hilliard diffuse interface model

Journal of Computational Physics 234 (2013) 140–171

Contents lists available at SciVerse ScienceDirect

Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp

On linear schemes for a Cahn–Hilliard diffuse interface model F. Guillén-González ⇑, G. Tierra Departamento de Ecuaciones Diferenciales y Análisis Numérico, Universidad de Sevilla, Sevilla, Spain

a r t i c l e

i n f o

Article history: Received 16 September 2011 Received in revised form 14 September 2012 Accepted 15 September 2012 Available online 5 October 2012 Keywords: Diffuse interface phase-field Cahn–Hilliard Time integration Mixed finite element Long time stability Equilibrium solutions

a b s t r a c t Numerical schemes to approximate the Cahn–Hilliard equation have been widely studied in recent times due to its connection with many physically motivated problems. In this work we propose two type of linear schemes based on different ways to approximate the double-well potential term. The first idea developed in the paper allows us to design a linear numerical scheme which is optimal from the numerical dissipation point of view meanwhile the second one allows us to design unconditionally energy-stable linear schemes (for a modified energy). We present first and second order in time linear schemes to approximate the CH problem, detailing their advantages over other linear schemes that have been previously introduced in the literature. Furthermore, we compare all the schemes through several computational experiments. Ó 2012 Elsevier Inc. All rights reserved.

1. Introduction The study of interfacial dynamics in the mixture of different fluids, solids, or gas is one of the fundamental issues in hydrodynamics and materials science due it plays an increasingly important role in many current scientific, engineering, and industrial applications. Traditional fluid dynamics treats these mixtures as sharp interfaces on which a set of interfacial balance conditions must be imposed. It is delicate and difficult to develop a proper notion of generalized solutions of these models and it becomes even more challenging to compute the generalized solutions when they are defined, leading to an almost intractable theoretical problem. A conceptually straightforward way of handling the moving interface is to employ a mesh that has grid points on the interfaces, and deforms in time according to the flow on both sides of the boundary. However, from a computational point of view, these problems are difficult to solve because during the evolution the fluid interface may experience topological changes such as self-intersection, pinch-off, splitting, and fattening, and in such a situation the sharp interface formulation breaks down. As an alternative, fixed-grid methods that regularize the interface have been highly successful in treating deforming interfaces. These include the volume-of-fluid (VOF) method, the front-tracking method and the level-set method (see [8] for a review of these type of schemes). Instead of formulating the flow of two domains separated by an interface, these methods represent the interfacial tension as a body force or bulk stress spread over a narrow region covering the interface. Then a single set of governing equations can be written over the entire domain, and solved on a fixed grid in a purely Eulerian framework. An alternative approach for solving interface problems is the diffuse interface theory, which was originally developed as methodology for modeling and approximating solid–liquid phase transitions in which the effects of surface tension and non-equilibrium thermodynamic behavior may be important at the surface. The diffuse interface model describes the interface by a mixing energy. This idea can be traced to van der Waals [29], and is the foundation for the phase-field theory for phase transition and critical phenomena. Thus, the structure of the interface is determined by molecular forces; the tenden⇑ Corresponding author. E-mail addresses: [email protected] (F. Guillén-González), [email protected] (G. Tierra). 0021-9991/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jcp.2012.09.020

F. Guillén-González, G. Tierra / Journal of Computational Physics 234 (2013) 140–171

141

cies for mixing and de-mixing are balanced through the non-local mixing energy. In the theory, the interface is represented as a layer of small thickness. The method uses an auxiliary function / (called phase-field function) to localize the phases. The phase-field function assumes distinct values in the bulk phases (for instance / ¼ 1 in a phase and / ¼ 1 in the other one) away from the interfacial regions over which the phase function varies smoothly. Moreover, the diffuse interface models converge to their corresponding sharp interface models as the width of the interfacial layer e tends to zero. The phase-field method is also a level-set method except that the surface motion can be viewed as due to the physical e ð/Þ energy dissipation /t ¼  dEd/ where Ee ð/Þ is the phase-field’s free energy functional, which depends on the interface thickness e. Although resolution of the interface for small e becomes difficult, the phase field motion is dictated by a bulk field over the whole domain. Therefore, it inherits all of the aforementioned qualities; ease of coupling with a fluid, indifference to morphological singularities in the interface, physical dissipation, and a spatially uniform discretization. In the context of the surface elasticity, the value Ee ð/Þ can represent different interfacial energies associated to the phase field. The most basic energy functional one may introduce is

Ee ð/Þ ¼

 Z  1 jr/j2 þ Fð/Þ dx; X 2

ð1:1Þ

where Fð/Þ is a double-well potential (competition between a convex and a concave part). It is known that the Allen–Cahn equation

/t ¼ D/  f ð/Þ is a gradient flow for the energy functional Ee ð/Þ where f ð/Þ ¼ F 0 ð/Þ and the elliptic operator LAC ð/Þ :¼ D/ þ f ð/Þ is the representation of the Fréchet derivative E0e ð/Þ in the space L2 ðXÞ. Another gradient flow for the same Liapunov energy functional is the Cahn–Hilliard equation

/t ¼ DðD/  f ð/ÞÞ; which was originally introduced by Cahn and Hilliard [7] to describe the phase separation and coarsening phenomena in a melted alloy that is quenched to a temperature at which only two different concentration phases can exist stably. For more physical background, derivation and discussion of the Cahn–Hilliard equation we refer to [4,27]. It is also known that the elliptic operator LCH ð/Þ :¼ DðD/  f ð/ÞÞ associated to the Cahn–Hilliard equation is the representation of the Fréchet derivative of Ee ð/Þ in the space H1 ðXÞ. The two boundary conditions that we will consider, the outward normal derivatives of / and D/  f ð/Þ vanish on @ X, imply that none of the mixture can pass through the walls of the container X which is the most natural way to ensure that the total ‘‘free energy’’ of the mixture decreases in time, which is required by thermodynamics, when there is no interaction between the alloy and the containing walls. The evolution of the concentration consists of two stages: the first stage (rapid in time) is known as phase separation and the second (slow in time) is known as phase coarsening. At the end of the first stage, fine-scaled phase regions are formed, which are separated by the interface. At the end of the second stage, the solution will generically tend to a equilibrium state, which minimizes the energy functional, at least for the trajectory of the solution. Due to the existence of nonlinearity, numerical approximations of the Cahn–Hilliard problem become a crucial mean for understanding the phase transition of the isothermal binary alloy. The primary numerical challenge for solving the Cahn– Hilliard equation results from the presence of the parameter e in the equation, which usually is small in phase transition applications. Then, an appropriate numerical resolution of the problem requires proper relation of numerical scales, that is, the (spatial) mesh size h and the (time) step size k have to properly relate to the interaction length e. In the last years, many authors have been working in the numerical approximation of Cahn–Hilliard problems, coupled Navier–Stokes–Cahn–Hilliard and in many other related problems. These models usually deal with two fluids, but it is possible to extend them to a case with three components as in [6]. Some authors focus their attention in the case of a logarithmic potential, where implicit approximations of the potential are considered [9,20,21]. In the polynomial potential case, several approximations have been considered. On one hand, implicit approximations of the potential are made in [11–13,18], where a Newton method is usually employed in order to compute the nonlinear scheme. There is also an implicit-explicit approximation of the potential that does not introduce any numerical dissipation which have been widely used in phase field models and in the liquid crystal context [10,19,23,26]. On the other hand, following the ideas introduced by Eyre in [16], many authors split the potential into a convex and a non-convex part in order to assure the existence of some numerical dissipation to obtain a unconditional stable scheme [5,17,22,30]. Another approach that introduces some numerical dissipation to derive unconditionally stable schemes has been considered in [25]. The resulting schemes are again nonlinear. Trying to circumvent the problem of the nonlinearity of the schemes and following Eyre’s ideas, some linear (conditionally energy-stable) schemes have been studied in [28] (which can be unconditionally energy-stable introducing large enough numerical dissipation and truncating the potential term) and the same type of schemes have been considered in [31], where a nonconforming finite element approximation is presented. In this work we study these already known linear schemes and we propose two new ways of dealing in a linear way with the potential term. One of them is optimal from the numerical dissipation point of view meanwhile with the second one it is possible to design unconditionally stable schemes. Moreover, we use these potential approximations to design two new first order in time schemes and two new second order ones.

142

F. Guillén-González, G. Tierra / Journal of Computational Physics 234 (2013) 140–171

The rest of this paper is organized as follows. First of all, we introduce the Cahn–Hilliard problem in Section 2. In Section 3 we present a generic time-discrete scheme, and in Sections 4 and 5 we develop new first and second order in time linear schemes, respectively. Section 6 is devoted to compare numerically the schemes presented through the paper by means of several numerical experiments. Finally, we summarize the main achievements of the paper in Section 7. 2. Cahn–Hilliard diffuse interface model The Cahn–Hilliard equation is a gradient flow of the total free energy functional

1 2

Eð/Þ ¼ Ephilic ð/Þ þ Ephobic ð/Þ :¼

Z

jr/j2 dx þ

Z

X

Fð/Þdx; X

where Ephilic and Ephobic are the philic and phobic energies, respectively. To derive the Cahn–Hilliard equation, the starting point is the following mass balance law,

/t þ cr  J ¼ 0; where c > 0 is a relaxation time constant and J is the phase flux defined as

  @Eð/Þ ; J ¼ Mð/Þr @/ with Mð/Þ the so-called mobility function. From now on, we are going to consider a constant value for the mobility M ¼ 1 and the Ginzburg–Landau double well potential for the phobic free energy,

Fð/Þ ¼

1 ð/2  1Þ2 4e2

with

f ð/Þ ¼ F 0 ð/Þ ¼

1

e2

ð/2  1Þ/:

Under these assumptions, the Cahn–Hilliard problem reads,

8 > < /t ¼ cDðD/ þ f ð/ÞÞ in X  ð0; TÞ; ð/ÞÞ @/ ¼ 0; @ðD/þf ¼0 on @ X  ð0; TÞ; @n @n > : /jt¼0 ¼ /0 in X: e ~tÞ ¼ /ðtÞ verifies the same probIn the following, we take c ¼ 1 because changing t by ~t ¼ tc, the corresponding unknown /ð lem with c ¼ 1. Our task is to develop linear numerical schemes for the Cahn–Hilliard problem, suitable to use a continuous finite element approximation and taking into account the Definition 3.4 about energy-stability. In order to avoid the problems related to the fourth order derivative, it is usual to rewrite the problem using an auxiliary unknown (see for instance [10,14])

8 in X  ð0; TÞ; > < /t ¼ Dw; w ¼ D/ þ f ð/Þ @/ @w ¼ 0; ¼ 0 on @ X  ð0; TÞ; @n @n > : /jt¼0 ¼ /0 in X:

ðPÞ

The weak formulation of problem (P) is defined as follows: find ð/; wÞ such that

/ 2 L1 ð0; T; H1 ðXÞÞ;

/t 2 L2 ð0; T; ðH1 ðXÞÞ0 Þ and w 2 L2 ð0; T; H1 ðXÞÞ

and satisfying the following variational formulation

( ðVPÞ

 þ ðrw; rwÞ  ¼ 0 8w  2 H1 ðXÞ; h/t ; wi    2 H1 ðXÞ:  ðw; /Þ ¼ ðr/; r/Þ þ ðf ð/Þ; /Þ 8/

R  the duality product between Hereafter, ðf ; gÞ :¼ X f ðxÞgðxÞdx denotes the scalar product in the L2 ðXÞ-space and h/t ; wi R 0 1 1 ðH ðXÞÞ and H ðXÞ. This problem is conservative, because the total mass X /ðtÞ remains constant in time, as we can realize  ¼ 1 in (VP), taking w

d dt

Z

/ ¼ 0;

i:e:

X

Z X

/ðtÞ ¼

Z /0

8t P 0:

ð2:2Þ

X

 ¼ / in (VP), we derive the energy law  ¼ w and / Moreover, taking w t

d Eð/ðtÞÞ þ dt

Z X

jrwðtÞj2 dx ¼ 0:

ð2:3Þ

143

F. Guillén-González, G. Tierra / Journal of Computational Physics 234 (2013) 140–171

That implies that the energy decreases in time with a rate given by the physical dissipation krwðtÞk2L2 . From these properties is possible to deduce the existence and uniqueness of global in time weak solution ð/; wÞ of (P) [1,3,15,18,24]. 3. A generic implicit-explicit scheme Let consider a fixed a partition of the time interval ½0; T given by t n ¼ nk where k > 0 denotes the time step (that we take constant for simplicity). We are going to introduce some first and second order (in time) linear schemes, using an implicit1 explicit approximation fk ð/ðt nþ1 Þ; /ðt n ÞÞ of the potential part, as follows: given /n 2 H1 ðXÞ, compute ð/nþ1 ; wnþa Þ 2 1 1 H ðXÞ  H ðXÞ solving

(

 þ ðrwnþa1 ; rwÞ  ¼ 0 8w  2 H1 ðXÞ; ðdt /nþ1 ; wÞ 1 1  þ ðfk ð/nþ1 ; /n Þ; /Þ  2 H1 ðXÞ;  ¼ ðr/nþa ; r/Þ  8/ ðwnþa ; /Þ 1

1

nþ1

1

ð3:4Þ

n

nþ1

n

where /nþa ¼ /nþ1 (if a ¼ 1) or /nþa ¼ /nþ2 :¼ / 2þ/ (if a ¼ 2). Hereafter, we denote dt /nþ1 ¼ / k/ . At the initial time, we fix /0 ¼ /ð0Þ ¼ /0 in X. We are going to introduce different linear schemes, although we can see all of them as different ways of dealing with the potential term fk ð/nþ1 ; /n Þ. Note that all the schemes defined using (3.4) are conservative, because the total mass remains  ¼ 1 in (3.4), constant in time, as we can check taking w

Z X

/nþ1 ¼

Z

/n ¼    ¼

X

Z

/0 :

X

3.1. Energy law, numerical dissipation and energy-stability In order to define the concept of energy-stable schemes, we consider a discrete energy law and we need to define some numerical residual terms. Definition 3.1. The philic numerical dissipation of scheme (3.4) is defined by

NDphilic ð/; wÞ ¼





r/ þ ða  1Þrw r/  rw

a

;

k



1 2k

Z

ðjr/j2  jrwj2 Þ;

8/; w 2 H1 ðXÞ:

X

The phobic numerical dissipation depends on the choice of the approximation fk ð/; wÞ and is defined as

NDphobic ð/; wÞ ¼

1 k

Z

fk ð/; wÞð/  wÞ 

X

1 k

Z

ðFð/Þ  FðwÞÞ:

X

Finally, the total numerical dissipation is defined as

NDð/; wÞ :¼ NDphilic ð/; wÞ þ NDphobic ð/; wÞ: After the previous definition, in order to derive a discrete energy law from the numerical scheme (3.4), we just test by  ¼ dt /nþ1 in (3.4), arriving at  ¼ wnþa1 , / w 1

dt Eð/nþ1 Þ þ krwnþa k2L2 ðXÞ þ NDð/nþ1 ; /n Þ ¼ 0:

ð3:5Þ

Remark 3.2. In the case NDð/nþ1 ; /n Þ P 0 we say that the scheme introduces numerical dissipation (in the time step n), meanwhile when NDð/nþ1 ; /n Þ 6 0 we say that a numerical source appears. Remark 3.3. The term NDphilic ð/; wÞ depends on the value of a. In fact, when a ¼ 1, NDphilic ð/; wÞ ¼ 12 krð/  wÞk2L2 P 0 and when a ¼ 2, NDphilic ð/; wÞ ¼ 0. 1

Definition 3.4. The numerical scheme (3.4) is energy-stable if for any ð/nþ1 ; wnþa Þ solution of (3.4), the following relation holds 1

dt Eð/nþ1 Þ þ krwnþa k2L2 ðXÞ 6 0;

8n:

In particular, energy-stable schemes satisfy the energy decreasing in time property,

Eð/nþ1 Þ 6 Eð/n Þ;

8n:

144

F. Guillén-González, G. Tierra / Journal of Computational Physics 234 (2013) 140–171

This energy-stability criterion is more restrictive than the energy decreasing property (presented in some previous works), although we are not going to be able to prove that the linear schemes presented in this paper satisfy any of these criteria. However, we consider the more restrictive form of energy-stability because it is important that the decay of the energy is related to krwk2L2 , splitting the physical dissipation from the numerical perturbations. 3.2. Consistency errors Let ð/; wÞ denote the unique solution of (P). In order to bound the consistency errors, we will assume ð/; wÞ regular enough. In fact, we will impose the following regularity assumptions (see Lemma 3.8 below): 2

 If a ¼ 1, /t 2 CðH1 ðXÞÞ, /tt 2 CðL2 ðXÞÞ and dtd 2 Eð/ðtÞÞ 2 C½0; T. 3  If a ¼ 2, /ttt 2 CðL2 ðXÞÞ and dtd 3 Eð/ðtÞÞ 2 C½0; T. Definition 3.5. The consistency errors associated to the generic scheme (3.4) are defined as:

nnþ1 :¼ dt /ðtnþ1 Þ  Dwðtnþ1a Þ; w nnþ1 /

ð3:6Þ

8 < D/ðt nþ1 Þ þ fk ð/ðt nþ1 Þ; /ðt n ÞÞ if a ¼ 1;   :¼ : D /ðtnþ1 Þþ/ðtn Þ þ f ð/ðt Þ; /ðt ÞÞ if a ¼ 2; nþ1 n k 2

ð3:7Þ

where dt /ðtnþ1 Þ ¼ /ðtnþ1 kÞ/ðtn Þ and tnþa1 ¼ t n þ 1a k. Furthermore, the consistency error with respect to the discrete energy law (3.5) will be defined as

E nþ1 :¼ dt Eð/ðt nþ1 ÞÞ þ krwðtnþ1a Þk2L2 ðXÞ : Remark 3.6. Testing (3.6) by wðtnþa1 Þ and (3.7) by dt /ðtnþ1 Þ, equality (3.5) for the exact solution can be rewritten as nþ1 E nþ1 þ NDphilic ð/ðt nþ1 Þ; /ðt n ÞÞ þ NDphobic ð/ðt nþ1 Þ; /ðt n ÞÞ ¼ ðnnþ1 w ; wðt nþ1a ÞÞ þ ðn/ ; dt /ðt nþ1 ÞÞ:

ð3:8Þ

Remark 3.7. From Remark 3.3, if a ¼ 1 then

NDphilic ð/ðtnþ1 Þ; /ðtn ÞÞ ¼

k 2 krdt /ðt nþ1 ÞkL2 ðXÞ : 2

ð3:9Þ

Taking into account the exact problem (P) at t ¼ t nþa1 , consistency errors can be rewritten as follows:

nnþ1 ¼ dt /ðtnþ1 Þ  /t ðtnþ1a Þ; w     /ðt nþ1 Þ þ /ðt n Þ nþ1 þ Dð/ðtnþ1 ÞÞ þ fk ð/ðt nþ1 Þ; /ðt n ÞÞ  f ð/ðt nþa1 ÞÞ: n/ ¼ ða  1Þ D 2 2

ð3:10Þ ð3:11Þ

We are going to see the order with respect to k of each term of (3.8). Lemma 3.8 (a) NDphilic ð/ðt nþ1 Þ; /ðt n ÞÞ



6 B1 k if a ¼ 1;

ð3:12Þ

¼ 0 if a ¼ 2;

where B1 ¼ B1 ðkr/t kCð½0;T;L2 ðXÞÞ Þ. (b) nþ1 a

knw kL2 ðXÞ 6 C a k

a ¼ 1; 2;

ð3:13Þ

with C 1 ¼ C 1 ðk/tt kCð½0;T;L2 ðXÞÞ Þ and C 2 ¼ C 2 ðk/ttt kCð½0;T;L2 ðXÞÞ Þ. (c) Assume fk ð/ðtnþ1 Þ; /ðtn ÞÞ defined such that a

kfk ð/ðt nþ1 Þ; /ðt n ÞÞ  f ð/ðtnþ1a ÞÞkL2 ðXÞ 6 Da k

ð3:14Þ

with D1 ; D2 > 0 (depending on the regularity of /). Then

8 < knnþ1 / kL2 ðXÞ 6 D1 k : knnþ1 k /

L2 ðXÞ

if a ¼ 1; 2

6 D2 k þ

k2 8

kDð/tt ÞkCð½0;T;L2 ðXÞÞ

if a ¼ 2;

ð3:15Þ

145

F. Guillén-González, G. Tierra / Journal of Computational Physics 234 (2013) 140–171

(d)

jE nþ1 j 6 Ea k

a

a ¼ 1; 2;

 2 with E1 ¼ E1 dtd 2 Eð/ðtÞÞ

C½0;T

 3 and E2 ¼ E2 dtd 3 Eð/ðtÞÞ

C½0;T

ð3:16Þ

 .

(e) Assuming (3.14), one has

jNDphobic ð/ðt nþ1 Þ; /ðt n ÞÞj 6 F a k

a

a ¼ 1; 2;

ð3:17Þ

where F a is a constant depending on Ba , C a , Da and Ea .

Proof. The estimates presented in (3.12), (3.13) and (3.15) can be deduced from the expressions given in (3.9), (3.10) and and nnþ1 respectively, and the following Taylor’s formulas: (3.11) of NDphilic ð/ðt nþ1 Þ; /ðt n ÞÞ, nnþ1 w /

dt /ðt nþ1 Þ ¼ /t ðtnþn Þ;

n 2 ð0; 1Þ; (1

dt /ðt nþ1 Þ  /t ðt nþa1 Þ ¼

2

if a ¼ 1;

n 2 ð0; 1Þ

/tt ðt nþn Þk;

1 ð/ttt ðtnþn Þ 6

2

 /ttt ðt nþg ÞÞk ;

 2 /ðt nþ1 Þ þ /ðtn Þ 1 k  /ðt nþ1 Þ ¼ /tt ðtnþn Þ ; 2 2 2 2

n; g 2 ð0; 1Þ if a ¼ 2;

n 2 ð0; 1Þ:

On the other hand, (3.16) follows from taking advantage of the energy law (2.3) (at t ¼ tnþ1 for a ¼ 1 or at t ¼ tnþ1=2 for a ¼ 2). Finally, we can bound NDphobic ð/ðtnþ1 Þ; /ðtn ÞÞ from (3.8), using (3.12), (3.13), (3.15) and (3.16):

8 < ðB1 þ E1 Þk þ C 1 kkwðt nþ1 ÞkL2 ðXÞ þ D1 kkdt /ðt nþ1 ÞkL2 ðXÞ if a ¼ 1;   jNDphobic ð/ðtnþ1 Þ; /ðt n ÞÞj 6 : E2 k2 þ C 2 k2 kwðt nþ1 ÞkL2 ðXÞ þ D1 þ 1 kD/tt kjCð½0;T;L2 ðXÞÞ k2 kdt /ðt nþ1 ÞkL2 ðXÞ 8 2

if a ¼ 2;

hence (3.17) is deduced because

kwðtnþ1a ÞkL2 ðXÞ 6 kwkCð½0;T;L2 ðXÞÞ

and kdt /ðtnþ1 ÞkL2 ðXÞ 6 k/t kCð½0;T;L2 ðXÞÞ :



Remark 3.9. In order to define first ða ¼ 1Þ and second ða ¼ 2Þ order in time schemes by using (3.4), it is sufficient to choose fk ð/; wÞ such that estimate (3.14) holds. 4. First order schemes We are going to consider first order semi-implicit schemes, which are implicit for the philic part (i.e. taking a ¼ 1 in (3.4)), meanwhile for the potential term we will use some semi-implicit approximations,

(

  þ ðrwnþ1 ; rwÞ  ¼ 0 8w  2 H1 ðXÞ; dt /nþ1 ; w  þ ðfk ð/nþ1 ; /n Þ; /Þ  2 H1 ðXÞ:  ¼ ðr/nþ1 ; r/Þ  8/ ðwnþ1 ; /Þ

Recalling Remark 3.3,

NDphilic ð/nþ1 ; /n Þ P 0 and NDphilic ð/ðt nþ1 Þ; /ðtn ÞÞ 6 B1 k: The first kind of schemes that we are going to derive are based on the ideas introduced by Eyre in his celebrated preprint [16], which consists in separating the potential term between the convex and non-convex part in order to assure the introduction of some numerical dissipation and hence the unconditional energy-stability of the scheme is deduced. Secondly, a linear scheme with a phobic numerical dissipation term NDphobic of second order in time will be considered. Finally, we will present an unconditionally energy-stable linear scheme based on a Lagrange multiplier. This idea has been previously introduced in [2] to deal with nematic liquid crystal flows. 4.1. Eyre’s decomposition (nonlinear scheme) The key point in Eyre’s work is to approximate the potential term introducing phobic numerical dissipation in the discrete energy law (3.5), splitting the potential term in the following way

Fð/Þ ¼ F c ð/Þ þ F e ð/Þ with F 00c ð/Þ P 0 and F 00e ð/Þ 6 0 8/ 2 R; 4

2

2

ð4:18Þ 2

for instance, one can take F c ð/Þ ¼ ð/ þ 1Þ=ð4e Þ and F e ð/Þ ¼ / =ð2e Þ.

146

F. Guillén-González, G. Tierra / Journal of Computational Physics 234 (2013) 140–171

Then, the (nonlinear) Eyre’s scheme take implicitly the convex term (or contractive) and implicitly the non-convex one (or expansive), i.e.

fk ð/nþ1 ; /n Þ ¼ fc ð/nþ1 Þ þ fe ð/n Þ; where fc ð/nþ1 Þ ¼ F 0c ð/nþ1 Þ and fe ð/n Þ ¼ F 0e ð/n Þ. By Taylor’s expansion, 00

nþn

F c ð/n Þ ¼ F c ð/nþ1 Þ þ F 0c ð/nþ1 Þð/n  /nþ1 Þ þ F c ð/2 F e ð/nþ1 Þ ¼ F e ð/n Þ þ F 0e ð/n Þð/nþ1  /n Þ þ

F 00e ð/nþg Þ 2

Þ

ð/n  /nþ1 Þ2 ;

ð/nþ1  /n Þ2 ;

for some n; g 2 ð0; 1Þ, where /nþn ¼ /n þ nð/nþ1  /n Þ ¼ n/nþ1 þ ð1  nÞ/n . Then, adding both equalities,

Fð/nþ1 Þ  Fð/n Þ ¼ ðfc ð/nþ1 Þ þ fe ð/n ÞÞð/nþ1  /n Þ 

F 00c ð/nþn Þ  F 00e ð/nþg Þ nþ1 ð/  /n Þ2 ; 2

hence

NDphobic ð/nþ1 ; /n Þ ¼

Z X

F 00c ð/nþn Þ  F 00e ð/nþg Þ nþ1 ð/  /n Þ2 ¼ k 2k

Z X

F 00c ð/nþn Þ  F 00e ð/nþg Þ ðdt /nþ1 Þ2 : 2

Therefore, NDphobic ð/nþ1 ; /n Þ P 0 (due to (4.18)) and

jNDphobic ð/ðt nþ1 Þ; /ðt n ÞÞj 6 F 1 k where F 1 ¼ F 1 ðk/kCð½0;T;L1 ðXÞÞ ; k/t kCð½0;T;L2 ðXÞÞ Þ: Therefore, the total numerical dissipation verifies

NDð/nþ1 ; /n Þ P 0 and NDð/ðt nþ1 Þ; /ðt n ÞÞ 6 ðB1 þ F 1 Þ k: Since this scheme is a particular case of (3.4) for a ¼ 1 and (3.14) holds for a ¼ 1, Lemma 3.8 yields to

jE nþ1 j 6 E1 k;

knnþ1 and knnþ1 / kL2 ðXÞ 6 D1 k w kL2 ðXÞ 6 C 1 k:

Consequently, Eyre’s scheme is nonlinear unconditionally energy-stable and first order accurate in time. 4.2. Eyre’s linear schemes: (E1), (E1mod) In order to obtain linear schemes following the ideas of Eyre [16], we rewrite the potential part of the free energy using a parameter b P 0, as follows:

Fð/Þ ¼

1 1 ð/2  1Þ2 ¼ 2 ð/4 þ 2b/2  2ðb þ 1Þ/2 þ 1Þ; 4e2 4e

and we consider the following separation of the potential,

F 1 ð/Þ ¼

1 b/2 2e2

and F 2 ð/Þ ¼

1 ð/4  2ðb þ 1Þ/2 þ 1Þ: 4e2

Thus, taking fk ð/nþ1 ; /n Þ ¼ F 01 ð/nþ1 Þ þ F 02 ð/n Þ, we arrive at the Eyre linear scheme

(

 þ ðrwnþ1 ; rwÞ  ¼ 0 8w  2 H1 ðXÞ; ðdt /nþ1 ; wÞ 3 nþ1 n nþ1  1  þ 2 ðð/ Þ  ðb þ 1Þ/n þ b/nþ1 ; /Þ  2 H1 ðXÞ:  8/ ðw ; /Þ ¼ ðr/ ; r/Þ e

ð4:19Þ

The phobic numerical dissipation can be bounded as

jNDphobic ð/ðt nþ1 Þ; /ðt n ÞÞj ¼

Z k 2 2 ð3/ðt Þ þ 2b þ 1Þðd /ðt ÞÞ nþn t nþ1 6 F 1 k; 2 e

ð4:20Þ

X

with n 2 ð0; 1Þ and F 1 ¼ F 1 ðk/kCð½0;T;L1 ðXÞÞ ; k/t kCð½0;T;L2 ðXÞÞ Þ. In particular, using the same arguments of Section 4.1

jE nþ1 j 6 E1 k;

knnþ1 and knnþ1 w kL2 ðXÞ 6 C 1 k / kL2 ðXÞ 6 D1 k:

The scheme considering b ¼ 2 (E1) was introduced in [16], and it satisfies NDphobic ð/nþ1 ; /n Þ P 0 if j/n j 6 5=3, whereas the scheme for b ¼ 1 (E1mod) verifies NDphobic ð/nþ1 ; /n Þ P 0 only if j/n j 6 1. But, it is not clear that always NDphobic ð/nþ1 ; /n Þ P 0 due the non-existence of a maximum principle for the Cahn–Hilliard model. Under this decomposition, condition (4.18) does not hold in general, only considering j/j2 6 ð2b þ 1Þ=3. In fact, the choice of b is a non-trivial task because it is related to the size of the phobic numerical dissipation. In the numerical simulations we are going to consider two cases: b ¼ 2 (E1) and b ¼ 1 (E1mod), because b ¼ 2 is the scheme introduced in [16] and taking b ¼ 1 the phobic numerical dissipation term is the lowest one whenever j/n j  1, due to (4.20).

F. Guillén-González, G. Tierra / Journal of Computational Physics 234 (2013) 140–171

147

Remark 4.1. In the last years, the introduction of an extra dissipative term eS2 ð/nþ1  /n Þ to Euler’s schemes and the idea of considering a truncated potential e F such that the second derivative of this truncated potential is bounded (j e F 00 ð/Þj 6 M for each / 2 R) have been widely used [28], obtaining energy-stable schemes for a choice of S depending on M. It is interesting to note that this kind of schemes could be seen as Eyre’s linear schemes. 4.3. Optimal dissipation scheme (OD1) Previously, we have seen that the phobic numerical dissipation terms (NDphobic ð/ðt nþ1 Þ; /ðt n ÞÞ) in the Eyre’s schemes are of order OðkÞ and non-negative (in the nonlinear case). Now, our idea is to define a linear scheme such that 2 NDphobic ð/ðtnþ1 Þ; /ðtn ÞÞbe of order Oðk Þ. The disadvantages of this approach will be mainly two: 1. It is not possible to control the sign of NDphobic ð/ðtnþ1 Þ; /ðt n ÞÞ, but at least formally, we are able to control this term for k small enough, because its order will be lower than the positive philic numerical dissipation one (since NDphilic ð/ðt nþ1 Þ; /ðt n ÞÞ P 0 and NDphilic ð/ðt nþ1 Þ; /ðt n ÞÞ ¼ OðkÞ). 2. The scheme will be only conditionally solvable, i.e., a condition arises to assure the existence and uniqueness of solution of the scheme. The optimal dissipation approximation of the potential can be derived by using the following Hermite quadrature formula (exact for P1 )

Z a

b

1 gðxÞdx ¼ ðb  aÞgðaÞ þ ðb  aÞ2 g 0 ðaÞ þ Cðb  aÞ3 g 00 ðnÞ 2

in the following sense:

Fð/nþ1 Þ  Fð/n Þ ¼

Z

/nþ1 n

/

1 f ð/Þd/ ¼ ð/nþ1  /n Þf ð/n Þ þ ð/nþ1  /n Þ2 f 0 ð/n Þ þ Cð/nþ1  /n Þ3 f 00 ð/nþf Þ; 2

where /nþf ¼ /n þ fð/nþ1  /n Þ with f 2 ð0; 1Þ. In this case, if we define

1 fk ð/nþ1 ; /n Þ ¼ f ð/n Þ þ ð/nþ1  /n Þf 0 ð/n Þ; 2

ð4:21Þ

then

NDphobic ð/nþ1 ; /n Þ ¼

C k

Z

ð/nþ1  /n Þ3 f 00 ð/nþf Þ ¼

X

C

e

2

k

2

Z

ðdt /nþ1 Þ3 /nþf ;

X

hence 2

jNDphobic ð/ðt nþ1 Þ; /ðt n ÞÞj 6 F 2 k ;

ð4:22Þ

where F 2 ¼ F 2 ðk/kCð½0;T;L1 ðXÞÞ ; k/t kCð½0;T;L3 ðXÞÞ Þ. Using this potential approximation we can not control the sign of NDphobic ð/nþ1 ; /n Þ even if we were able to assure that j/n j 6 1. On the other hand, one has (see Appendix B) 2

kfk ð/ðt nþ1 Þ; /ðtn ÞÞ  f ð/ðtnþ1 ÞÞkL2 ðXÞ 6 D2 k ;

ð4:23Þ

2

where D2 ¼ D2 ðk/kCð½0;T;L1 ðXÞÞ ; k/t kCð½0;T;L4 ðXÞÞ ; k/tt kCð½0;T;L2 ðXÞÞ Þ. Therefore, 2

2

kfk ð/ðt nþ1 Þ; /ðtn ÞÞ  f ð/ðtnþ1 ÞÞkL2 ðXÞ 6 D2 k þ kf ð/ðt nþ1 ÞÞ  f ð/ðt nþ1 ÞÞkL2 ðXÞ 6 D2 k þ D1 k; 2

where D1 ¼ D1 ðk/kCð½0;T;L1 ðXÞÞ ; k/t kCð½0;T;L2 ðXÞÞ Þ. Since we are using a backward Euler approximation for the philic part, from Lemma 3.8 one has

jE nþ1 j 6 E1 k;

2

knnþ1 and knnþ1 w kL2 ðXÞ 6 C 1 k / kL2 ðXÞ 6 D1 k þ D2 k :

By using f ð/Þ ¼ e12 ð/2  1Þ/ and f 0 ð/Þ ¼ e12 ð3/2  1Þ, we arrive at

fk ð/nþ1 ; /n Þ ¼

" # 1 3 n 2 nþ1 1 n 3 /nþ1 þ /n ð/ ð/ Þ /  Þ  : 2 e2 2 2

ð4:24Þ

Therefore, the linear scheme with second order phobic numerical dissipation is:

8  nþ1 n >  þ ðrwnþ1 ; rwÞ  ¼ 0; < / k/ ; w   >  þ 12 3 ð/n Þ2 /nþ1  1 ð/n Þ3  1 ð/nþ1 þ /n Þ; /  ¼ ðr/nþ1 ; r/Þ  ; : ðwnþ1 ; /Þ 2 2 e 2

ð4:25Þ

148

F. Guillén-González, G. Tierra / Journal of Computational Physics 234 (2013) 140–171

 2 H1 ðXÞ2 .  /Þ for each ðw; Remark 4.2. A similar idea has been recently used by Gomez and Hughes in [21] to obtain a second order unconditional stable (nonlinear) scheme using quadrature formulas of Hermite type, exact for P2 . Remark 4.3. As we can see in Appendix A, there is another way of deriving this linear second order in time approximation of the potential, by using a combination of quadrature formulas (exact for P0 ). In this way, we can observe that the approach of the potential (4.21) is optimal with respect to the order in time of the phobic numerical dissipation, because using quadrature formulas it is not possible to design a linear scheme with a higher order in time phobic numerical dissipation. Remark 4.4. Although the phobic numerical dissipation of the scheme (4.25) is second order in time, a rigorous proof of conditional energy-stability remains as an open problem. Lemma 4.5. There exists a unique ð/nþ1 ; wnþ1 Þ solution of (4.25) under the constraint k < 8e4 . nþ1 nþ1 nþ1 nþ1  /nþ1 and w ¼ wnþ1  wnþ1 Proof. Let ð/nþ1 1 ; w1 Þ and ð/2 ; w2 Þ be two solutions of (4.25) and denoting / ¼ /1 2 1 2 , we arrive at:

81  Þ þ ðrw; rwÞ  ¼ 0 8w  2 H1 ðXÞ; > < k ð/; w > :

ð4:26Þ

  þ 1 f 0 ð/n Þ/; /  2 H1 ðXÞ:   ¼ ðr/; r/Þ 8/ ðw; /Þ 2

 ¼ 1 /,  ¼ w and / Taking w k

1 1 krwk2L2 ðXÞ þ kr/k2L2 ðXÞ þ k 2k

Z

f 0 ð/n Þj/j2 dx ¼ 0;

X

hence

1 1 krwk2L2 ðXÞ þ kr/k2L2 ðXÞ þ k 2ke2

Z

3ð/n Þ2 j/j2 dx ¼ X

In order to control the right hand-side term

1 k

Z

/2 þ X

Z

1 2ke2

R X

1 2ke2

Z

j/j2 dx:

ð4:27Þ

X

 ¼ /, arriving at j/j2 dx, we take w

rw  r/ ¼ 0:

X

Therefore

1 2ke2

Z X

/2 ¼ 

1 2e2

Z X

1 2

rw  r/ 6 krwk2L2 ðXÞ þ

1 kr/k2L2 ðXÞ : 8e4

Using this bound in (4.27):

  1 1 1 krwk2L2 ðXÞ þ  4 kr/k2L2 ðXÞ 6 0: 2 k 8e Considering the assumption k < 8e4 , one has that / and w are constant functions. Finally, from (4.26)1 we derive / ¼ 0, and combining this with (4.26)2 we obtain w ¼ 0. h Remark 4.6. In order to assure the validity of Lemma 4.5 when a space discretization is considered (for instance, using a FE approximation), it is important to note that we can mimic the previous argument when the approximated discrete space for / is in included in the approximated discrete space for w. This constraint has been used to control k/k2L2 ðXÞ by means of  ¼ /). krwk2L2 ðXÞ and kr/k2L2 ðXÞ (taking w 4.4. Lagrange multiplier scheme (LM1) In this section, we present a new approach of the potential term where the basic idea consists on thinking on the double potential Fð/Þ as a penalization of the constraint j/j2 ¼ 1. Following the ideas presented in the nematic liquid crystal framework in [2], we can introduce qðx; tÞ 2 R that is not exactly a Lagrange multiplier but allows to enforce the sphere condition j/j ¼ 1, getting a linear unconditional stable and uniquely solvable scheme. Basically, the new formulation consists in taking advantage of the Ginzburg–Landau potential, writing it as

F. Guillén-González, G. Tierra / Journal of Computational Physics 234 (2013) 140–171

f ð/Þ ¼ q/ where q ¼

1

e2

149

ð/2  1Þ

taking the time derivative in the q-equation:

e2 qt ¼ 2//t jointly with the initial condition qð0Þ ¼ e12 ð/ð0Þ2  1Þ. The weak formulation of this problem is defined as follows: Find ð/; w; qÞ such that

/ 2 L1 ðð0; TÞ; H1 ðXÞ \ L1 ðXÞÞ; /t 2 L2 ð0; T; ðH1 ðXÞÞ0 Þ;

w 2 L2 ðð0; TÞ; H1 ðXÞÞ;

q 2 L2 ðð0; TÞ  XÞ;

qt 2 L2 ð0; T; ðH1 ðXÞ \ L1 ðXÞÞ0 Þ

and satisfying the following equations

8 1    > < h/t ; wi þ ðrw; rwÞ ¼ 0 8w 2 H ðXÞ;  þ ðq/; /Þ  8/  2 H1 ðXÞ;  ¼ ðr/; r/Þ ðw; /Þ > : e2 1  2 H ðXÞ \ L2 ðXÞ: i ¼ h/ /t ; q  i 8q hqt ; q 2 We present a semi-implicit scheme for this problem (considering the forward Euler approximation for the last equation): Given ð/n ; qn Þ 2 H1 ðXÞ  L2 ðXÞ, compute ð/nþ1 ; wnþ1 ; qnþ1 Þ 2 H1 ðXÞ  H1 ðXÞ  L2 ðXÞ solving

8  nþ1 n  / / >  þ ðrwnþ1 ; rwÞ  ¼ 0 8w  2 H1 ðXÞ; ;w > k > <  þ ðqnþ1 /n ; /Þ  2 H1 ðXÞ;  ¼ ðr/nþ1 ; r/Þ  8/ ðwnþ1 ; /Þ >     > nþ1 n > e2 qnþ1 qn :  ¼ /n / k/ ; q  ;q 8q 2 L2 ðXÞ; 2 k

ð4:28Þ

where we consider at the initial time q0 ¼ e12 ðð/0 Þ2  1Þ and /0 ¼ /ð0Þ in X. The scheme is conservative, because the total  ¼ 1, mass remains constant in time, as we can check taking w

Z

/nþ1 ¼

X

Z

/n ¼    ¼

X

Z

/0 :

X

Lemma 4.7. Scheme (4.28) is unconditionally energy-stable, in the sense that the following discrete energy law holds

e Eð/nþ1 ; qnþ1 Þ  e Eð/n ; qn Þ k ke2 þ krwnþ1 k2L2 ðXÞ þ krdt /nþ1 k2L2 ðXÞ þ kdt qnþ1 k2L2 ðXÞ ¼ 0; k 2 4

ð4:29Þ

where 2

1 e e Eð/; qÞ ¼ kr/k2L2 ðXÞ þ kqk2L2 ðXÞ ; 2 4 can be view as a modified version of the total free energy Eð/Þ.  ¼ dt /nþ1 and q  ¼ qnþ1 , we obtain unconditional long-time stability for  ¼ wnþ1 , / Proof. Taking in (4.28) as test function w this scheme in the sense that (4.29) holds. h Lemma 4.8. Scheme (4.28) has a unique solution. nþ1 nþ1 nþ1 nþ1 nþ1 nþ1 nþ1  /nþ1  wnþ1 Proof. Let ð/nþ1 1 ; w1 ; q1 Þ and ð/2 ; w2 ; q2 Þ be two solutions of (4.28) and denoting / ¼ /1 2 , w ¼ w1 2 nþ1 nþ1 and q ¼ q1  q2 , then:

Z

/¼0

X

and

8 1 1 >    > < k ð/; wÞ þ ðrw; rwÞ ¼ 0 8w 2 H ðXÞ;  þ ðq/n ; /Þ  2 H1 ðXÞ;  8/  ¼ ðr/; r/Þ ðw; /Þ > > : e2 ðq; q 2 n 1  2 L ðXÞ: Þ ¼ ð/ /; q  Þ 8q 2k

k

150

F. Guillén-González, G. Tierra / Journal of Computational Physics 234 (2013) 140–171

 ¼ 1 / and q  ¼ q,  ¼ w, / Taking w k

1 e2 krwk2L2 ðXÞ þ kr/k2L2 ðXÞ þ kqk2L2 ðXÞ ¼ 0; k 2k

hence the proof of the unicity of the scheme (4.28) holds. h Remark 4.9. In order to assure the validity of Lemma 4.8 when a space discretization is considered (for instance, using a FE approximation), it is important to note that we can mimic the previous argument with any conformed (finite-element) approximation. Remark 4.10. It is interesting to realize that in this scheme we can rewrite the last equation as:

Þ ¼ ðqnþ1 ; q

2

e2

Þ þ ðqn ; q Þ: ðð/nþ1 /n  ð/n Þ2 Þ; q

If we take into account the relation q0 ¼ e12 ðð/0 Þ2  1Þ, we can eliminate the qn unknown writing the scheme as

8  < /nþ1 /n ; w   ¼ 0 8w  2 H1 ðXÞ; þ ðrwnþ1 ; rwÞ k : nþ1   þ ðqnþ1 /n ; /Þ  2 H1 ðXÞ;  8/ ðw ; /Þ ¼ ðr/nþ1 ; r/Þ

where nþ1

q

¼

" n 2 X

e2

# i

iþ1

ð/ /

i 2

0

1

 ð/ Þ Þ þ / /



i¼1

1

e2

ðð/0 Þ2 þ 1Þ:

This formulation can be viewed as an unconditional stable and uniquely solvable linear scheme in the unknowns ð/n ; wn Þ. 5. Second order schemes The aim of this section is to modify the linear schemes that we have developed in the previous section in order to improve their accuracy in time preserving its stability properties. 5.1. Second order optimal dissipation scheme (OD2) The idea is to modify scheme (OD1) considering the approximation of the potential term presented in (4.21) and a Crank– Nicholson discretization of the philic term: 1 Given /n 2 H1 ðXÞ, compute ð/nþ1 ; wnþ2 Þ 2 H1 ðXÞ  H1 ðXÞ solving

8  nþ1 n  1 < / / ; w  ¼ ðrwnþ2 ; rwÞ   2 H1 ðXÞ; 8w k

 1 : nþ12   þ f ð/n Þ þ 1 f 0 ð/n Þð/nþ1  /n Þ; /  2 H1 ðXÞ;  8/ ðw ; /Þ ¼ ðr/nþ2 ; r/Þ 2 1

ð5:30Þ

1

where /nþ2 ¼ ð/nþ1 þ /n Þ=2. Notice that wnþ2 is computed directly as approximation of wðt nþ1 Þ. For this scheme, the philic 2

numerical dissipation vanish (NDphilic ð/nþ1 ; /n Þ ¼ 0) and the following discrete energy law holds, 1

dt Eð/nþ1 Þ þ krwnþ2 k2L2 ðXÞ þ NDphobic ð/nþ1 ; /n Þ ¼ 0; where we already know from (4.22) and (4.23)

jNDphobic ð/ðt nþ1 Þ; /ðt n ÞÞj 6 F 2 k

2

2

and kfk ð/ðt nþ1 Þ; /ðt n ÞÞ  f ð/ðt nþ1 ÞÞkL2 ðXÞ 6 D2 k ; 2

where D2 ¼ D2 ðk/kCð½0;T;L1 ðXÞÞ ; k/t kCð½0;T;L4 ðXÞÞ ; k/tt kCð½0;T;L2 ðXÞÞ Þ and F 2 ¼ F 2 ðk/kCð½0;T;L1 ðXÞÞ ; k/t kCð½0;T;L3 ðXÞÞ Þ. The drawback of this scheme is again that we can not control the sign of NDphobic ð/nþ1 ; /n Þ. Lemma 5.1. Scheme (5.30) has a unique solution under the constraint k < 4e4 . Proof. Following same arguments presented in Lemma 4.5. h 5.2. Variant of the second order optimal dissipation scheme (OD2mod) In this case, we want to introduce some philic numerical dissipation to avoid, at least partially, the possibility of computing increasing discrete energies. The idea consists in introducing a positive second order in time philic numerical dissipation such that

151

F. Guillén-González, G. Tierra / Journal of Computational Physics 234 (2013) 140–171

0 6 NDphilic ð/ðtnþ1 Þ; /ðtn ÞÞ 6 B2 k

2 2

2

nþ1 without losing the second order accuracy, that is knnþ1 w kL2 ðXÞ 6 C 2 k and kn/ kL2 ðXÞ 6 D2 k . For this purpose, we consider the diffusion term D/ at time tnþ1þh with h ¼ OðkÞ and h > 0: 2

(

 1  ¼ ðrwnþ2 ; rwÞ   2 H1 ðXÞ; ;w 8w

 1 1  þ f ð/n Þ þ 1 f 0 ð/n Þð/nþ1  /n Þ; /  2 H1 ðXÞ;  ¼ ðr/nþ2þh ; r/Þ  8/ ðwnþ2 ; /Þ 2 /nþ1 /n k

ð5:31Þ

where 1

/nþ2þh ¼ /n þ



 1 1 þ h ð/nþ1  /n Þ ¼ /nþ2 þ hð/nþ1  /n Þ: 2

Since h ¼ OðkÞ, we have 2

/ðt nþ1þh Þ ¼ /ðtnþ1 Þ þ OðhkÞ ¼ /ðt nþ1 Þ þ Oðk Þ: 2

2

2

Moreover, 1

ðr/nþ2þh ; rdt /nþ1 Þ ¼ dt

 Z  Z 1 jr/nþ1 j2 þ h k jrdt /nþ1 j2 ; 2 X X

hence

NDphilic ð/nþ1 ; /n Þ ¼ h k

Z

2

jrdt /nþ1 j2 P 0 and NDphilic ð/ðtnþ1 Þ; /ðtn ÞÞ 6 B2 k ;

X

where B2 ¼ B2 ðkr/t kCð½0;T;L2 ðXÞÞ Þ. In conclusion,

dt Eð/nþ1 Þ þ krwnþ1=2 k2L2 ðXÞ þ NDphilic ð/nþ1 ; /n Þ þ NDphobic ð/nþ1 ; /n Þ ¼ 0; 2

2

where jNDphobic ð/ðtnþ1 Þ; /ðt n ÞÞj 6 F 2 k and 0 6 NDphilic ð/ðtnþ1 Þ; /ðtn ÞÞ 6 B2 k , with

F 2 ¼ F 2 ðk/kCð½0;T;L1 ðXÞÞ ; kr/t kCð½0;T;L3 ðXÞÞ Þ and B2 ¼ B2 ðkr/t kCð½0;T;L2 ðXÞÞ Þ: 5.3. Second order Lagrange multiplier scheme (LM2) Now we want to design a second order scheme from the Lagrange multiplier scheme (LM1). The general form of this new scheme is

8 nþ1  þ ðrwnþ12 ; rwÞ  ¼ 0; > < ðdt / ; wÞ 1 nþ12   þ ð/ e nþ12 qnþ12 ; /Þ;  ðw ; /Þ ¼ ðr/nþ2 ; r/Þ > : e2 nþ1 nþ12 nþ1  e  ðd q ; q Þ ¼ ð / d / ; q Þ; t t 2

ð5:32Þ

e nþ2 is going to be critical for defining a linear or a nonlinear scheme and for conserving the second order where the choice of / in time accuracy. In the following Lemma we will see that the energy-stability of the scheme is independent of the choice of ~ nþ12 . / 1

Lemma 5.2. The scheme (5.32) is unconditionally energy-stable, in the sense that the following modified energy law holds

Eð/nþ1 ; qnþ1 Þ þ krwnþ2 k2L2 ðXÞ ¼ 0; dt e 1

ð5:33Þ

where e Eð/; qÞ is defined as in the scheme (LM1).  ¼ dt /nþ1 and q  ¼ qnþ12 , the discrete energy law (5.33) holds. h  ¼ wnþ1 , / Proof. Taking w e nþ12 ¼ /nþ12 , we obtain a nonlinear (stable) second order scheme. In order to avoid this nonlinearity preUnder the choice / serving the second order accuracy, we consider the following expression 2

2

/ðt nþ1 Þ ¼ /ðt n1 Þ þ k /t ðtn1 Þ þ Oðk Þ ¼ /ðt n1 Þ þ k Dwðtn1 Þ þ Oðk Þ; 2

2

2

2

2

e nþ12 as where we rewrite the term /t ðt n1 Þ by Dwðt n1 Þ. Hence, we can define / 2

2

e nþ12 ¼ /n12 þ k Dwn12 : / Taking into account the following relations 1 1   ¼ ðrwn12  rqnþ12 ; /Þ   ðqnþ12 ; rwn12  r/Þ; ðqnþ2 Dwn2 ; /Þ

152

F. Guillén-González, G. Tierra / Journal of Computational Physics 234 (2013) 140–171 1

1

1

Þ; Þ ¼ ðrwn2  rdt /nþ1 ; q Þ  ðdt /nþ1 ; rwn2  rq ðdt /nþ1 Dwn2 ; q we are going to consider the following two-step linear scheme (LM2): 1 1 Given ðwn2 ; /n1 ; /n ; qn Þ 2 H1 ðXÞ  H1 ðXÞ  H1 ðXÞ  L2 ðXÞ compute ðwnþ2 ; /nþ1 ; qnþ1 Þ 2 H1 ðXÞ  H1 ðXÞ  L2 ðXÞ such that:

8 nþ1 nþ12   > > < ðdt / ; wÞ þ ðrw ; rwÞ ¼ 0; 1 1 nþ  þ ð/n12 qnþ12 ; /Þ   ¼ ðr/ 2 ; r/Þ   kððrwn12  rqnþ12 ; /Þ  þ ðqnþ12 ; rwn12  r/ÞÞ; ðwnþ2 ; /Þ > > 1 1 : e2 n12 nþ1 nþ1 nþ1 n n nþ1  ÞÞ: Þ  kððrw 2  rdt / ; q Þ þ ðdt / ; rw 2  rq ðdt q ; qÞ ¼ ð/ dt / ; q 2

ð5:34Þ

Since we have developed a two-step scheme, it is necessary to include a detailed description of the first iteration: given /0 ¼ /ð0Þ in X, the first iteration consists in, 1. Compute q0 such that q0 ¼ e12 ðð/0 Þ2  1Þ.  ¼ ðr/0 ; rwÞ   þ ðfe ð/0 Þ; wÞ. 2. Compute w0 such that ðw0 ; wÞ 1 3. Compute ð/1 ; w2 ; q1 Þ such that

8 1 1 >  2  > < ðdt / ; wÞ þ ðrw ; rwÞ ¼ 0; 1 1  þ ð/ e 12 q12 ; /Þ;  ¼ ðr/2 ; r/Þ  ðw2 ; /Þ > > 1 2 :e 1 1  e  2 ðdt q ; qÞ ¼ ð / dt / ; qÞ; 2

where

e 12  /0 þ k Dw0 / 2 with the corresponding integration by parts to deal with Dw0 . Remark 5.3. In order to save computational time, it is possible to ‘‘rewrite’’ qnþ1 depending on /nþ1 and the previous steps in the same sense of Remark 4.10. Lemma 5.4. Scheme (5.34) is uniquely solvable. Proof. Following the same argument presented in Lemma 4.8. h 6. Numerical simulations The aim of this section is to compare the results of several simulations that we have carried out using the schemes derived through the paper. The section is organized as follows: In the first part we compare the results of the first order schemes considering different values of the time steps meanwhile in the second part we compare the behavior of the second order schemes. Finally, we compare together all the new schemes presented in the paper (first and second order). We are considering a finite element discretization in space for all the schemes, where the P 1 approximation is assumed for /h and wh . For the Lagrange multiplier approach, a P 1 discretization is also assumed for qh . We have chosen the domain 1 X ¼ ½0; 12 using a structured mesh with h ¼ 90 . In all the simulations we take c ¼ 0:0001, a fixed value for the interaction length parameter e ¼ 0:01 and we also take the same initial condition corresponding to a random initial data with values between 102 and 102 , in order to simulate a spinodal decomposition (see Fig. 1). All the simulations are carried out using FreeFem++ software, solving the corresponding algebraic systems by a LU direct method. It is interesting to note that in the simulations we have not had any solvability problem for the schemes OD1 and OD2 even when we consider a time step that does not satisfy the solvability constraint derived in Lemmas 4.5 and 5.1 (k < 8e4 ). The figures show the dynamic of the different schemes and the equilibrium solutions obtained in each case (we are considering that a scheme has achieved an equilibrium solution when the energy remains constant in time), the behavior of the nþ1 free energy Eð/nþ1 ; /n Þ and the total numerical dissipation NDð/nþ1 ; /n Þ of h Þ, the phobic numerical dissipation NDphobic ð/ nþ1 each scheme. In the simulations, we are calculating k  NDphobic ð/ ; /n Þ and by abuse of language we denote this term by NDphobic ð/nþ1 ; /n Þ. It is important to note that in the Lagrange multiplier schemes (LM1, LM2) we are also presenting the results obtained with respect to the energy Eð/nþ1 Þ and not the ones related to the energy e Eð/nþ1 Þ. 6.1. First order schemes: E1, E1mod, LM1 and OD1 6.1.1. Time step k ¼ 105 In Figs. 2 and 3 it is shown the evolution in time of the different schemes and how different equilibrium solutions are achieved. The reason of the different dynamics is due to the high numerical dissipation introduced in some of the schemes. Figs. 4 and 5 show the evolution of the free energy in time for the four schemes and how the equilibrium solutions achieved have different free energy values for each scheme. The phobic numerical dissipation behaves as we derived in Section 4.3, where the order of OD1 is the lowest one, and it is not always positive, as can be seen in Figs. 6 and 7. As we

F. Guillén-González, G. Tierra / Journal of Computational Physics 234 (2013) 140–171

Fig. 1. Random initial data.

Fig. 2. Comparison of the dynamic of the schemes in [0, 0.41].

153

154

F. Guillén-González, G. Tierra / Journal of Computational Physics 234 (2013) 140–171

also expected, the phobic numerical dissipation of E1mod is lower than the one related to E1. Meanwhile the phobic numerical dissipation of LM1 has a behavior like E1 and E1mod. The evolution in time of the total numerical dissipation (with respect to the energy Eð/nþ1 Þ) is shown in Figs. 8 and 9. In this case the dissipation related to OD1 is again the lowest one, the one related to E1mod is also lower than the one related to E1 and for LM1 we obtain alternating behavior (sometimes numerical dissipation and others numerical source). This fact

Fig. 3. Comparison of the dynamic of the schemes in [3.9, 25].

Spinodal k=1e−05 3000 E1mod E1 OD1 LM

2500

Mixing Energy

2000

1500

1000

500

0

0

0.05

0.1

0.15

0.2 Time

0.25

Fig. 4. Free energy in [0, 0.4].

0.3

0.35

0.4

155

F. Guillén-González, G. Tierra / Journal of Computational Physics 234 (2013) 140–171 Spinodal k=1e−05 550 E1mod E1 OD1 LM

500 450

Mixing Energy

400 350 300 250 200 150 100 50

2

4

6

8

10 Time

12

14

16

18

20

Fig. 5. Free energy in [0.4, 20].

Spinodal k=1e−05 0.35 E1mod E1 OD1 LM

0.25 0.2 0.15 0.1 0.05 0 −0.05 0

0.005

0.01

0.015

0.02 Time

0.025

0.03

0.035

0.04

Fig. 6. NDphobic in [0, 0.4].

Spinodal k=1e−05

−4

5

x 10

E1mod E1 OD1 LM

4 3 Phobic Numerical Dissipation

Phobic Numerical Dissipation

0.3

2 1 0 −1 −2 −3 −4 −5 2

3

4

5

6 Time

7

Fig. 7. Zoom of the NDphobic for the OD1.

8

9

10 −3

x 10

156

F. Guillén-González, G. Tierra / Journal of Computational Physics 234 (2013) 140–171 Spinodal k=1e−05

4

4

x 10

E1mod E1 OD1 LM

3.5

Numerical Dissipation

3 2.5 2 1.5 1 0.5 0 −0.5 0

0.005

0.01

0.015

0.02 Time

0.025

0.03

0.035

0.04

Fig. 8. ND in [0, 0.04].

Spinodal k=1e−05 250 E1mod E1 OD1 LM

200

Numerical Dissipation

150

100

50

0

−50

−100

0.05

0.1

0.15

0.2

0.25 Time

Fig. 9. ND in [0.04, 0.4].

Fig. 10. E1 equilibrium solution.

0.3

0.35

0.4

F. Guillén-González, G. Tierra / Journal of Computational Physics 234 (2013) 140–171

157

Fig. 11. E1mod equilibrium solution.

Fig. 12. OD1 equilibrium solution.

Fig. 13. LM1 equilibrium solution.

does not contradict with the results obtained in Section 4.4 because in that case we obtain the unconditional energy-stability of the scheme LM1 for a modified energy e Eð/nþ1 Þ but we do not have any result concerning the behavior of the total numerical dissipation with respect to the energy Eð/nþ1 Þ.

6.1.2. Time step k ¼ 107 In this case, the schemes OD1, E1 and E1mod reach the same equilibrium solution as can be viewed in Figs. 10–13. Furthermore this solution is exactly the same that we have obtained for OD1 taking k ¼ 105 . The reason of achieving the same solution with k ¼ 107 is that Eyre’s schemes introduces too much numerical dissipation if the time step is not small enough, so they need a very small time step in order to achieve a correct equilibrium solution. Only the LM1 scheme arrives to a slightly different equilibrium solution. As a conclusion of this subsection we claim that from the computational point of

158

F. Guillén-González, G. Tierra / Journal of Computational Physics 234 (2013) 140–171 Spinodal k=1e−07 2600 E1mod E1 OD1 LM

2400

Mixing Energy

2200

2000

1800

1600

1400

0

0.001

0.002

0.003

0.004

0.005 Time

0.006

0.007

0.008

0.009

0.01

Fig. 14. Free energy in [0, 0.01]. Spinodal k=1e−07 1600 E1mod E1 OD1 LM

1500 1400

Mixing Energy

1300 1200 1100 1000 900 800 700 0.01

0.02

0.03

0.04

0.05

0.06 Time

0.07

0.08

0.09

0.1

0.11

Fig. 15. Free energy in [0.01, 0.11].

Spinodal k=1e−07 760 E1mod E1 OD1 LM

740

Mixing Energy

720

700

680

660

640

620 0.11

0.12

0.13

0.14

0.15

0.16 Time

0.17

Fig. 16. Free energy in [0.11, 0.21].

0.18

0.19

0.2

0.21

159

F. Guillén-González, G. Tierra / Journal of Computational Physics 234 (2013) 140–171 Spinodal k=1e−07 630 E1mod E1 OD1 LM

620 610

Mixing Energy

600 590 580 570 560 550 540 530 0.21

0.22

0.23

0.24

0.25

0.26 Time

0.27

0.28

0.29

0.3

0.31

Fig. 17. Free energy in [0.21, 0.31].

Spinodal k=1e−07 550 E1mod E1 OD1 LM

540

Mixing Energy

530

520

510

500

490

480 0.31

0.32

0.33

0.34

0.35

0.36 Time

0.37

0.38

0.39

0.4

0.41

Fig. 18. Free energy in [0.31, 0.41]. Spinodal k=1e−07 520 E1mod E1 OD1 LM

510

Mixing Energy

500

490

480

470

460

450 0.41

0.42

0.43

0.44

0.45

0.46 Time

0.47

Fig. 19. Free energy in [0.41, 0.51].

0.48

0.49

0.5

0.51

160

F. Guillén-González, G. Tierra / Journal of Computational Physics 234 (2013) 140–171 Spinodal k=1e−07 500 E1mod E1 OD1 LM

450

Mixing Energy

400

350

300

250

200

150

1

2

3

4

5

6

7

Time

Fig. 20. Free energy in [0.51, 7]. Spinodal k=1e−07 197 E1mod E1 OD1 LM

196 195

Mixing Energy

194 193 192 191 190 189 188 7.5

8

8.5 Time

9

9.5

10

Fig. 21. Free energy in [7, 10].

Spinodal k=1e−07

−4

10

x 10

E1mod E1 OD1 LM

Phobic Numerical Dissipation

8

6

4

2

0

−2 0.01

0.02

0.03

0.04

0.05

0.06 Time

0.07

Fig. 22. NDphobic in [0.01, 0.11].

0.08

0.09

0.1

0.11

161

F. Guillén-González, G. Tierra / Journal of Computational Physics 234 (2013) 140–171 Spinodal k=1e−07

−5

3

x 10

E1mod E1 OD1 LM

Phobic Numerical Dissipation

2.5

2

1.5

1

0.5

0

−0.5 0.21

0.22

0.23

0.24

0.25

0.26 Time

0.27

0.28

0.29

0.3

0.31

Fig. 23. NDphobic in [0.21, 0.31].

Spinodal k=1e−07 1000 E1mod E1 OD1 LM

900

Numerical Dissipation

800 700 600 500 400 300 200 100 0 0.01

0.02

0.03

0.04

0.05

0.06 Time

0.07

0.08

0.09

0.1

0.11

Fig. 24. ND in [0.01, 0.11]. Spinodal k=1e−07 30 E1mod E1 OD1 LM

25

Numerical Dissipation

20

15

10

5

0

−5 0.21

0.22

0.23

0.24

0.25

0.26 Time

0.27

Fig. 25. ND in [0.21, 0.31].

0.28

0.29

0.3

0.31

162

F. Guillén-González, G. Tierra / Journal of Computational Physics 234 (2013) 140–171

Fig. 26. Dynamic of OD2 and OD2mod.

Fig. 27. LM2 equilibrium solution.

2600 OD2 OD2mod LM2

2400 2200

Mixing Energy

2000 1800 1600 1400 1200 1000 800 600

0

0.01

0.02

0.03

0.04

0.05 Time

0.06

Fig. 28. Free energy in [0, 0.1].

0.07

0.08

0.09

0.1

163

F. Guillén-González, G. Tierra / Journal of Computational Physics 234 (2013) 140–171 800 OD2 OD2mod LM2

750 700

Mixing Energy

650 600 550 500 450 400 350 0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Time

Fig. 29. Free energy in [0.1, 1].

420 OD2 OD2mod LM2

400 380

Mixing Energy

360 340 320 300 280 260 240 220

1

1.5

2

2.5

3 Time

3.5

4

4.5

5

Fig. 30. Free energy in [1, 5.1].

250 OD2 OD2mod LM2

240

Mixing Energy

230

220

210

200

190

5.5

6

6.5

7

7.5 Time

8

Fig. 31. Free energy in [5.1, 10].

8.5

9

9.5

10

164

F. Guillén-González, G. Tierra / Journal of Computational Physics 234 (2013) 140–171 0.06 OD2 OD2mod LM2

Phobic Numerical Dissipation

0.05

0.04

0.03

0.02

0.01

0

−0.01 0

0.01

0.02

0.03

0.04

0.05 Time

0.06

0.07

0.08

0.09

0.1

Fig. 32. NDphobic in [0, 0.1].

−5

3.5

x 10

3

Phobic Numerical Dissipation

2.5 2 1.5 1 0.5 0 OD2 OD2mod LM2

−0.5 −1 −1.5 3

3.1

3.2

3.3

3.4

3.5 Time

3.6

3.7

3.8

3.9

4

Fig. 33. NDphobic in [3, 4].

1000 OD2 OD2mod LM2

800 600

Total Dissipation

400 200 0 −200 −400 −600 −800 0

0.01

0.02

0.03

0.04

0.05 Time

0.06

Fig. 34. ND in [0, 0.1].

0.07

0.08

0.09

0.1

165

F. Guillén-González, G. Tierra / Journal of Computational Physics 234 (2013) 140–171

view, OD1 is the best of the first order schemes because we can take bigger time-steps (reducing the computational cost) and still obtain a good behavior of the scheme. Figs. 14–21 show the evolution of the free energy in time for the four schemes. The energies do not evolve in the same way for all the schemes but at last all of them reach the same equilibrium energy (except LM1, it achieves almost the same energy). The evolution of the phobic and total numerical dissipation is shown in Figs. 22–25, respectively. They behave as considering k ¼ 105 , i.e., the order of OD1 is the lowest one, the order of E1mod is lower than the one related to E1 and for LM1 we

0.25 OD2 OD2mod LM2

0.2 0.15

Total Dissipation

0.1 0.05 0 −0.05 −0.1 −0.15 −0.2 3

3.1

3.2

3.3

3.4

3.5 Time

3.6

Fig. 35. ND in [3, 4].

Fig. 36. OD1 equilibrium solution.

Fig. 37. OD2 equilibrium solution.

3.7

3.8

3.9

4

166

F. Guillén-González, G. Tierra / Journal of Computational Physics 234 (2013) 140–171

obtain again an oscillatory behavior. The difference with respect to the previous case is that the smaller time step allows to obtain smaller values of NDphobic ð/nþ1 ; /n Þ and NDð/nþ1 ; /n Þ for each scheme. 6.2. Second order schemes: OD2, OD2mod and LM2 with k ¼ 105 In this section we present the results concerning to the simulations using the second order schemes presented in the paper. We have observed that OD2 and OD2mod have exactly the same behavior and they achieve the same equilibrium solutions that we have considered valid in the first order schemes section as can be viewed in Fig. 26. Furthermore, the dynamic

Fig. 38. LM1 equilibrium solution.

Fig. 39. LM2 equilibrium solution.

2600 OD2 OD1 LM1 LM2

2400 2200

Mixing Energy

2000 1800 1600 1400 1200 1000 800 600

0

0.01

0.02

0.03

0.04

0.05 Time

0.06

Fig. 40. Free energy in [0, 0.1].

0.07

0.08

0.09

0.1

167

F. Guillén-González, G. Tierra / Journal of Computational Physics 234 (2013) 140–171 800 OD2 OD1 LM1 LM2

700

Mixing Energy

600

500

400

300

200

0.5

1

1.5

2

2.5 Time

3

3.5

4

4.5

5

Fig. 41. Free energy in [0.1, 5.1].

260 OD2 OD1 LM1 LM2

240

Mixing Energy

220

200

180

160

140

120

5.5

6

6.5

7

7.5 Time

8

8.5

9

9.5

10

Fig. 42. Free energy in [5.1, 10.1].

200 190 180

Mixing Energy

170 160 OD2 OD1 LM1 LM2

150 140 130 120 110 100

10.5

11

11.5

12

12.5 Time

13

Fig. 43. Free energy in [10.1, 15].

13.5

14

14.5

15

168

F. Guillén-González, G. Tierra / Journal of Computational Physics 234 (2013) 140–171 7 OD2 OD1 LM1 LM2

Phobic Numerical Dissipation

6 5 4 3 2 1 0 −1 0

0.01

0.02

0.03

0.04

0.05 Time

0.06

0.07

0.08

0.09

0.1

Fig. 44. NDphobic in [0, 0.1].

0.8 OD2 OD1 LM1 LM2

0.7

Phobic Numerical Dissipation

0.6 0.5 0.4 0.3 0.2 0.1 0 −0.1 −0.2 0.11

0.115

0.12

0.125 Time

0.13

0.135

0.14

Fig. 45. NDphobic in [0.11, 0.14].

4

9

x 10

OD2 OD1 LM1 LM2

8 7

Total Dissipation

6 5 4 3 2 1 0 −1 0

0.01

0.02

0.03

0.04

0.05 Time

0.06

Fig. 46. NDphobic in [0, 0.1].

0.07

0.08

0.09

0.1

169

F. Guillén-González, G. Tierra / Journal of Computational Physics 234 (2013) 140–171

of LM2 is almost the same of the schemes OD2 and OD2mod and at the end the scheme arrives at the same equilibrium solution as we can see in Fig. 27. The behavior of the numerical dissipation of LM2 has an oscillatory behavior as in the previous LM1 case. The free energies evolve exactly in the same way for OD2 and OD2mod meanwhile in the LM2 case, the energy evolves in another way arriving at the same equilibrium energy (Figs. 28–31). Figs. 32–35 show the evolution in time of the phobic and the total numerical dissipation, respectively. LM2 has an oscillatory behavior (like LM1) but the size of the dissipation terms is much lower than in the LM1 case (considering the same time step). The dissipation terms related to OD2 and OD2mod are very small (as we expected). 6.3. Comparison between first and second order: OD1, OD2, LM1 and LM2 Hitherto we have shown the simulations of the four new schemes presented in the paper when the considered time step is k ¼ 105 . Three of them (OD1, OD2 and LM2) have what we have considered the correct behavior and the other one (LM1) arrives to slightly different equilibrium solution. The next step in our simulations consists on taking a bigger time step and compare the results obtained for each scheme. 6.3.1. Time step k ¼ 104 In this case we obtain a different equilibrium solution for each of the schemes considered as can be observed in Figs. 36– 39. It is important to note that OD2 is the best of the schemes presented in the paper because it is the only one that even considering a bigger time step, arrives to the correct equilibrium solution. Figs. 40–43 show the different evolution of the free energy for each scheme and how they arrive to different equilibrium solutions. The evolution of the phobic and total numerical dissipation is shown in Figs. 44–47, respectively. The behavior of OD1 and OD2 is quite similar in the phobic numerical case (as we expected because both schemes use exactly the same approximation of the potential term f ð/nþ1 ; /n Þ) but in the total dissipation case we observe the difference between both schemes. For the Lagrange multiplier approaches we see how both of them introduce a lot of phobic and total numerical dissipation (in particular LM1 introduces much more than LM2). 6.4. Numerical computations refining the spatial grid and the time step In order to assure that the results obtained in the previous experiments do not depend on the size of the space mesh considered, the authors have carried out several simulations considering a refined space mesh. In particular, we have chosen a 1 structured mesh with h ¼ 150 in the domain X ¼ ½0; 12 . The results obtained agree with the results presented when a struc1 tured mesh with h ¼ 90 is considered, in the following sense:  OD1 arrives to the correct equilibrium solution considering a time step k 6 105 while E1 and E1mod arrives to the correct equilibrium solution when the time step considered is k 6 106 . Moreover, considering k 6 105 , LM1 arrives at an equilibrium solution that is slightly different from the correct one.  Using the second order schemes, we have obtained that OD2 and OD2mod arrives to the correct equilibrium solution considering a time step k 6 104 while LM2 arrives to the correct equilibrium solution when the time step considered is k 6 105 . 1400 OD2 OD1 LM1 LM2

1200

Total Dissipation

1000 800 600 400 200 0 −200 0.11

0.115

0.12

0.125 Time

Fig. 47. NDphobic in [0.11, 0.14].

0.13

0.135

0.14

170

F. Guillén-González, G. Tierra / Journal of Computational Physics 234 (2013) 140–171

7. Conclusions In this paper we have developed two new ways of handling the double-well (polynomial) potential term that allows us to design four new linear schemes for the Cahn–Hilliard problem, OD1, LM1, OD2 and LM2, the first two of first order in time and the other two of second order. OD1 and OD2 are optimal from the phobic numerical dissipation point of view and its behavior can be considered more effective in the computational experiments (compared with the other linear schemes that we have studied in this work) but the numerical analysis to assure the long time stability of these kind of schemes remains as an open problem for future works. On the other hand, LM1 and LM2 are designed as unconditionally energy-stable linear schemes (for a modified energy e Eð/nþ1 ; qnþ1 Þ, where qnþ1 play a role of Lagrange multiplier) and the numerical experiments also work well despite the oscillations of the numerical dissipation. Moreover, during the simulations is observed that the second order schemes are more accurate than the first order ones due to the possibility of considering larger time steps without loosing accuracy of the equilibrium solutions (saving computational cost). In particular, from this point of view it is shown in the simulations section that OD2 is the most accurate of all the computed schemes. Acknowledgments The authors thank the anonymous reviewers for their valuable comments and suggestions to improve the quality of this manuscript. This work has been partially supported by project MTM2009-12927, Spain. Appendix A In order to obtain a second order potential approximation, we are going to use the following quadrature formulas (exact for P0 ):

Z

gðxÞdx ¼ ðb  aÞgðbÞ  g 0 ðnÞ

ðb  aÞ2 ; 2

gðxÞdx ¼ ðb  aÞgðaÞ þ g 0 ðnÞ

ðb  aÞ2 : 2

b

a

Z

b

a

Applying these quadrature formulas to f ¼ f1 þ f2 in a separated form, where f1 ¼ F 01 , f2 ¼ F 02 with

F 1 ð/Þ ¼

1 b/2 2e2

and F 2 ð/Þ ¼ Fð/Þ  F 1 ð/Þ ¼

1 ð/4  2ðb þ 1Þ/2 þ 1Þ 4e2

b > 0:

Using the first one for F 1 ð/Þ and the second one for F 2 ð/Þ,

Fð/nþ1 Þ  Fð/n Þ ¼

Z

/nþ1 /n

dFð/Þ d/ ¼ d/

Z

/nþ1

/n

½f1 ð/Þ þ f2 ð/Þd/

¼ f1 ð/nþ1 Þð/nþ1  /n Þ  f10 ð/nþn Þ

ð/nþ1  /n Þ2 ð/nþ1  /n Þ2 þ f2 ð/n Þð/nþ1  /n Þ þ f20 ð/nþl Þ ; 2 2

ðA:35Þ

where /nþn ¼ /n þ nð/nþ1  /n Þ, /nþl ¼ /n þ lð/nþ1  /n Þ with n; l 2 ð0; 1Þ. Assuming /nþl ’ /n and taking into account that f10 is constant,

  1 Fð/nþ1 Þ  Fð/n Þ ’ f1 ð/nþ1 Þ þ f2 ð/n Þ þ ½f20 ð/n Þ  f10 ð/n Þð/nþ1  /n Þ ð/nþ1  /n Þ 2 2

3

being the error proportional to k ðdt /nþ1 Þ2 ð/nþl  /n Þ ¼ Oðk Þ. Therefore, the choice of the potential term in this case reads,

fk ð/nþ1 ; /n Þ ¼

  1 0 n nþ1 n n nþ1 n 0 ½f f ð/ Þ þ f ð/ Þ þ ð/ Þ  f ð/ Þð/  / Þ : 1 2 2 1 2 e2 1

ðA:36Þ

Since f1 is linear, this expression can be rewritten as the approximation of the potential used in scheme (OD1) (and (OD2)) and in particular, this potential approximation does not depend on b. Moreover, we can not assure that using this idea we will obtain a linear scheme for all kind of potentials. Appendix B Using Taylor’s formula

F. Guillén-González, G. Tierra / Journal of Computational Physics 234 (2013) 140–171

171

 2 2 k 1 d k f ð/ðtnþ1 ÞÞ ¼ f ð/ðtn ÞÞ þ f 0 ð/ðtn ÞÞ/t ðt n Þ þ f ð/ðtÞÞj t¼t nþs 2 2 2 dt 2 2

 2 2 k 1 d k ¼ fk ð/ðt nþ1 Þ; /ðt n ÞÞ þ f 0 ð/ðt n ÞÞ ð/t ðt n Þ  dt /ðt nþ1 ÞÞ þ f ð/ðtÞÞj : t¼t nþs 2 2 dt 2 2

for some s 2 ð0; 12Þ. Then

0

1 d2 A kf ð/ðtn ÞÞ/tt ðt nþn ÞkL2 ðXÞ þ 2 f ð/ðtÞÞjt¼tnþs 2 dt L ðXÞ

2@

kfk ð/ðt nþ1 ; /ðtn ÞÞÞ  f ð/ðtnþ1 ÞÞkL2 ðXÞ 6 C k 2

0

for some n 2 ð0; 1Þ. Since 2

d

dt

2

f ð/ðtÞÞ ¼

d 0 ðf ð/ðtÞÞ/t ðtÞÞ ¼ f 00 ð/ðtÞÞð/t ðtÞÞ2 þ f 0 ð/ðtÞÞ/tt ðtÞ; dt

we conclude that

kfk ð/ðt nþ1 ; /ðtn ÞÞÞ  f ð/ðtnþ1 ÞÞkL2 ðXÞ 6 D2 k

2

2

with D2 ¼ D2 ðk/kCð½0;T;L1 ðXÞÞ ; k/t kCð½0;T;L4 ðXÞÞ ; k/tt kCð½0;T;L2 ðXÞÞ Þ. References [1] H. Abels, Diffuse Interface Models for Two-Phase Flows of Viscous Incompressible Fluids, Habilitation Thesis, . [2] S. Badia, F. Guillén-González, J.V. Gutiérrez-Santacreu, Finite element approximation of nematic liquid crystal flows using a saddle-point structure, Journal of Computational Physics 230 (2011) 1686–1706. [3] J.W. Barrett, J.F. Blowey, H. Garcke, Finite element approximation of the Cahn–Hilliard equation with degenerate mobility, SIAM Journal on Numerical Analysis 37 (1999) 286–318. [4] P.W. Bates, P.C. Fife, The dynamics of nucleation for the Cahn–Hilliard equation, SIAM Journal on Applied Mathematics 53 (1993) 990–1008. [5] R. Becker, X. Feng, A. Prohl, Finite element approximations of the Ericksen–Leslie model for nematic liquid crystal flow, SIAM Journal on Numerical Analysis 46 (2008) 1704–1731. [6] F. Boyer, S. Minjeaud, Numerical schemes for a three component Cahn–Hilliard model, ESAIM: Mathematical Modelling and Numerical Analysis. 45 (2011) 697–738. [7] J.W. Cahn, J.E. Hilliard, Free energy of a non-uniform system. I. Interfacial free energy, The Journal of Chemical Physics 28 (1958) 258–267. [8] H.D. Ceniceros, Tracking fluid interfaces approaching singular events, Boletín de la Sociedad Española de Matemática Aplicada 48 (2009) 31–57. [9] M.I.M. Copetti, C.M. Elliott, Numerical analysis of the Cahn–Hilliard equation with a logarithmic free energy, Numerische Mathematik 63 (1992) 39–65. [10] Q. Du, R.A. Nicolaides, Numerical analysis of a continuum model of phase transition, SIAM Journal on Numerical Analysis 28 (1991) 1310–1322. [11] C.M. Elliott, D.A. French, A nonconforming finite-element method for the two-dimensional Cahn–Hilliard equation, SIAM Journal on Numerical Analysis 26 (1989) 884–903. [12] C.M. Elliott, D.A. French, Numerical studies of the Cahn–Hilliard equation for phase separation, IMA Journal of Applied Mathematics 38 (1987) 97–128. [13] C.M. Elliott, D.A. French, F.A. Milner, A second order splitting method for the Cahn–Hilliard equation, Numerische Mathematik 54 (1989) 575–590. [14] C.M. Elliott, H. Garcke, On the Cahn–Hilliard equation with degenerate mobility, SIAM Journal on Mathematical Analysis 27 (1996) 404–423. [15] C.M. Elliott, Z. Songmu, On the Cahn–Hilliard equation, Archive for Rational Mechanics and Analysis 96 (1986) 339–357. [16] David J. Eyre, An Unconditionally Stable One-Step Scheme for Gradient System, , unpublished. [17] X. Feng, Fully discrete finite element approximations of the Navier–Stokes–Cahn–Hilliard diffuse interface model for two-phase fluid flows, SIAM Journal on Numerical Analysis 44 (2006) 1049–1072. [18] X. Feng, A. Prohl, Error analysis of a mixed finite element method for the Cahn–Hilliard equation, Numerische Mathematik 99 (2004) 47–84. [19] D. Furihata, A stable and conservative finite difference scheme for the Cahn–Hilliard equation, Numerische Mathematik 87 (2001) 675–699. [20] H. Gomez, V.M. Calo, Y. Bazilevs, Thomas J.R. Hughes, Isogeometric analysis of the Cahn–Hilliard phase-field model, Computer Methods in Applied Mechanics and Engineering 197 (2008) 4333–4352. [21] H. Gomez, Thomas J.R. Hughes, Provably unconditionally stable, second-order time-accurate, Journal of Computational Physics 230 (2011) 5310–5327. [22] Z. Hu, S.M. Wise, C. Wang, J.S. Lowengrub, Stable and efficient finite-difference nonlinear-multigrid schemes for the phase field crystal equation, Journal of Computational Physics 228 (2009) 5323–5339. [23] J. Hua, P. Lin, C. Liu, Q. Wang, Energy law preserving C 0 finite element schemes for phase field models in two-phase flow computations, Journal of Computational Physics 230 (2011) 7115–7131. [24] D. Kay, V. Styles, E. Suli, Discontinuous Galerkin finite element approximation of the Cahn–Hilliard Equation with convection, SIAM Journal on Numerical Analysis 47 (2009) 2660–2685. [25] J. Kim, K. Kang, J. Lowengrub, Conservative multigrid methods for Cahn–Hilliard fluids, Journal of Computational Physics 193 (2004) 511–543. [26] P. Lin, C. Liu, H. Zhang, An energy law preserving C 0 finite element scheme for simulating the kinematic effects in liquid crystal dynamics, Journal of Computational Physics 227 (2007) 1411–1427. [27] A. Novick-Cohen, L.A. Segel, Nonlinear aspects of the Cahn–Hilliard equation, Physica D 10 (1984) 277–298. [28] J. Shen, X. Yang, Numerical approximations of Allen–Cahn and Cahn–Hilliard equations, Discrete and Continuous Dynamical Systems 28 (2010) 1669– 1691. [29] J.D. van der Waals, The thermodynamic theory of capillarity flow under the hypothesis of a continuous variation of density, Verhandelingen Koninklijke Nederlandse Akademie van Wetenschappen Amsterdam (Sec 1) 1 (1893) 1–56. [30] S.M. Wise, Unconditionally stable finite difference, nonlinear multigrid simulation of the Cahn–Hilliard–Hele–Shaw system of equations, Journal of Scientific Computing 44 (1) (2010). [31] S. Zhang, M. Wang, A nonconforming finite element method for the Cahn–Hilliard equation, Journal of Computational Physics 229 (2010) 7361–7372.