Error analysis for an ALE evolving surface finite element method

Error analysis for an ALE evolving surface finite element method

Charles M. Elliott Mathematics Institute, Zeeman Building, University of Warwick, Coventry, UK, CV4 7AL.  and  Chandrasekhar Venkataraman Department of Mathematics, University of Sussex, Falmer, UK, BN1 9QH.
12th July 2019
Abstract.

We consider an arbitrary-Lagrangian-Eulerian evolving surface finite element method for the numerical approximation of advection and diffusion of a conserved scalar quantity on a moving surface. We describe the method, prove optimal order error bounds and present numerical simulations that agree with the theoretical results.

C.M. Elliott]C.M.Elliott@warwick.ac.uk

C. Venkataraman]c.venkataraman@sussex.ac.uk

1. Introduction

For each let be a smooth connected hypersurface in , oriented by the normal vector field , with . We assume that there exists a diffeomorphism satisfying . We set with (the identity). Furthermore we assume that . The given velocity field may contain both normal and tangential components, i.e., and .

We focus on the following linear parabolic partial differential equation on ;

 (1.1) ∂∙vu+u∇Γ(t)⋅v−\UpdeltaΓ(t)u=0 on Γ(t),

where, denotes the surface gradient, the Laplace Beltrami operator and

 ∂∙vu=∂tu+v⋅∇u=∂tu+vν⋅∇u+vT⋅∇Γ(t)u

is the material derivative with respect to the velocity field . For simplicity we will assume that the boundary of is empty and hence no boundary conditions are needed. The method is easily adapted to surfaces with boundary. In the case that the surface has a boundary, under suitable assumptions (c.f., Remark 3.3), our analysis is valid for homogeneous Neumann boundary conditions, i.e.,

 (1.2) ∇Γ(t)u⋅μ=0 on ∂Γ(t),

where is the conormal to the boundary of the surface. The upshot is that the total mass is conserved i.e., . Note that the case being constant in space and time corresponds to the -dimensional hypersurface being flat. In this case (1.1) is a standard bulk PDE.

We expect that our results apply also to the case of Dirichlet boundary conditions however in this setting one must also estimate the error due to boundary approximation which we neglect in this work.

The following variational formulation of (1.1) was derived in [1] and makes use of the transport formula (A.1) and the integration by parts formula on the evolving surface [1, (2.2)] ;

 (1.3) ddt∫Γ(t)uφ+∫Γ(t)∇Γ(t)u⋅∇Γ(t)φ=∫Γ(t)u∂∙vφ,

where is a sufficiently smooth test function defined on the space-time surface

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

In [1], a finite element approximation was proposed for (1.3) using piecewise linear finite elements on a triangulated surface interpolating (at the nodes) , the vertices of the triangulated surface were moved with the material velocity (of points on ) . In this work we adopt a similar setup, in that we propose a finite element approximation using piecewise linear finite elements on a triangulated surface interpolating (at the nodes) , however we move the vertices of the triangulated surface with the velocity , where is an arbitrary tangential velocity field (). Furthermore we assume that satisfies the same smoothness assumptions as the material velocity , i.e., there exits a diffeomorphism satisfying with and with (the identity) and . We remark, that if has a boundary this assumption implies that the arbitrary velocity has zero conormal component on the boundary, i.e., for ,

 (1.4) (v−va)⋅μ=0 on ∂Γ(t).

For a sufficiently smooth function , we have that

 ∂∙vaf

Thus we may write the following equivalent variational formulation to (1.1), which will form the basis for our finite element approximation

 (1.5) ddt∫Γ(t)uφ+∫Γ(t)∇Γ(t)u⋅∇Γ(t)φ=∫Γ(t)(u∂∙vaφ−uaT⋅∇Γ(t)φ),

