• No results found

Space-time adaptive methods for phase-field models

N/A
N/A
Protected

Academic year: 2021

Share "Space-time adaptive methods for phase-field models"

Copied!
120
0
0

Bezig met laden.... (Bekijk nu de volledige tekst)

Hele tekst

(1)

Space-time adaptive methods for phase-field models

Citation for published version (APA):

Wu, X. (2017). Space-time adaptive methods for phase-field models. Technische Universiteit Eindhoven.

Document status and date: Published: 29/06/2017

Document Version:

Publisher’s PDF, also known as Version of Record (includes final page, issue and volume numbers)

Please check the document version of this publication:

• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website.

• The final author version and the galley proof are versions of the publication after peer review.

• The final published version features the final layout of the paper including the volume, issue and page numbers.

Link to publication

General rights

Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights.

• Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain

• You may freely distribute the URL identifying the publication in the public portal.

If the publication is distributed under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license above, please follow below link for the End User Agreement:

www.tue.nl/taverne

Take down policy

If you believe that this document breaches copyright please contact us at:

openaccess@tue.nl

providing details and we will investigate your claim.

(2)

S

PACE

-

TIME ADAPTIVE METHODS

(3)

639.031.033.

© 2017 by Xunxun Wu. All rights reserved.

A catalogue record is available from the Eindhoven University of Technology Library ISBN: 978-90-386-4308-3

Cover designed by Xunxun Wu Typeset with LATEX

(4)

S

PACE

-

TIME ADAPTIVE METHODS

FOR PHASE

-

FIELD MODELS

PROEFSCHRIFT

ter verkrijging van de graad van doctor aan de Technische Universiteit Eindhoven, op gezag van de rector magnificus prof.dr.ir. F.P.T. Baaijens, voor een commissie aangewezen

door het College voor Promoties, in het openbaar te verdedigen op donderdag 29 juni 2017 om 16:00 uur

door

Xunxun Wu

(5)

voorzitter: prof.dr. L. P. H. de Goey

promotor: prof.dr.ir. E. H. van Brummelen

copromotor: dr.ir. K. G. van der Zee (University of Nottingham) externe leden: prof.dr. O. Lakkis (Free University of Bozen-Bolzano)

prof.dr. S. Prudhomme (École Polytechnique de Montréal) overige leden: prof.dr.ir. P. D. Anderson

prof.dr.ir. B. Koren prof.dr. J. G. M. Kuerten

Het onderzoek of ontwerp dat in dit proefschrift wordt beschreven is uitgevoerd in overeen-stemming met de TU/e Gedragscode Wetenschapsbeoefening.

(6)

S

UMMARY

The phase field method, originally proposed for phase separation in material science, has rapidly gained popularity in various fields over the years. Owing to its diffuse-interface assumption, this technique presents a new modeling paradigm for predicting interfacial phenomena, meanwhile avoiding the problems associated with sharp interface models. The computational challenge of phase-field models is that the solutions exhibit very com-plex alternating fast and slow variations in space and time that would highly benefit from adaptivity.

This work is devoted to present and investigate solutions to issues regarding the accu-rate and efficient simulation using adaptive discretization methods for phase-field models. These include the derivation of stable discretization methods and suitable error estimates, as well as the development of efficient space-time adaptive algorithms.

We begin by developing unconditionally energy-stable mixed variational formulations for phase-field models, which rely on the use of an implicit-explicit (IMEX) time-stepping scheme and a mixed finite element method in space. Then, for assessing the accuracy of the numerical solutions, we propose a new methodology for duality-based a posteriori er-ror estimation for fully discrete approximations, which hinges on a specially-tailored IMEX discretization of the dual problem. The resultant error estimates can separate the effects of spatial and temporal discretization error, and can be used to drive adaptive mesh refine-ment and adaptive time-step selection. The performance of the presented error estimates and the adaptive algorithm are numerically demonstrated through a series of numerical experiments.

(7)
(8)

C

ONTENTS

1 Introduction 1

2 Phase-field models as gradient flows 5

2.1 Abstract gradient flows. . . 5

2.2 Allen-Cahn and Cahn-Hilliard equation . . . 7

2.3 A phase-field model of tumor growth . . . 10

3 Energy-stable time-stepping schemes 13 3.1 Motivation . . . 13

3.2 The Cahn–Hilliard equation. . . 15

3.2.1 Strong formulation. . . 15

3.2.2 First-order accurate scheme . . . 17

3.2.3 Second-order accurate scheme. . . 18

3.3 Diffuse-interface tumor-growth model . . . 22

3.3.1 Strong formulation. . . 22

3.3.2 First-order accurate scheme . . . 23

3.3.3 Second-order accurate scheme. . . 24

3.4 Numerical experiments . . . 27

3.4.1 The Cahn-Hilliard equation: Accuracy test. . . 28

3.4.2 The Cahn-Hilliard equation: Spinodal decomposition. . . 28

3.4.3 Tumor-growth model: Accuracy test. . . 30

3.4.4 Tumor-growth model: Single tumor. . . 33

3.4.5 Tumor-growth model: Three tumors . . . 35

4 A posteriori error estimation for IMEX–Galerkin discretization 41 4.1 Motivation . . . 41

4.2 Abstract setting . . . 43

4.2.1 Weak formulation and error representation . . . 43

4.2.2 IMEX–FEM Discretization. . . 45

4.3 Space-time decomposed a posteriori error estimate . . . 46

4.3.1 Time-discrete error representation . . . 47

(9)

4.4 Computable error estimate . . . 53

5 Adaptive methodology for space and time refinement 57 5.1 Space adaptive strategy. . . 57

5.2 Space-time adaptive strategy . . . 60

5.2.1 Error indicators . . . 60

5.2.2 The space-time adaptive algorithm . . . 62

6 Applications of the error estimates and adaptive algorithms 65 6.1 Heat equation . . . 65

6.1.1 Computable error indicator. . . 67

6.1.2 Numerical results. . . 68

6.2 Allen-Cahn equation . . . 75

6.2.1 Computable error indicator. . . 78

6.2.2 Numerical results. . . 80

6.3 Tumor growth model. . . 86

6.3.1 Computable error indicator. . . 91

6.3.2 Numerical results. . . 92 7 Conclusions 95 Appendix 99 Bibliography 101 Acknowledgements 109 Curriculum Vitae 111

(10)

1

I

NTRODUCTION

The phase-field method has proved to be a promising computational approach for mod-eling and predicting a wide range of interfacial phenomena in various areas of research. The success of the method has led to a rapidly growing list of topics ranging from materials science to fluid dynamics, see, e.g., [2,23, 15,30,63] and, very recently, to mechanobiol-ogy such as tumor growth [60,50]. Surfaces and interfaces in the phase-field method are implicitly described by a continuous scalar field, known as the phase field or order param-eter. The phase field takes constant values in the bulk phases and varies continuously but steeply from one bulk value to the other. Interfaces are identified by the transition of the phase field, which are diffuse and have a finite width. A distinct advantage of phase-field models is their flexibility in treating complex geometries and topological changes such that there is no need to explicitly track the evolution of interfaces, as in the case of sharp inter-face models.

The evolution of the phase field is driven by a dissipative minimization of the total free energy of the system that may include the chemical bulk free energy, interfacial energy and the energy contributions associated to external fields such as stress and temperature fields. The behavior of the dynamic variable is governed by a particular class of highly nonlinear parabolic equations, which may exhibit very complex solutions with alternating fast and slow variations in space and time which would highly benefit from adaptivity. Determining efficient and accurate discretization techniques for capturing the variations is a significant challenge.

The aim of this thesis is to present and investigate solutions to issues regarding efficient and accurate simulation of phase-field models using space-time adaptive methods. This encompasses both the derivation of stable discretizations and the development of suitable

(11)

error estimates. In particular, we identify three primary objectives:

• To investigate the essential characteristics of phase-field models: free-energy dissi-pation in the context of general gradient flows;

• To derive unconditionally gradient-stable time-stepping schemes such that the free energy is dissipated at the discrete level;

• To develop duality-based a posteriori error estimates for fully discrete approxima-tions which allows to separate the effects of spatial and temporal discretization error, and can be used to drive adaptive mesh refinement and adaptive time-step selection.

T

HE ESSENTIAL CHARACTERISTICS OF PHASE

-

FIELD MODELS

Employing the continuum theory of mixtures, phase-field models can be derived through the framework of irreversible thermodynamics introduced by Gurtin in [46] such that con-stitutive equations satisfy a dissipation inequality. The key feature of this approach is that constitutive equations are allowed to depend on the variational derivative of the free en-ergy, which directly links to energy dissipation. Therefore, phase-field models can also be interpreted as general gradient flows, in which the energy decreases along solutions as fast as possible, see e.g. [76,59,62,26,44]. This is the point of view that we introduce in Chapter

