Evolving surface finite element method for the Cahn-Hilliard equationThe work of C. M. Elliott was supported by the UK Engineering and Physical Sciences Research Council EPSRC Grant EP/G010404 and the work of T. Ranner was supported by a EPSRC Ph.D. studentship (Grant EP/P504333/1 and EP/P50516X/1) and the Warwick Impact Fund.

Evolving surface finite element method for the Cahn-Hilliard equation††thanks: The work of C. M. Elliott was supported by the UK Engineering and Physical Sciences Research Council EPSRC Grant EP/G010404 and the work of T. Ranner was supported by a EPSRC Ph.D. studentship (Grant EP/P504333/1 and EP/P50516X/1) and the Warwick Impact Fund.

Charles M. Elliott C.M.Elliott Mathematics Institute, Zeeman Building, University of Warwick, Coventry. CV4 7AL. UK
Tel.: +44 (0)24 7615 0773
22email: C.M.Elliott@warwick.ac.ukT. Ranner Mathematics Institute, Zeeman Building, University of Warwick, Coventry. CV4 7AL. UK
Present address: School of Computing, University of Leeds. LS2 9JT. UK
E-mail: T.Ranner@leeds.ac.uk
Thomas Ranner C.M.Elliott Mathematics Institute, Zeeman Building, University of Warwick, Coventry. CV4 7AL. UK
Tel.: +44 (0)24 7615 0773
22email: C.M.Elliott@warwick.ac.ukT. Ranner Mathematics Institute, Zeeman Building, University of Warwick, Coventry. CV4 7AL. UK
Present address: School of Computing, University of Leeds. LS2 9JT. UK
E-mail: T.Ranner@leeds.ac.uk
The final publication is available at springerlink.com. DOI: 10.1007/s00211-014-0644-y
Abstract

We use the evolving surface finite element method to solve a Cahn-Hilliard equation on an evolving surface with prescribed velocity. We start by deriving the equation using a conservation law and appropriate transport formulae and provide the necessary functional analytic setting. The finite element method relies on evolving an initial triangulation by moving the nodes according to the prescribed velocity. We go on to show a rigorous well-posedness result for the continuous equations by showing convergence, along a subsequence, of the finite element scheme. We conclude the paper by deriving error estimates and present various numerical examples.

Keywords:
Evolving surface finite element method Cahn-Hilliard equation triangulated surfaces error analysis
journal: Numerische Mathematik

1 Introduction

In this paper, we will study a Cahn-Hilliard equation posed on an evolving surface with prescribed velocity. The key methodology is to discretise the equations using the evolving surface finite element method DziEll07 () originally proposed for a surface heat equation. The idea is to take a triangulation of the initial surface and evolve the nodes along the velocity field. This leads to a family of discrete surfaces on which we can pose a variational form of the Cahn-Hilliard equation.

There are two key results in this paper: first, we show well posedness of the continuous scheme and, second, we show convergence of a finite element scheme. The well posedness result is proven by rigorously showing convergence, along a subsequence, of the discrete scheme. In contrast to the planar setting, there are extra difficulties in this work since the classical Bochner space set-up is unavailable to us. The finite element method is analysed under the assumption of higher regularity of the solution and shown to converge to the true solution quadratically with respect to the mesh size in an norm. The paper concludes with some numerical examples to show various properties of the methodology.

1.1 The Cahn-Hilliard equation

We assume we are given an evolving surface , for , which evolves according to a given underlying velocity field which can be decomposed into normal and tangential components so that . We seek a solution of

 ∂∙u+u∇Γ⋅v=ΔΓ(−εΔΓu+1εψ′(u)) on ⋃t∈(0,T)Γ(t)×{t} (1.1)

subject to the initial condition

 u(⋅,0)=u0 on Γ(0)=Γ0. (1.2)

Here denotes the material derivative of and the Laplace-Beltrami operator of . The function is a double well potential, which we will take to be given by

 ψ(z)=14(z2−1)2. (1.3)

The behaviour of the Cahn-Hilliard equation in the planar case is well studied Ell89 (). Extra effects such as spatial or concentration dependent mobilities or more physically realistic potentials could also be solved with similar methods to those suggested in this paper. Such considerations are left for future work.