where is a sufficiently smooth test function defined on . We note that (1.5) may be thought of as a weak formulation of an advection diffusion-equation on a surface with material velocity , in which the advection is governed by some external process other than material transport. Hence the results we present are also an analysis of a numerical scheme for an advection diffusion equation on an evolving surface with another source of advective transport other than that due to the material velocity.

The original (Lagrangian) evolving surface finite element method (ESFEM) was proposed and analysed in [1], where optimal error bounds were shown for the error in the energy norm in the semidiscrete (space discrete) case. Optimal error bounds for the semidiscrete case were shown in [2] and an optimal error bound for the full discretisation was shown in [3]. High order Runge-Kutta time discretisations and BDF timestepping schemes for the ESFEM were analysed in [4] and [5] respectively. There has also been recent work on the analysis of ESFEM approximations of the Cahn-Hilliard equation on an evolving surface [6], scalar conservation laws on evolving surfaces [7] and the wave equation on an evolving surface [8]. For an overview of finite element methods for PDEs on fixed and evolving surfaces see [9]. Although the analytical results have thus far focussed on the case where the discrete velocity is an interpolation of the continuous material velocity, the Lagrangian setting, in many applications it proves computationally efficient to consider a mesh velocity which is different to the interpolated material velocity. In particular it appears that the arbitrary tangential velocity, that we consider in this study can be chosen such that beneficial mesh properties are observed in practice. This provides the motivation for this work in which we analyse an ESFEM where the material velocity of the mesh is different to (the interpolant of) the material velocity of the surface, i.e., an arbitrary Lagrangian-Eulerian ESFEM (ALE-ESFEM). We refer to [10] for extensive computational investigations of the ALE-ESFEM that we analyse in this study. For examples in the numerical simulation of mathematical models for cell motility and biomembranes, where the ALE approach proves computationally more robust than the Lagrangian approach, we refer to [11, 12, 13, 14].

Our main results are Theorems 4.3 and 5.4 where we show optimal order error bounds for the semidiscrete (space discrete, time continuous) and fully discrete numerical schemes. The fully discrete bound is proved for a second order backward difference time discretisation. An optimal error bound is also stated for an implicit Euler time discretisation. While the fully discrete bound is proved independently of the bound on the semidiscretisation, we believe that the analysis of the semidiscrete scheme may prove a useful starting point for the analysis of other time discretisations. We also observe that, under suitable assumptions on the evolution, the analysis holds for smooth flat surfaces, i.e, bulk domains with smooth boundary. Thus the analysis we present is also an analysis of ALE schemes for PDEs in evolving bulk domains.

We report on numerical simulations of the fully discrete scheme that support our theoretical results and illustrate that the arbitrary tangential velocity may be chosen such that the meshes generated during the evolution are more suitable than in the Lagrangian case. Proposing a general recipe for choosing the tangential velocity is a challenging task that is beyond the scope of this article and we do not address this issue. Moreover the choice of the tangential velocity and what constitutes a good computational mesh is likely to depend heavily on the specific application. We also investigate numerically the long time behaviour of solutions to (1.1) with different initial data when the evolution of the surface is a periodic function of time. Our numerical results indicate that in the example we consider the solution converges to the same periodic solution for different initial data.

The original ESFEM was formulated for a surface with a smooth material velocity that had both normal and tangential components [1]. Hence many of the results from the literature are applicable in the present setting of a smooth arbitrary tangential velocity.

2. Setup

We start by introducing an abstract notation in which we formulate the problem.

2.1 Definition (Bilinear forms).