2. Using the approach, we illustrate phase-field models as gradient flows for three exam-ples: the classical Allen–Cahn and Cahn–Hilliard models and a phase-field model of tumor growth.

U

NCONDITIONALLY GRADIENT

-

STABLE TIME

-

STEPPING SCHEMES

When it comes to numerically solving phase-field models that have the structure of a gradi-ent flow, it is natural to ask for appropriate time-discrete schemes that respect the special structure and inherit energy dissipation at the discrete level. The major obstructions to the development of such schemes for phase-field models arise from the nonlinear term originating from the nonconvex (bulk) free energy, which locally produces ‘backwards’ diffusion. This nonlinearity can impose prohibitive stability restrictions for naive time-discrete schemes. For example, for the Cahn-Hilliard equation, the well-know backward Euler scheme is conditionally energy-stable (with time-step size restrictions) [27,38]. Here, energy-stability (or gradient-stability) of a time-discrete scheme means that the free energy is dissipated in time, analogous to the time-continuous evolution.

A fundamental unconditionally gradient-stable implicit-explicit (IMEX) first-order ac-curate scheme for Cahn-Hilliard equation was established by Elliott and Stuart [29,

(12)

3

Eq. (5.4)] and Eyre [37] in the 1990s, and has been utilized in many other gradient flow problems; see, e.g. [79,42,50] and references therein. The scheme is based on a convex-concave splitting of the free energy functional in which the convex part is treated implicitly and the concave part explicitly. In particular, for a splitting with a quadratic convex part, the resulting system is linear.

Similar ideas have been developed to stabilize the second-order trapezoidal or Crank– Nicolson schemes, where a secant form of the nonlinearity is employed to achieve uncon-ditional gradient stability, see e.g. [25, 27, 80, 5]. Other alternatives to the secant form can be found in the work of Gomez and Hughes [42,43] using a standard Crank–Nicolson scheme with additional stabilization, and the work of Hu, Wise, Wang and Lowengrub [52] and Shen, Wang, Wang and Wise [67] using second-order splitting-based schemes. Unfor-tunately, although these second-order time-accurate methods are provably energy-stable, the resulting schemes are inherently nonlinear.

In Chapter3, we resolve this by presenting new unconditionally energy-stable second-order linear schemes for phase-field models that combine a particular convex-concave splitting of the free energy with an artificial-diffusivity stabilization. To maintain second-order accuracy, the convex part is treated by an implicit Taylor approximation and the concave part by an explicit Taylor approximation. Then, an artificial-diffusivity stabiliza-tion is introduced to proof uncondistabiliza-tional energy-stability. The case of variable mobility is treated using extrapolation. Remarkably, for suitable convex-concave splitting, the result-ing schemes are indeed linear. Here, we apply our ideas to the Cahn–Hilliard equation, and a phase-field model of tumor-growth. Our idea can also be applied to many other models (e.g., phase-field crystal equation: see [75] which is based on our work).

D

UALITY

-

BASED A POSTERIORI ERROR ESTIMATES

To assess the accuracy of numerical simulations and to control the discretization error, there is an obvious need to develop a posteriori error estimates and corresponding space-time adaptive algorithms for fully discrete approximation methods for phase-field mod-els. In particular, we focus on time discretizations that rely on the use of unconditionally gradient-stable IMEX time-stepping schemes.

Over the years, rigorous theory of a posteriori error analysis has been established for the linear parabolic case, see for example the survey by Ainsworth and Oden [1] and Verfürth [74]. However, the theory is much less satisfactory for (strongly) nonlinear parabolic equations. The accuracy of error estimates in the nonlinear case, which per-tain to global bounds in the energy norm or L2norm, is affected by a problem-dependent pre-multiplication constant (the notorious stability constant). For strongly nonlinear

(13)

prob-lems, such as the Cahn–Hillard equation, it can be very hard to obtain quantitatively-accurate estimates because the invoked bounds typically consider worst-case scenarios, leading to huge stability constants [53,9,8,69].

In the late 90’s so-called goal-oriented a posteriori error estimates were developed to estimate the error with respect to user specified quantities of interest. The dual-weighted-residual method (see, [10, 12, 61, 64] for an overview) is a well-known framework for de-riving such error estimates. During the last decade, a variety of dual-based error estimates for parabolic equations were published, see for instance Eriksson and Johnson [32,33,34], Schmich and Vexler [66], Carey et al. [20], Bermejo and Carpio [13], Braack et al. [18], Besier and Rannacher[14], and Asner et al. [4]. Duality-based error estimates in general provide good indicators of the errors for nonlinear problems. However, the underlying dis-cretizations are mostly in the context of space-time (discontinuous) Galerkin finite element discretization. There is very little work on rigorous a posteriori estimates and adaptivity for nonlinear parabolic equations discretized using IMEX time-stepping schemes.

In Chapter4–6, we present a methodology for duality-based a posteriori error estima-tion for fully discretized semi-linear parabolic PDEs using conforming finite elements in space and first-order IMEX schemes in time, and subsequently applied to the linear heat equation, the Allen-Cahn equation and the diffuse-interface tumor-growth model. The main result in our work is a decomposition of the error estimate for which the temporal and the spatial discretization errors are addressed separately. The decomposition is based on adapting the residual decomposition by Verfürth [74] in the context of norm-based a posteriori error analysis to our duality-based error representation through a specifically-tailored IMEX time-discrete dual problem. We further employ the two-level methodology [69] to obtain computable error indicators. Based on the indicators, a general space-time adaptive algorithm is proposed for adaptive mesh refinement and adaptive time-step se-lection.

We have organized the contents of the thesis as follows: in Chapter2 we review the gradient-flow structure of phase-field models. The unconditionally energy-stable schemes for phase-field models are presented in Chapter3. In Chapter4we present a methodology for duality-based a posteriori error estimation for nonlinear parabolic PDEs. The associ-ated adaptive algorithm is proposed in Chapter5. In Chapter 6, the performance of the error estimates and the proposed adaptive algorithm is demonstrated on three different applications: the linear heat equation, the Allen–Cahn equation and a nonlinear parabolic system relevant to tumor growth. Concluding remarks are drawn in Chapter 7, where we also state some open questions. Additionally in the Appendix A we present a brief descrip-tion of our fully-discrete scheme for the diffuse-interface tumor model of Secdescrip-tion 3.3.

(14)

2

P

HASE

-

FIELD MODELS AS GRADIENT FLOWS

A gradient flow is a process that follows the path of steepest descent in an energy landscape. Starting from an initial state, the evolution process is driven by an energy in the direction for which the energy decreases along solutions as fast as possible until arriving at a local minimum. Such formalism provides a natural framework to describe various important phenomena of thermodynamic systems. A wide range of well-established evolution equa-tions can be interpreted as gradient flows, such as convection-diffusion equaequa-tions, reactive moving-boundary problems, phase field models and so on, see e.g. [76, 59, 62, 26] and references therein.

In this chapter, we first give an overview about how a gradient flow is determined by the driving energy and the dissipation mechanism. We then illustrate the gradient-flow structure of three phase-field models: the classical Allen–Cahn and Cahn–Hillard equations and a phase-field model of tumor growth. The derivation of this chapter mainly follows the work in [44].

2.1.

A

BSTRACT GRADIENT FLOWS

LetV be a Hilbert space with norm k·k and inner product (·,·) and consider a continuously differentiable functional E :V → R. The gradient flow associated to energy functional E is given by the evolution equation

∂tu = −KδE(u)

(15)

where dissipation K is a positive-definite operator, i.e. (K x, x) > 0 for all x 6= 0, andδE(u)

δu is

the variational derivative of E , defined through the Gâteaux derivative: µδE(u) δu , v= E0(u)(v) = lim s→0 E (u + s v) − E (u) s , (2.2)

for all u, v ∈ V satisfying suitable homogeneous boundary conditions. The evolution of u in (2.1) is driven by the energy functional E in the direction opposite to KδE(u)

δu , which

indicates the following dissipation mechanism: d dtE (u(t )) = µδE(u(t)) δu ,∂tu(t ) ¶ = − µδE(u(t)) δu , K δE(u(t)) δu ¶ ≤ 0. The most classical example is that of the heat equation

∂tu = ∆u (2.3)

subject to initial conditions u(0) = u0 and suitable boundary conditions, for example,

∇u · n = 0 on ∂Ω, which can be interpreted as the gradient flow in L2(Ω) for the quadratic functional E (u) :=1 2 Z Ωk∇uk 2dx. (2.4)

The Gâteaux derivative of E is