This Cahn-Hilliard equation is a simplification of the model for surface dissolution set out in EilEll08 (); ErlAziKar01 () arising from a conservation law. The model MerPtaKuh12 () takes a different approach and considers a gradient flow for an energy consisting of the sum of the Ginzburg-Landau functional and a Helfrich energy on a stationary surface. One could alternatively couple the evolution of the surface to the surface field and recover a gradient flow of the Ginzburg-Landau functional EllSti10b (); EllSti10a ().

The results in this work can be seen as a generalisation of the work of DuJuTia09 () to evolving surfaces. That work considers a fully discrete approximation of a Cahn-Hilliard equation posed on a two-dimensional stationary surface with boundary (with a zero Dirichlet boundary condition) under the assumption and for . Their method uses a triangulated surface for the spatial discretisation and a Crank-Nicolson scheme in time. They show an error estimate of the form

 maxm∥∥umh−u−ℓ(tm)∥∥L2(Γh)≤c(h2+τ2),

where is a partition of time with fixed time step and is the inverse lift (3.21) of the continuous solution .

1.2 Outline of paper

The paper is laid out as follows. In Section two, we will derive a Cahn-Hilliard equation on an evolving surface using a local conservation law. We introduce the notation for partial differential equations on evolving surfaces taken from DecDziEll05 (); DziEll13a () and state any assumptions on the smoothness of the surfaces and its evolution we require. The third section introduces a finite element discretisation of the continuous equations. We describe the process of triangulating an evolving surface and how we formulate the space discrete-time continuous problem as a system of ordinary differential equations. This section is completed by showing some domain perturbation results relating geometric quantities on the discrete and smooth surfaces. Well posedness of the continuous equations is addressed in the fourth section. An existence result is achieved by showing convergence, along a subsequence, of the discrete solutions as the mesh size tends to zero. In Section five, we analyse the errors introduced by our finite element scheme and go on to show an optimal order error estimate. Some numerical experiments are shown in the sixth section backing up the analytical results.

We will use a Gronwall inequality as a standard tool in the analysis which leads to exponential dependence on in most bounds. We are not interested in taking in this work so will simply write for a generic constant which depends on .

2 Derivation of continuous equations

In this section, we will derive a Cahn-Hilliard equation on an evolving surface as a conservative advection-diffusion equation. We will also introduce functional analytic setting and definition of solution that will be used.

2.1 Assumptions on the evolving surface

Given a final time , for each time , we write for a compact, smooth, connected -dimensional hypersurface in for or and . We assume that is the boundary of an open, bounded domain . It follows that admits a description as the zero level set of a signed distance function so that in and in . We denote by for the space-time domain given by

 GT=⋃t∈[0,T]Γ(t)×{t}. (2.1)

For our analysis, it is sufficient to consider locally to . We restrict our considerations to , an open neighbourhood of . We choose so that for and assume that

 d,dt,dxi,dxixj∈C2(NT) for i,j=1,…,n+1;

here . The orientation of is fixed by choosing as the outward pointing normal, so that . For , we denote the projection operator onto the tangent space , given by and by the (extended) Weingarten map (or shape operator),

 Hij(x,t)=(νi(x,t))xj=dxixj(x,t).

We will use the fact that . Finally, we denote by the mean curvature of

 H(x,t)=traceH(x,t)=n+1∑i=1Hii(x,t).

For a function , we define its tangential gradient by

 ∇Γη=∇˜η−∇˜η⋅νν=P∇˜η,

where is a smooth extension of away from . It can be shown that this definition is independent of the choice of extension. We denote the components of by

 ∇Γη=(D––1η,…,D––n+1η).

The Laplace-Beltrami operator is given by

 ΔΓη=∇Γ⋅∇Γη=n+1∑j=1D––jD––jη.

We will denote by the surface measure on which admits the following formula for partial integration for a portion (DziEll13a, , Theorem 2.10):

 ∫M(t)∇Γηdσ=∫M(t)ηHνdσ+∫∂M(t)ημdσ, (2.2)