For we define the following bilinear forms

 a(φ(⋅,t),ψ(⋅,t)) =∫Γ(t)∇Γ(t)φ(⋅,t)⋅∇Γ(t)ψ(⋅,t) m(φ(⋅,t),ψ(⋅,t)) =∫Γ(t)φ(⋅,t)ψ(⋅,t) g(φ(⋅,t),ψ(⋅,t);w(⋅,t)) =∫Γ(t)φ(⋅,t)ψ(⋅,t)∇Γ(t)⋅w(⋅,t) b(φ(⋅,t),ψ(⋅,t);w(⋅,t)) =∫Γ(t)φ(⋅,t)∇Γ(t)ψ(⋅,t)⋅w(⋅,t) ~a(φ(⋅,t),φ(⋅,t);va(⋅,t)) =∫Γ(t)(∇Γ(t)⋅va(⋅,t)−2DΓ(t)(va(⋅,t)))∇Γ(t)φ(⋅,t)⋅∇Γ(t)ψ(⋅,t) ~b(φ(⋅,t),ψ(⋅,t);w(⋅,t);va(⋅,t)) =∫Γ(t)∇Γ(t)⋅va(⋅,t)(φ(⋅,t)w(⋅,t)⋅∇Γ(t)ψ(⋅,t)) −∫Γh(t)φ(⋅,t)w(⋅,t)⋅BΓ(t)(va(⋅,t))∇Γ(t)ψ(⋅,t),

with the deformation tensors and as defined in Lemma A.1.

We may now write the equation (1.5) as

 (2.1) ddtm(u,φ)+a(u,φ)=m(u,∂∙vaφ)−b(u,φ;aT).

In [1] the authors showed existence of a weak solution to (1.3) and hence a weak solution exists to the (reformulated) problem (1.5), furthermore for sufficiently smooth initial data the authors proved the following estimate for the solution of (1.3) and hence of (1.5)

 (2.2) (2.3) ∫T0∥∂∙vu∥2L2(Γ(t))dt+supt∈(0,T)∥∇Γu∥L2(Γ)≤c∥u0∥2H1(Γ0).

We immediately conclude that as the bound (2.3) holds with the material derivative with respect to the material velocity replaced with the material derivative with respect to the ALE-velocity. See [15, 16, 17] for further discussion on the well-posedness of the weak formulation of the continuous problem.

3. Surface finite element discretisation

3.1. Surface discretisation

The smooth surface is interpolated at nodes () by a discrete evolving surface . These nodes move with velocity and hence the nodes of the discrete surface remain on the surface for all . The discrete surface,

 Γh(t)=⋃K(t)∈Th(t)K(t)

is the union of -dimensional simplices that is assumed to form an admissible triangulation ; see [9, §4.1] for details. We assume that the maximum diameter of the simplices is bounded uniformly in time and we denote this bound by which we refer to as the mesh size.

We assume that for each point on the exists a unique point on such that for (see [9, Lemma 2.8] for sufficient conditions such that this assumption holds)

 (3.1) x=p(x,t)+d(x,t)ν(p(x,t),t),

where is the oriented distance function to (see [9, §2.2] for details).

For a continuous function defined on we define its lift onto by extending constantly in the normal direction (to the continuous surface) as follows

 (3.2) ηlh(p,t)=ηh(x(p,t),t)for p∈Γ(t),

where is defined by (3.1).

We assume that the triangulated and continuous surfaces are such that for each simplex there is a unique , whose edges are the unique projections of the edges of onto . The union of the induces an exact triangulation of with curved edges. We refer, for example, to [9, §4.1] for further details.

We also find it convenient to introduce the discrete space-time surface

 Gh,T:=⋃t∈[0,T]Γh(t)×{t}.
3.2 Definition (Surface finite element spaces).

For each we define the finite element spaces together with their associated lifted finite element spaces

 Sh(t) ={Φ∈C0(Γh(t))|Φ|K is linear affine for each K∈Th(t)}, Slh(t) ={φ=Φl|Φ∈Sh(t)}.

Let () be the nodal basis of , so that, denoting by the vertices of , . The discrete surface moves with the piecewise linear velocity and by we denote the interpolant of the arbitrary tangential velocity

 (3.3) Vah(x,t) =J∑j=1va(Xj(t),t)χj(x,t), (3.4) Tah(x,t) =J∑j=1aT(Xj(t),t)χj(x,t).