E0(u)(v) = lim s→0 E (u + s v) − E (u) s = Z Ωlims→0 1 2|∇(u + s v)| 21 2|∇u| 2 s dx = (∇u, ∇v), ∀v ∈ H1(Ω). Integration by parts, and using∂nu = 0 gives

E0(u)(v) =

µδE(u)

δu , v

= (−∆u, v), ∀v ∈ H1(Ω),

in other words, the variational derivative isδE(u)

δu = −∆u, from which we see that the heat

equation (2.3) is equivalent to the gradient flow (2.1) for K = 1. In particular, if u is the solution to (2.3), then we have the following energy dissipation:

d dtE (u) = E 0(u)( tu) = − ¡ ∆u,∂tu¢ = −k∆uk2≤ 0.

Note that the definition of the gradient flow for (2.3) is not unique. The choice of the energy as well as the dissipation mechanism allows for various formulations [26,62]. For example, the heat equation (2.3) can also fit into the gradient-flow structure (2.1) by setting K := −∆(·) for the energy functional defined by E(u) :=1

2 Z

u 2

dx with the variational derivative

δE(u)

(16)

2.2.ALLEN-CAHN ANDCAHN-HILLIARD EQUATION 7

2.2.

A

LLEN

-C

AHN AND

C

AHN

-H

ILLIARD EQUATION

The principal characteristic of phase-field models is the diffuse-interface assumption. In-stead of sharp interfaces, the interface is described by a steep, but smooth, phase transition. A second characteristic of phase-field models is the dependence on the Helmholtz free en-ergy. The evolution of the phase field follows a relaxation process which minimizes the total free energy of the system. This links directly to the gradient-flow structure.

The Allen-Cahn equation and the Cahn-Hilliard equation are the key examples of phase-field models. For their derivation, we start by assuming a material continuum body B composed of two constituent species BαandBβ, which completely fill an open region

Ω ⊂ Rd. Any arbitrary partP ⊂ Ω is allowed to be occupied by the constituents

simultane-ously. For the sake of simplicity, we suppose that the material body is not undergoing any deformation, and the density of the mixtureρ in Ω is constant and equal to 1, i.e. ρ = 1.

Let m be the total mass of the mixture, i.e. m =R

ρ dx. And let mα, mβbe the masses

of each species such that m = mα+ mβ. We denote byρα,ρβthe mass of each constituent

per unit volume of the mixture, namely

mα= Z

ραdx, =

Z

ρβdx.

The local mass fraction or concentration of each constituent is then defined by

cα=ρα

ρ , =

ρβ

ρ .

For a general mixture, each of the constituents must satisfy its own mass balance law as follows d dt Z Pρcidx = Z ∂P−hi· n dS + Z P γidx, i = α,β, (2.5)

where hi is the mass flux andγi is the mass supply. By recalling thatρ is assumed to be

equal to 1 (and the material velocities are assumed to be negligible for simplicity), we con-clude from (2.5)

∂tci+ div hi = γi, i = α,β. (2.6)

To distinguish one phase from the other one, an additional scalar functionϕ, known as phase field or order parameter, is introduced as a label for each constituent. Such a func-tion allows a smooth transifunc-tion between phases in the interfacial regions, which may be interpreted as the difference of local mass fraction of the two constituents, namely

ϕ = cα− cβ=        1, in phaseα, −1, in phaseβ,

(17)

We thus have the following conservation equation:

∂tϕ + divh = γ (2.7)

where h = hα− hβandγ = γα− γβ.

Then, we consider another fundamental property of phase-field models: energy. The Allen–Cahn and Cahn–Hilliard equations are driven by the Ginzburg–Landau energy func-tional E which is defined as

E (ϕ) :=Z Ω ³ Ψ(ϕ) +²2 2|∇ϕ|dx,

where² > 0 is the equilibrium interface thickness and Ψ(ϕ) is a double-well potential func-tion with wells atϕ = 1 and ϕ = −1. The Gâteaux derivative of E(ϕ) is

E0(ϕ)(v)= lim s→0 E¡ ϕ + s v¢ − E ¡ϕ¢ s = Z Ωlims→0 Ψ(ϕ + s v) +²22|∇(ϕ + s v)|2− Ψ(ϕ) −²22|∇ϕ|2 s dx = Z Ω ³ Ψ0(ϕ)v + ²2 ∇ϕ · ∇v ´ dx, ∀v ∈ H1(Ω). After applying integration by parts, it follows

E0(ϕ)(v) = ¡Ψ0(ϕ) − ²2∆ϕ,v¢ + Z

∂Ω²

2

∇ϕ · n v dS, ∀v ∈ H1(Ω) Thus, the variational derivative of E (ϕ) is identified as

δE(u)

δϕ = Ψ0(ϕ) − ²

2∆ϕ.

The so-called chemical potential is defined as the variational derivative

µ :=δE(u)

δϕ = Ψ0(ϕ) − ²

2∆ϕ.

Following the thermodynamical framework of phase-field models introduced by Gurtin (see [46] for more details), the system is closed by appropriately choosing the constitutive equations in equation (2.7), which must be consistent with the dissipation inequality such that the following constraint holds:

(18)

2.2.ALLEN-CAHN ANDCAHN-HILLIARD EQUATION 9

THEALLEN-CAHN EQUATION

To satisfy the required constraints, the constitutive equations of (2.7) can be set to be h = 0 andγ = −µ, which leads to the classical Allen-Cahn equation

∂tϕ = ²2∆ϕ − Ψ0(ϕ). (2.8)

We see that the Allen-Cahn equation (2.8) is indeed a gradient flow for the Ginzburg– Landau energy with K = 1, i.e.

∂tϕ = −µ

with the following energy dissipation: d dtE (ϕ(t)) = E 0(ϕ)(∂ tϕ) = (µ,∂tϕ) + Z ∂Ω² 2 ∇ϕ · n ∂tϕ dS = −kµk2+ Z ∂Ω² 2 ∇ϕ · n ∂tϕ dS

Note that the second termR

∂Ω²2∇ϕ · n ∂tϕ dS is associated with a nonstandard working

caused by internal forcesξ = ²2∇ϕ conjugate to changes in ϕ [46,44]. If we take homoge-neous Neumann boundary conditions, i.e. ∇ϕ · n = 0 on ∂Ω, equation (2.8) guarantees the energy decay

d

dtE (ϕ(t)) = −kµk

2

≤ 0.

THECAHN–HILLIARD EQUATION

Taking the consistent constitutive equations h = −M(ϕ)∇µ with mobility M(ϕ) > 0 and

γ = 0 in (2.7) gives the Cahn–Hilliard equation

∂tϕ = div

³

M (ϕ)∇¡ − ²2∆ϕ + Ψ0(ϕ)¢´

which is equivalent to the gradient flow formulation

∂tϕ = −K µ

for K := −div(M(ϕ)∇(·)). The energy functional E decreases along ϕ as d dtE (ϕ(t)) = E 0(ϕ)(∂ tϕ) = (µ, ∂tϕ) + Z ∂Ω² 2 ∇ϕ · n ∂tϕ dS = ³ µ,div¡M(ϕ)∇µ¢´+ Z ∂Ω² 2 ∇ϕ · n ∂tϕ dS = −¡∇µ, M(ϕ)∇µ¢ + Z ∂Ω n M (ϕ)∇µ · nµ + ²2∇ϕ · n ∂ o dS

(19)

We note that the second termR∂ΩM (ϕ)∇µ · nµ dS = −R∂Ωh · nµ dS can be understood as a free-energy flux through∂Ω. If we take homogeneous Neumann boundary conditions for bothϕ and µ, i.e. ∇ϕ · n = 0 and M(ϕ)∇µ · n = 0 on ∂Ω, energy dissipation is guaranteed:

d

dtE (ϕ(t)) = −¡∇µ,M(ϕ)∇µ¢ ≤ 0.

2.3.

A

PHASE

-

FIELD MODEL OF TUMOR GROW TH

In this section, we introduce a more complicated phases-field model, viz. a diffuse-interface avascular-tumor-growth model proposed by Hawkins-Daarud et al [47,50]. The derivation of the governing equations is based on the continuum theory of mixtures, which starts from the assumption that the system under consideration is an isothermal mixture of four constituents (i = 1,2,3,4) with zero velocity:

• B1: tumor cells • B2: healthy cells

• B3: nutrient-rich extracellular water • B4: nutrient-poor extracellular water

Each constituent is assigned a mass densityρi which is the mass of constituent i per unit

volume of the mixture. The mass fraction or concentration of constituent i is given by

ci=ρi

ρ

whereρ is the density of the mixture. For each constituent, the mass balance of the form

∂t

¡

ρci¢ + divhi= γi, i = 1,2,3,4 (2.9)

holds, where hi andγi represent the mass flux and the mass supply, respectively, which

