Numerical Analysis for a System Coupling Curve Evolution to Reaction-Diffusion on the Curve

# Numerical Analysis for a System Coupling Curve Evolution to Reaction-Diffusion on the Curve

John W. Barrett222Department of Mathematics, Imperial College London, London, SW7 2AZ, UK    Klaus Deckelnick 333Institut für Analysis und Numerik, Otto-von-Guericke-Universität, 39106 Magdeburg, Germany     Vanessa Styles 444Department of Mathematics, University of Sussex, Brighton, BN1 9RF, UK
###### Abstract

We consider a finite element approximation for a system consisting of the evolution of a closed planar curve by forced curve shortening flow coupled to a reaction-diffusion equation on the evolving curve. The scheme for the curve evolution is based on a parametric description allowing for tangential motion, whereas the discretisation for the PDE on the curve uses an idea from [6]. We prove optimal error bounds for the resulting fully discrete approximation and present numerical experiments. These confirm our estimates and also illustrate the advantage of the tangential motion of the mesh points in practice.

s

urface PDE, forced curve shortening flow, diffusion induced grain boundary motion, parametric finite elements, tangential motion, error analysis

{AMS}

65M60, 65M15, 35K55, 53C44, 74N20

## 1 Introduction

The aim of this paper is to analyse a fully discrete numerical scheme for approximating a solution of the following system: find a family of planar, closed curves and a function such that

 (1a) v =κ+f(w) on Γ(t),t∈(0,T], (1b) ∂∙tw =dwss+κvw+g(v,w) on Γ(t),t∈(0,T],

subject to the initial conditions

 (2) Γ(0)=Γ0,w(⋅,0)=w0 on Γ0.

Here, and are the normal velocity and the curvature of corresponding to the choice of a unit normal, while is the arclength parameter on . Furthermore, denotes the material derivative of , i.e. . Finally, , , , the closed curve and are all given. The system (1a,b) couples the evolution of the curves by forced curve shortening flow to a parabolic PDE on the moving curves. It occurs for example as a sharp–interface model for diffusion induced grain boundary motion: in this setting represents a grain boundary separating the crystals of a thin polycrystalline film of metal that is placed in a vapour containing another metal; atoms from the vapour diffuse into the film along the grain boundaries causing them to move. A thorough description of the physical set-up can be found in [11], while an existence and uniqueness result has been obtained in [12].

In what follows we shall describe the evolving curves with the help of a parametrisation , where is the periodic unit interval. Then

 (3) →τ=→xs=→xρ|→xρ|,→ν=→τ⊥

are a unit tangent and unit normal to respectively, where denotes counter-clockwise rotation by . The normal and tangential velocities of are then

 (4) v=→xt⋅→ν,ψ=→xt⋅→τ.

Note that (1a) only prescribes and with it the shape of , so that there is a certain freedom in choosing the tangential velocity. Since one may consider

 (5) →xt=→xss+f(˜w)→ν=1|→xρ|(→xρ|→xρ|)ρ+f(˜w)→ν,

where . Clearly, (1a) holds with the additional property that the velocity vector points in the normal direction. Coupling (5) to the PDE satisfied by , Pozzi and Stinner have derived and analysed in [13] a finite element scheme for (1a,b) (with ) based on continuous piecewise linears. They are able to prove the following error bounds in the spatially discrete case:

 sup[0,T]∫I(|→xρ−→xhρ|2+|˜w−˜wh|2)dρ+∫T0∫I(|→xt−→xht|2+|˜wρ−˜whρ|2)dρdt≤Ch2.

A major difficulty in the analysis arises from the fact that (5) is only weakly parabolic. Following [5], this problem is solved in [13] by deriving additional equations for the continuous and discrete length elements and by splitting the error into a tangent part and a length element part. Apart from these analytical difficulties, the motion in purely the normal direction also may lead to the accumulation of mesh points in numerical simulations. A natural way to handle the above mentioned difficulties is to introduce a tangential part in the velocity, which can be seen as a reparametrisation, an approach that has recently been explored in a systematic way by Elliott and Fritz in [7]. The underlying idea uses the DeTurck trick in coupling the motion of the curve to the harmonic map heat flow, which results in the following equation replacing (5) (cf. (3.1) in [7]):

 (6) α→xt+(1−α)(→xt⋅→ν)→ν=→xρρ|→xρ|2+f(˜w)→ν.