where is the co-normal to which is normal to but tangent to . If and has no boundary, the boundary term vanishes. Furthermore, we have a Green’s formula on (DziEll13a, , Theorem 2.14):

 ∫Γ(t)∇Γη⋅∇Γφdσ=−∫Γ(t)φΔΓηdσ. (2.3)

These formulae allow the definition of weak derivatives and Sobolev spaces. We define the space by

 W1,q(Γ(t)):={η∈Lq(Γ(t)):D––jη∈Lq(Γ(t)) for j=1,…,n+1},

with norm

 ∥η∥W1,q(Γ(t))=(∥η∥qLq(Γ(t))+∥∇Γη∥qLq(Γ(t)))1q.

This can be easily extended to higher order spaces. See DziEll13a () for details. We will use the notation for .

We will make use of the following Sobolev embeddings:

Lemma 2.1 ((Heb00, , Theorems 2.5 and 2.6))

For as above, we have

 W1,q(Γ(t))⊂{Lnq/(n−q)(Γ(t)) for qn. (2.4)

Furthermore there exists a constant , independent of , such that for any ,

 ∥η∥Lnq/(n−q)(Γ(t)) ≤c∥η∥W1,q(Γ(t)) for qn. (2.5b)

In particular, this allows us to embed in for all dimensions so that .

Further, we assume that for each there exists a unique , such that

 x=p(x,t)+d(x,t)ν(p(x,t),t). (2.6)

See (GilTru01, , Chapter 14) for a proof. We extend and to functions on by setting

 ν(x,t)=ν(p(x,t),t)=∇d(x,t),

and similarly and for .

Although it is sufficient to describe the evolution of the surface through a normal velocity, we wish to consider material surfaces for which a material particle, at on , has a material velocity not necessarily only in the normal direction. The normal velocity of the surface can be calculated to be . We say is a tangential velocity field if in . Given a tangential velocity field , we call

 v:=vτ+vν

a material velocity field. We assume that we are given a global velocity field so that points evolve with the velocity . We will assume that .

2.2 Material derivative and transport formulae

Given a family of surfaces evolving in time with normal velocity field , we define the normal time derivative of a function by

 ∂∘η:=∂˜η∂t+vν⋅∇˜η. (2.7)

Here, denotes a smooth extension of to . This derivative describes how a quantity evolves in time with respect to the evolution of . It can be shown that this definition is an intrinsic surface derivative, independent of the choice of extension.

Given a tangential vector field , we define the material derivative of a scalar function , by

 ∂∙η:=∂∘η+vτ⋅∇Γη=∂˜η∂t+v⋅∇˜η.

The following formula shows the significance of the material derivative. The result is a generalisation of the classical Reynolds’ Transport Formula to curved domains.

Lemma 2.2 (Transport formula (DziEll13, , Lemma 2.1))

Let be an evolving surface with normal velocity . Let be a tangential velocity field on . Let the boundary evolve with velocity . Assume that are functions such that all the following quantities exist. Then, we obtain the identity

 ddt∫M(t)ηdσ=∫M(t)∂∙η+η∇Γ⋅vdσ. (2.8)

Furthermore, we have

 ddt∫M(t)ηφdσ=∫M(t)∂∙ηφ+η∂∙φ+ηφ∇Γ⋅vdσ. (2.9)

Let be a matrix which is positive definite on the tangent space to . Denote by the rate of deformation tensor given by

 D(v)ij=12n+1∑k=1(AikD––kvj+AjkD––kvi) for i,j=1,…,n+1, (2.10)

and by the tensor

 B(v):=∂∙A+∇Γ⋅vA−2D(v). (2.11)

Then we have the formula

 ddt∫M(t)A∇Γη⋅∇Γφdσ =∫M(t)A∇Γ∂∙η⋅∇Γφ+A∇Γη⋅∇Γ∂∙φdσ (2.12) +∫M(t)B(v)∇Γη⋅∇Γφdσ.

We conclude this subsection with a result allowing us to extend functions defined on one surface to the whole space-time domain.

Lemma 2.3

Fix and let , respectively . Then there exists an extension such that and , resp. , for all times and .

Proof

The ordinary differential equation:

 ddsX(s)=v(X(s),s) for s∈[0,T],X(t)=x,

determines a flow on for such that

 ϕs(x)∈Γ(s) for all s∈[0,T] and ϕt(x)=x.

Our assumptions on imply that and are both mappings (Har64, , Theorem 3.1).

We define the extension by

 ˜η(x,s):=η((ϕs)−1(x)) for (x,s)∈GT.

It is clear that since , we have (resp. ) for all times .

Finally, we can calculate for ,

 ∂∙˜η(x,s)=dds˜η(ϕs(y),s)=ddsη(y)=0 for (x,s)∈GT,

which shows the result. ∎

2.3 Derivation of Cahn-Hilliard equations

We will consider a conservation law on an evolving surface with a diffusive flux driven by a chemical potential. This is the approach taken by ErlAziKar01 (). In general, the Ginzburg-Landau functional on will not decrease along the trajectory of solutions.

Let represent a density of a scalar quantity on . Following DziEll13 (), we arrive at the pointwise conservation law

 ∂∘u+u∇Γ⋅vν+∇Γ⋅q=0. (2.13)

Here represents the tangential flux of on .

We will assume that the flux is the sum of a diffusive flux and an advective flux :

 qd=−∇Γw and qa=uvτ.

The diffusive flux is driven by the gradient of chemical potential gives us the split system EllFreMil89 ()

 ∂∙u+u∇Γ⋅v−ΔΓw =0 (2.14a) −εΔΓu+1εψ′(u)−w =0. (2.14b)

This leads to the fourth order Cahn-Hilliard equation on :

 ∂∙u+u∇Γ⋅v=ΔΓ(−εΔΓu+1εψ′(u)). (2.15)

We close the system with the initial condition

 u(⋅,0)=u0 on Γ0. (2.16)

There are no boundary conditions since the boundary of is empty.

Remark 2.1

One can derive the Cahn-Hilliard equations posed in a Cartesian domain as an gradient flow of the Ginzburg-Landau functional. To obtain a gradient flow on an evolving surface, there would need to be a model for and which would lead to a coupled system for and . In terms of modelling, we feel these extra terms are geometric terms determining an evolution equation for the surface, which we assume is given. Therefore, we do not consider such terms in this work.

2.4 Solution spaces

In standard parabolic theory one looks for solutions in Bochner spaces. Considering our Cahn-Hilliard equation on a Cartesian domain Ell89 (), one would expect solutions to live in the spaces

 u∈L∞(0,T;H1(Ω)),u′∈L2(0,T;H−1(Ω)),w∈L2(0,T;H1(Ω)).

These spaces are constructed by considering as a function from into the Hilbert space . We would like to extend this definition so that is in the now time-dependent Hilbert space . We consider Sobolev spaces over the space-time domain . We will write for the space-time gradient and for the space-time measure on . This approach is similar to the Eulerian formulation of OlsReuXu13-pp (). We contrast our approach with that of Vie11-pp (), who proposed using an equivalent formulation using a reference domain.

We start by presenting the space-time domains and defined by

 L2(GT) :={η∈L1loc(GT):∫GTη2dσT<+∞} H1(GT) :={η∈L2(GT):∇GTη∈L2(GT)}.

with norms

 ∥η∥L2(GT) :=(∫GTη2dσT)12 ∥η∥H1(GT) :=(∥η∥2L2(GT)+∥∥∇GTη∥∥2L2(GT))12.
Proposition 2.1 ((Heb00, , Theorem 2.9))

The space is compactly embedded into .

Using the identities,

and

our assumptions on imply that the space-time norms can be replaced with the equivalent norms

 ∥η∥′L2(GT) :=(∫T0∫Γ(t)η2dσdt)12 ∥η∥′H1(GT) :=(∫T0∫Γ(t)η2+|∇Γη|2+(∂∘η)2dσdt)12.

We will use the equivalent primed norms (dropping the prime) on and in the following.

We define the space by

 L2L2:={η∈L1loc(GT):∫T0∫Γ(t)ηdσdt<+∞},

with the inner product

 (η,ξ)L2L2:=∫T0∫Γ(t)ηξdσdt.

It is clear that is equivalent to and hence is a Hilbert space.

Next, we define the space as

 L2H1:={η:GT→R:η∈L2L2 and ∇Γη∈(L2L2)n+1},

with the inner product

 (η,ξ)L2H1:=∫T0∫Γ(t)∇Γη⋅∇Γξ+ηξdσdt,

where should be interpreted in the weak sense. Notice that elements of this space are weakly differentiable at almost every time.

Lemma 2.4

The space is a Hilbert space.

Proof

It is clear that is an inner product space and we are left to show completeness. Let be a Cauchy sequence in . This implies that and are Cauchy sequences in and . This means that there exists such that

 ∥ηk−η∥L2(GT)+∥∇Γηk−ξ∥L2(GT)→0 as k→∞.

Fix and let and . Using Lemma 2.3, we can construct such that and for each time . Then, for , we obtain

 ∫T0∫Γ(t)ηD––j(α˜φ)+ξj(α˜φ)dσdt =∫T0∫Γ(t)(η−ηk)D––j(α˜φ)+(ηkD––j(α˜φ)+ξj(α˜φ))dσdt =∫T0∫Γ(t)(η−ηk)D––j(α˜φ)+(−D––jηk+ξj)(α˜φ)dσdt,

where we have used the fact that is weakly differentiable at almost every time. Taking the limit , we infer

 ∫T0α(∫Γ(t)ηD––j˜φ+ξj˜φdσ)dt=0.

Since this holds for all , by the Fundamental Lemma of the Calculus of Variations, at , we have

 ∫Γ(t∗)ηD––jφ+ξjφdσ=0 for all φ∈C1(Γ(t∗)).

Since the choice of was arbitrary, we infer that is the weak gradient of for almost every time and the proof is complete. ∎

The equivalence of norms implies that with if, and only if, .

For , we will define the space by

 LqH1:={η∈Lq(GT):∥η∥LqH1<+∞},

with norm

 ∥η∥LqH1:=⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩(∫T0∥η∥qH1(Γ(t))dt)1q for q<∞,esssupt∈(0,T)∥η∥H1(Γ(t)) for q=∞.

It is clear that and that

 ∥η∥L2H1≤√T∥η∥L∞H1 for all η∈L∞H1.

Finally, we define and by

 L∞H2 :={η∈L2(GT):esssupt∈(0,T)∥η∥H2(Γ(t))<+∞} L2H2 :={η∈L2(GT):∫T0∥η∥2H2(Γ(t))dt<+∞}.
Remark 2.2

As a restriction on our analysis we will only consider as a function in since we do not wish to consider a weak material derivative. Such considerations are left to future work.

We conclude this section with a result which will take an integral in time equality into an almost everywhere in time equality. The proof is the generalisation of a similar result given in (Rob01, , Lemma 7.4) for planar domains.

Lemma 2.5

Let with

 ∫T0∫Γ(t)∇Γη⋅∇Γξ+ηξdσdt=0 for all ξ∈L2H1. (2.17)

Then for almost all times ,

 ∫Γ(t)∇Γη⋅∇Γφ+ηφdσ=0 for all φ∈L2H1. (2.18)
Proof

Fix and , then choosing and

 0=∫T0∫Γ(t)∇Γη⋅∇Γξ+ηξdσdt=∫T0α(∫Γ(t)∇Γη⋅∇Γφ+ηφdσ)dt.

Since the choice of was arbitrary, the Fundamental Lemma of the Calculus of Variations implies the result. ∎

2.5 Weak and variational form

We start by multiplying (2.14a, 2.14b) by a test function and apply integration by parts to the Laplacian terms to give the weak form. This will be the definition of solution used throughout this paper. Existence and uniqueness of solutions will be shown Section 4.

Definition 2.1 (Weak solution)

We say that the pair , with and , are a weak solution of the Cahn-Hilliard equation (2.15) if, for almost every time ,

 ∫Γ(t)∂∙uφ+uφ∇Γ⋅v+∇Γw⋅∇Γφdσ =0 (2.19a) ∫Γ(t)ε∇Γu⋅∇Γφ+1εψ′(u)φ−wφdσ =0, (2.19b)

for all ,

and pointwise almost everywhere in .

Restricting our thoughts to , applying the transport formula to the first two terms in (2.19a) gives the variational formulation:

 ddt(∫Γ(t)uφdσ)+∫Γ(t)∇Γw⋅∇Γφdσ =∫Γ(t)u∂∙φdσ (2.20a) ∫Γ(t)ε∇Γu⋅∇Γφ+1εψ′(u)φdσ =∫Γ(t)wφdσ. (2.20b)

We remark that this formulation has no explicit mention of the velocity field and will be the basis of our finite element calculations.

It will be useful to write these equations using abstract bilinear forms. We define the following three to describe the above equations for :

 m(η,φ)=∫Γ(t)ηφdσa(η,φ)=∫Γ(t)∇Γη⋅∇Γφdσg(v;η,φ)=∫Γ(t)ηφ∇Γ⋅vdσ.

This lets us write (2.19) as

 m(∂∙u,φ)+g(v;u,φ)+a(w,φ) =0 (2.21) εa(u,φ)+1εm(ψ′(u),φ)−m(w,φ) =0,

and (2.20) as

 ddtm(u,φ)+a(w,φ) =m(u,∂∙φ) (2.22) εa(u,φ)+1εm(ψ′(u),φ) =m(w,φ).

We may also write the results of Lemma 2.2 in this form:

 ddtm(η,φ) =m(∂∙η,φ)+m(η,∂∙φ)+g(v;η,φ) ddta(η,φ) =a(∂∙η,φ)+a(η,∂∙φ)+b(v;η,φ),

 b(v;η,φ)=∫Γ(t)B(v)∇Γη⋅∇Γφdσ,

using in the definition of .

3 Finite element approximation

In this section, we propose a finite element method for approximating solutions of the Cahn-Hilliard equation (2.15) based on the evolving surface finite element method DziEll07 ().

3.1 Evolving triangulation and discrete material derivative

Let be a polyhedral approximation of the initial surface with the restriction that the nodes of lie on . We evolve the nodes by the smooth surface velocity:

 ˙Xj(t)=v(Xj(t),t),Xj(0)=X0j, for j=1,…,N.

Linearly interpolating between these nodes defines a family of discrete surfaces . At each time, we assume that we have a triangulation of , with the maximum diameter of elements in uniformly in time:

 h:=supt∈(0,T)maxE(t)∈Th(t)diamE(t). (3.1)

We assume this triangulation is quasi-uniform BreSco02 () uniformly in time.

Remark 3.1

In practical situations, assuming a uniformly regular mesh may not be feasible. Large surface deformations can lead to poor quality triangulations with deformed elements. In such cases, re-meshing may be required ClaDieDzi04 (); EilEll08 (). Alternatively, one may use an arbitrary Lagrangian-Eulerian formulation by allowing extra tangential mesh motions EllSti13 (); EllSty12 ().

We define element-wise as the unit outward pointing normal to and denote by the tangential gradient on defined element-wise by

 ∇Γhηh:=∇˜ηh−(∇˜ηh⋅νh)νh=(Id−νh⊗νh)∇˜ηh=:Ph∇˜ηh.

This is a vector-valued quantity and we will denote its components by

 ∇Γhηh=(D––h,1ηh,…,D––h,n+1ηh)

We define the finite element space of piecewise linear functions on by

 Sh(t):={ϕh∈C(Γh(t)):ϕh|E(t) is affine % linear, for each E(t)∈Th(t)}. (3.2)

We will write for the nodal basis of given by .

The definition of a basis of allows us to characterise the velocity of the surface . An arbitrary point on evolves according to the discrete velocity given by

 ˙X(t)=Vh(X(t),t):=N∑j=1˙Xj(t)ϕNj(X(t),t)=N∑j=1v(Xj(t),t)ϕNj(X(t),t). (3.3)

We will write as the discrete equivalent to :

 Gh,T:=⋃t∈(0,T)Γh(t)×{t}. (3.4)

The discrete velocity induces a discrete material derivative. For a scalar quantity on , we define the discrete material derivative by

 ∂∙hηh:=∂t˜ηh+∇˜ηh