can be interpreted as the mass exchange between constituent i and the rest of the mixture due to diffusion and chemical reactions.

Let us assume for simplicity that the density of mixtureρ is constant and equal to 1, i.e.

ρ = 1. Hence, the conservation equation (2.10) simplifies to

∂tci+ div hi= γi, i = 1,2,3,4 (2.10)

We further assume that the mixture is saturatedP

ici = 1, and the concentration of

cell-phases remains a fixed fraction ¯ccel l∈ (0, 1), i.e.

(20)

2.3.APHASE-FIELD MODEL OF TUMOR GROWTH 11

which imply that

∂tc1+ ∂tc2= 0, ∂tc3+ ∂tc4= 0.

Thus, the equations of the mass balance (2.10) for each phase can be reduced to only two equations: the equations for the tumor phaseϕ := c1and nutrient concentration n := c3:

∂tϕ + divh1= γ1, (2.11a)

∂tn + divh3= γ3. (2.11b)

As in the Allen–Cahn and Cahn–Hilliard equation, the evolution of the system is also based on the total free energy, which can be written as follows:

E (ϕ,n) := Z Ω µ f (ϕ) +² 2 2|∇ϕ| 2 + 1 2σn 2 ¶ dx, (2.12)

where f (ϕ) is a double-well potential function with wells at ϕ = 0 and ϕ = ¯ccel l andσ > 0 is

the reversibility parameter that controls the relative strength of the interaction among the cell and nutrient species. The Gâteaux derivative of E (ϕ,n) with respect to ϕ and n is

E0ϕ(ϕ,n)(v) = ¡f0(ϕ) − ²2∆ϕ,v¢ + Z ∂Ω² 2 ∇ϕ · n v dS, ∀v ∈ H1(Ω), E0n(ϕ,n)(v) = ³n σ, v ´ , ∀v ∈ L2(Ω).

The relevant chemical potentialsµ1andµ3, which are associated withϕ and n, respectively,

are defined as the corresponding variational derivatives:

µ1:=δE(ϕ,n) δϕ = f0(ϕ) − ² 2∆ϕ, µ 3:= δE(ϕ,n) δn = n σ. (2.13)

To appropriately close this system, the constitutive equations in (2.11) must be con-sistent with both the classical balance laws for the entire mixture and the second law of thermodynamics for the mixture. In particular, they have to comply with the following constraint (see Ref. [47,50]):

¡

µ1γ1+ µ3γ3¢ + ¡∇µ1· h1+ ∇µ3· h3¢ ≤ 0.

Choosing the constitutive equations

h1= −M(ϕ)∇µ1

h3= −σD(n)∇µ3

γ1= σP (ϕ)(µ3− µ1)

(21)

for a mobility function M (ϕ) > 0, a diffusivity function D(n) > 0 and a nonnegative

prolifer-ation function P (ϕ), leads to a Cahn–Hilliard reaction-diffusion system for tumor-growth:

∂tϕ = div¡M(ϕ)∇µ1¢ + P(ϕ)(n − σµ1) (2.14a)

∂tn = div¡D(n)∇n¢ − P(ϕ)(n − σµ1) (2.14b)

Let c := [c1, c3]T = [ϕ, n]T denote the variables, and letµ := [µ1,µ3]T denote the chemical

potentials. The PDE system in (2.14) fits into the generalized gradient-flow structure (see Ref. [48]): ∂tc = −K µ for K := " − div¡M(ϕ)∇(·)¢ 0 0 − div¡σD(n)∇(·)¢ # + σP (ϕ) " 1 −1 −1 1 #

which satisfies the dissipation law: d dtE (ϕ,n) = E 0 ϕ(ϕ,n)(∂tϕ) + En0(ϕ,n)(∂tn)µ1,∂tϕ¢ + Z ∂Ω² 2 ∇ϕ · n ∂tϕ dS + ¡µ3,∂tn ¢ = ³ µ1, div¡M(ϕ)∇µ1¢ + P(ϕ)(n − σµ1) ´ + Z ∂Ω² 2 ∇ϕ · n ∂tϕ dS + ³ µ3, div¡D(n)∇n¢ − P(ϕ)(n − σµ1) ´ = −¡∇µ1, M (ϕ)∇µ1¢ − 1 σ¡∇n,D(n)∇n¢ − Z ΩσP(ϕ)(µ3− µ1) 2dx + Z ∂Ω n − h1· n µ1− h3· n µ3+ ²2∇ϕ · n ∂ o dS. Note that If we take homogeneous Neumann boundary conditions forϕ, µ1 andµ3, i.e.

∇ϕ · n = 0, M(ϕ)∇µ1· n = 0 and σD(n)∇µ3· n = 0 on ∂Ω, energy dissipation is guaranteed:

d dtE (ϕ,n) = −¡∇µ1, M (ϕ)∇µ1¢ − 1 σ¡∇n,D(n)∇n¢ − Z ΩσP(ϕ)(µ3− µ1) 2 dx ≤ 0.

(22)

3

E

NERGY

-

STABLE TIME

-

STEPPING SCHEMES

Phase-field models have a structure of gradient flows for which the systems dissipate the total free energy. For the stability of the computations, it is essential to develop appro-priate time-stepping schemes for phase-field models that inherit energy dissipation at the discrete level.

In this chapter, we present unconditionally energy-stable second-order time-accurate linear schemes for phase-field models; in particular, in Section3.2we consider the Cahn– Hilliard equation, while in Section3.3a phase-field tumor-growth model consisting of a reactive Cahn–Hilliard equation and a reaction-diffusion equation. Apart from the new second-order schemes, we also present well-known first-order variants for the sake of com-pleteness. In Section3.4we present numerical results for fully discrete schemes employ-ing mixed finite-elements for spatial discretization. The numerical results fully verify the second-order accuracy and stability of the proposed schemes.1

3.1.

M

OTIVATION

Diffuse-interface models are characterized by a non-convex (bulk) free energyΨ(u) (for a binary system with phase u), and a distributed interface energy²22|∇u|2, both of which are dissipated during evolution in time (given by a typically fourth-order, nonlinear parabolic system).

Time-discretization methods for diffuse-interface models have to deal with the nonlin-ear termΨ0(u) originating from the nonconvex (double-well in case of a binary system) free energyΨ(u). The nonlinearity allows for “backwards” diffusion, which is only mildly 1This chapter is based on [82].

(23)

regularized by the higher-order operator coming from the interfacial energy, so that naive time-discrete schemes can be unstable. This stability issue is the motivation for a large amount of literature on the development of schemes that are provably energy-stable. Here energy-stability (or simply stability) of a discrete scheme means that the total energy is dis-sipated in time, analogous to the continuous evolution. Note that energy-stability is also referred to as gradient-stability in the case of a gradient-flow system [71].

The nontrivial stability of first-order time-accurate schemes has been studied since the late 1980s, mostly in the context of the Cahn–Hilliard equation. For example, Elliott [27, Eq. (4.10b)] established the remarkable conditional stability (time-step restriction) of the backward Euler scheme. A fundamental unconditionally-stable semi-implicit scheme, that is often attributed to Eyre [36,37] but actually goes back to Elliott and Stuart [29, Eq. (5.4)], involves the convex–concave splitting of the nonconvex free energy:Ψ(u) = Ψc(u) − Ψe(u).

In the scheme the convex partΨ0c(u) is treated implicitly while the concave partΨ0e(u) is treated explicitly. The concise proof of its stability is based on an elegant convexity ar-gument which can be found in [29, proof of Proposition 5.9] and has been re-discovered in [79]. An additional advantage of the splitting scheme is that if the convex part is quadratic, then the resulting scheme is linear, i.e., it requires the solution of a linear sys-tem at each time-step.2An alternative to splitting is to add stabilization, similar to classical artificial-viscosity stabilization. For an explicitly-treated nonlinearity this has been studied in [51,68].

Second-order time-accurate methods have mostly been studied for variants of the

Crank–Nicolson scheme, as standard treatment of the nonlinear termΨ0(u) is not suffi-cient for stability. The earliest unconditionally-stable second-order scheme was studied by Du and Nicolaides [25] and employs a secant form of the nonlinearity; see also [27]. The secant form is very powerful for stability, however the resulting scheme is inherently nonlinear. Applying an implicit Taylor approximation of the secant form also results in an unconditionally-stable (and nonlinear) scheme; see Kim, Kang and Lowengrub [54]. Unfortunately, an explicit Taylor approximation of the secant form, although leading to a linear scheme, is only conditionally stable; see the recent work by Guillén-González and Tierra [45]. We note that nonlinear solvers are described in, e.g., [80,5], while extensions of the secant-based schemes to systems with more than two phases can be found in [55,17]. An alternative to the secant form has recently been proposed by Gomez and Hughes [42]; see also [43]. Their unconditionally-stable nonlinear scheme can be interpreted as a stan-dard Crank–Nicolson scheme with additional stabilization.

