FINITE ELEMENT METHODS FOR SURFACE DIFFUSION AND APPLICATIONS TO STRESSED EPITAXIAL FILMS Ricardo H. Nochetto Department of Mathematics and Institute for Physical Science and Technology University of Maryland, College Park, USA http://www.math.umd.edu/~rhn joint work with Eberhard B¨ ansch and Pedro Morin WIAS - Berlin, Germany Santa Fe, Argentina
Outline 1. Problem Description and Challenges 2. FEM for Surface Diffusion • Time Discretization • Variational Formulation • Space Discretization • Mesh regularization and adaptivity • Simulations • Graphs: Formulation, Estimates and Simulations 3. Stressed Epitaxial Films • Coupling: 1st Approach • Coupling: 2nd Approach • Simulations • Related Issues and Open Problems Ricardo H. Nochetto - NIST, May 4, 2004
1. Problem Description Physical problem: morphological changes in epitaxial thin films Γ σ σ film substrate Missfit between crystalline structures (linear) elasticity in bulk plus surface diffusion on free boundary ⇒ large deformations of Γ( t ) = morphological instabilities ⇒ crack formation and fracture ⇒ Ricardo H. Nochetto - NIST, May 4, 2004
Simplest Model � � Dynamics of free surface Γ( t ) V = − ∆ Γ κ − ε � V = normal velocity ∆ Γ = surface Laplacian κ = mean curvature ε = elastic energy density • First step: Understand the purely geometric PDE V = − ∆ Γ κ ( ε = 0 or given ) Surface diffusion � Related work: U.F. Mayer; Falk et al.; Deckelnick/Dziuk/Elliott; Sethian; Smereka. • Second step: Couple with elasticity ( ε = solution of a problem in the bulk) Ricardo H. Nochetto - NIST, May 4, 2004
Basic Properties for Closed Surfaces • Volume conservation � � � d dt | Ω( t ) | = V = − ∆ Γ ( κ + ε ) = ∇ Γ ( κ + ε ) · ∇ Γ 1 = 0 . Γ( t ) Γ( t ) Γ( t ) • Area decrease (for ε = 0 ) d � � |∇ Γ κ | 2 ≤ 0 . dt | Γ( t ) | = − V κ = − Γ( t ) Γ( t ) • A surface that starts as a graph may cease to be such in finite time. • A closed embedded hypersurface may selfintersect in finite time. Ricardo H. Nochetto - NIST, May 4, 2004
Numerical Challenges • Definition of curvature κ for a discrete surface • Definition of ∆ Γ κ : surface laplacian of a discrete variable • 4th order problem • Lack of maximum principle • Volume conservation • Area decrease • Stability • Error Analysis Ricardo H. Nochetto - NIST, May 4, 2004
2. General (closed) Surfaces Issue: How to deal with V = − ∆ Γ κ ν = ∆ Γ � Basic identity for Γ given: κ := κ� � X KEY IDEA: write the problem in the scalar and vector quantities � � ⇒ κ = � κ · � V = V � � κ, κ, V , V ν, ν κ = ∆ Γ � � X κ = � κ · � ν ⇒ (Mixed Method) V = − ∆ Γ κ � V = V � ν Ricardo H. Nochetto - NIST, May 4, 2004
Time Discretization: Semi-Implicit Given Γ n , describe Γ n +1 as the image of a mapping defined on Γ n : Γ n − → Γ n +1 , � → � X + τ � V n +1 X − Semi-Implicit Discretization: Take Γ n as a fixed domain ν on Γ n • Compute ∆ Γ and � = ⇒ • Take � X implicitly in the curvature equation: κ n +1 := ∆ Γ � X n +1 = ∆ Γ ( � X n + τ � V n +1 ) Ricardo H. Nochetto - NIST, May 4, 2004
Time Discretization: Semi-Implicit κ n +1 = ∆ Γ n ( � X n + τ � V n +1 ) κ = ∆ Γ � � � X κ n +1 = � κ n +1 · � ν n κ = � κ · � ν ⇒ V n +1 = − ∆ Γ n κ n +1 V = − ∆ Γ κ V n +1 = V n +1 � � � ν n V = V � ν κ n +1 − τ ∆ Γ n � V n +1 = ∆ Γ n � X n � κ n +1 − � κ n +1 · � ν n = 0 V n +1 + ∆ Γ n κ n +1 = 0 V n +1 − V n +1 � ν n = 0 � Ricardo H. Nochetto - NIST, May 4, 2004
Variational Formulation � Γ := Γ n , V (Γ) := H 1 (Γ) , V (Γ) := V (Γ) d , κ n +1 ∈ � V n +1 , κ n +1 ∈ V (Γ) � V n +1 ,� Seek V (Γ) , s.t. � � � � � � κ n +1 , � ∇ Γ � V n +1 , ∇ Γ � ∇ Γ � X n , ∇ Γ � ∀ � φ ∈ � � φ + τ φ = − φ V (Γ) κ n +1 · � κ n +1 , φ � � � � − = 0 ∀ φ ∈ V (Γ) � ν, φ V n +1 , φ ∇ Γ κ n +1 , ∇ Γ φ � � � � − = 0 ∀ φ ∈ V (Γ) � � � � V n +1 , � ν · � ∀ � � V n +1 , � φ ∈ � − = 0 V (Γ) φ φ � Γ n V n +1 = 0 V n +1 , 1 � � = = ⇒ discrete volume conservation Ricardo H. Nochetto - NIST, May 4, 2004
Finite Element Discretization V h (Γ) ⊆ � � Γ = Γ n h , V h (Γ) ⊆ V (Γ) , V (Γ) . κ n +1 ∈ � V n +1 , κ n +1 ∈ V h (Γ) � V n +1 ,� Seek V h (Γ) , s.t. � � � � � � κ n +1 , � ∇ Γ � V n +1 , ∇ Γ � ∇ Γ � X n , ∇ Γ � ∀ � φ h ∈ � � φ h + τ φ h = − φ h V h (Γ) κ n +1 · � κ n +1 , φ h � � � � − = 0 ∀ φ h ∈ V h (Γ) � ν, φ h V n +1 , φ h ∇ Γ κ n +1 , ∇ Γ φ h � � � � − = 0 ∀ φ h ∈ V h (Γ) � � � � V n +1 , � � V n +1 , � ν · � ∀ � φ h ∈ � − = 0 V h (Γ) φ h φ h Γ n V n +1 = 0 V n +1 , 1 � � � = = ⇒ discrete volume conservation • S κ n +1 | 2 ≤ | Γ n | | Γ n +1 | + τ n � Γ n |∇ = ⇒ area decrease + stability • Ricardo H. Nochetto - NIST, May 4, 2004
Nodal Representation and Schur Complement τ � � � − � A� 0 0 A M X n V 0 − A 0 M K 0 = � − � � 0 0 M N 0 K − � N T 0 0 0 V M Schur complement for V : � � N T � M − 1 � M − 1 � N T � M − 1 � τ � A � Q V = − Q � A X n N + MSM Q on ker( A ) ⊥ S is the inverse of A | ker( A ) ⊥ : AS = I = SA Q is the L 2 (Γ) projection onto X h (Γ) = { φ ∈ V h (Γ) : � Γ φ = 0 } The system is symmetric and positive definite ⇒ Solvability Ricardo H. Nochetto - NIST, May 4, 2004
(Basic) Final Procedure 1. Let T be the initial triangulation of Γ with nodes � X . 2. Build the matrices A , � A , M , � M , � N . 3. Solve for V the system � N T � M − 1 � M − 1 � � N T � M − 1 � τ � A � Q V = − Q � N + MSM Q A X . 4. Solve for � M � � V = � V the system: N V . 5. Update � X ← � X + τ � V . 6. Go to step 2. Ricardo H. Nochetto - NIST, May 4, 2004
Mesh Regularization Regularization sweep 1. For each node z of the mesh do the following: (a) Compute a normal � ν z to the node z . The area of the shaded triangle (b) Compute a weighted average ˆ z of coincides with that of the triangle marked with thick lines. all the vertices that belong to Then the area of the whole bulk remains unchanged. the star centered at z . (c) Consider the line that passes midpoint of through ˆ the element z in the direction of midpoint of \hat{z} the element the normal � ν z . Replace the node z by the only point belonging to this line that keeps the volume new vertex \tilde{z} vertex to update z enclosed by the surface unchanged. direction of normal \nu_z Ricardo H. Nochetto - NIST, May 4, 2004
Timestep Control Two goals: 1. Prevent large timesteps for which the position change of a node, is larger than the element size (to avoid crossing). 2. Allow large timesteps when the normal velocity does not exhibit large variations. Relative position change = τ | � V ( z 0 ) − � V ( z ) | ≈ τh T |∇ Γ � V |≈ ǫ t h T ǫ t Compute ρ = and try to use τ ≈ ρ max |∇ Γ V | Ricardo H. Nochetto - NIST, May 4, 2004
Space Adaptivity Goal: Have an accurate representation of Γ in the sense that the density of nodes should correlate with the local variation (regularity) of Γ . We achieve this by enforcing h S | ∠ ( � ν 2 ) | ≈ α ν 1 ,� on every side S of the mesh. Angle Width Control Split those elements with an angle wider than a certain α max . Ricardo H. Nochetto - NIST, May 4, 2004
(Natural) Boundary Conditions t = 0 . 113 × 10 − 5 t = 0 . 932 × 10 − 5 t = 0 t = 0 . 4300 × 10 − 4 t = 0 . 35039 × 10 − 3 t = 0 . 31211 × 10 − 2 t = 0 . 02545 t = 0 . 07545 t = 0 . 12545 Ricardo H. Nochetto - NIST, May 4, 2004
Features of Final Procedure • Consistent approximation, no smoothing of normals etc. needed • Only C 0 regularity for the finite element spaces • Arbitrary polynomial degree for the finite element spaces • Nearly volume conservative (exact volume conservation in the graph case) • Area decrease / stability • Time/Space Adaptation and volume conservative Mesh Regularization • Simulations using ALBERT with P 1 elements (A. Schmidt and K. Siebert) and GEOMVIEW (Geomety Center-Minneapolis) Ricardo H. Nochetto - NIST, May 4, 2004
Volume Conservation and Area Decrease 1.014 1 1.012 Cube 0.95 Prism-16-1 1.01 Prism-8-1 0.9 Prism-4-1 1.008 0.85 1.006 0.8 1.004 Prism-4-1 0.75 1.002 Cube Prism-8-1 0.7 1 Prism-16-1 0.998 0.65 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 Volume Area Relative volume and surface area with respect to the initial values vs. time. The computations were performed with the full adaptive algorithm. Ricardo H. Nochetto - NIST, May 4, 2004
Recommend
More recommend