In the above, is a parameter so that corresponds to the diffusion coefficient in the harmonic map heat flow. Note that we obtain (1a) by taking the scalar product of (6) with . It turns out that gets closer to a parametrisation proportional to arclength as gets small. At the numerical level this means that mesh points along the curve become more and more equidistributed. Setting formally one recovers an approach introduced by Barrett, Garcke and Nürnberg, see [1], [2]. A nice feature of (6) is that the problem now is strictly parabolic allowing for a more straightforward error analysis; for the spatially discrete case with , see [3], [4] for , and [7] for .

It remains to derive the equation for , which in contrast to [13] will also involve the tangential velocity . It is easily seen that

 ˜wt(ρ,t)=ddt[w(→x(ρ,t),t)]=∂∙tw(→x(ρ,t),t)+(→xt(ρ,t)⋅→τ(→x(ρ,t))ws(→x(ρ,t),t).

Inserting this relation into (1b) and recalling (4), we obtain

 (7) ˜wt−(→xt⋅→τ)1|→xρ|˜wρ−d|→xρ|(˜wρ|→xρ|)ρ−κv˜w=g(v,˜w).

The initial conditions for (6), (7) now are: , , , where is a parametrization of . We shall derive in Section 2 a weak formulation for (6), (7) which forms the basis for a discretisation by continuous piecewise linear finite elements in space and a backward Euler scheme in time. As the main result of our paper we shall present an error analysis for the resulting fully discrete scheme. The precise result will be formulated at the end of Section 2, while the proof is carried out in detail in Section 3. In Section 4 we present a series of test calculations that confirm our theoretical results and demonstrate the above mentioned improvement of the mesh quality for small values of .

Finally, we end this section with a few comments about notation. We adopt the standard notation for Sobolev spaces, denoting the norm of (, and a bounded interval in ) by and the semi-norm by . For , will be denoted by with the associated norm and semi-norm written, as respectively, and . For ease of notation, in the common case when the subscript “” will be dropped on the above norms and semi-norms. The above are naturally extended to vector functions, and we will write for a vector function with two components. In addition, we adopt the standard notation (, , an interval in , a Banach space) for time dependent spaces with norm . Once again, we write if . Furthermore, denotes a generic constant independent of the mesh parameter and the time step , see below.

## 2 Weak formulation and finite element approximation

We shall assume that the data and the solution of (6), (7) (writing again instead of for ease of notation) are such that

 (1a) f∈C1(R,R),g∈C1(R×R,R), (1b) ∀ L>0, ∃ CL≥0: |g(v1,w)−g(v2,w)|≤CL|v1−v2|∀ v1,v2, |w|≤L, (1c) →x∈W1,∞(0,T;[H2(I)]2)∩H2(0,T;[H1(I)]2)∩W2,∞(0,T;[L2(I)]2), (1d) w∈C([0,T];H2(I))∩W1,∞(0,T;H1(I))∩H2(0,T;L2(I)), (1e) 00.

Since our error analysis will also include the discretisation in time our regularity assumptions are slightly stronger than those made in [13], see Assumption 2.2 there. Let us fix . For a test function , we obtain on multiplying (6) by that

 (2) ∫I|→xρ|2[α→xt+(1−α)(→xt⋅→ν)→ν]⋅→ξdρ+∫I→xρ⋅→ξρdρ=∫I|→xρ|2f(w)→ν⋅→ξdρ.

In order to derive a weak formulation for (7) we employ an idea from [6], [9], and calculate for with the help of integration by parts and (3)

 (3) = ∫I|→xρ|wtηdρ−∫I(→xt⋅→τ)wρηdρ−∫I(→xt⋅→τ)wηρdρ−∫Iκ(→xt⋅→ν)wη|→xρ|dρ =−d∫Iwρηρ|→xρ|dρ−∫Iψwηρdρ+∫I|→xρ|g(v,w)ηdρ,

where we have also used (4), (7) and the fact that .

We now use (2), (3) in order to discretise our system and begin by introducing the decomposition , where . We set , where and assume that

 (4) h≤Chj,j=1,…,J.

Let

 (5) Vh1:={χ∈C(I):χ∣σj is affine, j=1,…,J}⊂H1(I)

and be the standard Lagrange interpolation operator such that , . We require also the local interpolation operator , and recall for , , and that

 (6a) h1pj|ηh|0,∞,σj+hj|ηh|1,p,σj ≤C|ηh|0,p,σj ∀ ηh∈Vh1, (6b) |(I−Ihj)η|k,p,σj ≤Chℓ−kj|η|ℓ,p,σj ∀ η∈Wℓ,p(σj), (6c) |(I−Ihj)η|ℓ−1,∞,σj ≤Ch12j|η|ℓ,σj ∀ η∈Hℓ(σj).

As well as the standard inner product , we introduce the discrete inner product defined by

 (7) (η1,η2)h:=J∑j=1∫σjIhj(η1η2),

where are piecewise continuous functions on the partition of . We note for and for all that

 (8a) ∫σj|ηh|2dρ ≤∫σjIhj[|ηh|2]dρ≤3∫σj|ηh|2dρ, (8b) ∣∣∣∫σj(I−Ihj)(ηhχh)dρ∣∣∣ ≤Ch2j|ηh|1,σj|χh|1,σj≤Chj|ηh|1,σj|χh|0,σj.

The result (8b) follows immediately from (6a,b). The inner products and are naturally extended to vector functions. In addition to the above spatial discretisation, let be a partitioning of with time steps , , and . Before we define our scheme we assign to an element (the upper index referring to the time level ) a piecewise constant discrete unit tangent and normal by

 (9) →Tn=→Xnρ|→Xnρ|,→Vn=(→Tn)⊥ on% σj,j=1,…,J.

Our discretisation of (2) now reads: given and , find such that

 (10) =(|→Xn−1ρ|2f(Wn−1)→Vn−1,→ξh)h∀ →ξh∈[Vh1]2.

Here, and in what follows, we abbreviate . We next use the solution of (10) in order to discretise (3). To do so, we define approximations , of the normal and tangential velocities by

 (11) Vn=Dt→Xn⋅→VnandΨn=Dt→Xn⋅→Tn on σj,j=1,…,J.

Then find such that

 (12) =(|→Xnρ|g(Vn,Wn−1),ηh)h∀ ηh∈Vh1.

Let us formulate the main result of this paper, which will be proved in Section 3. {theorem} Let and . There exists such that for all and the discrete problem (10), (12) has a unique solution , and the following error bounds hold:

 supn=0,…,N[|→xn−→Xn|21+|wn−Wn|20] (13) +N∑n=1Δtn[∣∣→xnt−Dt→Xn∣∣20+|wn−Wn|21]≤Ch2,

where .

## 3 Error analysis

To begin, it follows from (1c,d) and (4) that for

 (1) ∥→xn∥2+∥→xnt∥1+∥vn∥1+∥ψn∥1+∥wn∥2≤C,

where . We abbreviate for

 (2) →En:=Ih→xn−→Xn∈[Vh1]2andZn:=Ihwn−Wn∈Vh1.

Our proof of Theorem 2 will be based on induction. We assume for an that

 (3) |→En−1|21+β2(|→Xn−1ρ|Zn−1,Zn−1)h ≤h2eγ∫tn−10ζ(t)dt,for  h∈(0,h⋆],

where the function is defined by

 (4) ζ(t):=1+∥→xtt(t)∥21+|wtt(t)|20

and is chosen so small that

 (5) (h⋆)12eγK≤βwithK:=∫T0ζ(t)dt.

Here, and are independent of and , and will be chosen a posteriori. Clearly, (3) holds for in view of our choice of initial data for and .
It follows from (4), (6c), (3), (1) and (5) that

 |→xn−1−→Xn−1|1,∞ ≤|→En−1|1,∞+|(I−Ih)→xn−1|1,∞≤Ch−12|→En−1|1+Ch12|→xn−1|2 (6) ≤Ch12(eγK2+1)≤C(h⋆)12eγK≤Cβ≤min{m2,M},

provided that is chosen small enough. Combining this inequality with (1e) we infer that

 (7) 0

If we use (7) and argue similarly as in (6) we further obtain for any

 (8a) |→τn−1−→Tn−1|0+|→νn−1−→Vn−1|0 ≤C|→xn−1−→Xn−1|1≤C[|→En−1|1+h], |→τn−1−→Tn−1|0,∞+|→νn−1−→Vn−1|0,∞ ≤C|→xn−1−→Xn−1|1,∞ (8b) ≤Ch−12[|→En−1|1+h].

In the same way as (6), we obtain from (8a), (7), (3) and (5) that

 (9) |Zn−1|20≤C(|→Xn−1ρ|Zn−1,Zn−1)h≤Cβ2h2eγK≤Cβ2hh⋆e2γK≤Ch,

which combined with (6a,c), (4) and (1) yields for

 (10) |Wn−1|0,∞≤|Zn−1|0,∞+|Ihwn−1|0,∞≤Ch−12|Zn−1|0+C≤C.

### 3.1 The curve equation

We assume throughout that . Since (10) forms a linear problem it is easily seen that (7) implies the existence of a unique solution to (10). We deduce from (10) and (2) with that

 LHS +Δtn(→Enρ,(Dt→En)ρ) +Δtn(→xnρ,(Dt→En)ρ)−Δtn(|→Xn−1ρ|2f(Wn−1)→Vn−1,Dt→En)h =Δtn[(|→Xn−1ρ|2[αDt→xn+(1−α)(Dt→xn⋅→Vn−1)→Vn−1],Dt→En)h −(|→xnρ|2[α→xt(tn)+(1−α)(→xt(tn)⋅→νn)→νn],Dt→En)] +Δtn[(|→xnρ|2f(wn)→νn,Dt→En)−(|→Xn−1ρ|2f(Wn−1)→Vn−1,Dt→En)h] (11) =:A1+A2.

Using (7) and (8a) we find, with the help of an elementary calculation, that

 (12) LHS≥Δtnαm24∣∣Dt→En∣∣20+12[|→En|21+|→En−→En−1|21−|→En−1|21].

Let us analyse the term defined in (11) and note that

 A1 =Δtn[(|→Xn−1ρ|2[αDt(Ih→xn)+(1−α)(Dt(Ih→xn)⋅→Vn−1)→Vn−1],Dt→En)h −(|→Xn−1ρ|2[αDt(Ih→xn)+(1−α)(Dt(Ih→xn)⋅→Vn−1)→Vn−1],Dt→En)] +Δtn[α(|→Xn−1ρ|2[Dt(Ih→xn)−→xnt],Dt→En) +(1−α)(|→Xn−1ρ|2([Dt(Ih→xn)−→xnt]⋅→Vn−1)→Vn−1,Dt→En)] +Δtn[(1−α)(|→Xn−1ρ|2[(→xnt⋅→Vn−1)→Vn−1−(→xnt⋅→νn)→νn],Dt→En) (13) +([|→Xn−1ρ|2−|→xnρ|2][α→xnt+(1−α)(→xnt⋅→νn)→νn],Dt→En)]=:3∑i=1A1,i.

We now bound the terms defined in (13) on recalling (1c,e), (7), (3) and (9). It follows from (8b) and (6b) that

 (14) |A1,1| ≤ChΔtn|Dt(Ih→xn)|1|Dt→En|0≤ChΔtn|Dt→En|0,

since . Similarly,

 |A1,2| ≤CΔtn[|(I−Ih)Dt→xn|0+|Dt→xn−→xnt|0]|Dt→En|0 (15) ≤CΔtn⎡⎢⎣h+(Δtn∫tntn−1|→xtt|20dt)12⎤⎥⎦|Dt→En|0.

Next, we have from Sobolev embedding, (6b) and (1) that

 |A1,3| ≤CΔtn[|→νn−→Vn−1|0+∣∣|→xnρ|2−|→Xn−1ρ|2∣∣0]|Dt→En|0 ≤CΔtn[Δtn|Dt→xn|1+|→xn−1−→Xn−1|1]|Dt→En|0 (16) ≤CΔtn[h+|→En−1|1]|Dt→En|0,

since . We now turn our attention to the term in (11) and note that

 A2 =Δtn(|→xnρ|2f(wn)(→νn−→Vn−1)+(|→xnρ|2−|→Xn−1ρ|2)f(wn)→Vn−1,Dt→En) <