2A quadratic convex part is always possible if the double-well free energy has at most quadratic growth. Of course, for so-called regularized or truncated free energies this can be satisfied by construction.

(24)

3.2.THECAHN–HILLIARD EQUATION 15

Another seemingly natural alternative would be to search for second-order splitting-based schemes. This idea has been pursued by Hu, Wise, Wang and Lowengrub [52] and by Shen, Wang, Wang and Wise [67]. To maintain second-order accuracy, their idea is to treat the convex part by a secant form, and the concave part by extrapolation (e.g., a second-order backward difference). Unfortunately, although the energy for this scheme is always smaller than the initial energy, it need not decay at each time-step; a property which was referred to as weak energy-stability.

In the current work, we consider new second-order schemes that combine a particular second-order splitting with stabilization. We apply our ideas to the Cahn–Hilliard equa-tion, and the diffuse-interface tumor-growth system proposed in [50], which consists of a reactive Cahn–Hilliard equation and a reaction–diffusion equation. The particular split-ting maintains second-order accuracy for the nonlinear termΨ0(u) by treating the convex part by an implicit Taylor approximation and the concave part by an explicit Taylor ap-proximation. Then, to be able to proof unconditional energy-stability, artificial-diffusivity stabilization is crucial. Remarkably, because of splitting, the resulting scheme is also lin-ear (for a quadratic convex part). Furthermore, all properties are maintained for variable mobility if treated through extrapolation; a trick that was also recognized in [52]. For the tumor-growth system, we are also able to prove unconditional energy-stability, but only by adding additional (artificial-convexity) stabilization. The reactive terms are treated semi-implicitly, combining the mid-point method with extrapolation. In this manner, the result-ing scheme for the tumor-growth system is also linear.

3.2.

T

HE

C

AHN

–H

ILLIARD EQUATION

3.2.1.

S

TRONG FORMUL ATION

The key equation in diffuse-interface modeling is the Cahn–Hilliard equation, given by: find u :Ω × [0,T ] → R satisfying

ut = ∇ ·£M(u)∇¡Ψ0(u) − ²2∆u¢¤, for (x, t ) ∈ Ω × (0,T ]

∇u · n = M(u)∇(Ψ0(u) − ²2∆u) · n = 0, for (x,t) ∈ ∂Ω × (0,T ]

u(x, 0) = u0, for x ∈ Ω

where u represents the phase and at the same time serves to model the interface,² > 0 is the equilibrium interface thickness,Ψ(u) is the (bulk) free-energy density function, M(u) ≥ 0 is the mobility function, and u0is the initial condition. It is assumed that the domainΩ

(25)

simplicity we consider natural boundary conditions; n denotes the exterior unit normal. Note that we employed the notation (·)t= ∂(·)/∂t .

The nonlinear free-energy density functionΨ(u) is a double well potential. We consider the following C2,1-continuousΨ:

Ψ(u) :=              (u + 1)2 u < −1 1 4(u 2 − 1)2 u ∈ [−1,1] (u − 1)2 u > 1 (3.1)

One can view this function as a regularized (or truncated) version of the classical quartic function 14(u2− 1)2. For constant mobility M (u) = ˆM > 0 while for variable mobility, the

function M (u) is typically defined as:

M (u) :=    ˆ M (1 − u2) u ∈ [−1,1] 0 otherwise (3.2)

Let us briefly recall some well-known results concerning the Cahn–Hilliard equation. Since it is a fourth-order nonlinear parabolic equation, it can be formulated as a sys-tem of two coupled second-order equations. In spatial discretizations, this reformula-tion avoids the discretizareformula-tion of higher-order derivatives. The system is given by: find (u,µ) : Ω × [0,T ] → R such that

ut = ∇ ·¡M(u)∇µ¢, for (x, t) ∈ Ω × (0,T ]

µ = Ψ0(u) − ²2∆u, for (x,t) ∈ Ω × (0,T ]

∇u · n = ∇µ · n = 0, for (x, t ) ∈ ∂Ω × (0,T ]

u(x, 0) = u0, for x ∈ Ω

We note thatµ is called the chemical potential.

In the above Cahn–Hilliard model, there are two important global properties. One is the conservation of mass. Applying the no flux boundary condition, the total amount of phase inΩ must be always equal to the given original amount, i.e.

d dt Z Ωu d x = Z Ωut d x = Z Ω∇ ·¡M(u)∇µ¢dx = Z ∂ΩM (u)∇µ · ndS = 0 (3.3)

The second property is the dissipation of energy. The total free energy is defined by

E (u) := Z Ω µ²2 2|∇u| 2 + Ψ(u)d x. (3.4)

(26)

3.2.THECAHN–HILLIARD EQUATION 17

The Cahn–Hilliard equation is a gradient flow of this functional. The relaxation of u is driven by the local minimization of the total free energy subject to the above conservation constraint [40]. Therefore, for all times, the total energy is always non-increasing, i.e.

dE (u(t )) dt = E 0(u)(u t) = Z Ω ¡ Ψ0(u)u t+ ²2∇u · ∇ut ¢ = Z Ωµut = Z Ωµ∇ · ¡M(u)∇µ¢ = − Z ΩM (u)|∇µ| 2 ≤ 0

3.2.2.

F

IRST

-

ORDER ACCURATE SCHEME

Due to both the bilaplacian operator and the nonlinearity, the Cahn–Hilliard equation is stiff. A fine discretization in space is needed to meet the high resolution requirement for diffuse interfaces, while a nontrivial time-discretization is required to have a stable scheme that is independent of the time-step size. In this paper, we focus on the time-discretization. Therefore we present our ideas without discretizing in space. Of course, one may obtain fully discrete schemes by using finite-difference of finite-element methods in space; see for example [28,24,7,78]. We furthermore assume sufficient spatial and temporal regularity.

Elliott and Stuart [29, Eq. (5.4)] and Eyre [37] proposed a first-order accurate uncon-ditional energy-stable scheme for gradient-flow systems based on the splitting of E into a convex (contractive) and concave (expansive) part, i.e.

E = Ec− Ee

where both Ec and Eeare convex. In the scheme the convex part is treated implicitly and

the concave part explicitly, i.e.

uk+1− uk

τ = ∆ ˜µ (3.5a)

˜

µ = DuEc(uk+1) − DuEe(uk) = −²2∆uk+1+ Ψ0c,k+1− Ψ0e,k (3.5b)

withτ the time-step size, and M(u) has been set to 1 (see Sec.3.2.3for the variable mobility case). Note that the splitting of E induces a splitting ofΨ. For convenience, throughout this work, a subscript onΨ shall refer to its argument, for example, Ψ0

c,k+1≡ Ψ0c(uk+1).

The properties of the above scheme are analyzed in [29, 36, 79]. It is shown that the scheme is energy-stable, that is, it inherits the dissipation property of the Cahn–Hilliard equation in that the discrete total free energy always dissipates analogous to the time-continuous case:

(27)

It is also shown that the splitting ofΨ always exists, but it is not unique. In particular, for a splitting with a quadratic convex part, the resulting system is linear. For ourΨ, we may choose: Ψ = Ψc− Ψe=                  µ u2+1 4 ¶ − µ −2u −3 4 ¶ u < −1 µ u2+1 4 ¶ −µ 3 2u 2 −1 4uu ∈ [−1,1] µ u2+1 4 ¶ − µ 2u −3 4 ¶ u > 1

3.2.3.

S

ECOND

-

ORDER ACCURATE SCHEME

Following the idea of splitting, we propose a new second-order accurate, unconditionally energy-stable, linear scheme:

uk+1− uk τ = ∇ · ¡ ˜ Mk+1/2∇ ˜µ¢ (3.6a) ˜ µ = ˜Ψ0(u k, uk+1) − ²2∆ uk+1+ uk 2 − α1τ∆(uk+1− uk) (3.6b) where ˜ Mk+1/2:= Mµ 3 2uk− 1 2uk−1 ¶ ˜ Ψ0(u k, uk+1) := Ψ0c,k+1uk+1− uk 2 Ψ 00 c,k+1− Ψ0e,kuk+1− uk 2 Ψ 00 e,k

This scheme is a modification of the Crank–Nicolson method which includes splitting, stabilization and extrapolation. In particular ˜Ψ0(uk, uk+1) is a novel second-order

accu-rate splitting employing an implicit Taylor expansion of the convex part Ψ0c and an ex-plicit Taylor expansion of the concave partΨ0e. Stabilization is introduced in the form of

α1τ∆(uk+1− uk), withα1≥ 0. The stabilization can be thought of as artificial diffusivity.