The discrete surface gradient is defined piecewise on each surface simplex as

 ∇Γhg=∇g−∇g⋅νhνh,

where denotes the normal to the discrete surface defined element wise.

We now relate the material velocity of the triangulated surface to the material velocity of the smooth triangulated surface. For each on there is a unique with

 (3.5) ddtY(t)=∂tp(X(t),t)+Vah(X(t),t)⋅∇p(X(t),t)=:vah(p(X(t),t),t),

where is as in (3.1). We note that is not the interpolant of the velocity into the space (c.f., [2, Remark 4.4]). We denote by the lift of the velocity to the smooth surface.

3.3 Remark (Surfaces with boundary).

While the method is directly applicable to surfaces with boundary, for the analysis we require the lift of the triangulated surface to be the smooth surface i.e., . Thus in general we must allow the faces of elements on the boundary of the triangulated surface to be curved. For the natural boundary conditions we consider it is possible to define a conforming piecewise linear finite element space on a triangulation with curved boundary elements, see for example [19], assuming this setup and neglecting the variational crime committed in integrating over curved faces the analysis we present in the subsequent sections remains valid. However as the surface is evolving a further requirement is that the continuous material velocity and the material velocity of the smooth triangulated (lifted) surface must satisfy

 (3.6) (v−vah)⋅μ=0 on% ∂Γ(t),

where is the conormal to the surface. If (3.6) does not hold, the additional error due to domain approximation must also be estimated. We remark that this issue is not specific to the ALE scheme we consider in this study and also arises if we take , i.e., the Lagrangian setting.

We note that although restrictive the above assumptions are satisfied for some nontrivial examples that are of interest in applications. In §6 we present two such examples. In Example 6.1, we present an example of a moving surface with boundary where the lift of the polyhedral surface (with straight boundary faces) is the smooth surface. In Example 6.4, we present an example where the surface is the graph of a time dependent function over the unit disc. Here the boundary curve is fixed and the boundary edges of elements on the boundary of the triangulated surface are curved such that the triangulation of the boundary is exact.

3.4 Remark (ALE schemes for PDEs posed in bulk domains).

As a special case of a surface with boundary, the method is applicable to a moving domain in dimensions, that is when is a flat (i.e., the normal to the surface is constant) dimensional hypersurface in with smooth boundary. We note that the formulae for the discrete schemes (4.1), (5.1) and (5.63) are the same as in the case of a curved hypersurface. Under suitable assumptions on the velocity at the boundary, the analysis we present is valid in this setting. Specifically the analysis we present is valid for a flat three dimensional surface with zero normal velocity but nonzero tangential (and conormal) velocity (subject to (3.6)). In this case, as the domain is flat the geometric errors we estimate in the subsequent sections are zero (as since the lift is in the normal direction only). We note that this assumption necessitates curved boundary elements in this case. Therefore, as a consequence of our analysis we get an error estimate for an ALE scheme for a linear parabolic equation on an evolving three dimensional bulk domain in which all of the analysis is all carried out on the physical domain.

3.5. Material derivatives

We introduce the normal time derivative on a surface moving with material velocity defined by

 ∂∘η:=∂tη+v⋅νν⋅∇η,

and define the space

 H1(GT):={η∈L2(GT)|∇Γη∈L2(GT)|∂∘η∈L2(GT)}.

We are now in a position to define material derivatives on the triangulated surfaces. Given the velocity field and the associated velocity on we define discrete material derivatives on and element wise as follows, for with and ,

 (3.7) ∂∙h,VahΦh|K(t) =(∂tΦh+Vah⋅∇Φh)|K(t), (3.8) ∂∙h,vahφ|k(t) =(∂tφ+vah⋅∇φ)|k(t).

We find it convenient to introduce the spaces

 STh:={Φh and ∂∙VahΦh∈C0(Gh,T)|Φh(⋅,t)∈Sh(t),t∈[0,T]}

