Finite Element Methods for Surface Diffusion

✂✁✄✆☎✞✝✠✟✡☎☞☛✍✌✏✎✒✑✔✓✕✝✠✟☞✁✖✟☞✗✘✟ ✙✛✚✢✜ ✣ ✤✢✥✧✦✩★ ✪✫✤✭✬✂✮✯✦ ✣ ✤✢✪✧✰✲✱✂✳✵✴✶✳ ✷✸✤✭✬ ✹✺✮✯✻✆✼✾✽✿✪✫✳❀✮✯✴❂❁ ❛❝❜❚❞✶❡❢❜❤❣❥✐❧❦ ❃❅❄❇❆❉❈❂❊●❋■❍❑❏▼▲❖◆❖P❂❋❘◗❂❙❚❊■❯❖▲❖◆❖❱❳❲❀❙❚❊■❨❩❃❩◆❬❙❪❭❴❫❵❭ ♠♦♥♣♥❢qsr✖t✈✉❧✇②①④③✖✇✛⑤✖⑤ ⑥⑧⑦⑩⑨❶⑦❸❷❺❹ ❻✿❼⑩❹❾❽ ❹❾⑨❿❷➁➀ ❹➂❷➄➃❶➅✠➆❬➇➁➈➉➅✒➊✕➋➍➌➎➊♣➈✲➏✒➐❀❹➒➑➓⑦❸➔✧➌❬➇→⑦⑩➅➣⑨ ↔❝↕➙❞✶❜❤➛➝➜✈❜❤➞➠➟✔➡✈✐❢➢➥➤❪➛➧➦❪➨➩❛❾❞✶➞♣❜❚➫✂➭✢➫✛❜❤❣➯✐♣➲➳➨➧➵✠❣❥➤✶➜✖❜❤➞➄➫✆➸✞➺➧q✠➫❺➤✲➛➄❞➻❦❤❦➥➫✖➼ ➽➚➾➻➪➻➶✾➹➴➘➷➘➷➬♦➮❖➱❺✃➥❐❤❒➉➾▼❐❤❮➷❰❬Ï❑Ð❤➘➷ÑÓÒ❤Ð➥ÐÕÔ Ö✭×❵Ø■ÙÚØ■ÛÝÜ❥Þ➯Ûàß❸ÜÝÜ❧áãâ➉Ü❥Þ➯ÙÚÞ➯ä➥ÞÝØ➝å➴æ❸Û❺ç❺è➉èÕéÚÙêØ■ëìç➄â➉ß♦éêí➥Ü❥ÙîÜ→ß♦â❪ëðï❚ÞÝæ➥ñàò➉ß❸Ü❥Þ➯Ùîñ❘Ü■óÕôìæ❸ò➉Û➯Ø■â➉ÜõÞÝÛÝß⑩Ü➯Ü➯Ø♣ö⑩÷ÕóÕø❳ù❺ú●ûÕú⑩úýü➙þ❧Ø■Û➯éêÙÚâ ß♦â➉ë❝ÿ➉Û➯Ø■ÙÚ✁Ø ❺â➉✄Ù ✂⑩Ø❘ÛàÜ❥ÙÚ✆Þ ☎♦Þ→þ❧Ø■Û➯éêÙÚâ✶ó✞✝❢Ø❘✆Û ✟➂ß❸â❤í ✠Ø ✟✵ß♦Ùê☛é ✡✌☞➉ß♦Ø■â➉Ü➯ñà✎ò ✍✑✏→Ùêß⑩Üõù✒☞❪Ø■Û➯éêÙê✔â ✓ ëÕØ ✕ ø➄Ø■è➉ß♦Û➯ÞÝ✖ß ✟➂Ø❘â❚Þ➯æ✵ë➥Ø➙ô✾ß♦Þ➯Ø❘✆Û ✟✘✗❑ÞÝÙêñ■ß❀ÿ❪ß❸ñ❘äÕé❴Þàß❸ëðë➥Ø➙áã✚â ✙⑩Ø❘âÕÙêØ❘ÛÝÙê✘ ß ✛♣äÕ✜Ù ✟➂Ùêñ■ßÕ✚ó ❺â➉✄Ù ⑩✂ Ø❘ÛàÜ❥ÙîëÕß❸ë✣➄✢ ß⑩ñ➚Ùêæ❸â➉ß❸é➻ë➥Ø❘✥é ✤✶Ù❴ÞÝæ❸Ûàß♦é ßØ➚♦✒ùâ➉✟✵ëðß♦áãÙêâ➉✱é ✡❖Ü❥Þ➯✚è ÙÚÞ➯✟➂äÕæ❸Þ➯æ❾ÛÝÙê✩â ë➥✍✑Ø➙✵✟ ô✾ß❑ß♦Þ➯Þ➯✥ò ✠Ø ✓✟✘äÕâ✶✗♦Þ➯ú✖Ùî✓ ñ❘Ø■ß➂ë➥✔ä ç❺✓ ß♦èÕÛ éêÙîñ❘ß❸ë➉ß❀ëÕØ❘✥é ✤✶Ù❴ÞÝæ❸Ûàß♦✌é ✦➴áõô✾✑ç ✤★✧➚✩ó ✝✫Õ✪ ✬Ø ✟❀Ø●Ü→✮ö ✭✰✯♦ûÕó➥ö⑩û❸û❸û➂ï➥ß❸â⑩Þàß➂ÿÕØ❸ó❪ç❺✆Û ✙❸Ø❘â❚ÞÝÙÚâ➉ß ✲ ø➄Ø■è➉ß♦Û➯✳Þ ✟➂Ø■â⑩Þ❺æ❸å ô✾ß♦Þ➯òÕ✬Ø ✟✵ß❑Þ➯Ùîñ❘Ü→ß❸â➉ëðáãâ➉Ü❥Þ➯ÙÚÞ➯ä➥ÞÝØ➝å➴æ❸Û✵✴✈ò❤í➥Ü❥Ùîñ❘ß❸é✶ï➥ñ➚ÙêØ❘â❪ñ➚Ø➩ß♦â❪ë✷✶ Ø■ñàòÕâ➉æ❸éê✖æ ✙❸í⑩✸ó ❺âÕ✜Ù ✂❸Ø■ÛÝÜ➯Ù❴Þõí æ♦å ô✒ß♦ÛÝí❚éîß♦â❪ë Ø➚✹ ✒ù æ❸✟✵éêéÚß♦✬Ø Ùê✙❸✱é ✺Ø✡❖ÛÝ✴❉òÕß♦✩â ✆Û ✍✑✻❂ó❤✟✵ô✒ß❑ÞÝ✽ø ✔ò ✓ ✼♦✞ä û❤✟➂✾ü ✔ë ✭✰✓✼➥Ø●✰ó ë➥♣ïä ➥ç ✿ ❈❖❭❁❀✩❂❄❃ ❲❀❙❚❊■❨❩❃❩◆❆❅❇❂✎❂✎❈ ✼♦û❸û⑩û✣❉❋❊❍●✱■✰❏✠❑▲❊✮●✱▼❖◆◗P✁❘❇❙✞❚❖❯✠❏❱◆◗●✑❲❨❳✜❊✾P❱P❩▼ ❬❭◆❱❊❍●☛▼❖❪❍❫✎❴✺ö❵✯✖❛❁✯✖✯Õó❝❜❝✯♦ô❿ú✾✼➥ó✸❜❵✯❸ô❿ú❞✯Õó✰❜❝✯♦ô❡❜⑩ûÕó✸❜❵✯❝❢❂û❵✯✚✓ ❏✠❤❥✐❦❪❍❧✳♠❍P✘❊✮❫✩♠♦♥✩■❵❧✳❊❍P✠❏❩P✠❾ ❴ ï❤ä➉Û❥å ß⑩ñ➚Ø➂ë➥qÙ ❂♣ ä➉Ü➯ÙÚæ⑩â▼ó➻å➴æ⑩äÕÛ❥ÞÝò➥ù æ❸Ûàë➥Ø❘Û❢è➉ß♦Ûàß✮✲☞ æ❸éêÙîñ➩èÕÛÝæ✖Õ☞ éê✠Ø ❝ ✟ ó❄❪r âÕÙ❴ÞÝØ➂Ø❘éêØ✠✟➂Ø❘â❚ÞÝÜ■ó▼ß èÕÛÝÙÚæ⑩Û➯❣ Ù❂Ø❘ÛÝÛÝæ❸Û✄Ø■Ü❥Þ➯✜Ù ✵✟ ß❑Þ➯Ø●Ü❘ó➉ïÕñàò❚ä➉Û❺ñ➚✖æ ➂✟ èÕéêØ✠➂✟ Ø■â⑩Þ●óÕÜs➂✟ æ❤æ♦ÞÝòÕÙêâ✚➂✙ ◗Ø ♣❂Ø■ñ➷Þ●óÕèÕÙêâ➉ñàò➥ù æ✮♣t✓ ✂✁☎✄✝✆✟✞✠✁☛✡✌☞ ✍✎✞✏✄✑✞✏✒✔✓✔✆✟✒✖✕✘✗✚✙✜✛✣✢☎✓✔✆✖✄✑✆✟✤✥✆✧✦✩★☎✒✧✪✫✢☎✬✘✞✠✭✮✕✘✢☎✁✯✆✖✞✰✪✱✢✲✕✘✳✑☞✴✓✔✄✵✓✶✤☎✢☎✁✸✷✴✆✟✹✴✺✼✻✲✕✘✓✽✆✟✄✑✾❀✿❁✍❂✛✣✪❃✷☎❄ ❅☛✹❆✻☎✒✖✞✠✢☎✓✔✆✖✒✟✕✘✗✘✞❈❇❆❉ ❊●❋ ❍❏■✯❍❆❍❏❑▼▲✮✞✏✒✔✳✵✄✵✢ ◆ ✞✏✒✔❖P✕✚✢◗☞ ❘❙✕❯❚❲❱ ❳❩❨❆❉❬❇✘■❪❭✚■✘❨❆❨✌❉✌❑✘❫ ❴✙✣❅❀✕✘✄✵✳❵❱ ❛☎✒✔✞✏❛☎✒✔✄✵✢✌✆❝❜✫✭✶✄❞✕✚✓✽✙✣✡❡✞✏✒✔✳✵✄✑✢❣❢❤✁✯✞ ✍✎✹❆✒✔✳✵✁✐✍❂✄✵✁☎✞❃✍✎✞✏✡❥❱❦✻✌✆✖✆✖❛❣❱♠❧✘❧❯✭✧✭✧✭♥❢♦✭✶✄✵✕✘✓✽✙✜✡❲✞✠✒✖✳✑✄✵✢❣❢♦✁☎✞❏❧ Finite Element Methods for Surface Diffusion Eberhard Bänsch, Pedro Morin, and Ricardo H. Nochetto Abstract. Surface diffusion is a (4th order highly nonlinear) geometric driven motion of a surface with normal velocity proportional to the surface Laplacian of mean curvature. We present a novel variational formulation for the parametric case, develop a finite element method, and propose a Schur complement approach to solve the resulting linear systems. We also introduce a new graph formulation and state an optimal a priori error estimate. We conclude with several significant simulations, some with pinch-off in finite time. Keywords: Surface diffusion, fourth-order parabolic problem, finite elements, a priori error estimates, Schur complement, smoothing effect, pinch-off. AMS subject classification: 35K55, 65M12, 65M15, 65M60, 65Z05. 1. Surface Diffusion and Mixed Formulation The overall goal of this work is to devise efficient numerical tools for simulating morphological changes in stressed epitaxial films and thereby study their complicated nonlinear dynamics. To model the misfit between the crystalline structure of the substrate and epitaxial film, the film may be thought of as subjected to mechanical stresses. This causes a plastic deformation of the free surface of the film. This morphological instability of the free surface may eventually lead to crack formation and fracture, an issue of paramount importance in Materials Science; see for instance [1, 4, 14] and also the list of references in [2]. The dynamics of the free surface Γ in Rd is governed by V = −∆S (κ + ε), (1) which is a 4th order highly nonlinear PDE. Hereafter, d ∈ {2, 3}, V and κ are the (scalar) normal velocity and mean curvature of Γ, respectively, ∇S is the surface gradient, ∆S = divS ∇S the Laplace-Beltrami operator and ε is the elastic energy density of the bulk Ω(t) enclosed by Γ(t). In this paper, we take ε = 0 throughout. A number of issues arise, from existence, well posedness and regularity to effective algorithms for (1). The chief mathematical and numerical difficulties arise from the 4th order nonlinear operator ∆S κ and the fact that one cannot work directly with the curvature vector ~κ as in [8]. In contrast to [5, 10], we have been able to derive, from a suitable semi-implicit time discretization, a variational finite element formulation for surfaces in Rd which involves the 4 unknowns scalar 2 E. Bänsch, P. Morin, and R. H. Nochetto ~ , and (scalar) normal velocity curvature κ, curvature vector ~κ, normal velocity V V as follows: ~ ~ = V ~ν , κ = ~κ · ~ν , V = −∆S κ, V (2) ~κ = ∆S X, where ~ν denotes the unit normal vector to Γ, pointing outward of the bulk enclosed by Γ. In view of this mixed formulation, which involves only second order operators, our finite element method solely requires continuity of the discrete functions; in particular we do not need C 1 elements to handle curvature but can accommodate any polynomial degree. A Schur complement approach, with a symmetric and positive definite operator, is used to reduce the system (2) to the single unknown V . This allows for an efficient solution technique via preconditioned CG. This paper is organized as follows. We introduce the time discretization and natural variational formulation in §2, and present its finite element discretization in §3. We discuss the ensuing algebraic problem along with a Schur complement approach to its solution in §4. We present in §5 the graph formulation together with an error estimate. We document the performance of parametric and graph methods in §6 via several simulations which exhibit singularities in finite time. 2. Time Discretization and Variational Formulation We now consider a semi-implicit time discretization of first order with time-step τ , representing the next surface Γn+1 in terms of the current surface Γn as follows: ~ n+1 = X ~ n + τ V~ . X (3) This time discretization implies that all geometric quantities such as ~ν are evaluated on the current surface Γn . Consequently, if ∆S denotes now the LaplaceBeltrami operator on Γ = Γn , then the only modification to (2) is as follows: ~ n. ~κ − τ ∆S V~ = ∆S X (4) To derive a weak formulation to (2) and (4), we simply multiply them by test functions and integrate. If V(Γ) := H 1 (Γ) and X (Γ) is the subspace of V(Γ) of ~ functions with mean value zero, we thus seek V~ ∈ V(Γ), κ ∈ V(Γ), ~κ ∈ X~ (Γ) and V ∈ X (Γ) such that E E D D ~ ∇S ϕ ∀ ϕ ∈ X~ (Γ), (5) h~κ, ϕi + τ ∇S V~ , ∇ϕ = − ∇S X, hκ, φi − h~κ · ~ν , φi = 0 ∀ φ ∈ V(Γ), (6) hV, φi − h∇S κ, ∇S φi = 0 D E ~ , ϕ − hV, ϕ · ~ν i = 0 V ∀ φ ∈ X (Γ), (7) ~ ∀ ϕ ∈ V(Γ). (8) Hereafter the symbol h·, ·i stands for the L2 scalar product over the current surface ~ =X ~ n ; we thus suppress the Γ = Γn , which is described by the vector function X superscripts n and n + 1 since no confusion arises. Since R V ∈ X (Γ) has mean value zero, we realize that volume is preserved in the sense Γ V = 0. Finite Element Methods for Surface Diffusion 3 3. Finite Element Discretization Let T be a regular but possibly graded mesh of triangular finite elements over the surface Γ, which from now on is assumed to be polyhedral. Let h denote the local meshsize of T . Let T ∈ T be a typical triangle and let ~νT = (νTk )dk=1 be the unit normal to T pointing outwards. We denote by ~ν the outward unit normal to Γ, which satisfies ~ν |T = ~νT for all T ∈ T , and is thus discontinuous across ∂T . Let {φi }Ii=1 be the set of canonical basis functions of the finite element space Vh (Γ) of piecewise linear functions over Γ; this is a conforming approximation of ~ ∈ Rd×d is the identity matrix, consider the following matrix entries V(Γ). If Id Mij := hφi , φj i , Aij := h∇S φi , ∇S φj i , ~ ~ ij := Mij Id, M ~ ij := φi , φj ν k N d k=1 , (9) (10) and corresponding mass and stiffness matrices M := (Mij )Ii.j=1 , ~ := (M ~ ij )I M i.j=1 , A := (Aij )Ii.j=1 , ~ := (N ~ ij )I N i.j=1 , ~ := (A ~ ij )Ii.j=1 . A (11) (12) ~ ,A ~ and N ~ possess matrix-valued entries and therefore the We point out that M matrix-vector product is understood in the following sense ~V ~ = M I X ~ ij V~j M j=1 I , i=1 ~V ~ , is itself a vector in Rd . ~ = (V ~i )I , as well as of M because each component of V i=1 The vector of nodal values of a finite element functionPis written in bold I face, namely V = (Vi )Ii=1 ∈ V := RI is equivalent to V = i=1 Vi φi ∈ Vh (Γ). Let Xh (Γ) be the subspace of Vh (Γ) of functions with mean value zero, and let X = {V ∈ V : V · M 1 = 0} be the corresponding subspace of V with 1 := (1)Ii=1 : V = I X Vi φi ∈ Xh (Γ) ⇔ V = (Vi )Ii=1 ∈ X. (13) i=1 Upon expanding the functions V , κ, V~ , ~κ in terms of the basis functions and testing against the latter, a discretization of system (5)–(8) can be ~ ∈ V, ~ ~κ ∈ V, ~ K ~ ∈ X, ~ V ∈ X such that written in matrix form as follows: find V     ~  ~ ~ ~X ~ τA 0 M 0 V −A     0 −A  0 M   K  =  0  .  (14)     M ~ ~ ~ ~ 0 0 −N 0  K ~T 0 V 0 M −N 0 {φi }Ii=1 We discuss the solvability of (14) and propose an algorithm for its solution in §4. 4 E. Bänsch, P. Morin, and R. H. Nochetto 4. Schur Complement Approach Consider the following vector equation with a (possibly singular) square block A:      A B U F = . C D Q G Let A be symmetric with (nontrivial) kernel ker(A). Then the range Y of A is the orthogonal complement of ker(A). Let S : Y → Y be the right inverse of A: AS = Id on Y. If P denotes the orthogonal projection onto ker(A), we have SAV = V − P V = (Id − P )V ∀ V ∈ RI . (15) The Schur complement equation for Q then reads (−CSB + D)Q + CP U = G − CSF . (16) Solvability of this system depends on the structure of the two terms on the left hand-side of (16). We intend to apply this splitting to (14), which involves dealing ~ and A on the diagonal. with the block containing A Since the kernel of A in (12) is Z = span{1}, with 1 = (1)Ii=1 , the range of A is Y = Z⊥ . Then the spaces X, defined in (13), and Y are related as follows: V ∈X ⇔ M V ∈ Y. (17) Let S : Y → Y be the right inverse of A and P : V → Z be the orthogonal projection into Z, thereby satisfying (15) with P = 1⊗1 1·1 . ~ , K]T and Q = [K, ~ V ]T , we Applying (16) to (14) with vectors U = [V obtain the symmetric system  1    1 ~~~ ~ ~S ~A ~X ~ +M ~ P~ V ~ ~ N K −τ M τ M SM . (18) = ~T MPK N −M SM V ~ this could be viewed as a ~A ~X ~ makes sense because A ~X ~ ∈ Y; We observe that S compatibility condition. To investigate the solvability of (18), we note that both ~ ∈X ~ and V ∈ X or, in view of (17), M ~K ~ ∈ Y, ~ and components of Q satisfy K ~ ~ ~ ~ ~ ~ M V ∈ Y. Since the upper left block of (18), M S M : X → M Y, is nonsingular ~ −1 A ~M ~ −1 , we can apply (16) again to arrive at with inverse M   ~TM ~ −1 A ~M ~ −1 N ~ + M SM V + M P K = −N ~ TM ~ −1 A ~ X. ~ τN (19) Let Π := Id − M1⊗M1 M1·M1 be the orthogonal projection onto X. Since M P K ∈ span {M 1} = X⊥ , applying Π to (19) yields the final form of the Schur complement :   ~ ~TM ~ −1 A ~M ~ −1 N ~ + M SM ΠV = −ΠN ~TM ~ −1 A ~ X, Π τN (20) because ΠV = V . The ensuing matrix of this reduced system is nonsingular because ΠM SM Π : X → X is symmetric and positive definite, and the first term in (20) is symmetric and positive semi-definite. The method actually implemented consists of first solving for V using (20), ~ =N ~ and finally updating X ~ via X ~ + τV ~. ~V ~ V for V next solving M Finite Element Methods for Surface Diffusion 5 5. Graph Formulation and Error Analysis We now discuss the case of Γ being a graph. To this end let Ω ⊆ IR d−1 with t ∈ [0, T ] forpsome T > 0, d ∈ {2, 3} and Γ(t) := {(x, u(t, x)) | x ∈ Ω} ⊆ IR d . If Q(u) := 1 + |∇u|2 denotes the elementary surface area, then we have the following formulas for the geometric quantities ν, κ and V :  ∇u  ∂t u 1 , V = (−∇u, 1)T , κ = ∇ · , (21) ν= Q(u) Q(u) Q(u) with ∇ = (∂x1 , . . . , ∂xd−1 )T . According to (21), (1) can be written as the system  ∇u  ∂t u . (22) = −∆S κ, κ=∇· Q(u) Q This system has to be supplied with suitable initial and boundary conditions. We restrict ourselves to periodic boundary conditions for ease of presentation, but Neumann and Dirichlet boundary conditions can be handled as well [2]. In this vein, let Ω be a parallelogram and X be the subspace of H 1 (Ω) with periodic boundary values. Then, (22) admits the following variational formulation: find u and κ with u(t), κ(t) ∈ X for a.e. t and u(0) = u0 in a suitable sense, fulfilling Z h∂t u, ψi − ∇S κ · ∇S ψ = 0 ∀ ψ ∈ X, (23) ZΓ ∇u · ∇ϕ hκ, ϕi + =0 ∀ φ ∈ X, (24) Q(u) Ω where h·, ·i denotes the L2 (Ω) inner product; compare with §2. Note that, in contrast to the mean curvature flow of graphs [6, 7], (23) does not have the weight Q(u)−1 because the equation is written on Γ. If Xh is a finite element subspace of X , then a semi space discrete scheme is obvious from (23) and (24). RIn particular, upon taking ψ = 1 ∈ Xh we deduce d exact volume conservation dt u = 0. We can also prove stability as well as Ω h coercivity properties, which are crucial in deriving the following error estimate [2]. Theorem 5.1 (A Priori Error Estimate for the Semidiscrete Scheme). Let u, κ be sufficiently smooth in [0, T ], and let k ≥ 1 be the polynomial degree of Xh . Then, there exists a constant C > 0, solely depending on the regularity of u, κ and T , such that Z   2 sup ||(u − uh )(t)||L2 (Ω) + |∇S u − ∇S uh )|2 ≤ C h2k , t∈[0,T ] Z T 0 ||(κ − κh )(t)||2L2 (Ω) + Γh (t) Z Γh (t)  |∇S κ − ∇S κh |2 dt ≤ C h2k . (25) Since (25) involve the tangential gradients ∇S u and ∇S κ, the order of convergence O(hk ) is optimal. This is corroborated by the numerical experiments of [2], which are obtained via a semi-implicit time discretization similar to [7]. It consists of writing the equations for step n + 1 on the current surface Γn , which linearizes the system and allows for a Schur complement solver similar to that in §4. 6 E. Bänsch, P. Morin, and R. H. Nochetto 6. Implementation and Simulations In this section we explain briefly the implementation of both the parametric and graph formulations within ALBERT [12, 13], and document their performance. 6.1. Mesh Regularization Figure 1. Pathological ear formation in the evolution of a 4 × 1 × 1 prism by surface diffusion. Surface at time t = 0 and t = 0.07. Figure 2. Steps towards the pathological formation of ears. Zoom into a vertex of the initial prism with surfaces at t = k × 10−3 for 0 ≤ k ≤ 9. After 6 time steps some triangles collapse into segments, thereby making the mesh degenerate and producing numerical artifacts. The geometric flow by parametric surface diffusion is not as gentle as the corresponding mean curvature flow [8], and leads to severe mesh distortions. Even if the parametric formulation of §3 allows corners and edges, which are rather singular for surface diffusion, they give rise to fast node motion and mesh distortion. This is illustrated by the creation of ears during the evolution of a 4 × 1 × 1 prism in Figs. 1 and 2 towards a ball. This is clearly a numerical artifact and cannot be cured by mesh refinement and/or coarsening. We use mesh regularization instead, which consists of replacing each node by the projection onto the surface of Finite Element Methods for Surface Diffusion 7 a weighted average of the nodes belonging to its finite element star, see also [11]. This regularization is equivalent to one time step of the discretized (tangential) heat equation. Its beneficial effect is reflected in the subsequent simulations of Figs. 3-6, all computed with forcing term ǫ = 0. 6.2. Example 1: Smoothing Effect We illustrate the evolution of a 4×1×1 prism towards a ball with the same volume in Fig. 3, thereby showing the smoothing effect of surface diffusion. t=0 t = 0.08 t = 0.02 t = 0.04 t = 0.16 t = 0.32 Figure 3. Smoothing effect of surface diffusion. Evolution of a 4 × 1 × 1 prism towards a ball with equal volume at various time instants. 6.3. Example 2: Pinch-Off in Finite Time t=0 t = 0.32 t = 0.08 t = 0.36 t = 0.16 t = 0.38 Figure 4. Pinch-off in finite time. Evolution of an 8 × 1 × 1 prism at various time instants leading to a dumbbell and cusp formation. 8 E. Bänsch, P. Morin, and R. H. Nochetto Surface diffusion is not always a regularizing flow and may in fact create singularities in finite time depending on the initial configuration. This is depicted in Figs. 4-6 which correspond to prisms with larger aspect ratios than before. Figs. 4-5 display the evolution of an 8 × 1 × 1 prism towards a dumbbell and cusp formation in finite time. Increasing the aspect ratio of the prism seems to be an effective mechanism to produce several simultaneous cusps as reported in Fig. 6. t = 0.390 t = 0.394 t = 0.398 t = 0.402 t = 0.404 Figure 5. Detailed view of the pinch-off for the 8 × 1 × 1 prism. Mesh regularization appears to cure mesh distortion until the moment of pinch-off when the elements are rather elongated but not degenerate. t=0 t = 0.48 t = 0.16 t = 0.64 t = 0.32 t = 0.72 Figure 6. Evolution of a 12 × 1 × 1 prism towards two simultaneous cusps revealing that the number of singularities depends on the aspect ratio of the initial prism. 6.4. Example 3: Mushroom Formation We consider the graph formulation of §5 of a curve in 2d starting from the initial condition u0 (x) = 1 + δ(x), where δ is a steep perturbation with zero meanvalue (see Figure 7). Since the slope seems to become vertical around t = 4.8 × 10 −5 in Fig. 7, the classical solution might cease to exist in finite time. To investigate the formation of singularities in finite time, we use the parametric formulation with the same initial data upon embedding the graph of u0 into a closed curve (see Figure 8 top left). For the time scale of Fig. 7, the effect of Finite Element Methods for Surface Diffusion 9 this extension is negligible. Fig. 8 displays a sequence of solutions obtained for the same eight time instants of Fig. 7. It is worth noticing the striking similarity of the solutions obtained with both methods. Even though the parametric solution develops a mushroom shape at t = 9.6 × 10−5 , and thus the solution to the graph formulation is questionable thereon, they still exhibit an excellent quantitative agreement for t > 9.6 × 10−5 (compare the last two plots of Figs. 8 and 7). t=0 t = 8 × 10−6 t = 16 × 10−6 t = 24 × 10−6 t = 4.8 × 10−5 t = 9.6 × 10−5 t = 19.2 × 10−5 t = 38.4 × 10−5 Figure 7. Evolution of a graph in 2d starting from u0 (x) = 1 + δ(x) with a steep perturbation δ(x). In all plots for various times t, the x-axis ranges from −1 to 1, and the y-axis ranges from 0 to 1.5. t=0 t = 8 × 10−6 t = 16 × 10−6 t = 24 × 10−6 t = 4.8 × 10−5 t = 9.6 × 10−5 t = 19.2 × 10−5 t = 38.4 × 10−5 Figure 8. Mushroom formation. Parametric evolution of a curve in 2d with the same initial condition and the same times of Fig. 7. The rectangles in thin lines are [−1, 1] × [0, 1.5]. The effect of embedding the graph of Fig. 7 into a close curve (see top-left picture) is negligible for the time scale of this evolution. 10 E. Bänsch, P. Morin, and R. H. Nochetto Acknowledgments We would like to thank K. G. Siebert for his assistance in incorporating new data structures to handle surfaces in R3 and curves in R2 within ALBERT [12, 13]. We were all partially supported by the NSF-DAAD international cooperation grant INT-9910086. References [1] R. J. Asaro and W. A. Tiller, Surface morphology development during stress corrosion cracking: Part I: via surface diffusion, Metall. Trans., 3 (1972), 1789-1796. [2] E. Bänsch, P. Morin, and R.H. Nochetto, Surface diffusion of graphs: variational formulation, error analysis and simulations, SIAM J. Numer. Anal. (submitted). [3] E. Bänsch, P. Morin, and R.H. Nochetto, A finite element method for surface diffusion: The parametric case, (to appear). [4] J. Cahn and J. Taylor, Surface motion by surface diffusion, Acta Metall. Mater., 42 (1994), 1045-1063. [5] B.D. Coleman, R.S. Falk, and M. Moakher, Space-time finite element methods for surface diffusion with applications to the theory of stability of cylinders, SIAM J. Sci. Comp., 17 (1996), 1434-1448. [6] K. Deckelnick and G. Dziuk, Convergence of a finite element method for the non-parametric mean curvature flow, Numer. Math. 72 (1995), 197-222. [7] K. Deckelnick and G. Dziuk, Error estimates for a semi-implicit fully discrete finite element scheme for the mean curvature flow of graphs, Interfaces and Free Boundaries 2, (2000), 341-359. [8] G. Dziuk, An algorithm for evolutionary surfaces, Numer. Math., 58 (1991), 603-611. [9] K. Deckelnick, G. Dziuk, and C. Elliott, Semidiscrete finite elements for axially symmetric surface diffusion, CMAIA Research Report No 2002/05 (2002). [10] U.F. Mayer, Numerical solutions for the surface diffusion flow in three space dimensions, Comp. Appl. Math. (to appear). [11] A. Schmidt, Computation of three dimensional dendrites with finite elements, J. Comput. Physics, 125 (1996), 293-312. [12] A. Schmidt and K.G. Siebert, ALBERT: An Adaptive Hierarchical Finite Element Toolbox, Documentation, Preprint 06/2000 Universität Freiburg, 244 p; avaiulable online from http://www.mathematik.uni-freiburg.de/IAM/ALBERT. [13] A. Schmidt and K.G. Siebert, ALBERT–Software for scientific computations and applications, Acta Math. Univ. Comenian. (N.S.), 70 (2001), 105-122. [14] B.J. Spencer, S.H. Davis, and P.W. Voorhees, Morphological instability in epitaxially-strained dislocation-free solid films: nonlinear evolution, Phys. Rev. B, 47 (1993), 9760-9777. Finite Element Methods for Surface Diffusion 11 Weierstrass Institute for Applied Analysis and Stochastics Mohrenstrasse 39, 10117 Berlin, and Freie Universität Berlin, GERMANY E-mail address: [email protected] Instituto de Matemática Aplicada del Litoral, Güemes 3450, 3000 Santa Fe, and Departamento de Matemática, Facultad de Ingenierı́a Quı́mica, Universidad Nacional del Litoral, Santa Fe, ARGENTINA Partially supported by CONICET of Argentina and NSF Grant DMS-9971450. E-mail address: [email protected] Department of Mathematics and Institute for Physical Science and Technology, University of Maryland, College Park, MD 20742, USA Partially supported by NSF Grants DMS-9971450 and DMS-0204670. E-mail address: [email protected]