Finally, to allow for a linear scheme, a second-order extrapolation ˜Mk+1/2is employed for the mobility M (u). For the time step, we use a first-order extrapolation, or equivalently, set u−1= u0.

We summarize the properties of the above scheme in the following theorem.

Theorem 3.2.A Let the free energy densityΨ be C2,1-continuous and have a convex split-tingΨ = Ψc− Ψewith C2,1-continuousΨcandΨeand finite second derivatives, i.e., |Ψ00c| ≤

Lc and |Ψ00e| ≤ Lefor some constants Lc, Le≥ 0. Let the mobility M be C0,1-continuous and

satisfy 0 ≤ M ≤ ¯M . If the stabilization is large enough, i.e. α1≥ ¯M (Lc+ Le)2/16, then the

time-stepping scheme (3.6) has the following properties: 1. Unconditional energy-stability: E (uk+1) ≤ E(uk)

(28)

3.2.THECAHN–HILLIARD EQUATION 19 2. Mass conservation: Z Ωuk+1d x = Z Ωu0d x 3. Second-order accuracy. 

Proof 1. Since both functionΨc andΨe are convex, the property of the convex function

leads to the inequalities

Ψc(φ) − Ψc(ϕ) ≥ Ψ0c(ϕ)(φ − ϕ)

Ψe(ϕ) − Ψe(φ) ≥ Ψ0e(φ)(ϕ − φ)

Subtracting the two equations yields

Ψ(ϕ) − Ψ(φ) = Ψc(ϕ) − Ψe(ϕ) − (Ψc(φ) − Ψe(φ))

≤ −Ψ0c(ϕ)(φ − ϕ) − Ψ0e(φ)(ϕ − φ)

=¡Ψ0c(ϕ) − Ψ0e(φ)¢(ϕ − φ)

(3.7)

Base on the above inequality (3.7), we have

E (uk+1) − E(uk) = Z Ω ½ Ψ(uk+1) − Ψ(uk) +² 2 2(|∇uk+1| 2 − |∇uk|2) ¾ ≤ Z Ω ½ ³ Ψ0 c,k+1− Ψ0e,k ´ (uk+1− uk) +² 2 2(|∇uk+1| 2 − |∇uk|2) ¾

Next, using the time-stepping scheme (3.6b) and integration by parts, it follows that

E (uk+1) − E(uk) ≤ Z Ω n ³ ˜ µ +uk+1− uk 2 (Ψ 00 c,k+1+ Ψ00e,k) ´ (uk+1− uk) −α1τ∇(uk+1− uk) · ∇(uk+1− uk) −² 2 2∇(uk+1+ uk) · ∇(uk+1− uk) +² 2 2(|∇uk+1| 2 − |∇uk|2) o = Z Ω n³ ˜ µ +uk+1− uk 2 (Ψ 00 c,k+1+ Ψ00e,k) ´ (uk+1− uk) o − α1τk∇(uk+1− uk)k2L2 (3.8) Using that Z Ω| f (x)g (x)| ≤ supx∈Ω|g (x)| Z

| f (x)|, the following inequality holds:

Z Ω(Ψ 00 c,k+1+ Ψ00e,k)(uk+1− uk) 2 ≤ (sup R Ψ 00 c,k+1+ sup R Ψ 00 e,k) Z Ω(uk+1− uk) 2 (3.9) Moreover, according to the definition (3.2), the mobility function M is non-negative and bounded, i.e. 0 < ˜Mk+1/2≤ sup R M µ 3 2uk− 1 2uk−1 ¶ = ¯M (3.10)

(29)

From (3.6a), (3.9) and integrating by parts, and then using Cauchy-Schwarz inequality, we obtain Z Ω n³ ˜ µ +uk+1− uk 2 (Ψ 00 c,k+1+ Ψ00e,k) ´ (uk+1− uk) o ≤ Z Ωµ(u˜ k+1− uk) + 1 2(supR Ψ 00 c,k+1+ sup R Ψ 00 e,k) Z Ω(uk+1− uk) 2 = −τ Z Ω ˜ Mk+1/2∇ ˜µ · ∇ ˜µ −τ 2(supR Ψ 00 c,k+1+ sup R Ψ 00 e,k) Z Ω ˜ Mk+1/2∇ ˜µ · ∇(uk+1− uk) ≤ −τ Z Ω ˜ Mk+1/2|∇ ˜µ|2+τ 2(supR Ψ 00 c,k+1+ sup R Ψ 00 e,k) s Z Ω ˜ Mk+1/2|∇ ˜µ|2 s Z Ω ˜ Mk+1/2|∇(uk+1− uk)|2 (3.11) Finally, combining the above inequalities (3.8), (3.11), (3.10) and applying Young’s inequal-ity, 2 f g ≤ β f2+ β−1g2, forβ > 0, we have

E (uk+1) − E(uk) ≤τ 4(supR Ψ 00 c,k+1+ sup R Ψ 00 e,k) Z Ω µ ¯ Mβ|∇(uk+1− uk)|2+ ˜ Mk+1/2 β |∇ ˜µ| 2 ¶ − α1τk∇(uk+1− uk)k2L2− τ Z Ω ˜ Mk+1/2|∇ ˜µ|2 ≤ −τ µ 1 − 1 4β(supR Ψ 00 c,k+1+ sup R Ψ 00 e,k) ¶ Z Ω ˜ Mk+1/2|∇ ˜µ|2 − τ µ α1− ¯ 4 (supR Ψ 00 c,k+1+ supR Ψ00e,k) ¶ k∇(uk+1− uk)k2L2

In order to have E (uk+1) ≤ E(uk) for any time-step sizeτ > 0, it follows that

β ≥1 4(supR Ψ 00 c,k+1+ sup R Ψ 00 e,k) α1≥ ¯ 4 (supR Ψ 00 c,k+1+ sup R Ψ 00 e,k), (3.12)

which means that the stabilization constantα1need to satisfy

α1≥ ¯ M 16(supR Ψ 00 c,k+1+ supR Ψ00e,k) 2 (3.13)

2. Integrating (3.6a) over the domainΩ, and applying the divergence Theorem and the homogeneous Neumann boundary conditions, we get

Z Ω uk+1− uk τ d x = Z Ω∇ · ¡ ˜ Mk+1/2∇ ˜µ¢dx = Z ∂Ω ˜ Mk+1/2∇ ˜µ · ndS = 0 Then, it follows by induction that

Z

uk+1d x =

Z

(30)

3.2.THECAHN–HILLIARD EQUATION 21

3. To compute the local truncation error, we assume that uk+1is determined given exact values for ukand uk−1. Then

uk+1= u(tk) + τ∇ ·

¡ ˜

Mk+1/2∇ ˜µ(u(tk+1), u(tk))

¢

(3.14) with ˜Mk+1/2= M¡32u(tk) −12u(tk−1)¢. We next Taylor-expand the second term of the

equa-tion. SinceΨ,Ψc andΨeare C2,1-continuous, it follows that

Ψ0

c(u(tk+1)) = Ψ0c(u(tk+1/2)) − Ψc00(u(tk+1))(u(tk+1/2) − u(tk+1)) + O(τ2)

= Ψ0c(u(tk+1/2)) − Ψ00c(u(tk+1))( u(tk+1) + u(tk) 2 − u(tk+1)) + O(τ 2 ) = Ψ0c(u(tk+1/2)) − u(tk) − u(tk+1) 2 Ψ 00 c(u(tk+1)) + O(τ2) (3.15) Ψ0

e(u(tk)) = Ψ0e(u(tk+1/2)) − Ψe00(u(tk))(u(tk+1/2) − u(tk)) + O(τ2)

= Ψ0e(u(tk+1/2)) − u(tk+1) − u(tk) 2 Ψ 00 e(u(tk)) + O(τ2) (3.16)

Taking the difference of (3.15) and (3.16), it follows Ψ0(u(t k+1/2)) = Ψ0c(u(tk+1/2)) − Ψ0e(u(tk+1/2)) = Ψ0c(u(tk+1)) − u(tk+1) − u(tk) 2 Ψ 00 c(u(tk+1)) −Ψ0e(u(tk)) − u(tk+1) − u(tk) 2 Ψ 00 e(u(tk)) + O(τ2) (3.17)

The interface energy term and the mobility function are approximated using the midpoint rule and extrapolation, respectively, both of which are well-known second-order accurate, thus M (u(tk+1/2)) = Mµ 3 2u(tk) − 1 2u(tk−1) ¶ + O(τ2) ²2∆u(t k+1/2) = ²2∆ u(tk+1) + u(tk) 2 + O(τ 2 ) (3.18)

The stabilization terms are O(τ2) because

α1τ∆(uk+1− uk) = α1τ2∆(

uk+1− uk

τ ) = α1τ