and

 (STh)l:={φh and ∂∙vahφh∈C0(GT)|φh(⋅,t)∈Slh(t),t∈[0,T]}.

The following transport property of the finite element basis functions was shown in [1, §5.2]

 (3.9) ∂∙h,Vahχj=0,∂∙h,vahχlj=0,

which implies that for with

 (3.10) ∂∙h,VahΦh(⋅,t)=J∑j=1(ddtΦj(t))χj(⋅,t),∂∙h,vahφh(⋅,t)=J∑j=1(ddtΦj(t))χj(⋅,t)l.

We now introduce the notation we need to formulate and analyse the fully discrete scheme. Let be a positive integer, we define the uniform timestep . For each we set . We also occasionally use the same shorthand for time dependent objects, e.g., and , etc. For a discrete time sequence , we introduce the notation

 (3.11) ∂τfn=1τ(fn+1−fn).

For we denote by and by . For , we set

 (3.12) χnj=χj(⋅,tn),χn,lj=χlj(⋅,tn)

and employ the notation

 (3.13) Φnh=J∑j=1Φnjχnj∈Snh,φnh=Φn,l∈Sn,lh.

Following [3] we find it convenient to define for and

 (3.14) Φ––n+αh(⋅,t)=J∑j=1Φn+αjχj(⋅,t)∈Sh(t), (3.15) φ––n+αh(⋅,t)=(Φ––n+αh(⋅,t))l∈Slh(t).

We now introduce a concept of material derivatives for time discrete functions as defined in [3]. Given and we define the time discrete material derivative as follows

 (3.16) ∂∙,τh,VahΦnh=J∑j=1∂τΦnjχnj∈Snh,∂∙,τh,vahφnh=J∑j=1∂τΦnjχn,lj∈Sn,lh.

The following observations are taken from [3, §2.2.3], for

 (3.17) ∂∙,τh,Vahχnj=0,∂∙,τh,vahχn,lj=0.

On , for

 (3.18) ∂∙h,VahΦ––n+αh=0,∂∙h,vahφ––n+αh=0,

which implies

 (3.19) Φ––n+1h(⋅,tn)=Φnh+τ∂∙,τh,VahΦnh,φ––n+1h(⋅,tn)=φnh+τ∂∙,τh,vahφnh.

We will also make use of the following notation, for given and , with lifts and , we define and , such that for

 (3.20) ΦLh(⋅,t)=tn+1−tτΦ––nh(⋅,t)+t−tnτΦ––n+1h(⋅,t), (3.21) φLh(⋅,t)=tn+1−tτφ––nh(⋅,t)+t−tnτφ––n+1h(⋅,t).

We note that (3.18) implies

 (3.22) ∂∙h,VahΦLh(⋅,t)=1τ(Φ––n+1h(⋅,t)−Φ––nh(⋅,t)), (3.23) ∂∙h,vahφLh(⋅,t)=1τ(φ––n+1h(⋅,t)−–φnh(⋅,t)).
3.6 Definition (Discrete bilinear forms).

We define the analogous bilinear forms to those defined in Definition 2.1 as follows, for , and

 ah(Φh(⋅,t),Ψh(⋅,t)) =∫Γh(t)∇Γh(t)Φh(⋅,t)⋅∇Γh(t)Ψh(⋅,t) mh(Φh(⋅,t),Ψh(⋅,t)) =∫Γh(t)Φh(⋅,t)Ψh(⋅,t) gh(Φh(⋅,t),Ψh(⋅,t);Wh(⋅,t)) =∫Γh(t)Φh(⋅,t)Ψh(⋅,t)∇Γh(t)⋅Wh(⋅,t) bh(Φh(⋅,t),Ψh(⋅,t);Wh(⋅,t)) =∫Γh(t)Φh(⋅,t)Wh(⋅,t)⋅∇Γh(t)Ψh(⋅,t) ~ah(Φh(⋅,t),Ψh(⋅,t);Vah(⋅,t)) = ∫Γh(t)(∇Γh⋅Vah(⋅,t)−2 DΓh(Vah(⋅,t)))∇Γh(t)Φh(⋅,t)⋅∇Γh(t)Ψh(⋅,t) ~bh(Φh(⋅,t),Ψh(⋅,t);Wh(⋅,t);Vah(⋅,t)) =∫Γh(t)∇Γh⋅Vah(⋅,t)(ΦWh(⋅,t)⋅∇ΓhΨh(⋅,t)) −∫Γh(t)Φ(⋅,t)Wh(⋅,t)⋅BΓh(Vah(⋅,t))∇ΓhΨh(⋅,t),