2(∆u

t) + O(τ3) = O(τ2) (3.19)

Bringing (3.17), (3.18) and (3.19) together with (3.14) we have

uk+1= u(tk) + τ∇ ·¡M(u(tk+1/2))∇

¡ Ψ0(u(t

k+1/2)) − ²2∆u(tk+1/2)¢¢ +O(τ3) (3.20)

Further, the second order Taylor expansion of the exact value u(tk+1/2) is

u(tk+1/2) = u(tk+1) − τ 2u(t˙ k+1/2) − 1 2 ³τ 2 ´2 ¨ u(tk+1/2) + O(τ3) u(tk+1/2) = u(tk) + τ 2u(t˙ k+1/2) − 1 2 ³τ 2 ´2 ¨ u(tk+1/2) + O(τ3) (3.21)

(31)

Taking the difference of the two equation gives

u(tk+1) − u(tk) = τ ˙u(tk+1/2) + O(τ3)

= τ∇ ·¡M(u(tk+1/2))∇

¡ Ψ0(u(t

k+1/2)) − ²2∆u(tk+1/2)¢¢ +O(τ3)

(3.22)

Finally, subtraction (3.20) from this equation, we immediately obtain that the local trunca-tion error is

u(tk+1) − uk+1= O(τ3) (3.23)

which completes the proof.

3.3.

D

IFFUSE

-

INTERFACE TUMOR

-

GROW TH MODEL

3.3.1.

S

TRONG FORMUL ATION

Let the domainΩ be as in Sec.3.2. We consider the diffuse-interface tumor-growth model proposed in Hawkins-Daarud, van der Zee and Oden [50]: find (u, n) such that

ut = ∆ · (−²2∆u + Ψ0(u)) + γu for (x, t ) ∈ Ω × (0,T ]

nt = ∆ · (

n

δ) + γn for (x, t ) ∈ Ω × (0,T ]

∇u · n = ∇n · n = ∇(−²2∆u + Ψ0(u)) · n = ∇(nδ) · n = 0 for x ∈ ∂Ω × (0,T ]

u(x, 0) = u0, n(x, 0) = n0 for x ∈ Ω

To avoids the discretization of higher-order derivatives, it is chosen here to split the fourth-order differential equation into two second-fourth-order differential equations: find (u,µu, n,µn)

such that

ut= ∆µu+ γu for (x, t ) ∈ Ω × (0,T ]

µu= −²2∆u + Ψ0(u) for (x, t ) ∈ Ω × (0,T ]

nt= ∆µn+ γn for (x, t ) ∈ Ω × (0,T ] µn= n δ for (x, t ) ∈ Ω × (0,T ] ∇u · n = ∇n · n = ∇µu· n = ∇µn· n = 0 for x ∈ ∂Ω × (0,T ] u(x, 0) = u0, n(x, 0) = n0 for x ∈ Ω

whereΨ(u) is the free-energy density function as in Sec.3.2, n denotes the concentration of nutrients,µuandµnare the chemical potentials corresponding to u and n, respectively.

(32)

3.3.DIFFUSE-INTERFACE TUMOR-GROWTH MODEL 23

γuandγnare the chemical reactions defined as following

γu= P (u)(µn− µu), (3.24)

γn= −γu (3.25)

where P (u) ≥ 0 is a nonnegative proliferation function defined as

P (u) :=    δ ˆP (1 − u2) u ∈ [−1,1] 0 otherwise (3.26)

with δ > 0 and ˆP ≥ 0. For the sake of simplicity, we consider homogeneous Neumann

boundary conditions, and constant mobility and diffusion (equal to 1). Furthermore, we do not consider chemotaxis; see [50] for details on chemotaxis.

The total free energy of the tumor-growth model is defined as

E (u, n) := Z Ω µ²2 2|∇u| 2 + Ψ(u) + 1 2δn 2 ¶ (3.27) Similar to the Cahn–Hilliard equation, the tumor-growth model dissipates the total free energy: d dtE (u(t ), n(t )) = −k∇µuk 2 − k∇µnk2− Z ΩP (u)(µn− µu) 2 ≤ 0 which is proved in [50]. Furthermore, the total “mass” is conserved, i.e.

d dt Z Ω(u + n) d x = Z Ω(∆µu+ ∆µn)d x = Z ∂Ω(∇µu· n + ∇µn· n)dS = 0

where we employed the homogeneous Neumann boundary conditions.

3.3.2.

F

IRST

-

ORDER ACCURATE SCHEME

The tumor-growth model is a gradient-flow system. Hence, we can apply the same energy splitting idea:

E (u, n) = Ec(u, n) − Ee(u, n)

= Z Ω µ²2 2|∇u| 2 + Ψc(u) + 1 2δn 2 ¶ − Z ΩΨe(u)

A first-order accurate time-stepping scheme can then be written as

uk+1− uk

τ = ∆ ˜µu+ P (uk)( ˜µn− ˜µu) (3.28a) nk+1− nk

(33)

where ˜µuand ˜µnare discretized based on the splitting of the energy as, ˜ µu= DuEc(uk+1, nk+1) − DuEe(uk, nk) = −²2∆uk+1+ Ψ0c,k+1− Ψ0e,k (3.28c) ˜ µn= DnEc(uk+1, nk+1) − DnEe(uk, nk) = nk+1 δ (3.28d)

Using that Ec and Ee are convex, and that P is a nonnegative function, it can be

straight-forwardly proven, see [50], that the above scheme is unconditionally energy-stable, i.e.,

E (uk+1, nk+1) ≤ E(uk, nk). Note that for quadraticΨc(u), the scheme is linear.

3.3.3.

S

ECOND

-

ORDER ACCURATE SCHEME

Our second-order accurate scheme for the Cahn–Hilliard equation can be extended to the tumor-growth model as follows:

uk+1− uk τ = ∆ ˜µu+ ˜Pk+1/2( ˜µn− ˜µu) (3.29a) ˜ µu= ˜Ψ0(uk, uk+1) − ²2∆ uk+1+ uk 2 − α1τ∆(uk+1− uk) + α2τ(uk+1− uk) (3.29b) nk+1− nk τ = ∆ ˜µn− ˜Pk+1/2( ˜µn− ˜µu) (3.29c) ˜ µn= nk+1+ nk 2δ (3.29d) where ˜ Pk+1/2:= Pµ 3 2uk− 1 2uk−1 ¶ ˜ Ψ0(u k, uk+1) = Ψ0c,k+1uk+1− uk 2 Ψ 00 c,k+1− Ψ0e,kuk+1− uk 2 Ψ 00 e,k

Similar as before, the above scheme is a modification of the Crank–Nicolson method. It involves the second-order accurate convex splitting of Ψ introduced in Sec. 3.2.3, and two types of stabilizations. Theα1-stabilization is an artificial diffusivitiy, while theα2

-stabilization can be thought of as artificial convexity. Furthermore, we apply the extrap-olation technique for the treatment of P (u), which allows the scheme to be linear (for quadraticΨc(u)). The initial step of the scheme done by setting u−1= u0.

We summarize the properties of the above scheme in the following theorem.

Theorem 3.3.A Let the free energy densityΨ be C2,1-continuous and have a convex split-tingΨ = Ψc− Ψewith C2,1-continuousΨcandΨeand finite second derivatives, i.e., |Ψ00c| ≤

Lc and |Ψ00e| ≤ Le for some constants Lc, Le ≥ 0. Let the proliferation function P be C0,1

-continuous and satisfy 0 ≤ P ≤ ¯P . If the stabilization is large enough, i.e.α1≥ (Lc+ Le)2/16

(34)

3.3.DIFFUSE-INTERFACE TUMOR-GROWTH MODEL 25

1. Unconditional energy-stability: E (uk+1, nk+1) ≤ E(uk, nk)

2. Total mass conservation: Z

(uk+1+ nk+1)d x =

Z

(u0+ n0)d x

3. Second-order accuracy. 

Remark 3.3.1 Variable mobilities can also be considered in the tumor-growth model. If

they satisfy the same mobility conditions as in Theorem3.2.A, and if they are treated by extrapolation, then Theorem3.3.Aalso holds, but with the constraint:

α1≥ ¯M (Lc+ Le)2/16 

Proof 1. According to (3.7), we have, by convexity

E (uk+1, nk+1) − E(uk, nk) = Z ΩΨ(uk+1) − Ψ(uk) + ²2 2(|∇uk+1| 2 − |∇uk|2) + (nk+1)2− (nk)2 2δ ≤ Z Ω ³ Ψ0 c,k+1− Ψ0e,k ´ (uk+1− uk) +² 2 2(|∇uk+1| 2 − |∇uk|2) + (nk+1− nk)(nk+1+ nk) 2δ

Using the time-stepping scheme (3.29b), (3.29d) and integration by parts, it follows that

E (uk+1, nk+1) − E(uk, nk) ≤ Z Ω n³ ˜ µu+ uk+1− uk 2 (Ψ 00 c,k+1+ Ψ00e,k) ´ (uk+1− uk) + ˜µn(nk+1− nk) o −α1τk∇(uk+1− uk)kL22− α2τkuk+1− ukk 2 L2 (3.30) Since 0 ≤ P ≤ ¯P , we have 0 ≤ ˜Pk+1/2≤ ¯P (3.31)

From (3.9), (3.29a), (3.29c), (3.31) and integration by parts, we have Z Ω n³ ˜ µ +uk+1− uk 2 (Ψ 00 c,k+1+ Ψ00e,k) ´ (uk+1− uk) + ˜µn(nk+1− nk) o ≤τ 2(supR Ψ 00 c,k+1+ supR Ψ00e,k) Z Ω(uk+1− uk) ¡ ∆ ˜µu+ ˜Pk+1/2( ˜µn− ˜µu) ¢ + Z Ω n τ ˜µu ¡ ∆ ˜µu+ ˜Pk+1/2( ˜µn− ˜µu)¢ + τ ˜µn ¡ ∆ ˜µn− ˜Pk+1/2( ˜µn− ˜µu) ¢o ≤ −τ 2(supR Ψ 00 c,k+1+ sup R Ψ 00 e,k) Z Ω∇ ˜µ · ∇(uk+1− uk) +τ 2(supR Ψ 00 c,k+1+ sup R Ψ 00 e,k) Z Ω ¯ ¯ ˜P1/2 k+1/2( ˜µn− ˜µu) ¯ ¯ ¯ ¯ ˜P1/2 k+1/2(uk+1− uk) ¯ ¯ −τk∇ ˜µuk2L2− τk∇ ˜µnk2L2− τ Z Ω ˜ Pk+1/2( ˜µn− ˜µu)2

(35)

Applying Young’s inequality to the above result, and substituting this into (3.30), we have E (uk+1, nk+1) − E(uk, nk) ≤τ 4(supR Ψ 00 c,k+1+ sup R Ψ 00 e,k) Z Ω µ β1|∇(uk+1− uk)|2+ 1 β1 |∇ ˜µu|2 ¶ − τk∇ ˜µuk2L2 +τ 4(supR Ψ 00 c,k+1+ sup R Ψ 00 e,k) Z Ω µ ˜ Pk+1/2 β2 ( ˜µn− ˜µu)2+ ¯2(uk+1− uk)2 ¶ − τk∇ ˜µnk2L2 − α1τk∇(uk+1− uk)k2L2− α2τkuk+1− ukk2L2− τ Z Ω ˜ Pk+1/2( ˜µn− ˜µu)2 ≤ −τ µ 1 − 1 4β1 (sup R Ψ 00 c,k+1+ sup R Ψ 00 e,k) ¶ k∇ ˜µuk2L2 − τ µ 1 − 1 4β2 (sup R Ψ 00 c,k+1+ sup R Ψ 00 e,k) ¶ Z Ω ˜ Pk+1/2( ˜µn− ˜µu)2 − τ µ α1−β1 4 (supR Ψ 00 c,k+1+ sup R Ψ 00 e,k) ¶ k∇(uk+1− uk)k2L2 − τ µ α2− β2P¯ 4 (supR Ψ 00 c,k+1+ sup R Ψ 00 e,k) ¶ kuk+1− ukk2L2− τk∇ ˜µnk2L2

Arguing as in the proof of Theorem3.2.A, we will have E (uk+1, nk+1) ≤ E(uk, nk)

indepen-dent ofτ > 0, if the stabilization constants satisfy

α1≥ (sup R Ψ 00 c,k+1+ sup R Ψ 00 e,k) 2/16, α 2≥ ¯P (sup R Ψ 00 c,k+1+ sup R Ψ 00 e,k) 2/16

2. Integrating (3.29a) and (3.29c) over the domainΩ and applying the divergence The-orem and homogeneous Neumann boundary conditions, we get

Z Ω ³u k+1− uk τ + nk+1− nk τ ´ d x = Z Ω ¡ ∆ ˜µu+ ∆ ˜µn¢ d x = Z ∂Ω(∇ ˜µu· n + ∇ ˜µn· n)dS = 0

Then, it follows by induction that Z

(uk+1+ nk+1)d x =

Z

(u0+ n0)d x, ∀k = 1,..., N

3. To compute the local truncation error, we assume that uk+1and nk+1are determined from exact values for uk, uk−1and nk. Then

uk+1= u(tk) + τ∆ ˜µu(u(tk+1), u(tk))

+τ ˜Pk+1/2(u(tk), u(tk−1))( ˜µn(n(tk+1), n(tk)) − ˜µu(u(tk+1), u(tk)))

nk+1= n(tk) + τ∆ ˜µu(u(tk+1), u(tk))

−τ ˜Pk+1/2(u(tk), u(tk−1))( ˜µn(n(tk+1), n(tk)) − ˜µu(u(tk+1), u(tk)))

(36)

3.4.NUMERICAL EXPERIMENTS 27

According to (3.17),(3.18) and (3.19), we know that ˜

µu(u(tk+1), u(tk)) = µu(u(tk+1/2)) + O(τ2)

˜

µn(n(tk+1), n(tk)) = µn(n(tk+1/2)) + O(τ2)

˜

Pk+1/2(u(tk), u(tk−1)) = P(u(tk+1/2)) + O(τ2)

(3.33)

Inserting (3.33) into (3.32), it follows that

uk+1= u(tk) + τ∆µu(u(tk+1/2))

+τ P (u(tk+1/2))(µn(n(tk+1/2)) − µu(u(tk+1/2))) + O(τ3)

nk+1= n(tk) + τ∆µn(n(tk+1/2))

−τ P (u(tk+1/2))(µn(n(tk+1/2)) − µu(u(tk+1/2))) + O(τ3)

(3.34)

Then, making use of a Taylor expansion to expand the exact value of u and n about tk+1/2 similar to (3.21), we obtain

u(tk+1) − u(tk) = τ ˙u(tk+1/2) + O(τ3)

n(tk+1) − n(tk) = τ ˙n(tk+1/2) + O(τ3)

Combining the above results, the local truncation error of u and n is

u(tk+1) − uk+1 = O(τ3)

n(tk+1) − nk+1 = O(τ3) ■

3.4.

N

UMERICAL EXPERIMENTS

In this section, we investigate numerically the robustness, stability and accuracy of our second-order time-stepping schemes for both the Cahn–Hilliard equation and the diffuse-interface tumor-growth model. We verify the time accuracy of our schemes with constant mobility M (u) = 1 and variable mobility M(u) = ˆM (1 − u2), and demonstrate the discrete energy dissipation and the discrete mass conservation properties of the scheme. We also make a comparison with the first-order accurate schemes. Notice that because of our choice ofΨ, see (3.1), all schemes are linear.

For the discretization in space, we consider linear finite element approximations on fixed uniform meshes based on mixed weak formulations; see for such discretizations, e.g. [73, Sec. II.C] for the Cahn–Hilliard equation and [50, Sec. 4.3] or [49, Sec. A.2] for the tumor-growth system. We have included a brief description of our fully-discrete scheme in the Appendix. For the time discretization, we use constant time steps.

Referenties

GERELATEERDE DOCUMENTEN

Van materialen op basis van stro of hout is bekend dat ze stikstof vastleggen (imobilisatie) waardoor er minder stikstof beschikbaar komt voor het gewas.. In Topsoil+ zijn enkele

We kunnen echter niet stellen dat goed uitgelopen vier keer zo groot is als bijna dood, of dat het verschil tussen dood en bijna dood twee keer zo groot is als het verschil

Namelijk, aantrekkende prijzen voor de meeste snijbloemen en kamerplanten, een stijging van de exportwaarde ten opzichte van vergelijkbare preioden vorig jaar en hogere kosten

Door te maaien in twee of drie blokken en tijdig met de koeien naar het etgroen te gaan (eerste perceel van blok met halve weidesnede) kunnen alle etgroenpercelen beweid worden..

oorzaken met een relatief hogere frequentie van voorkomen en negatieve impact zijn: fouten in het ontwerp van de geotechnische constructie, een onvoldoende robuust ontwerp en

The contingency approach illustrates which factors influence transfer success (i.e., hypothesis testing approach). Yet it does not provide rich description on the process

interpretatie van het onderschrift, waarbij de protagonisten de ogen van de vogels met bladeren bedekken, kan de hand van Loplop richting het oog van de vogel gelezen worden als

In this study, no difference was found in the association of parental expressed anxiety and children’s fear and avoidance between mothers and fathers, therefore it makes no