with the deformation tensors and as defined in Lemma A.1.

3.7. Transport formula

We recall some results proved in [2] and [3] that state (time continuous) transport formulas on the triangulated surfaces and define an adequate notion of discrete in time transport formulas and certain corollaries. The proofs of the transport formulas on the lifted surface (i.e., the smooth surface) follow from the formula given in Lemma A.1, the corresponding proofs on the triangulated surface follow once we note that we may apply the same transport formula stated in Lemma A.1 (with the velocity of replacing the velocity of ) element by element.

We note the transport formula are shown for a triangulated surface with a material velocity that is the interpolant of a velocity that has both normal and tangential components. Hence the formula may be applied directly to the present setting where the velocity of the triangulated surface is the interpolant of the velocity .

3.8 Lemma (Triangulated surface transport formula).

Let be an evolving admissible triangulated surface with material velocity . Then for ,

 (3.24) ddtmh(Φh,Ψh) =mh(∂∙h,VahΦh,Ψh)+mh(∂∙h,VahΨh,Φh)+gh(Φh,Ψh;Vah) (3.25) ddtah(Φh,Ψh) =ah(∂∙h,VahΦh,Ψh)+ah(∂∙h,VahΨh,Φh)+~ah(Φh,Ψh;Vah) (3.26) ddtbh(Φh,Ψh;Wh) =bh(∂∙h,VahΦh,Ψh;Wh)+bh(Φh,∂∙h,VahΨh;Wh) +bh(Φh,Ψh;∂∙VahWh)+~bh(Φh,Ψh;Wh;Vah).

Let be an evolving surface made up of curved elements whose edges move with velocity . Then for ,

 (3.27) ddtm(φ,ψ) =m(∂∙h,vahφ,ψ)+m(φ,∂∙h,vahψ)+g(φ,ψ;vah) (3.28) ddta(φ,ψ) =a(∂∙h,vahφ,ψ)+a(φ,∂∙h,vahψ)+~a(φ,ψ;vah) (3.29) ddtb(φ,ψ;w) =b(∂∙h,vahφ,ψ;w)+b(φ,∂∙h,vahψ;w) +b(φ,ψ;∂∙h,vahw)+~b(φ,ψ;w;vah).

We find it convenient to introduce the following notation for and and for a given and corresponding lift

 (3.30) L2,h(Wh,Φn+1)= 32τ(mh(Wh(⋅,tn+1),Φ––n+1h(⋅,tn+1))−mh(Wh(⋅,tn),Φ––n+1h(⋅,tn))) −12τ(mh(Wh(⋅,tn),Φ––n+1h(⋅,tn))−mh(Wh(⋅,tn−1),Φ––n+1h(⋅,tn−1))), (3.31) L2(wh,φn+1)= 32τ(m(wh(⋅,tn+1),φ––n+1h(⋅,tn+1))−m(wh(⋅,tn),φ––n+1h(⋅,tn))) −12τ(m(wh(⋅,tn),φ––n+1h(⋅,tn))−m(wh(⋅,tn−1),φ––n+1h(⋅,tn−1))).

The following Lemma defines an adequate notion of discrete in time transport and follows easily from the transport formula (3.24) and (3.18).

3.9 Lemma (Discrete in time transport formula).

For and and for a given and corresponding lift