Simulation and Analysis of Coupled Surface and Grain Boundary Motion by Zhenguo Pan B.Sc., Shandong University, 1999 M.Sc., Shandong University, 2002 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF Doctor of Philosophy in The Faculty of Graduate Studies (Mathematics) The University of British Columbia (Vancouver) October, 2008 c© Zhenguo Pan 2008 Abstract At the microscopic level, many materials are made of smaller and randomly oriented grains. These grains are separated by grain boundaries which tend to decrease the electrical and thermal conductivity of the material. The motion of grain boundaries is an important phenomenon controlling the grain growth in materials processing and synthesis. Mathematical modeling and simulation is a powerful tool for studying the motion of grain boundaries. The research reported in this thesis is focused on the numerical simulation and analysis of a coupled surface and grain boundary motion which models the evolution of grain boundary and the diffusion of the free surface during the process of grain growth. The “quarter loop” geometry provides a convenient model for the study of this coupled motion. Two types of normal curve velocities are involved in this model: motion by mean curvature and motion by surface diffusion. They are coupled together at a triple junction. A front tracking method is used to simulate the migration. To describe the problem, different formulations are presented and discussed. A new formulation that comprises partial differen- tial equations and algebraic equations is proposed. It preserves arc length parametrization up to scaling and exhibits good numerical performance. This formulation is shown to be well-posed in a reduced, linear setting. Numeri- cal simulations are implemented and compared for all formulations. The new formulation is also applied to some other related problems. We investigate numerically the linear stability of the travelling wave so- lutions for the quarter loop problem and a simple grain boundary motion problem for both curves in two dimensions and surfaces in three dimensions. The numerical results give evidence that they are convectively stable. A class of high order three-phase boundary motion problems are also studied. We consider a region where three phase boundaries meet at a triple junction and evolve with specified normal velocities. A system of partial differential algebraic equations (PDAE) is proposed to describe this class of problems by extending the discussion for the coupled surface and grain ii Abstract boundary motion. The linear well-posedness of the system is analyzed and numerical simulations are performed. iii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . xii 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Numerical Methods . . . . . . . . . . . . . . . . . . . . . . . 10 1.3 Linear Well-posedness . . . . . . . . . . . . . . . . . . . . . . 13 1.4 Linear Stability of Travelling Waves . . . . . . . . . . . . . . 13 1.5 Thesis Overview . . . . . . . . . . . . . . . . . . . . . . . . . 15 2 Mathematical Formulation . . . . . . . . . . . . . . . . . . . . 17 2.1 Cartesian Formulation . . . . . . . . . . . . . . . . . . . . . . 17 2.2 Parabolic Formulation . . . . . . . . . . . . . . . . . . . . . . 21 2.2.1 Motion by Mean Curvature . . . . . . . . . . . . . . . 21 2.2.2 Motion by Surface Diffusion . . . . . . . . . . . . . . . 23 2.2.3 The Parabolic System . . . . . . . . . . . . . . . . . . 26 2.2.4 Boundary Conditions for the Parabolic System . . . . 27 2.2.5 Artificial Tangential Conditions . . . . . . . . . . . . . 28 2.2.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . 28 2.3 A Formulation of Partial Differential Algebraic Equations . . 29 2.3.1 Motion by Mean Curvature . . . . . . . . . . . . . . . 29 2.3.2 Motion by Surface Diffusion . . . . . . . . . . . . . . . 30 2.3.3 The PDAE System . . . . . . . . . . . . . . . . . . . . 31 2.3.4 Boundary Conditions for the PDAE System . . . . . . 32 iv Table of Contents 2.4 Exact Travelling Wave Solutions . . . . . . . . . . . . . . . . 33 2.4.1 Travelling Wave Solution for the Grain Boundary Mo- tion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.4.2 Travelling Wave Solution for the Coupled Motion . . . 34 2.5 Scaling the Motion Laws . . . . . . . . . . . . . . . . . . . . 38 3 Linear Well-posedness . . . . . . . . . . . . . . . . . . . . . . . 40 3.1 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.2 Well-posedness of the Parabolic System . . . . . . . . . . . . 42 3.2.1 Well-posedness . . . . . . . . . . . . . . . . . . . . . . 42 3.2.2 Analysis of Artificial Tangential Conditions . . . . . . 45 3.3 Well-posedness of the PDAE System . . . . . . . . . . . . . . 47 4 Numerical Simulation . . . . . . . . . . . . . . . . . . . . . . . 51 4.1 Cartesian Formulation . . . . . . . . . . . . . . . . . . . . . . 51 4.2 Parabolic Formulation . . . . . . . . . . . . . . . . . . . . . . 53 4.2.1 Discretization of the Motions Laws . . . . . . . . . . . 54 4.2.2 Discretization of Boundary Conditions . . . . . . . . . 55 4.2.3 Time Stepping . . . . . . . . . . . . . . . . . . . . . . 57 4.2.4 Numerical Results and Tangential Adjustment . . . . 58 4.3 PDAE Formulation . . . . . . . . . . . . . . . . . . . . . . . 61 4.3.1 Discretization of Boundary Conditions . . . . . . . . . 61 4.3.2 Numerical Results . . . . . . . . . . . . . . . . . . . . 62 4.4 More Examples . . . . . . . . . . . . . . . . . . . . . . . . . . 64 4.4.1 Examples of Surface Diffusion . . . . . . . . . . . . . . 65 4.4.2 More Examples of Coupled Surface and Grain Bound- ary Motion . . . . . . . . . . . . . . . . . . . . . . . . 67 4.4.3 Grain Boundary Motion with Triple Junction . . . . . 71 5 Stability of the Travelling Wave Solutions . . . . . . . . . . . 73 5.1 Theoretical Background . . . . . . . . . . . . . . . . . . . . . 74 5.1.1 Stability Criteria for Travelling Waves . . . . . . . . . 74 5.1.2 Spectrum on Unbounded Domains . . . . . . . . . . . 76 5.1.3 Spectrum on Large Bounded Domains . . . . . . . . . 77 5.1.4 A Simple Example . . . . . . . . . . . . . . . . . . . . 78 5.2 Linear Stability . . . . . . . . . . . . . . . . . . . . . . . . . . 79 5.3 Numerical Results on Stability . . . . . . . . . . . . . . . . . 81 5.3.1 Stability for the Grain Boundary Motion . . . . . . . 81 v Table of Contents 5.3.2 Stability for the Coupled Surface and Grain Boundary Motion . . . . . . . . . . . . . . . . . . . . . . . . . . 84 5.4 3D Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 5.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 6 Extension to a Class of High Order Problems . . . . . . . . 94 6.1 Mathematical Formulation . . . . . . . . . . . . . . . . . . . 94 6.2 Linear Well-posedness . . . . . . . . . . . . . . . . . . . . . . 97 6.2.1 Linearization of the System . . . . . . . . . . . . . . . 97 6.2.2 Equal Order Problems . . . . . . . . . . . . . . . . . . 99 6.2.3 Mixed Order Problems . . . . . . . . . . . . . . . . . 101 6.3 Numerical Simulation . . . . . . . . . . . . . . . . . . . . . . 103 6.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 7 Thesis Summary and Future Work . . . . . . . . . . . . . . . 110 7.1 Thesis Summary . . . . . . . . . . . . . . . . . . . . . . . . . 110 7.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 vi List of Tables 4.1 Estimated errors and convergence rates for the adjusted parabolic formulation with m = 0.5. Errors are evaluated at time t = 0.02. 60 4.2 Estimated errors and convergence rates for the PDAE formu- lation with m = 0.5. Errors are evaluated at t = 0.02. . . . . . 63 vii List of Figures 1.1 Bicrystal geometries. (a) quarter loop geometry; (b) Sun- Bauer geometry. . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2 Experimentally observed grain boundary shape using the Sun- Bauer technique . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.3 A thermal groove forms when grain boundary is exposed to a free surface. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.4 3D AFM (atomic force microscope) image of a groove devel- oped on the surface of tungsten upon 128h of annealing in vacuum at 1350◦C. This picture is published by Zhang and Gladwell in [67]. . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.5 The quarter loop bicrystal geometry. . . . . . . . . . . . . . . 7 1.6 A grain boundary motion with triple junction. All three curves undergo mean curvature motion. . . . . . . . . . . . . . . . . . 8 1.7 Video frames for the motion of grain boundaries with a triple junction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.1 A normal direction motion in cartesian coordinates with nor- mal velocity Vn. . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.2 Sketch of the grain boundary groove with angle description. . 19 2.3 An example when the right branch of the free surface is not single-valued. . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.4 Sketch of the grain boundary groove showing the curve num- bering. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 4.1 Sketch of the discretization at the triple junction showing a staggered grid and five ghost values. The triple junction is on the y-axis and the ghost values are shown by the five points on the dashed lines. . . . . . . . . . . . . . . . . . . . . . . . . 52 viii List of Figures 4.2 Numerical result for system (4.1.3) with m = 0.5,∆t = 1e−6. The initial spatial size is ∆s = 0.2 for the top free surface and ∆s = 0.004 for the grain boundary. Dotted line: numerical result; Solid line: travelling wave solution. . . . . . . . . . . . 53 4.3 Sketch of the ghost points at the triple junction for the dis- cretization of the parabolic formulation. . . . . . . . . . . . . 55 4.4 Plot of results for parabolic formulation using Backward Euler method with m = 0.5. Left: initial status with grid points; Right: result at t = 0.5 zoomed in near triple junction. Dotted line: numerical result; Solid line: travelling wave solution; Time step size: ∆t = 0.01; Initial spatial size: ∆s = 0.05. . . . 58 4.5 Plot of results for the parabolic formulation adjusted by equa- tion (4.2.16) with m = 0.5, α = −100,∆t = 0.01,∆s = 0.05. Left: result at t=0.5; Right: result at time t=0.5 zoomed in near triple junction. Dotted line: numerical result; Solid line: travelling wave solution. . . . . . . . . . . . . . . . . . . . . . 59 4.6 Comparison of grid distribution for the parabolic formulation and the adjusted parabolic formulation. Top: parabolic for- mulation; Bottom: adjusted parabolic formulation. Both pic- tures are zoomed in near triple junction. Dotted line: numer- ical result; Solid line: travelling wave solution; Time step size: ∆t = 0.01; Initial spatial size: ∆s = 0.05. . . . . . . . . . . . . 60 4.7 Sketch of ghost values at the triple junction for the discretiza- tion of the PDAE formulation. . . . . . . . . . . . . . . . . . . 62 4.8 Plot of result for the PDAE formulation withm = 0.5, zoomed in near triple junction. Dotted line: numerical result; Solid line: travelling wave solution; Time step size: ∆t = 0.01; Initial spatial size: ∆s = 0.05. . . . . . . . . . . . . . . . . . 63 4.9 Comparison of grid distribution for the parabolic formulation and the PDAE formulation. Top: parabolic formulation; Bot- tom: PDAE formulation. Both pictures are zoomed in near triple junction. Dotted line: numerical result; Solid line: trav- elling wave solution; Time step size: ∆t = 0.01; Initial spatial size: ∆s = 0.05. . . . . . . . . . . . . . . . . . . . . . . . . . . 64 4.10 Plot of result for the PDAE formulation with m = 1.96, which has non-single-valued free surface. Dotted line: numerical re- sult; Solid line: travelling wave solution. . . . . . . . . . . . . 64 ix List of Figures 4.11 A computational example that involves only the motion by surface diffusion solved by a PDAE formulation. ∆t = 5e − 7;∆s = 0.05; number of node points: 160. The enclosed area changes by 0.032%. . . . . . . . . . . . . . . . . . . . . . . . 66 4.12 A computational example that involves only the motion by surface diffusion solved by a PDAE formulation. ∆t = 1e − 5;∆s = 0.04; number of node points: 100. The enclosed area changes by 0.07%. . . . . . . . . . . . . . . . . . . . . . . . . 67 4.13 The Sun-Bauer geometry. . . . . . . . . . . . . . . . . . . . . 68 4.14 Sun-Bauer geometry simulation with β < γ/6 (PDAE formu- lation). Top: the groove and grain boundary profile; Middle: groove velocity versus time; Bottom: inclination angle versus time. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 4.15 Sun-Bauer geometry simulation with β > γ/6 (PDAE formu- lation). Top: the groove and grain boundary profile; Middle: groove velocity versus time; Bottom: inclination angle versus time. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 4.16 The evolution of film morphology during breakup. . . . . . . . 71 4.17 A computational example that involves only the motion by mean curvature solved by PDAE formulation. . . . . . . . . . 72 5.1 Eigenvalue distribution for the problem of grain boundary mo- tion. Left: distribution of eigenvalues; Right: Close-up near zero. The dots are numerical eigenvalues. The solid line and the star point are predictions. The solid parabola is given by x+ (y/c)2 = 0 and the star point is given by (−c2/4, 0). Note that both the middle line segment and the circle in the left picture are composed of numerical eigenvalues. . . . . . . . . . 82 5.2 Typical eigenvectors are plotted after being added to the trav- elling wave solution. Left: when eigenvalue is real number; Right: when eigenvalue is complex number. . . . . . . . . . . . 84 5.3 Eigenvalue distribution for the coupled surface and grain bound- ary motion. Left: distribution of eigenvalues; Right: Close-up near zero. The dots are numerical eigenvalues. The solid line and the star point are predictions. The solid parabola is given by x+ (y/c)4 = 0 and the star point is given by (−c2/4, 0). . . 85 x List of Figures 5.4 Typical eigenvectors are plotted after being added to the trav- elling wave solution. Left: when eigenvalue is real number; Right: when eigenvalue is complex number. . . . . . . . . . . . 86 5.5 Results for initial data that is a perturbation of a travelling wave. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 5.6 Extension of the 2D travelling wave solution to 3D. . . . . . . 90 5.7 Eigenvalue distribution for the coupled surface and grain bound- ary motion in 3D with wave number ω = 1. Left: distribution of eigenvalues; Right: Close-up near zero. The dots are numer- ical eigenvalues. The solid line and the star point are predic- tions. The solid parabola is given by x+(y/c)4/2+ω2(y/c)2+ ω4/2 = 0 and the star point is given by (−c2/2− ω2/2, 0). . . 91 6.1 A three-phase geometry. . . . . . . . . . . . . . . . . . . . . . 95 6.2 Sketch of the ghost values at the triple junction for the dis- cretization of the sixth order problem with m = 2. Two ghost points (X0 and X−1) and one ghost value (κ−1) are needed for each curve. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 6.3 Simulation of equal order problem with m = 2 (sixth order). Angles between curves at the triple junction are 2pi/3. . . . . 105 6.4 Simulation of mixed order problem with m1 = 0,m2 = m3 = 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 6.5 Plot of the junction speed versus time for the mixed order problem with m1 = 0,m2 = m3 = 2. The junction speed converges approximately to 1.78. . . . . . . . . . . . . . . . . 107 6.6 Simulation of mixed order problem with m1 = 1,m2 = m3 = 2. 109 xi Acknowledgements I would like to sincerely thank my supervisor, Brian Wetton, for his guidance, patience and encouragement throughout my studies at UBC. Without his devoted help this thesis would not have been possible. I am also grateful for the constant support and insightful comments of my supervisory committee, Uri Ascher, Anthony Peirce and Michael Ward. Special thanks to Amy Novick-Cohen whose research motivated this thesis work. She gave me insight into the research in the literature in this field and encouraged my own work. I would also like to thank Lloyd Bridge, Paul Chang, Yujin Guo and Tai-Peng Tsai who have been willing to answer my questions and give me help. Thanks also to the people in the Institute of Applied Mathematics and Department of Mathematics for the helpful discussions on math problems and the great time spent together. I gratefully acknowledge the financial support from NSERC and the De- partment of Mathematics for this work. I also wish to thank my family for their firm and continuous love and support. I am forever grateful to them all. Finally, thanks to all my friends and everyone who made this work pos- sible. xii Chapter 1 Introduction We study a coupled surface and grain boundary motion which arises in ma- terials science. It is an important phenomenon controlling the grain growth in materials processing and synthesis. In this thesis, we focus on the anal- ysis and numerical simulations of simplified mathematical models. A com- prehensive introduction to the physical background, including experiments, mechanisms and applications of the grain boundary motion is available in [26]. In this introduction, we shall first briefly introduce some background on the motion of grain boundaries, including the simplified models designed for both experimental investigations and numerical simulations. We then discuss available numerical methods for solving the problems. Other aspects that arise with the problems, such as the well-posedness of the formulations and the stability of travelling wave solutions are also described. 1.1 Background Materials science is an interdisciplinary field in which the properties of solid materials are studied. There is interest in how those properties are deter- mined by the materials’ composition and structure, both macroscopically and microscopically. In addition, the field includes attempts to improve the physical and chemical properties of the materials for best performance. With the increasing demand of novel materials this is becoming more and more important in many engineering fields, including electronics, aerospace, telecommunications, nuclear power and many others. With significant me- dia attention on nanoscience and nanotechnology in recent years, materials science has been propelled to the forefront of scientific research and evolved to the point at which scientists from physics, engineering and applied math- ematics are all working on problems of common interest. Mathematical modeling and simulation is an increasingly powerful tool being used, both for prediction and for obtaining a more thorough under- 1 Chapter 1. Introduction standing of material properties. It is increasingly being accepted as a third methodology in materials science, complementing the traditional approaches of theory and experiment. The Ginzburg-Landau model for superconduc- tivity is an example. As a theory introduced more than half a century ago it successfully predicted the existence of the so-called type II superconduc- tors before they were observed experimentally. This theory is now widely accepted in the study of both low-temperature and high-temperature super- conductors. However, due to the equations’ highly nonlinear nature, analytic solutions are impossible. Thus numerical simulations have become valuable tools in order to better understand the properties of the Ginzburg-Landau model and to provide further theoretical insight. It has attracted the inter- est of many computational scientists in recent years (see [14, 18, 19] and the references therein). Today, with the spectacular development in computer hardware and the new concepts and algorithms in computational mathematics, we can predict what properties a material will have and how it will behave under certain conditions based only on our knowledge of a single atom. For example, in [63], the properties of a novel class of semiconductors in the optoelec- tronics industry were addressed based on computational simulations before the materials are actually synthesized. In this way, one can learn how to improve the properties of known materials and design new materials. Nu- merical simulations are also used for determining the material properties if experimental conditions are not cost-efficient or practical. For instance, the study of materials under conditions existing only near the earth’s core ([33]). The macroscopic properties of nanophotonic and composite materials can also be predicted and improved through numerical simulation of microscopic quantities ([1, 13, 27, 64]). The work reported in this thesis covers a particular topic in materials science: grain growth in materials processing and synthesis. Our purpose is to provide an efficient mathematical approach to simulate this process. Therefore a better understanding and control of the resulting material are possible, though this is beyond the scope of this work. It was recognized that the properties of a metallic material are not pri- marily determined by the overall chemical composition, but rather by the distribution of chemical elements and crystal defects, i.e. by its microstruc- ture ([26]). The major goal of modern materials science is the optimization of the material microstructure to obtain the best properties. An important aspect of microstructure is the polycrystalline structure. Almost all common 2 Chapter 1. Introduction metals, and many ceramics are polycrystalline. A polycrystal is a material that is made of many smaller and randomly oriented crystallites which are often referred to as grains. The grains usually have the same physical and chemical properties but have different crystalline orientations. They are separated by grain boundaries, across which orien- tation change discontinuously. Grain boundaries are defects in the crystal structure. They are regions of high internal energy and tend to decrease the electrical and thermal conductivity of the material. Under certain con- ditions, e.g. when annealed (a heat treatment that alters the microstructure of a material), the grain boundaries migrate to reduce the interfacial en- ergy and heal the orientation mismatch and hence improve the physical and chemical properties of the material. The motion of grain boundaries is the dominant process of microstructural evolution during recrystallization and grain growth. There are various sources of driving forces for the migration of grain boundaries, such as stored deformation energy, grain boundary energy, mag- netic field, temperature gradient and so on. Thus grain boundary motion is an extremely complicated process. In order to produce a relatively simple driving force to study the mobility of grain boundaries, several bicrystal ge- ometries are designed, e.g. the quarter loop geometry introduced by Dunn et al. in [20] as shown in Figure 1.1(a) and the Sun-Bauer geometry, also called reversed-capillary technique, introduced by Sun and Bauer in [59, 60] as shown in Figure 1.1(b). In the quarter loop geometry, the grain boundary runs parallel to a free surface before it turns up and attaches to the upper free surface. In the Sun-Bauer geometry, an initially straight grain boundary attaches to the free surface with an inclination angle β. After a short time of migration, the inclination angle changes to α at the junction. Figure 1.2 shows the experimentally observed grain boundary shape using the Sun-Bauer technique. Note that, the free surface is located at the bottom of the picture instead of the top as shown in Figure 1.1(b). The material shown in the picture is electrical steel (specialty steel tailored to produce certain magnetic properties) containing about 3.5 wt% of silicon. The grain boundary (the interface between the black and the grey region) started from the solid line (only a part of it is shown) and migrated to the right and formed an angle different from the initial inclination angle. This picture was originally published in [23] and we have added the solid line to show the initial position of the grain boundary. 3 Chapter 1. Introduction Grain 2 Grain 1θ Grain 2 Grain 1 β α (a) (b) Figure 1.1: Bicrystal geometries. (a) quarter loop geometry; (b) Sun-Bauer geometry. In the quarter loop geometry, the grain boundary runs parallel to a free surface before it turns up and attaches to the upper free surface. In the Sun-Bauer geometry, an initially straight grain boundary attaches to the free surface with an inclination angle β. Figure 1.2: Experimentally observed grain boundary shape using the Sun- Bauer technique. This picture is published by Furtkamp et al. in [23]. The solid line is added to the original experimental picture to show the initial position of the grain boundary. However, if a grain boundary is exposed to a free surface, which is true 4 Chapter 1. Introduction Grain 1 Grain 2 grain boundary free surface Figure 1.3: A thermal groove forms when grain boundary is exposed to a free surface. for almost all bicrystal geometries, a new situation arises. The interaction between the grain boundary and the free surface yields a new surface defect – a thermal groove. Figure 1.3 shows a more realistic junction profile in which the groove root is lower than the asymptotic level of the free surface. Figure 1.4 shows an experimentally observed 3D groove root on the surface of tungsten ([67]). This defect can be interpreted as a balance between the surface tension of the free surface and the grain boundary. It serves as a dragging force to slow down the motion of the grain boundary. In the ideal case that grain boundaries do not undergo any progressive change in their structure or mode of motion, and thus the free energy per unit area and the mobility constant remain unchanged, the grain boundary migration is controlled by the curvature of the interface. It was first consid- ered by Mullins in [41] where he proposed that the boundary speed at any point will be proportional to the mean curvature at that point, i.e. Vn = Aκ (1.1.1) where A = mbγb is the reduced mobility (mb is the grain boundary mobility and γb is the grain boundary surface energy) and κ is mean curvature. This motion is called mean curvature motion. Not only in the area of materials science, mean curvature motion also plays an important role in image processing, such as noise removal and edge detection. The basic idea for noise removal is to take the noisy image as an intensity function which describes a surface and the height of the surface is 5 Chapter 1. Introduction Figure 1.4: 3D AFM (atomic force microscope) image of a groove developed on the surface of tungsten upon 128h of annealing in vacuum at 1350◦C. This picture is published by Zhang and Gladwell in [67]. the pixel value at that point. If we let each contour undergo motion by cur- vature, then very small contours, corresponding to spikes of noise, will disap- pear quickly. Real features will remain sharp because of the relatively lower curvatures at the boundaries. Edge detection can also be implemented in these schemes. These image processing applications are usually implemented by level set methods and detailed descriptions can be found in [45, 53]. Many mechanisms are involved in the transport of matter at the thermal groove, such as surface diffusion, volume diffusion, evaporation and conden- sation. However, only surface diffusion dominates for temperatures far below the melting temperature. Thus for most metal materials, surface diffusion is the most important mechanism for grain boundary grooving. During the process of surface diffusion, atoms move from areas of high mean curvature to areas of low mean curvature. Neglecting other effects and assuming that the properties of the interface do not depend on their orientation, Mullins formulated the equation to describe the motion for the free surface in [42, 43]. He stated that the normal velocity is proportional to the second derivative 6 Chapter 1. Introduction of mean curvature with respect to arc length, i.e. V = −Bκss (1.1.2) where B is a material constant and s is arc length. Note that in the case of surfaces in three dimensions (3D), κss in the equation above is replaced by the surface Laplacian of mean curvature ∇2sκ. A motion that follows this evolution law is called surface diffusion. Surface diffusion is one of the most important growth mechanisms in nanoscale surface growth. Strong interest in the theory and simulation has been initiated recently, driven by applications to electromigration of voids, shape transition in alloys, and in particular, by the growth of nanostructures (cf. [11]). As the major model discussed in this thesis, the quarter loop geometry has attracted special attention in both the past and present. In this model. the grain boundary runs parallel to a free surface and attaches to the free surface at a groove root as shown in Figure 1.5. When the grain boundary migrates to reduce the internal energy, the shape of the groove is maintained after the decay of initial transients, assuming the driving force is constant. Therefore a travelling wave solution exists for this problem. Figure 1.5: The quarter loop bicrystal geometry. In the literature, it has been assumed that the bicrystal is uniform along the cross-section direction. Thus it is reasonable for us to consider only 2D geometry. One can describe both the grain boundary motion and the diffusion of the free surface with the fundamental formulations proposed by 7 Chapter 1. Introduction Mullins in [41, 42, 43]. However, Mullins limited his consideration to two decoupled motions by assuming that the angle between the two free surfaces at the groove remains constant and is determined by the equilibrium values of the surface tension of the grain boundary and the free surface. This limitation is removed by Kanel et al. in [29] where a fully coupled nonlinear formulation was proposed to treat both the exterior free surface and the grain boundary simultaneously. The Sun-Bauer geometry can be treated analogously due to its similarity to the quarter loop geometry. Therefore we shall focus only on the quarter loop problem and a numerical example for the Sun-Bauer model will be shown. We also consider a system that involves only grain boundary motion as shown in Figure 1.6. In this model, three grain boundaries meet at a triple junction and evolve by mean curvature motion. Far away from the triple junction they run parallel to the horizontal direction. This system also admits a travelling wave solution which will be derived later. In the real situation of grain growth, triple junctions are unlikely to experience steady state motion. However, it is a useful, simple model we consider as an exercise for our linear travelling wave stability analysis. Grain 1 Grain 2 Grain 3θ θ Figure 1.6: A grain boundary motion with triple junction. All three curves undergo mean curvature motion. Figure 1.7 shows an experimental observation for the grain boundary model mentioned above. The material for this experiment is high purity Zn. The straight line grain boundary between Grain 1 and Grain 2 (cf. Figure 8 Chapter 1. Introduction Figure 1.7: Video frames for the motion of grain boundaries with a triple junction. One of the grain boundary is invisible, and its position is indicated by the straight solid line in the right bottom frame. The curved solid line is the numerical approximation of the grain boundary location. This picture is published by Czubayko et al. in [16]. 9 Chapter 1. Introduction 1.6) is not visible in the pictures. For convenience, we highlight its location by the straight solid line in the right bottom frame. The curved solid line is the numerical approximation of the grain boundary location. This picture is originally published by Czubayko et al. in [16]. 1.2 Numerical Methods In addition to experiments, mathematical modeling and numerical simula- tions are also crucial to the study of the motion of the grain boundaries. Many approaches have been developed to simulate the motion by mean curvature and the motion by surface diffusion, e.g., front tracking methods ([10, 35, 39, 68]), level set methods ([15, 34, 46, 54]) and the threshold dynam- ics algorithms (MBO algorithms developed by Merriman, Bence and Osher) ([22, 37, 36, 38, 50]). For an overview of numerical methods on curvature driven interface evolution problems we refer to [17]. Brief descriptions of these three approaches are given below. Conceptually, the front tracking methods are the simplest and most ef- ficient since they explicitly approximate the motion of the interface. More precisely, the location of the interface is tracked by moving a set of particles on the interface. However, it is usually impossible, or extremely difficult, to track topological changes. The authors in [10] attempted to track the topological changes for curve networks moving with curvature motion and reasonable results were achieved. In this case, topological changes are easy to detect. For general curve and surface tracking, this is not true. Level set methods are widely used for front propagation problems espe- cially when there are topological changes during the evolution. The interface is embedded as the zero level set of a signed distance function (level set func- tion). The level set function is evolved by an initial value partial differential equation in a higher dimension and the interface is matched with the zero level set of the level set function. In this setting, topological changes are captured in a natural manner. Level set methods provide a stable method for problems in both two and three dimensions. Detailed descriptions for level set methods can be found in the monographs by Sethian [53] and Os- her and Fedkiw [44]. Level set methods can handle smooth interfaces and also isolated (in time) events where smoothness is lost due to topological changes. The underlying formulation is that the interface separates two dis- tinct regions. The interfaces of junction problems that separate three regions 10 Chapter 1. Introduction and are not smooth at the junction, are not handled naturally by level set methods. The MBO algorithms uses a 1(inside)-0(outside) characteristic function to describe regions separated by the interface. Unlike level set methods, the characteristic function is evolved by diffusion and the resulting function is thresholded to locate the new position of the interface. The mean curvature motion, for example, is produced by the 1/2 level set of the resulting func- tion. These methods exhibit excellent stability properties and they can also be used for multiphase problems with triple junctions [50, 22]. For triple junction problems, a separate characteristic function is assigned to each re- gion and it is evolved independently according to the given velocity. The profile of the junction is captured by an appropriate coupling of all char- acteristic functions after each time step. However, this approach does not allow manual implementation of junction conditions and there are difficult unanswered questions on how to implement these techniques for mixed order problems. Similar approaches for the handling of triple junctions were also discussed in [69, 70]. When compared to mean curvature motion, surface diffusion is an intrinsi- cally difficult problem to solve numerically even in two dimensions [12, 15, 54]. It is stiff due to the fourth order derivatives and an explicit time-stepping strategy requires very small time steps. Moreover, owing to the lack of a maximum principle, a curve may self-intersect during its evolution. Smereka describes a level set method for curves that move with surface diffusion in [54]. Esedoglu et al. developed a MBO method to simulate the motion of surface diffusion with triple junctions in [22]. However, as noted above, both methods can not handle the quarter loop problem. Therefore, we have turned to a method that involves curve tracking. A parameterized representation for the curves is more helpful in this case. Normal direction motions specify only the normal velocity. That is enough to determine the location of the interface at any time uniquely . However, when a parametric representation is employed for front tracking methods, two unknowns (in 2D) are associated with each marker particle. Thus one equation is not enough to determine the position of each particle. To track a 2D curve, one also needs to specify the tangential velocity. Though the tangential velocity has no impact on the resulting shape of the curves, in- appropriate choices may lead to serious numerical difficulties, such as the overlap of grid points. In [4, 21], the authors presented formulations al- lowing only normal direction evolution and therefore redistribution of grid 11 Chapter 1. Introduction points in the tangential direction was needed to avoid coalescence. The im- portance of nontrivial tangential velocity has been emphasized and utilized in [10, 39]. In [10], the authors developed a parabolic formulation and utilized its smoothing property to adjust the grid points along the curve. Another approach for uniform grid distribution is proposed by Mikula and Sevcovic in [39]. They explicitly introduced a tangential velocity function into the system and derived an evolutionary partial differential equation to move the points tangentially such that a uniform grid was maintained. Other front tracking methods have also been developed based on vari- ational formulations and solved by the finite element method, for example [57, 58, 61] and more recently [5, 6]. The analysis and simulation of the mixed order quarter loop problem has been considered in different contexts. In many cases, linearized problems are considered as in [62, 65, 66]. In other cases, steady state or travelling wave solutions are considered [31, 29, 30, 32, 40]. Despite the continued experimental and mathematical interest in the mixed order problem (see also [20, 42, 43]), numerical methods that describe the general, nonlinear, time-dependent evolution of this phenomenon have not been developed. In this thesis, we contribute a robust numerical method using a finite difference discretization. Our method maintains a uniform distribution of grid points in arc length without redistribution by interpolation or using an evolution equation for the tangential velocity. Instead, we use an algebraic equation to adjust the position of grid points. This approach does not depend on the type of the motion and therefore the same equation for both the curvature motion and surface diffusion can be used. Other authors [5] recently and independently considered a similar numer- ical approach to the type of curve networks we consider here. Their method is based on a variational formulation of the problem and they used a finite el- ement method to approximate it. Careful consideration of their formulation shows that in the semi-discrete case (time continuous) their method leads to uniform distribution of grid points in a manner similar to our approach. However, this exact discrete property is lost in their fully discrete (space and time) scheme but retained in ours. In addition, our method is much easier to implement, allowing modifications to the junction conditions and normal motion in a straightforward way. Neither their method nor ours is particu- larly amenable to computations in three dimensions (3D) in the situations of interest where complex topological changes occur. 12 Chapter 1. Introduction 1.3 Linear Well-posedness In [9], a derivation for the mean curvature motion with a triple junction was obtained as a singular limit of a system of Ginzburg-Landau equations. The short time existence and uniqueness of the limit problem were investigated therein. A similar derivation and well-posedness analysis for the more com- plicated case when all three curves undergo surface diffusion were presented in [25]. The idea for the existence proofs in both studies is the linearization of the problem about the initial data and the verification of the existence of the resulting linear system using the fundamental theory for parabolic system in, e.g. [55]. The existence for the full nonlinear problem is then obtained by means of a fixed-point argument. However, there is a gap in the theory for the quarter loop problem, which is mixed order and the underlying theoretical estimates in [55] for the equal order problems are not available. We investigate the well-posedness of the mixed order problem by linearizing around fixed straight line solutions which gives a system that has the same highest order behaviour as the original prob- lem near the junction. The analysis gives the conditions that match those that more complicated nonlinear analysis gives as in [9, 25]. Therefore, we believe the results of the analysis should apply to the full nonlinear problem although the technical details are beyond the scope of this work. Such investigation is extended to a class of three-phase problems, in which three interfaces meet at a triple junction and evolve with specified normal velocities. Problems that are higher than the order of surface diffusion are considered for both equal order case and mixed order case. 1.4 Linear Stability of Travelling Waves Travelling waves are solutions to partial differential equations that move with constant speed while maintaining their shapes. Travelling waves arise in many applied problems. For the quarter loop geometry, a travelling wave solution was found in [30, 31] using an angle formulation for a range of physical parameter m which measures the relative surface tensions between the grain boundary and the free surface. The simpler problem that involves only grain boundary motion (see Figure 1.6) also admits a travelling solution as shown in section 2.4. We are interested in the stability of these travelling waves, that is the 13 Chapter 1. Introduction response of these solutions to small perturbations. If any such solution stays close to the translates of the travelling wave, then it is said to be stable. For a detailed classification of the type of travelling waves and the approaches for stability investigation, we refer to the review article [51] and the references therein. A frequently used approach for stability analysis is to linearize the partial differential equation about the solution and investigate the spectrum of the resulting linear differential operator. The spectrum is the union of point spectra, defined as the set of isolated eigenvalues with finite multiplicity, and its complement, the essential spectrum. Roughly speaking, if the spectrum of the linear operator lies in the left half plane, then small perturbations will not grow with time and therefore the travelling wave is stable. Note that the zero eigenvalue corresponding to the translates of the travelling wave will always be present. Unfortunately, for many applications, including the problems considered in this thesis, it is impossible, or very difficult to compute the spectrum analytically. In such a situation, numerical computations are often the only way to investigate stability. To investigate the stability of a travelling wave numerically, we linearize the system together with the junction conditions and discretize it using finite difference schemes. This yields a generalized eigenvalue problem. The eigen- values are then found numerically. However, numerical computations can only be done on a finite domain with artificial boundary conditions applied. Therefore another issue arises: how accurately can the numerical eigenvalues approximate the spectrum of the original differential operator on the infinite domain? This issue has been discussed in, for example, [7, 51, 52]. A com- mon conclusion is that the behaviour of the spectrum depends on the type of the artificial boundary conditions used. Some of the details will be discussed later in Chapter 5. However, the general theory does not fit our problem because of the pres- ence of the triple junctions. The numerical results give us some confidence that the travelling waves for both the grain boundary motion and the coupled motion are at least convectively stable, which means any unstable eigenmode can be stabilized by an appropriately chosen exponentially weighted norm. We do point out that there are no concrete theories to support our numerical evidence of stability. A tentative interpretation of our numerical results is given in Chapter 5. 14 Chapter 1. Introduction 1.5 Thesis Overview In this thesis, we mainly consider a coupled surface and grain boundary mo- tion which models the grain growth in materials processing and synthesis. Mathematical formulations are proposed to describe this coupled motion and other related problems. The emphasis throughout this thesis is the devel- opment of an efficient numerical approach for solving a class of three-phase boundary motion problems. We also analyze the well-posedness of the prob- lems and the stability of travelling waves, in a simple, linear setting. As discussed above, our results on this subject are not mathematically com- plete. The remainder of this thesis is organized as follows. In Chapter 2, three formulations that describe the coupled surface and grain boundary motion are developed, i.e. cartesian formulation, parabolic formulation and a formu- lation of partial differential algebraic equations (PDAE). The PDAE formu- lation preserves arc length and outperforms the other methods numerically. Travelling wave solutions are also presented for both the coupled motion and the simpler grain boundary motion. In Chapter 3, the well-posedness of the parabolic formulation and the PDAE formulation are analyzed in a simple, linear setting. Both formulations are shown to be well-posed with appropriate conditions. The rationality of the artificial conditions appearing in the parabolic formulation is also analyzed. Numerical simulations are presented for the quarter loop problem and other related problems in Chapter 4. The performance for each formulation is discussed and compared with others. Equations and junction conditions are approximated by finite difference methods on a staggered grid. Numer- ical convergence to exact travelling wave solutions is shown. The results of Chapter 2-4 are summarized in our paper [47]. The stability of the travelling waves for both the quarter loop problem and the simple grain boundary motion problem is discussed in Chapter 5. The analysis is based on a linearized system whose spectrum is approximated numerically on a finite computational domain. Numerical results give evi- dence that the travelling waves are convectively stable for both curves in 2D surfaces in 3D. A class of high order three-phase boundary motion problems are discussed in Chapter 6. The PDAE formulation is extended naturally to describe this class of problems. The well-posedness of the PDAE formulations is analyzed 15 Chapter 1. Introduction and numerical simulations are performed. The results of Chapter 5 and 6 are being prepared for journal publications. We are also trying to involve other collaborators who have expertise in the rigorous analysis of the stability of travelling waves. A summary of the thesis is given in Chapter 7. 16 Chapter 2 Mathematical Formulation As introduced in previous chapter, we mainly consider the coupled surface and grain boundary motion (see Figure 1.5) in this thesis. The simpler prob- lem that involves only curvature motion (see Figure 1.6) is also considered. The first model is a mixed fourth and second order problem. This problem has been studied in, for example, [29, 30, 31, 32, 40], both numerically and analytically. The second problem is a second order problem. Mathematical formulations and numerical simulations have been studied in [10, 24] and the well-posedness of this problem is discussed in [9]. In this chapter, we consider different formulations to describe both prob- lems: the cartesian formulation, the parabolic formulation and a new formu- lation containing both partial differential equations and algebraic equations. The problem of grain boundary motion involves only curvature motion. Thus it can be treated easily using the results for the coupled problem. There- fore we will focus on the coupled problem in this chapter and throughout the thesis. Relevant conclusions and results for the simpler grain boundary mo- tion problem will be given when necessary. 2.1 Cartesian Formulation We first consider the cartesian formulation presented in [29]. The grain boundary and the free exterior surface are denoted by u = u(x, t) and y = y(x, t) respectively (see Figure 1.5). Recall that the motion of the grain boundary is modeled by mean curvature motion, i.e., the normal velocity is proportional to its curvature. This motion can be represented in cartesian coordinates as (see Figure 2.1) ut = Aκ · √ 1 + u2x = A uxx 1 + u2x (2.1.1) where the curvature is given by κ = uxx (1 + u2x) 3/2 . (2.1.2) 17 Chapter 2. Mathematical Formulation θ tan θ = ux Vn · dt du Figure 2.1: A normal direction motion in cartesian coordinates with normal velocity Vn. The evolution of the free surface is modeled by surface diffusion. To find an expression for the normal velocity, i.e. the second derivative of curvature with respect to arc length s, we define the arc length in terms of x, s(x) = ∫ x x0 √ 1 + y2ξdξ from which we have sx = √ 1 + y2x. The second derivative of curvature is derived using the chain rule: κs = κx · dx ds = κx · 1√ 1 + y2x , κss = [ κx · 1√ 1 + y2x ] x · 1√ 1 + y2x . The evolution of the surface y(x) is then described by yt = −Bκss · √ 1 + y2x = −B [ 1 (1 + y2x) 1/2 [ yxx (1 + y2x) 3/2 ] x ] x . (2.1.3) The full system that describes the coupled motion is then written as ut = A uxx 1 + u2x , t > 0, x > j(t), yt = −B [ 1 (1 + y2x) 1/2 [ yxx (1 + y2x) 3/2 ] x ] x , t > 0, x ∈ (−∞, j(t)−) ∪ (j(t)+,∞), (2.1.4) 18 Chapter 2. Mathematical Formulation Free Surface Grain Boundary θ α + α − Figure 2.2: Sketch of the grain boundary groove with angle description. where j(t) denotes the x value of the junction where the three curves meet. Since the free surface on the right side of the groove root is different from that on the left, variables are separated into two branches: right (+) and left (–), where necessary. The conditions at the triple junction are given as follows. y+ = y− = u, α+ + α− = 2arcsin m 2 , θ = 1 2 [α+ − α−], (2.1.5) κ+ = κ−, κ+s = κ − s , where (see Figure 2.2) α+ = arctan(y+x (j(t), t)) α− = − arctan(y−x (j(t), t)) θ = pi 2 + arctan(ux(j(t), t)) m = γg γs . The constants γg and γs represent the surface energy of the grain boundary and the exterior free surface respectively. For the free surface, the surface 19 Chapter 2. Mathematical Formulation energy is assumed to be the same on both sides of the groove. Thus the parameter m is a constant measuring the relative surface energy. Of the five equations listed in (2.1.5), the first guarantees the triple junc- tion does not pull apart. The balance of surface tensions (Young’s law) gives the next two angle conditions. The fourth condition reflects the continuity of surface chemical potentials and the last condition demonstrates the balance of mass flux. Note that the two angle conditions can also be expressed as α+ − θ + pi 2 = arcsin m 2 + pi 2 , α− + θ + pi 2 = arcsin m 2 + pi 2 . (2.1.6) It is easy to see that the angles at the left hand side of the equations above represent the angles between the grain boundary and the two branches of the free surface, respectively. Since m is constant, the two angles are equal. We will find later that this form is more convenient when we consider other types of formulations. The far field conditions can be written as y(+∞) = y(−∞) = 0, u(+∞) = −1. (2.1.7) We will show later that the proportionality constants A and B in system (2.1.4) can both be taken as one by rescaling time and space. Even though, these scalings leave no degree of freedom to normalize the distance between the free surface and the grain boundary, we still take the distance as one. This fixes the speed of the travelling wave used to validate the numerical methods in Chapter 4 and used in the stability analysis in Chapter 5. However, similar results are obtained when different distances are used. The derivation of this cartesian formulation is straightforward. However, this formulation can not handle the case when the interface is not single- valued. As mentioned in [31, 32], it is physically reasonable that the exterior free surface may be non-single-valued as shown in Figure 2.3. Furthermore, even if the initial surface is single-valued, that property might be lost over the course of evolution. Thus we shall seek other more general formulations. 20 Chapter 2. Mathematical Formulation Free Surface Grain Boundary Figure 2.3: An example when the right branch of the free surface is not single-valued. 2.2 Parabolic Formulation In this section we derive a parabolic formulation to model the coupled sur- face and grain boundary motion. For the sake of possible non-single-valued curves, we use a parametric representation for all interfaces. Here and throughout this thesis we use X(·) = (u(·), v(·)) to represent a parameterized curve with u(·) and v(·) being the coordinates. The arc length parameter is denoted by s. We also use the parameter σ which is in the interval [0,∞) with σ = 0 at the junction. The derivative of X is the derivative of each component with respect to s or σ and it is denoted by Xs or Xσ. Even though the final system is meant to be given in σ, parameter s will be used in many places to keep the expressions simple. The vectors ~τ and ~n represent the unit tangential direction and unit normal direction respectively. 2.2.1 Motion by Mean Curvature We first consider the motion by mean curvature for the grain boundary. With the notation introduced above one has the following two fundamen- tal equalities: Xs = ~τ , Xss = κ~n. As mentioned before, the subscript s represents the derivative of X with 21 Chapter 2. Mathematical Formulation respect to arc length s. The first equality is obvious since ~τ = ( du√ du2 + dv2 , dv√ du2 + dv2 ) = ( du ds , dv ds ) = Xs. The second equality comes from the fact that Xss ·Xs = 1 2 (Xs ·Xs)s = 1 2 (|Xs|2)s = 0, and |Xss| = |dτ ds | = κ. To derive a relation between Xs and Xσ, we define the arc length function in terms of σ s(σ) = ∫ σ σ0 √ u2σ + v 2 σdσ (2.2.1) which gives the length of the curve from point X(σ0) to X(σ). Direct com- putation shows that Xσ = Xs ds(σ) dσ = Xs|Xσ| (2.2.2) where |Xσ| is defined as |Xσ| = √ u2σ + v 2 σ. Differentiating equation (2.2.2) with respect to σ leads to Xσσ = Xss|Xσ|2 +Xs|Xσ|σ. (2.2.3) Dividing through the equation above by |Xσ|2 gives Xσσ |Xσ|2 = Xss +Xs |Xσ|σ |Xσ|2 . Now we see that the term Xσσ|Xσ|2 is expressed as two parts, the normal com- ponent and the tangential component. One can compute the magnitude of the normal component by Xσσ |Xσ|2 · ~n = Xss · Xss κ + |Xσ|s |Xσ| Xs · Xss κ = κ. 22 Chapter 2. Mathematical Formulation Thus, if we set up a formulation, Xt = Xσσ |Xσ|2 , then it describes a motion with normal velocity being equal to its curva- ture. This gives us an option to describe the mean curvature motion with proportionality constant A, i.e. Xt = A Xσσ |Xσ|2 . (2.2.4) Equation (2.2.4) is fully parabolic (its linearization is parabolic in both components). This equation has the property that it smooths out the distor- tions in the grid when discretized and solved numerically. There are other equations that can also describe mean curvature motion, for example, Xt = κ~n. (2.2.5) But this system is not fully parabolic. A linearization shows that it is parabolic in the normal component and hyperbolic in the tangential com- ponent. A discretization of such a system will not have good numerical properties as those of a fully parabolic system due to the lack of regularity in the parametrization as shown in [10]. 2.2.2 Motion by Surface Diffusion In this section, we derive a parabolic equation for the motion by surface diffusion. The proportionality constant B is temporally taken as one. By analogy with the parabolic equation for the mean curvature motion described above, we look for an equation in the form of Xt = −Xσσσσ|Xσ|4 + L(Xσσσ, Xσσ, Xσ) (2.2.6) where L(Xσσσ, Xσσ, Xσ) denotes the combination of all lower order terms. The exact expression for L will be determined such that the normal compo- nent of the right hand side of the equation above is equal to the expected normal velocity −κss, i.e.(− Xσσσσ|Xσ|4 + L(Xσσσ, Xσσ, Xσ)) · ~n = −κss (2.2.7) 23 Chapter 2. Mathematical Formulation The rest of this section is focused on determining an expression for L, which will not be unique since it has an arbitrary tangential component. Note first the following expression (−Xssss − κ2Xss) · ~n = −κss. (2.2.8) To keep the argument compact, the proof of this expression is left to the end of this section. Comparing the left hand side of equation (2.2.7) and (2.2.8) one may set −Xσσσσ|Xσ|4 + L(Xσσσ, Xσσ, Xσ) = −Xssss − κ 2Xss which yields L(Xσσσ, Xσσ, Xσ) = Xσσσσ |Xσ|4 −Xssss − κ 2Xss. (2.2.9) This expression gives an option for L if the two fourth order derivatives cancel. To show this point, we shall first express the derivative Xssss in terms of σ. We start with equation (2.2.2) and differentiate it three times with respect to σ to obtain the following relations, Xσσ = Xss|Xσ|2 +Xs|Xσ|σ, Xσσσ = Xsss|Xσ|3 + 3Xss|Xσ||Xσ|σ +Xs|Xσ|σσ, (2.2.10) Xσσσσ = Xssss|Xσ|4 + 6Xsss|Xσ|2|Xσ|σ + 4Xss|Xσ||Xσ|σσ + 3Xss|Xσ|2σ +Xs|Xσ|σσσ. Dividing through the last equation by |Xσ|4 and reordering the equation, one obtains Xssss = Xσσσσ |Xσ|4 −6 |Xσ|σ |Xσ|2Xsss−4 |Xσ|σσ |Xσ|3 Xss−3 |Xσ|2σ |Xσ|4Xss− |Xσ|σσσ |Xσ|4 Xs (2.2.11) Similarly, one has Xss = Xσσ |Xσ|2 −Xs |Xσ|σ |Xσ|2 , Xsss = Xσσσ |Xσ|3 − 3Xss |Xσ|σ |Xσ|2 −Xs |Xσ|σσ |Xσ|3 . 24 Chapter 2. Mathematical Formulation Substitute these expressions for Xssss, Xsss, Xss into equation (2.2.9) and rewrite all arc length expressions in terms of σ. We ignore the Xs terms since they have no contribution to the normal direction. An expression for L is then given by L(Xσσσ, Xσσ, Xσ) = 6 |Xσ|σXσσσ |Xσ|5 − 15 |Xσ|2σXσσ |Xσ|6 + 4 |Xσ|σσXσσ |Xσ|5 − κ 2 Xσσ |Xσ|2 . (2.2.12) We now substitute equation (2.2.12) back into (2.2.6) and reorganize to ob- tain Xt = −Xσσσσ|Xσ|4 + 6 |Xσ|σXσσσ |Xσ|5 − (15 |Xσ|2σ |Xσ|4 − 4 |Xσ|σσ |Xσ|3 + κ 2) Xσσ |Xσ|2 . (2.2.13) The equation above is parabolic after linearization and it describes motion by surface diffusion since the normal component of the right hand side ex- pression is equal to −κss. Again, we shall point out that the choice of L is not unique because one can choose arbitrary tangential component. A simi- lar expression for L is given by Garcke et al. in [25] with different tangential velocity. Proof of Expression (2.2.8): Recall that Xss = κ~n, (2.2.14) which gives Xss ·Xss = κ2. (2.2.15) Differentiating twice with respect to s leads to Xssss ·Xss +Xsss ·Xsss = κκss + κ2s. Therefore, κss = Xssss ·Xss +Xsss ·Xsss − κ2s κ . (2.2.16) On the other hand, one may differentiate equation (2.2.14) with respect to s to get Xsss = κs~n+ κ~ns = κs~n− κ2~t, and thus Xsss ·Xsss = κ2s + κ4. (2.2.17) 25 Chapter 2. Mathematical Formulation Now we substitute equation (2.2.17) into (2.2.16) which yields κss = Xssss · ~n+ κ3. (2.2.18) Noting the fact that κ = Xss · ~n, the expression above can be rewritten as κss = (Xssss + κ 2Xss) · ~n. This completes the proof. 2.2.3 The Parabolic System X 1 (grain boundary) θ θ X 2 (free surface) X 3 (free surface) Figure 2.4: Sketch of the grain boundary groove showing the curve number- ing. We label the grain boundary, the left branch and the right branch of the free surface as X1, X2 and X3 respectively (see Figure 2.4). To avoid the confusion with the power index, all superscripts that followX directly specify the index of the curves unless otherwise noted. The parabolic system that describes the coupled surface and grain boundary motion is given by X1t = A X1σσ |X1σ|2 , X2t = −B( X2σσσσ |X2σ|4 − 6 |X 2 σ|σX2σσσ |X2σ|5 + (15 |X2σ|2σ |X2σ|4 − 4 |X 2 σ|σσ |X2σ|3 + |X2σσ|2 |X2σ|4 ) X2σσ |X2σ|2 ), X3t = −B( X3σσσσ |X3σ|4 − 6 |X 3 σ|σX3σσσ |X3σ|5 + (15 |X3σ|2σ |X3σ|4 − 4 |X 3 σ|σσ |X3σ|3 + |X3σσ|2 |X3σ|4 ) X3σσ |X3σ|2 ), (2.2.19) 26 Chapter 2. Mathematical Formulation for t > 0 and σ ∈ [0,∞). 2.2.4 Boundary Conditions for the Parabolic System The conditions at the triple junction (σ = 0) can be transformed from those discussed for the cartesian formulation. Referring to Figure 2.4, the conditions are given as follows: first, the common coordinates at the junction gives X1(0) = X2(0) = X3(0). (2.2.20) The two angle conditions that come from Young’s law now read as X1σ |X1σ| · X 2 σ |X2σ| = cos θ = cos( pi 2 + arcsin m 2 ), (2.2.21) X1σ |X1σ| · X 3 σ |X3σ| = cos θ = cos( pi 2 + arcsin m 2 ). (2.2.22) Here we use the equivalent form of the angle conditions described in equation (2.1.6). This guarantees that the grain boundary makes the same angle with both branches of the free surface. The continuity of the surface chemical potentials gives κ2 = −κ3 (2.2.23) where the curvature can be calculated by κ = Xσσ |Xσ|3 ·X ⊥ σ . (2.2.24) Here the superscripts of κ are curve indices. The superscript “⊥” stands for the perpendicular of that vector. Last, the balance of mass flux implies κ2s = κ 3 s (2.2.25) where κs is given by κs = Xσσσ ·X⊥σ |Xσ|4 − 3 |Xσ|σ(Xσσ ·X⊥σ ) |Xσ|5 . (2.2.26) 27 Chapter 2. Mathematical Formulation Again, the superscripts are curve indices. The expression for κs is obtained by taking the derivative of κ directly. The reader may have noticed that there is a minus sign in equation (2.2.23). The continuity of the surface chemical potentials requires that the two free surfaces have the same convexity. However, opposite directions for σ are used in the parametrization of the free surfaces. The positive direction goes to the left for the left branch while for the right branch it goes to the right. Thus the normal direction changes sign between the branches and the minus sign in (2.2.23) is needed. Any further derivatives will change the sign for the left branch. Therefore no modification is needed for condition (2.2.25) since another derivative is taken. Similar conditions are derived in [25], where the authors considered a triple junction problem with all three interfaces moving by surface diffusion. The conditions as σ →∞ are simply expressed as v1 = −1, v2 = v3 = 0 where vi represents the y coordinate of curve X i. 2.2.5 Artificial Tangential Conditions The whole system contains two second order equations and four fourth order equations and therefore needs ten conditions at the triple junction. Recall that there are only eight physical junction conditions as addressed above. Thus, two more conditions are needed. There are several options to impose the extra conditions. We will show later that our conditions only change the parametrization of the curves and do not change the shape of the curves. We refer to the additional conditions as artificial tangential conditions. As one of the options, the following two conditions are applied into the system: X iσσ ·X iσ = 0 for i = 2, 3. (2.2.27) The meaning of these conditions will be shown in next section. 2.2.6 Conclusions This parabolic formulation leads to tangential adjustment of the grid points and it is possible to avoid singularities in the parametrization. However, we 28 Chapter 2. Mathematical Formulation will show later that this approach leads to a scheme with parameters in the tangential parabolicity for the fourth order equations that need to be chosen carefully to keep a monotone parametrization in a numerical approximation. 2.3 A Formulation of Partial Differential Algebraic Equations We already pointed out that a successful numerical simulation requires an efficient approach to adjust the distribution of grid points. Allowing only normal direction motion is known to be a bad idea. The parabolic formulation described in last section can partly improve the situation but still cannot completely solve the problem. In this section, we propose a new formulation to model the coupled sur- face and grain boundary motion. This formulation contains not only partial differential equations (PDE) but also algebraic equations (AE). Therefore it is called a PDAE formulation. It avoids possible loss of tangential mono- tonicity in the parametrization due to the fourth order PDE and maintains uniform distribution of grid points in arc length. We use the same notation as introduced in previous section. The reader will find that the discussion in this section is much shorter and simpler than that for the parabolic formulation. Although the discussion is shortened in part by the use of previously defined notation, the new approach is also particularly concise. 2.3.1 Motion by Mean Curvature To model the motion by mean curvature, the basic evolution law (1.1.1) should be satisfied, which can be expressed as Xt · ~n− Aκ = 0. (2.3.1) This equation describes a motion with normal velocity being proportional to its curvature. However, only the normal velocity is specified and it leaves the tangential velocity to be arbitrary. We need to choose another equation to specify the tangential velocity. Comparing with the results of other techniques for adjusting the tan- gential velocity, much better numerical performance is achieved using the 29 Chapter 2. Mathematical Formulation following approach. Equation (2.3.1) is augmented by a condition that di- rectly imposes uniform grid spacing: |Xσ|σ = 0. (2.3.2) Noting the fact that |Xσ|σ = ( √ Xσ ·Xσ)σ = Xσ ·Xσσ|Xσ| , equation (2.3.2) can be replaced by Xσ ·Xσσ = 0. (2.3.3) The κ in equation (2.3.1) is computed by κ = Xσσ |Xσ|3 ·X ⊥ σ . (2.3.4) The idea of maintaining equally spaced grid points on a curve has been considered in several other works ([28, 39]). Our problem is simpler than some described in [28] in that it is not necessary to keep track of material points on the curves. The idea of applying (2.3.3) to enforce an equi-spaced parametrization appears to be a new idea. In previous work, particular choices of tangential velocities were used that caused the parametrization to move towards an equi-spaced parametrization. However, this approach has been applied only to the mean curvature motion and it requires extra work for the fourth order problem. Our formulation does not depend on the type of the motion thus it can be applied directly to the motion by surface diffusion. Moreover, it keeps a strict equi-distributed grid which leads to a much easier expression for the motion by surface diffusion, i.e. κss, as shown below. 2.3.2 Motion by Surface Diffusion The PDAEs for the motion by surface diffusion are given analogously as Xt · ~n+Bκss = 0 (2.3.5) Xσ ·Xσσ = 0. 30 Chapter 2. Mathematical Formulation Using the fact that σ is a scaled arc length parameter, this leads to the expression κss = κσσ |Xσ|2 . (2.3.6) where κ is computed by (2.3.4). 2.3.3 The PDAE System We use the same notation for the three curves as shown in Figure 2.4. The full PDAE system for the coupled surface and grain boundary motion is given by X1t · ~n− Aκ = 0 X it · ~n+Bκss = 0 i = 2, 3 X iσ ·X iσσ = 0 i = 1, 2, 3 . (2.3.7) This is an index-1 PDAE system. The index of a DAE system is the minimum number of differentiations of the system which would be required to obtain an ODE system for all differential variables and algebraic variables. To show the index of this formulation, we consider the curvature motion: Xt · ~n− Aκ = 0 Xσ ·Xσσ = 0. (2.3.8) We treat the second equation as a constraint and differentiate with respect to time t: (Xσ ·Xσσ)t = 0. Further investigation gives (Xσ ·Xσσ)t = (Xt ·Xσ)σσ − (Xt ·Xσσ)σ = 0. (2.3.9) The first equation in system (2.3.8) gives Xt · ~n = Xt · Xσσ κ|Xσ|2 = Aκ. (2.3.10) Using equation (2.3.9) and (2.3.10), one obtains (Xt ·Xσ)σσ = (Aκ2|Xσ|2)σ. 31 Chapter 2. Mathematical Formulation Noticing Xσ = ~t, we obtain a system for the normal component and the tangential component of Xt: Xt · ~n = Aκ (Xt · ~t)σσ = (Aκ2|Xσ|2)σ. (2.3.11) A similar derivation can be applied to the original system. Therefore, the PDAE system is of index-1 since we only need to differentiate once. Usually an index-1 DAE can be handled without any numerical difficulties ([3]), and that is our experience with the discretization of this system. 2.3.4 Boundary Conditions for the PDAE System Since we use the same parametric representations as those for the parabolic system, the junction conditions are exactly the same as before. We repeat all conditions below for completeness. The conditions at the triple junction are X1(0) = X2(0) = X3(0), X1σ |X1σ| · X 2 σ |X2σ| = cos( pi 2 + arcsin m 2 ), X1σ |X1σ| · X 3 σ |X3σ| = cos( pi 2 + arcsin m 2 ), κ2 = −κ3, κ2s = κ 3 s. The conditions as σ →∞ are given by v1 = −1, v2 = v3 = 0. Though these conditions are the same as those for the parabolic formu- lation, the ways to discretize them are quite different when they are solved numerically. Because of the lower order of the constraint algebraic equations in the PDAE system, it is possible to avoid the use of artificial tangential conditions at the triple junction and the domain boundaries. More details are described in the numerical implementation in Chapter 4. 32 Chapter 2. Mathematical Formulation 2.4 Exact Travelling Wave Solutions In this section, we give the exact travelling wave solutions for both the equal order grain boundary motion and the mixed order coupled motion. The stability of these solutions are investigated in Chapter 5 and they are also used to compare with the numerical results in Chapter 4. 2.4.1 Travelling Wave Solution for the Grain Boundary Motion We first derive the travelling wave solution for the simple grain boundary motion problem that involves only curvature motion (see Figure 1.6). Let ui, i = 1, 2, 3 represent the left curve, the top right and the bottom right curve in Figure 1.6 respectively. The system is given by uit = uixx 1 + (uix) 2 , i = 1, 2, 3 (2.4.1) with the following triple junction conditions, u1 = u2 = u3, arctan(u1x)− arctan(u2x) = −φ, arctan(u1x)− arctan(u3x) = φ, (2.4.2) u1(−∞, t) = 0, u2(+∞, t) = 1, u3(+∞, t) = −1, where φ = pi − θ and θ is the angle between u1, u2 and u1, u3 as shown in Figure 1.6. We first consider the travelling wave solution for the top right branch (u2). We assume there is a travelling wave solution in the form of u(ξ) = u(x−ct), where u(ξ) satisfies 0 = uξξ 1 + u2ξ + cuξ. (2.4.3) Integrating the equation above and using the fact that u(ξ)→ 1 and uξ(ξ)→ 0 as ξ →∞, we obtain, arctanuξ + cu− c = 0. (2.4.4) 33 Chapter 2. Mathematical Formulation Evaluating at ξ = 0 and noticing u(0) = 0, uξ(0) = tanφ, it gives c = φ. From equation (2.4.4) and using c = φ, one has uξ = − tan(φ(u− 1)), and it is integrated from 0 to ξ to obtain u = 1− 1 φ arcsin(sinφe−φξ). It is obvious that u1 = 0 is the solution for the left branch. And since u3 is symmetric to u2 about the x-axis, we have the following travelling wave solution for the system, u1 = 0, u2 = 1− 1 φ arcsin(sinφe−φξ), (2.4.5) u3 = −1 + 1 φ arcsin(sinφe−φξ), where ξ = x− ct and c = φ is the travelling wave speed. 2.4.2 Travelling Wave Solution for the Coupled Motion A travelling wave solution for the coupled surface and grain boundary motion was derived by Kanel et al. in [30]. Without loss of generality, we take both proportionality constants A and B as one and introduce briefly the results found in by Kanel et al. in [30]. To be quite clear, this travelling wave is not new work. Its derivation is included here for completeness. Since we are looking for a travelling wave solution, i.e. a solution in the form of f(ξ) = f(x− ct), the cartesian formulation (2.1.4) can be rewritten in terms of ξ, 0 = uξξ 1 + u2ξ + cuξ, t > 0, ξ > 0, 0 = − [ 1 (1 + y2ξ ) 1/2 [ yξξ (1 + y2ξ ) 3/2 ] ξ ] ξ + cyξ, t > 0, ξ ∈ (−∞, 0−) ∪ (0+,∞). (2.4.6) 34 Chapter 2. Mathematical Formulation We define Φ = arctan uξ and Ψ = arctan yξ, and require that −pi 2 ≤ Φ ≤ 0 and − pi 2 ≤ Ψ ≤ pi 2 . We introduce the arc length parametrization s1 = ∫ ξ 0 √ 1 + u2ξdξ, ξ ∈ (0,∞), and s2 = ∫ ξ 0 √ 1 + y2ξdξ, ξ ∈ (−∞,∞). Note that s1 = s2 = 0 is at the triple junction and when s2 < 0 represents the left branch of the free surface. We shall now rewrite formulation (2.4.6) into an angle formulation: Φs1 = −c sinΦ, s1 ∈ (0,∞), Ψs2s2s2 = c sinΨ, s2 ∈ (−∞,∞). (2.4.7) In the original problem, m is the only physical parameter that determines the profile of the solution. But now let us treat ψ = Ψ(0−) as a parameter and m as a function of ψ. The arc length s2 is rescaled to s = 3 √ cs2 where c is the travelling wave speed yet to be determined. Then the solution for the left branch of the free surface is given in terms of s by Ψ(s) = ψes + 1 3 ∫ s −∞ g(Ψ(s̃))G(s− s̃)ds̃ −1 3 es ∫ 0 −∞ g(Ψ(s̃))G(−s̃)ds̃, −∞ < s < 0, (2.4.8) where g(u) := sin u− u, 35 Chapter 2. Mathematical Formulation and G(x) := ex − 2e− 12x sin( √ 3 2 x+ pi 6 ). This is an integral equation and it is solved by iteration with Ψ0(s) = 0, Ψ1(s) = ψe s, and Ψn+1(s) = Ψ1(s) + 1 3 ∫ s −∞ g(Ψn(s̃))G(s− s̃)ds̃ −1 3 es ∫ 0 −∞ g(Ψn(s̃))G(−s̃)ds̃, n = 1, 2, · · · . Once we have Ψ(s) for s < 0, the solution for the right branch is then given by Ψ(s) = Ψ1(s)− 1 3 ∫ ∞ s g(Ψ(s̃))G(s− s̃)ds̃ −2 3 e− 1 2 s ∫ ∞ 0 g(Ψ(s̃))H0(s, s̃)ds̃, 0 < s <∞, where Ψ1(s) = e − 1 2 s [− 2ψ cos(√3 2 s) + 2 3 ∫ 0 −∞ g(Ψ(s̃))H1(s, s̃)ds̃ ] , and H0(s, s̃) := e −s̃ cos √ 3 2 s− e 12 s̃ sin( √ 3 2 (s− s̃) + pi 6 ), H1(s, s̃) := e 1 2 s̃ [ 2 cos( √ 3 2 s) sin( √ 3 2 s̃− pi 6 )− sin( √ 3 2 (s− s̃) + pi 6 ) ] . This integral equation is again solved by iteration. We are still following the derivation in [30] and we shall now derive the solution for the grain boundary. We first define Ψm := Ψ(0 +)−Ψ(0−) (2.4.9) 36 Chapter 2. Mathematical Formulation where Ψm is given by Ψm = −3ψ + 2 ∫ 0 −∞ g(Ψ(s̃))e 1 2 s̃ sin( √ 3 2 s̃− pi 6 )ds̃ − ∫ ∞ 0 g(Ψ(s̃))e−s̃ds̃. (2.4.10) Recall that ψ = Ψ(0−) is the parameter whose value is given, therefore we can compute Ψ(0+) using equation (2.4.9) and (2.4.10). The solution for the grain boundary is then given by Φ(s1) = 2 arctan(tan[ 1 2 Φ(0+)]e−cs1), (2.4.11) where Φ(0+) = pi 2 + 1 2 [Ψ(0+) + Ψ(0−)]. (2.4.12) Now there is only one unknown, the travelling wave speed c, yet to be deter- mined. We define f(λ) := λ3 + λ2 ∂2 ∂s2 Ψ(0+)− [pi 2 − 1 2 [Ψ(0+) + Ψ(0−)] ] = 0, (2.4.13) then c = λ31, where λ1 is the unique positive root of equation (2.4.13). Now we can scale s back to s2 for the free surface Ψ and express the solu- tion in cartesian coordinate system, i.e. ( x(s1, t), u(s1, t) ) , ( x(s2, t), y(s2, t) ) as x(s1, t) = x(0, t) + ∫ s1 0 cosΦ(s, t)ds, s1 ∈ (0,∞), u(s1, t) = −1− ∫ ∞ s1 sinΦ(s, t)ds, s1 ∈ (0,∞), x(s2, t) = x(0, t) + ∫ s2 0 cosΨ(s, t)ds, s2 ∈ (−∞,∞), (2.4.14) y(s2, t) = − ∫ ∞ s2 sinΨ(s, t)ds, s2 ∈ (0,∞), y(s2, t) = ∫ s2 −∞ sinΨ(s, t)ds, s2 ∈ (−∞, 0). 37 Chapter 2. Mathematical Formulation All previous derivation of the travelling wave solution for the coupled problem was introduced by Kanel et al. in [30]. The existence and uniqueness of the travelling wave solution above was proved for 0 < m . .92 in [30] and was extended to m ∈ [0, 2) in [31]. 2.5 Scaling the Motion Laws Before continuing to the next chapter, we shall now do the promised scaling about the proportionality constants A,B for the motion laws. After this deduction, A,B will always be taken as one in following chapters. Without loss generality, we rescale only the PDAE formulation. The original curvature motion and surface diffusion are given by Xt · ~n = Aκ Yt · ~n = −Bκss } . (2.5.1) We show here that both A and B can be made unity by rescaling time and space. We scale time by t̃ = Tt. The spatial variable X,Y are scaled by X̃ = RX, Ỹ = RY. We can easily verify that s̃ = Rs where s̃ is the arc length in the new system. It also follows that X̃s̃ = Xs, Ỹs̃ = Ys, κ̃ = 1 R κ, κ̃s̃s̃ = 1 R3 κss. Now system (2.5.1) becomes X̃t̃ · ~n = AR2T κ̃ Ỹt̃ · ~n = −BR4T κ̃s̃s̃ } . (2.5.2) 38 Chapter 2. Mathematical Formulation Here we use the same notation ~n since the normal direction does not change. If we choose R = √ A B , T = A2 B then equation (2.5.2) becomes X̃t̃ · n = κ̃ Ỹt̃ · n = −κ̃s̃s̃. This completes the normalization of coefficients A and B. Note that in [30], the authors retain the ratio A/B in their equations. They use the remaining degree of freedom to normalize the distance between the grain boundary and the upper surface far from the junction in their travelling wave solutions. 39 Chapter 3 Linear Well-posedness In this chapter, we analyze the well-posedness of the parabolic system and the PDAE system for the coupled surface and grain boundary motion we proposed in Chapter 2. In [9], a problem of three curves moving with curvature motion and meet- ing at a junction is analyzed. In [25], a problem of three curves undergoing surface diffusion and meeting at a junction is analyzed. Both problems con- sidered in [9] and [25] are equal order. There is a gap in the theory for the quarter loop problem, which is of mixed order, for which the underlying theoretical estimates for the equal order problems in [55] are not available. We gain some insight by linearizing around fixed straight line solutions and obtaining a system that has the same highest order behaviour as the original problem near the junction. The resulting linear system is then analyzed. The analysis here gives results that match those that more complicated nonlinear analysis gives in [9, 25]. Therefore, we believe the results of the analysis should apply to the full nonlinear problem, although some technical details are left as open problems by our analysis. 3.1 Preliminaries To prepare for the well-posedness analysis, we first use a simple example to illustrate the basic idea. We consider the well-posedness of the following initial-boundary value problem on the half line: ut = uxx − u+ f(x, t), x > 0, u(x, 0) = u0(x), ux(0, t)− u(0, t) = g(t). (3.1.1) Extending functions f and u0 smoothly to the whole real line leads to the following cauchy problem: vt = vxx − v + f(x, t), v(x, 0) = u0(x). (3.1.2) 40 Chapter 3. Linear Well-posedness This problem can be solved by Fourier transform and the solution is given by v = ∫ ∞ −∞ eiαxṽdα where ṽ = e−(α 2+1)tũ0 + ∫ t 0 e−(α 2+1)(t−τ)f̃(α, τ)dτ and ũ0, f̃ are the Fourier transform of u0 and f respectively. Now we see that the solution v is smooth as long as u0 and f are smooth enough. If we introduce function w = u − v defined for x > 0, then we have a system for w given by wt = wxx − w, x > 0 w(x, 0) = 0, wx(0, t)− w(0, t) = g(t) + v(0, t)− vx(0, t). (3.1.3) To solve for w, we apply Laplace transform to system (3.1.3) which gives sŵ = ŵxx − ŵ ŵx − ŵ = ĝ(s)− v̂x(s) + v̂(s) and therefore ŵ has solution in the form of ŵ = A(s)e− √ s+1x. Here we consider Re(s) > 0 and the square root of s + 1 is chosen such that Re( √ s+ 1) > 0. Substituting the expression above into the boundary condition and evaluating at x = 0, we obtain (−√s+ 1− 1)A(s) = ĝ(s)− v̂x(s) + v̂(s), and therefore, A(s) = ĝ(s)− v̂x(s) + v̂(s) (−√s+ 1− 1) , if √ s+ 1 + 1 6= 0. Thus the existence of solution ŵ, and consequently, w and u is reduced to determining the value of the coefficient of unknown A. If the coefficient is not 41 Chapter 3. Linear Well-posedness zero for any s satisfying Re(s) > 0, the problem is well-posed with at most algebraic growth. Similarly, if the system contains more than one variable, one can obtain a coefficient matrix for the unknowns at the boundary and it must be invertible for a well-posed problem. A similar argument can be applied to our mixed order problem. Thus we only need to investigate the property at the triple junction for the well- posedness of the problem. Linearizing around the initial data and freezing the coefficient at the junction, one can obtain a constant-coefficient linear system. The linear system is then solved by Laplace transform and the solutions are substituted into the junction conditions leading to a coefficient matrix for the unknowns. The well-posedness of the original system is obtained by investigating the determinant of the matrix. We remark that one may linearize around the tangential directions at the triple junction to obtain the constant-coefficient linear system. In following sections, we use this approach to investigate the well-posedness for both the parabolic formulation and the PDAE formulation. 3.2 Well-posedness of the Parabolic System We first consider the well-posedness of the parabolic system (2.2.19). 3.2.1 Well-posedness To linearize the system we consider a perturbation expansion around the tangential direction at the triple junction for each curve, i.e., X1 = d1σ + ²X̄ 1 +O(²2) X2 = d2σ + ²X̄ 2 +O(²2) X3 = d3σ + ²X̄ 3 +O(²2) where di = (di1, di2) are constant vectors standing for the unit tangential directions. Substituting the expressions above into (2.2.19), keeping only the leading order terms (O(²)) and dropping the bar superscript yield the following linear system: X1t = X 1 σσ X2t = −X2σσσσ (3.2.1) X3t = −X3σσσσ. 42 Chapter 3. Linear Well-posedness The linearization of the triple junction conditions including the artificial tangential conditions are given below: • Common point at σ = 0: X1 = X2 = X3 • Angle conditions: d1 ·X2σ + d2 ·X1σ − (d1 · d2)(d1 ·X1σ + d2 ·X2σ) = 0 d1 ·X3σ + d3 ·X1σ − (d1 · d3)(d1 ·X1σ + d3 ·X3σ) = 0 • Continuity of surface chemical potentials: X2σσ · d⊥2 = −X3σσ · d⊥3 • Balance of mass flux: X2σσσ · d⊥2 = X3σσσ · d⊥3 • Artificial tangential conditions: X2σσ · d2 = 0 X3σσ · d3 = 0. After taking the Laplace transform, the linearized system (3.2.1) is solved to have solutions in the form u1 = A11e −√sσ v1 = A12e −√sσ u2 = A21e λ1σ +B21e λ2σ v2 = A22e λ1σ +B22e λ2σ u3 = A31e λ1σ +B31e λ2σ v3 = A32e λ1σ +B32e λ2σ . (3.2.2) The variables ui and vi stand for their Laplace transforms and s here stands for the Laplace transform variable. The unknowns Aij and Bij will be deter- mined from boundary conditions and λ1 = (− √ 2 2 + √ 2 2 i) 4 √ s, λ2 = (− √ 2 2 − √ 2 2 i) 4 √ s. 43 Chapter 3. Linear Well-posedness Here we consider Re(s) > 0 and the roots of s are chosen such that | arg(√s)| < pi 2 , | arg( 4√s)| < pi 4 . (3.2.3) Without changing the well-posedness of the problem we specify one of the tangential directions, say d1 = (0,−1)T and thus d1 = 0 −1 d2 = − sin θ − cos θ d3 = sin θ − cos θ . Solution (3.2.2) is substituted into boundary conditions with the values of di given above. Simplification yields a linear system for the ten unknowns Aij and Bij. The transpose of the 10× 10 coefficient matrix M is given by 1 0 0 0 sin θ √ s − sin θ√s 0 0 0 0 0 0 1 0 0 0 0 0 0 0 −1 1 0 0 12 sin 2θλ1 0 cos θλ12 cos θλ13 − sin θλ12 0 −1 1 0 0 12 sin 2θλ2 0 cos θλ22 cos θλ23 − sin θλ22 0 0 0 −1 1 − sin2 θλ1 0 − sin θλ12 − sin θλ13 − cos θλ12 0 0 0 −1 1 − sin2 θλ2 0 − sin θλ22 − sin θλ23 − cos θλ22 0 0 −1 0 0 0 −12 sin 2θλ1 cos θλ12 − cos θλ13 0 sin θλ12 0 −1 0 0 0 −12 sin 2θλ2 cos θλ22 − cos θλ23 0 sin θλ22 0 0 0 −1 0 − sin2 θλ1 sin θλ12 − sin θλ13 0 − cos θλ12 0 0 0 −1 0 − sin2 θλ2 sin θλ22 − sin θλ23 0 − cos θλ22 . (3.2.4) As discussed in previous section, the problem is well-posed with at most algebraic growth in time provided that the determinant of matrix M is not zero for any s satisfying Re(s) > 0. One can compute the determinant of M directly which gives |M | = 64(sin3 θ)s3 − 16 √ 2(sin2 θ sin 2θ)s11/4. This determinant is zero only when s = 0, or sin θ = 0, or 4 √ s = √ 2 2 cos θ. (3.2.5) 44 Chapter 3. Linear Well-posedness The first equation violates the condition that Re(s) > 0. The second and the third equation are also impossible if pi/2 ≤ θ < pi. In summary, the problem is well-posed provided θ ∈ [pi/2, pi). We remark that, physically any angle cannot be equal or larger than pi which is corresponding to zero or negative surface energy. However, the problem is still well-posed if θ = pi/2. 3.2.2 Analysis of Artificial Tangential Conditions Artificial tangential conditions are required for both well-posedness analy- sis and numerical simulations of the parabolic formulation. There are many options for these conditions. One would wonder if all these choices lead to the same solution: maybe different parametrization but at least the same geometric shape. It is true that the solution does not depend on the artifi- cial conditions provided they are compatible. To prove this statement, it is sufficient to show that the position of the junction and the three tangential directions do not depend on the artificial tangential conditions. If we can show one of the three curves, say X1, is uniquely determined up to arbi- trary tangential parametrization, then the location of the triple junction is uniquely determined. Further, the angle conditions will guarantee that the tangential directions of the other two curves are also uniquely determined. In summary, we only need to show that X1 or more specifically, the coefficients of X1, i.e., A11, A12 do not depend on the artificial tangential conditions. Let C = [A11, A12, · · · , A32, B32] represent the vector of unknown coef- ficients in expression (3.2.2) and let P = [p1, p2, · · · , p9, p10] be a constant vector indicating the initial values of the system at the junction. The coeffi- cients Aij, Bij can be solved from C =M−1 · P (3.2.6) where M is the coefficient matrix (3.2.4). If we keep the same order of the junction conditions as given before, then only p9, p10 and the last two rows of M are related to artificial tangential conditions. Thus, to prove A11, A12 do not depend on the artificial tangential conditions we shall prove A11, A12 do not depend on the last two rows of matrix M and p9, p10. For convenience, we rewrite M into a block form M = M1(8× 2) M2(8× 8) M3(2× 2) M4(2× 8) (3.2.7) 45 Chapter 3. Linear Well-posedness The Gauss elimination for the (8 × 8) block M2 shows that the rank of M2 is six for any angle conditions. Therefore, it is possible to make the last two rows of M2 zeros by row reductions and meanwhile making the last two rows ofM1 into a full rank (2×2) matrix. For example, when θ = 2pi/3 the matrix after such an operation is given by s 11 4 1 0 −1 −1 0 0 0 0 0 0 i 2 0 0 −i i √ 3 2 − i √ 3 2 i 2 − i 2 − i √ 3 2 i √ 3 2 0 1 0 0 −1 −1 0 0 0 0 1 0 0 0 0 0 −1 −1 0 0 i √ 2 2 − √ 6 4 0 0 0 0 0 −i√2− (1+2i) √ 6 4 − (1−2i) √ 6 4 0 1 0 0 0 0 0 0 −1 −1 −4 √ 3 4 √ s+ √ 6 8 3 √ 2 4 0 0 0 0 0 0 0 0 0 3 √ 2 2 0 0 0 0 0 0 0 0 0 0 i √ 3 2 − i √ 3 2 − i 2 i 2 0 0 0 0 0 0 0 0 0 0 − i √ 3 2 i √ 3 2 − i 2 i 2 . The new matrix after row reductions will still be denoted byM . LetM−1 represent the inverse of M . One has M ×M−1 = M1(8× 2) M2(8× 8) M3(2× 2) M4(2× 8) × M̄1(2× 8) M̄2(2× 2) M̄3(8× 8) M̄4(8× 2) = I(8× 8) 0 0 I(2× 2) . (3.2.8) Here M̄i represent the four blocks of M −1 and I is identity matrix. Expand (3.2.8) directly to get M1 × M̄1 +M2 × M̄3 = I(8× 8) (3.2.9) M1 × M̄2 +M2 × M̄4 = 0(8× 2) (3.2.10) M3 × M̄1 +M4 × M̄3 = 0(2× 8) (3.2.11) M3 × M̄2 +M3 × M̄4 = I(2× 2) (3.2.12) 46 Chapter 3. Linear Well-posedness Note that equation (3.2.9)-(3.2.10) do not involve M3,M4 which means they do not depend on the artificial tangential conditions. If we can solve M̄1, M̄2 from equation (3.2.9)-(3.2.10) then M̄1, M̄2 are not affected by the artificial tangential conditions. In fact, just equation (3.2.9) is enough to determine M̄1 uniquely. Recall that the last two rows of M2 are all zeros from which one can obtain sixteen equations involving only the sixteen entries of M̄1. For the same reason we can solve for M̄2 by equation (3.2.10). Further, since the last two rows of M1 is a full rank (2× 2) matrix, M̄2 must be zero. In summary, we have showed that M̄1, M̄2 do not depend on the artificial tangential conditions, and M̄2 = 0. Now, to compute A11, A12 we have( A11 A12 ) = ( M̄1 M̄2 )× P. Clearly, A11, A12 depend only on the original physical conditions since the zero matrix M̄2 guarantees that p9, p10 have no impact on the results. This completes the proof that the shape of the three curves do not depend on the artificial tangential conditions. Garcke et al. [25] also pointed out that the artificial conditions do not influence the solutions, although the problem there is a bit different. In [25], the authors investigated a three-phase problem in which all three interfaces evolve by minus the surface Laplacian of mean curvature and meet at a triple junction. 3.3 Well-posedness of the PDAE System In this section, we analyze the well-posedness of the PDAE system. The linearization of the system is done in a similar way as before. A perturbation expansion around the tangential directions at the triple junction is considered for each curve, i.e., X1 = d1σ + ²X̄ 1 +O(²2) X2 = d2σ + ²X̄ 2 +O(²2) X3 = d3σ + ²X̄ 3 +O(²2). 47 Chapter 3. Linear Well-posedness The resulting linear system is given by X1t · d⊥1 = X1σσ · d⊥1 X it · d⊥i = −X iσσσσ · d⊥i , i = 2, 3 di ·X iσσ = 0, i = 1, 2, 3. (3.3.1) Although, this linear system is different from the one for the parabolic sys- tem, the linearized junction conditions are the same, except that no artificial tangential conditions are needed. They are given by X1 = X2 = X3, d1 ·X iσ + di ·X1σ − (d1 · di)(d1 ·X1σ + di ·X iσ) = 0, i = 2, 3, X2σσ · d⊥2 = −X3σσ · d⊥3 , X2σσσ · d⊥2 = X3σσσ · d⊥3 . We again specify one of the tangential directions, say d1 = (0,−1)T . Further, we assume θ ∈ (pi 2 , pi) (3.3.2) so that di,1, di,2 6= 0 for i = 2, 3. The linearized system (3.3.1) is solved by Laplace transform and the solutions (after Laplace transform) are given by u1 = A11e −√sσ v1 = B11 u2 = A21e λ1σ +B21e λ2σ + C21 v2 = −k2(A22eλ1σ +B22eλ2σ) + 1k2C21 u3 = A31e λ1σ +B31e λ2σ + C31 v3 = −k3(A32eλ1σ +B32eλ2σ) + 1k3C31 (3.3.3) where ki = di1/di2 is constant and λ1 = (− √ 2 2 + √ 2 2 i) 4 √ s, λ2 = (− √ 2 2 − √ 2 2 i) 4 √ s. In the expressions above, we consider the same constraints for s, i.e. Re(s) > 0 and the roots of s are chosen such that | arg(√s)| < pi 2 , | arg( 4√s)| < pi 4 . (3.3.4) 48 Chapter 3. Linear Well-posedness These solutions are substituted into junction conditions and an (8× 8) coef- ficient matrix M is obtained: 1 0 −1 −1 −1 0 0 0 0 1 sin θ cos θ sin θ cos θ − cos θ sin θ 0 0 0 1 0 0 0 0 −1 −1 −1 0 1 0 0 0 − sin θ cos θ − sin θ cos θ cos θ sin θ sin θ √ s 0 sin θλ1 cos θ sin θλ2 cos θ 0 0 0 0 − sin θ√s 0 0 0 0 − sin θλ1 cos θ − sin θλ2 cos θ 0 0 0 λ1 2 cos θ λ2 2 cos θ 0 λ1 2 cos θ λ2 2 cos θ 0 0 0 λ1 3 cos θ λ2 3 cos θ 0 − λ13 cos θ − λ23 cos θ 0 The determinant of M can be computed directly which gives |M | = s 7/4(16s1/4 sin θ − 4√2 sin(2θ)) cos4 θ . (3.3.5) With the assumption of (3.3.2), we see that (3.3.5) is zero only when s1/4 is negative real, violating (3.3.4). Thus the problem is well-posed if θ ∈ (pi/2, pi). We now consider the marginal case when θ = pi/2. Under this situation, both k2 and k3 do not exist. The solutions are given in the form u1 = A11e −√sσ v1 = B11 u2 = C21 v2 = A22e λ1σ +B22e λ2σ u3 = C31 v3 = A32e λ1σ +B32e λ2σ . (3.3.6) The corresponding determinant of the matrix M is |M | = 16s2. Therefore |M | is not zero for any Re(s) > 0 and the problem is well-posed for θ = pi/2. 49 Chapter 3. Linear Well-posedness In summary, the analysis above gave the same well-posedness condition as the parabolic formulation did. The problem is well-posed provided all angles at the junction are less than pi. The marginal case θ = pi/2 is also well-posed, though physically unreasonable. 50 Chapter 4 Numerical Simulation In this chapter, we present in detail the discretization procedure for the three formulations discussed in Chapter 2. All formulations are applied to the cou- pled surface and grain boundary motion, i.e. the quarter loop problem and the numerical results are compared with the exact travelling wave solution. The most successful approach, the PDAE formulation, is also applied to the Sun-Bauer model and other examples. 4.1 Cartesian Formulation Recall that the cartesian formulation is given by ut = uxx 1 + u2x , t > 0, x > j(t), yt = − [ 1 (1 + y2x) 1/2 [ yxx (1 + y2x) 3/2 ] x ] x , t > 0, x ∈ (−∞, j(t)−) ∪ (j(t)+,∞). (4.1.1) Since the junction is moving it is not straightforward to approximate the system numerically. We shall first fix the junction by making the following transform, x̄ = x− j(t). (4.1.2) We then consider y(x̄, t) and u(x̄, t). Derivatives in the two coordinate sys- tems are related by yx(x, t) = yx̄(x̄, t), yt(x, t) = yt(x̄, t)− yx̄(x̄, t)jt, ux(x, t) = ux̄(x̄, t), ut(x, t) = ut(x̄, t)− ux̄(x̄, t)jt. 51 Chapter 4. Numerical Simulation Junction Free Surface Grain Boundary Free Surface p0 p −1 p1 p0 p −1 p0 p1 Figure 4.1: Sketch of the discretization at the triple junction showing a staggered grid and five ghost values. The triple junction is on the y-axis and the ghost values are shown by the five points on the dashed lines. The system is written in the new coordinate system as ut = ux̄x̄(1 + u 2 x̄) −1 + ux̄(x̄, t)jt, t > 0, x̄ > 0, (4.1.3) yt = − [ 1 (1 + y2x̄)1/2 [ yx̄x̄ (1 + y2x̄)3/2 ] x̄ ] x̄ + yx̄(x̄, t)jt, t > 0, x̄ ∈ (−∞, 0) ∪ (0,∞). The boundary conditions are the same as (2.1.5) except that they are now evaluated at x̄ = 0 instead of x = j(t). Note that the scalar unknown jt(t) appears non-locally in equation (4.1.3). The system in the fixed coordinates can be discretized using standard finite difference schemes on a fixed staggered grid. In the staggered grid, the origin is located at the center of two grid points and thus the triple junction is not a node point (see Figure 4.1). The discretization of the six junction conditions requires five ghost values as shown in Figure 4.1. The curvature κ for the free surface at the junction is approximated by the average of the curvature of the closest two node points p0, p1. The derivative of curvature κs is approximated by the difference of the two curvature divided by the arc length between point p0 and p1. The five ghost values plus the unknown jt(t) make it a solvable system. A numerical result is shown in Figure 4.2. 52 Chapter 4. Numerical Simulation −5 −4 −3 −2 −1 0 1 2 3 4 5 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 Result at time 0.1 Figure 4.2: Numerical result for system (4.1.3) with m = 0.5,∆t = 1e − 6. The initial spatial size is ∆s = 0.2 for the top free surface and ∆s = 0.004 for the grain boundary. Dotted line: numerical result; Solid line: travelling wave solution. The advantage of this method is that it is solved on a fixed grid and the discretization is rather standard except for the junction. Regarding the computational cost, the number of unknowns is only about half of that of the parametric formulations. The disadvantage of using cartesian formulation is that it is not applicable to the non-single-valued case (when 1.81 . m ≤ 2 for travelling waves derived in [31]). Moreover, since the grain boundary is nearly vertical at the junction, it requires a very small grid size for accuracy. 4.2 Parabolic Formulation In this section, we present the discretization of the parabolic scheme (2.2.19), junction conditions (2.2.20)-(2.2.25) and the artificial tangential conditions (2.2.27). The basic approach is to use a staggered grid in σ which is in the fixed interval [0,1] with σ = 0 at the junction and σ = 1 at the boundary of the computational domain. We shall denote the approximations by capital letters with subscripts, i.e., Xj(t) ' X((j−1/2)h, t) = (u((j−1/2)h, t), v((j− 1/2)h, t)) where h is grid spacing and N = 1/h is the number of interior grid points for σ ∈ [0, 1]. 53 Chapter 4. Numerical Simulation We introduce some additional finite difference notations. Let Dk denote the second order centered approximation of the kth derivative, i.e., D1Xj = (Xj+1 −Xj−1)/2h D2Xj = (Xj+1 +Xj−1 − 2Xj)/h2 and let D+ and F denote forward differencing and forward averaging, respec- tively, D+Xj = (Xj+1 −Xj)/h FXj = (Xj+1 +Xj)/2. We shall now discretize each motion separately. 4.2.1 Discretization of the Motions Laws The motion of the grain boundary is a curvature motion and it is described by Xt = Xσσ |Xσ|2 . (4.2.1) This equation is approximated at all grid points by standard differences, Ẋ ij = D2X i j |D1X ij|2 i = 1, j = 1, 2, · · · , N (4.2.2) where Ẋj stands for time derivative. Formally, these discrete equations re- quire values of X0 and XN+1 outside the computation domain. We shall use the boundary conditions to extrapolate the interior values of X1 and XN to the unknown exterior values of X0 and XN+1. We shall give the details of the extrapolation procedure later. The differential equation that describes the motion by surface diffusion is given by Xt = −Xσσσσ|Xσ|4 − 6 |Xσ|σXσσσ |Xσ|5 + (15 |Xσ|2σ |Xσ|4 − 4 |Xσ|σσ |Xσ|3 + |Xσσ|2 |Xσ|4 ) Xσσ |Xσ|2 . (4.2.3) 54 Chapter 4. Numerical Simulation Figure 4.3: Sketch of the ghost points at the triple junction for the discretiza- tion of the parabolic formulation. The higher order derivatives are approximated by (Xσσσ)j ' D3Xj = D2Xj+1 −D2Xj−1 2h (4.2.4) (Xσσσσ)j ' D4Xj = D2Xj−1 +D2Xj+1 − 2D2Xj h2 . (4.2.5) Other terms that involve the derivative of |Xσ| can be approximated similarly using the fact that |Xσ| = √ Xσ ·Xσ. 4.2.2 Discretization of Boundary Conditions The discretization at the triple junction is the most difficult part of the problem. Since there are fourth order derivatives for the surface diffusion we shall need two ghost points for each surface curve and one ghost point for the grain boundary. These ghost points are denoted by X10 , X 2 −1, X 2 0 , X 3 −1, X 3 0 respectively (see Figure 4.3). Because we are using staggered grid, the triple junction is not a grid point, instead it is the center between point X i0 and X i1. The junction conditions (2.2.20)-(2.2.25) are approximated as follows. Condition (2.2.20) is approximated by FX10 = FX20 = FX30 . (4.2.6) 55 Chapter 4. Numerical Simulation The angle conditions (2.2.21)-(2.2.22) are approximated by D+X 1 0 |D+X10 | · D+X 2 0 |D+X20 | = cos θ (4.2.7) D+X 1 0 |D+X10 | · D+X 3 0 |D+X30 | = cos θ. (4.2.8) Condition (2.2.23) is approximated by Fκ20 = −Fκ30, (4.2.9) where Fκi0 = κi1 + κ i 0 2 . (4.2.10) Here κij represents the curvature at point X i j which is approximated by κij = D2X i j |D1X ij|3 · (D1X ij)⊥. (4.2.11) Similarly, condition (2.2.25) is approximated by D+κ 2 0 |D+X20 | = D+κ 3 0 |D+X30 | . (4.2.12) We now consider the discretization of the artificial tangential conditions. Though these conditions are supposed to be applied at the triple junction, it is more convenient to discretize at point X i0. Since these conditions impose a constraint on the grid for uniformity, it is more reasonable to approximate at a grid point instead of a middle point. Therefore, the artificial tangential conditions are approximated by D2X i 0 ·D1X i0 = 0 for i = 2, 3. (4.2.13) We note that since we use staggered grid, both operator F and D+ give second order accuracy at the junction. The numerical simulation has to be carried out on a finite domain, there- fore we need to specify the conditions at the boundary of the computation domain, i.e. at σ = 1. The idea is rather simple. First, we shall keep the 56 Chapter 4. Numerical Simulation size of computational domain fixed and second, we keep all three curves flat. Therefore, the following conditions are used: X i(1, t) = far field given value for i = 1, 2, 3, X iσ(1) = 0 for i = 2, 3. The discretization of these conditions is straightforward. But again we need one ghost point (XN+1) for the grain boundary and two ghost points (XN+1, XN+2) for each of the free surface. The second condition is discretized at point XN+1 instead of the domain boundary. The discretized conditions are given by FX iN = given value for i = 1, 2, 3, D+X i N+1 = 0 for i = 2, 3. 4.2.3 Time Stepping The easiest time discretization is using explicit schemes, for example, the Forward Euler scheme given by X(tn+1) = X(tn) + ∆tF (X(tn)). (4.2.14) Here F (X(tn)) denotes the right hand side expressions in formulation (2.2.19) evaluated at time level tn. This scheme is easy to implement. Given the results at time tn we update the values of the interior grid points by Forward Euler method to time level tn+1 for the three curves respectively. Next, solve for the ghost points, and the triple junction by the boundary conditions. Then go on to the next time level. Time step size ∆t is chosen so that the full discrete scheme is stable and O(h4) is required. This time step is excessively small due to the stiffness of the fourth order equation as noted before. To avoid the excessively small time steps of an explicit scheme, implicit techniques seem more practical for this problem. The most popular implicit scheme is the Backward Euler (BE) method. Given the values at time tn we update the values at time level tn+1 by solving the nonlinear system X(tn+1)−∆tF (X(tn+1))−X(tn) = 0. (4.2.15) 57 Chapter 4. Numerical Simulation 4.2.4 Numerical Results and Tangential Adjustment Since the three curves are strongly coupled by the junction, we solve all unknown points simultaneously including the ghost points. This leads to a nonlinear system which is solved by Newton’s method. Numerical results with time step ∆t = 0.01 are shown in Figure 4.4. All later numerical simulations for the quarter loop problem start with the same position and the same arc length parametrization as shown in Figure 4.4. All results are compared with the travelling wave solution (solid line) that we described before. −6 −4 −2 0 2 4 6 8 10 12 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 initial status −1 −0.5 0 0.5 1 1.5 −0.25 −0.2 −0.15 −0.1 −0.05 0 0.05 0.1 Result at time 0.5 Figure 4.4: Plot of results for parabolic formulation using Backward Euler method with m = 0.5. Left: initial status with grid points; Right: result at t = 0.5 zoomed in near triple junction. Dotted line: numerical result; Solid line: travelling wave solution; Time step size: ∆t = 0.01; Initial spatial size: ∆s = 0.05. The parabolic formulation suffers from numerical difficulties. After some time of computation, the distribution of the grid points often loses its uni- formity, and becomes dense at one end and sparse at the other end. This situation becomes even worse as the computation continues. Grid refinement can not avoid the loss of monotonicity: it is present in the continuum problem due to the lack of maximum principle for the fourth order equations. To overcome this problem, we add an artificial term to adjust the tan- gential velocity of each grid point for the motion of surface diffusion. We 58 Chapter 4. Numerical Simulation consider rescaling the tangential component of the fourth order term in sys- tem (2.2.19). They are replaced by Xt = F (X) + α( Xσσσσ |Xσ|4 · ~t)~t (4.2.16) where F (X) represents the right hand side of the original equation and α is a numerical parameter. The newly added term in (4.2.16) does not influence the normal velocity but it does change the tangential velocity. The value of α will depend on the problem itself and has to be chosen carefully. By trial and error, α = −100 seems to work well for our problem. A numerical result is shown in Figure 4.5. As shown in Figure 4.6, the grid points are more uniform than those in Figure 4.4. The time step size for Figure 4.5 is ∆t = 0.01. A numerical convergence study for this adjusted formulation is shown in Table 4.1. −6 −4 −2 0 2 4 6 8 10 12 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 Result at time 0.5 −1 −0.5 0 0.5 1 1.5 −0.25 −0.2 −0.15 −0.1 −0.05 0 0.05 0.1 Result at time 0.5 Figure 4.5: Plot of results for the parabolic formulation adjusted by equation (4.2.16) with m = 0.5, α = −100,∆t = 0.01,∆s = 0.05. Left: result at t=0.5; Right: result at time t=0.5 zoomed in near triple junction. Dotted line: numerical result; Solid line: travelling wave solution. Although the newly added tangential term improves the numerical per- formance, it requires a good choice for the numerical parameter α. The artificial tangential conditions at the junction and domain boundaries also make the problem more complicated for numerical simulations. This scheme 59 Chapter 4. Numerical Simulation −1 −0.5 0 0.5 1 1.5 −0.2 −0.1 0 0.1 Result at time 0.5 −1 −0.5 0 0.5 1 1.5 −0.2 −0.1 0 0.1 Result at time 0.5 Figure 4.6: Comparison of grid distribution for the parabolic formulation and the adjusted parabolic formulation. Top: parabolic formulation; Bottom: adjusted parabolic formulation. Both pictures are zoomed in near triple junction. Dotted line: numerical result; Solid line: travelling wave solution; Time step size: ∆t = 0.01; Initial spatial size: ∆s = 0.05. ∆t ∆s L2 Norm Rate L∞ Norm Rate 0.2 3.1494e-04 2.0241e-03 ∆t = 0.01∆s2 0.1 7.9775e-05 1.9811 5.4797e-04 1.8852 0.05 2.1445e-05 1.8953 1.4530e-05 1.9151 Table 4.1: Estimated errors and convergence rates for the adjusted parabolic formulation with m = 0.5. Errors are evaluated at time t = 0.02. will never be completely robust: As noted previously, fourth order parabolic problems and other higher order problems as considered in Chapter 6 do not obey the maximum principle. Thus the high order tangential parabolicity does not guarantee a monotone parametrization. In next section, we con- sider the PDAE formulation for which all these numerical issues are resolved naturally. Remark: Amore efficient approach of tangential adjustment for this parabolic formulation is possible. In [39], the authors proposed a formulation to control the tangential velocity which yields a uniform grid for the evolution of curves driven by curvature and external forces. A similar approach is possible for the motion by surface diffusion. However, the PDAE formulation described below seems to be an even better solution. 60 Chapter 4. Numerical Simulation 4.3 PDAE Formulation In this section, we present the discretization of the PDAE system (2.3.7) with A,B = 1 and junction conditions (2.2.20)-(2.2.25). The basic approach is to use a staggered grid in σ and we shall use the same notations as we introduced in previous section. Noticing equation (2.3.4), the curvature is approximated by κij = D2X i j |D1X ij|2 · (D1X ij)⊥, (4.3.1) and consequently the expression (2.3.6) which is used to describe the normal velocity of surface diffusion is approximated by (κss) i j = D2κ i j |D1X ij|2 . (4.3.2) Note that the ghost values X10 , X 2 0 , X 3 0 , κ 2 0 and κ 3 0 are needed to complete the stencils of these discrete equations (see Figure 4.7). These eight values can be related to interior values using the eight junction conditions (2.2.20)-(2.2.25) as shown below. 4.3.1 Discretization of Boundary Conditions The discretization at the junction is subtle. Instead of using two ghost points for each fourth order surface diffusion problem we use only one ghost point (X i0) plus another ghost value (κ i 0) indicating the curvature at the ghost point. Other approaches, such as the parabolic formulation require additional tangential conditions at the triple junction as shown in previous section. The junction conditions are discretized in a similar way as in section 4.2.2, except that when discretizing the following two conditions, Fκ20 = −Fκ30, D+κ 2 0 |D+X20 | = D+κ 3 0 |D+X30 | , the ghost values κ20 and κ 3 0 are applied directly into to the system. Similar to the situation at the triple condition, we need only three ghost points plus two ghost values in total at the boundary of the computational 61 Chapter 4. Numerical Simulation 1 0 X 3 1 X 2 1 X ) ( 2 0 2 0 κ X 1 1 X ) ( 3 0 3 0 κ X t r i p l e j u n c t i o n g r a i n b o u n d a r y Figure 4.7: Sketch of ghost values at the triple junction for the discretization of the PDAE formulation. domain and the following conditions are used: X i(1, t) = far field given value for i = 1, 2, 3, κiσ(1) = 0 for i = 2, 3. The discretization of these conditions is straightforward and it is given by FX iN = given value for i = 1, 2, 3, D+κ i N = 0 for i = 2, 3. For time discretization, we use the Backward Euler method. 4.3.2 Numerical Results The simulation of the quarter loop problem is shown in Figure 4.8. Notice the uniformity of the grid points, enforced by our formulation. A numerical convergence study is shown in Table 4.2. The convergence rates are close to 2 as expected. 62 Chapter 4. Numerical Simulation −1 −0.5 0 0.5 1 1.5 −0.25 −0.2 −0.15 −0.1 −0.05 0 0.05 0.1 Result at time 0.5 Figure 4.8: Plot of result for the PDAE formulation with m = 0.5, zoomed in near triple junction. Dotted line: numerical result; Solid line: travelling wave solution; Time step size: ∆t = 0.01; Initial spatial size: ∆s = 0.05. ∆t ∆s L2 Norm Rate L∞ Norm Rate 0.2 2.7837e-04 1.8996e-03 ∆t = 0.01∆s2 0.1 7.2717e-05 1.9366 5.4444e-04 1.8029 0.05 1.8732e-05 1.9568 1.4732e-04 1.8858 Table 4.2: Estimated errors and convergence rates for the PDAE formulation with m = 0.5. Errors are evaluated at t = 0.02. A comparison of the numerical results for the PDAE formulation and the parabolic formulation without tangential adjustment is shown in Figure 4.9. Both plots are drawn with an even aspect ratio. The grid distribution for the PDAE formulation is much more uniform. Without difficulty we can apply this scheme to the case when the surface curve is not a single-valued function (when 1.81 . m ≤ 2 for travelling waves). An example is shown in Figure 4.10. 63 Chapter 4. Numerical Simulation −1 −0.5 0 0.5 1 1.5 −0.2 −0.1 0 0.1 Result at time 0.5 −1 −0.5 0 0.5 1 1.5 −0.2 −0.1 0 0.1 Result at time 0.5 Figure 4.9: Comparison of grid distribution for the parabolic formulation and the PDAE formulation. Top: parabolic formulation; Bottom: PDAE formulation. Both pictures are zoomed in near triple junction. Dotted line: numerical result; Solid line: travelling wave solution; Time step size: ∆t = 0.01; Initial spatial size: ∆s = 0.05. −0.5 0 0.5 1 1.5 2 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 Result at time 0.2 Figure 4.10: Plot of result for the PDAE formulation with m = 1.96, which has non-single-valued free surface. Dotted line: numerical result; Solid line: travelling wave solution. 4.4 More Examples The PDAE formulation can also efficiently simulate the motion of a curve that involves only mean curvature motion or surface diffusion, as well as 64 Chapter 4. Numerical Simulation the evolution of curve networks if no topological changes occur. In this section, we show several computational examples using the idea of the PDAE formulation. 4.4.1 Examples of Surface Diffusion We first consider examples of normal direction motion that involves only motion by surface diffusion. The first simulation starts with a closed star shaped curve (see Figure 4.11) evolving by surface diffusion. Despite its stiffness, the motion by surface diffusion has some interesting geometrical features. Let S denote the total length of a closed 2D curve, then dS dt = − ∫ V κds = ∫ κssκds = − ∫ κ2sds ≤ 0. Here V stands for the normal velocity and it is equal to −κss in this case. Thus the total length of a curve driven by surface diffusion is decreasing. Let A denote the area enclosed by the curve, then dA dt = − ∫ V ds = ∫ κssds = 0. Therefore, the motion by surface diffusion is also area preserving. According to these properties the star shape curve will tend to evolve into a circle and preserve the area enclosed by itself. The numerical simulation is shown in Figure 4.11. The method conserves the area quite well and the change is about 0.032%. Similar examples have been investigated using level set meth- ods in [15, 54]. Level set methods have unbeatable superiority for interface motion especially when there are topological changes. However, for this sim- ple problem (with no topological changes), our experience is that the PDAE formulation is more efficient and accurate. 65 Chapter 4. Numerical Simulation −0.5 0 0.5 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 Time 0 −0.5 0 0.5 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 Time 7.5e−005 −0.5 0 0.5 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 Time 0.0001745 −0.5 0 0.5 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 Time 0.0003575 Figure 4.11: A computational example that involves only the motion by surface diffusion solved by a PDAE formulation. ∆t = 5e − 7;∆s = 0.05; number of node points: 160. The enclosed area changes by 0.032%. The second simulation starts with a unit square and the results are shown in Figure 4.12. 66 Chapter 4. Numerical Simulation 0 0.5 1 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 Time 0 0 0.5 1 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 Time 5e−005 0 0.5 1 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 Time 0.0005 0 0.5 1 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 Time 0.001 Figure 4.12: A computational example that involves only the motion by surface diffusion solved by a PDAE formulation. ∆t = 1e − 5;∆s = 0.04; number of node points: 100. The enclosed area changes by 0.07%. 4.4.2 More Examples of Coupled Surface and Grain Boundary Motion In this section, we consider two more examples, the Sun-Bauer geometry [40, 65, 66] and the grain growth in a thin film array [8, 62]. Both examples preserve the complexity of the coupled motion but are considered in differ- ent contexts in which travelling wave solutions are not available. Previous studies considered linearized or approximate problems. Our method is able to accurately compute the full nonlinear problem. In the Sun-Bauer geometry, an initially straight grain boundary attaches to the exterior surface with an inclination angle β as shown in Figure 4.13. After annealing, the grain boundary migrates to the right with inclination 67 Chapter 4. Numerical Simulation Grain 2 Grain 1 αβ Figure 4.13: The Sun-Bauer geometry. angle α at the groove root. The same junction conditions as in the quarter loop geometry are applied. In [40, 65, 66], the authors studied the mobility of grain boundary under different inclination angle β. They assumed the grain boundary migrates at a constant speed which is measured by experiments and then the problem is solved by iterating linearized surface diffusion and grain boundary motion alternately. They considered the asymptotic behaviour of the grain boundary when β is small and drew the conclusion that the coupled motion can be separated into two time regimes. In Regime I, the grain boundary turns vertically at the groove root. In Regime II, the turning relaxes following two different paths depending on the value of γ/β where γ = 2arcsin(m 2 ) is the supplementary dihedral angle between the exterior surfaces. If β < γ/6, the grain boundary root moves continuously away from the initial position, but the grain boundary tries to remain straight. For β > γ/6, the grain boundary bends towards the free surface to be inclined with angle γ/6 at the groove root. Figure 4.14-4.15 show the numerical results solved from the PDAE for- mulation with different inclination angle β. The simulations start with a straight line attached to a horizontal surface. The relative surface tension is m = 0.5 and therefore γ/6 = 0.083. Figure 4.14 shows the results with β = 0.05 and Figure 4.15 uses β = 0.9. The grain boundary turns almost perpendicular (α is small) at the beginning and then follows two different paths as described in [40, 65, 66]. For β = 0.05 < γ/6, the inclination angle α remains small. For β = 0.9 > γ/6, α begins to increase and approaches to a value (roughly 0.078) that is smaller than but roughly comparable to the value γ/6 in the linear, approximate theory [40, 65, 66]. The velocity of groove root is also shown as a function of time. It is not constant but the acceleration is small especially after some time of evolution. Our accu- rate nonlinear simulations of the full model show that the predictions of the earlier work are approximately valid. 68 Chapter 4. Numerical Simulation −10 −5 0 5 10 15 20 25 −30 −25 −20 −15 −10 −5 0 5 initial inclination angle β=0.05, γ=0.5, β < γ/6 x y 0 20 40 60 80 100 120 0 0.01 0.02 0.03 0.04 groove root velocity t v 0 20 40 60 80 100 120 2 3 4 5 6 7 8 x 10−3 inclination angle α t α Figure 4.14: Sun-Bauer geometry simulation with β < γ/6 (PDAE formula- tion). Top: the groove and grain boundary profile; Middle: groove velocity versus time; Bottom: inclination angle versus time. 69 Chapter 4. Numerical Simulation 0 5 10 15 20 25 30 35 −25 −20 −15 −10 −5 0 5 initial inclination angle β=0.9, γ=0.5, β > γ/6 = 0.083 x y 0 20 40 60 80 100 120 0 0.1 0.2 0.3 0.4 0.5 groove root velocity t v 0 20 40 60 80 100 120 0.04 0.05 0.06 0.07 0.08 inclination angle α t α Figure 4.15: Sun-Bauer geometry simulation with β > γ/6 (PDAE formula- tion). Top: the groove and grain boundary profile; Middle: groove velocity versus time; Bottom: inclination angle versus time. The second example is the grain growth in a thin film array. In [62], it is predicted that film breakup occurs when grain sizes are larger than the film thickness. They simulated the case when breakup occurs using an approximate linear system, assuming the surface and the boundary deviate slightly from their inclinations. We employed our formulation to compute the evolution of a periodic array in a 2D film where breakup occurs. The computation starts with grain boundaries perpendicular to exterior surface. The qualitative behaviour in our simulation is similar to the results in [62]. However, there are certain geometric anomalies in the earlier work (angles at triple junctions are not matched accurately) from approximations made in their model that are computed correctly by our approach. 70 Chapter 4. Numerical Simulation 0 10 20 30 40 50 60 70 80 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 Result at time 1 x y 0 10 20 30 40 50 60 70 80 −1 −0.5 0 0.5 Result at time 20 x y 0 10 20 30 40 50 60 70 80 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 Result at time 50 x y 0 10 20 30 40 50 60 70 80 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 Result at time 70 x y Figure 4.16: The evolution of film morphology during breakup. 4.4.3 Grain Boundary Motion with Triple Junction Due to its simplicity compared to the mixed order problem, we did not pay much attention to the grain boundary motion shown in Figure 1.6. In this section, we shall consider an example that involves only mean curvature motion. In this model all three curves are enclosed in a unit circle and meet at a triple junction as shown in the top left frame of Figure 4.17. The angles between the curves are 2pi/3 at the junction. The other end of each curve is attached perpendicularly to the circumference and is allowed to move freely along the circle. This model can be used to simulate the evolution of a simple grain boundary network in a polycrystalline structure. The numerical simulation is shown in Figure 4.17. The results are comparable to those in [10] and the PDAE approach is much easier to implement. 71 Chapter 4. Numerical Simulation −1 −0.5 0 0.5 1 −0.5 0 0.5 x y time = 0 −1 −0.5 0 0.5 1 −0.5 0 0.5 x y time = 0.1 −1 −0.5 0 0.5 1 −0.5 0 0.5 x y time = 0.5 −1 −0.5 0 0.5 1 −0.5 0 0.5 x y time = 1 Figure 4.17: A computational example that involves only the motion by mean curvature solved by PDAE formulation. 72 Chapter 5 Stability of the Travelling Wave Solutions In this chapter, we investigate the linear stability of the travelling waves for both the coupled quarter loop problem (see Figure 1.5) and the equal order grain boundary motion problem (see Figure 1.6). The basic approach is to linearize the partial differential equation about the travelling wave solutions and investigate the spectrum of the resulting linear differential operator. The infinite domain is truncated to a bounded but large enough domain and artificial boundary conditions are chosen appropriately. The resulting linear systems together with the linearized boundary conditions are discretized on a staggered grid. Two generalized eigenvalue problems are obtained and solved numerically for the coupled problem and the equal order problem, respectively A necessary condition for stability is that the spectrum of the linear operator lies in the left half plane. The zero eigenvalue corresponding to the translates of travelling waves will always present. However, the truncated problems are not completely equivalent to the original problems which are proposed on an unbounded domain and their spectra may be very different. Much research has been devoted to the study of the underlying connection between the spectrum on an unbounded domain and the spectrum on a bounded domain, e.g., [7, 51], and theories are derived for some cases. Because of the presence of the triple junction, our problems do not fit any existing theoretical framework. However, it turns out that our numer- ical results are similar in appearance to the relevant spectrum theory. Our numerical results suggest that both the coupled problem and the simple equal order problem are stable on any bounded domain. Therefore we have confidence that the travelling wave solutions for the original problems are convectively stable. That is, unstable eigenmodes, if any, can be stabilized by an exponentially weighted norm ([51, 52]). As we have pointed out, it is still an open question on how the numerical eigenvalues on a finite domain 73 Chapter 5. Stability of the Travelling Wave Solutions and the spectrum on an infinite domain are related when triple junctions are present. 5.1 Theoretical Background We first review briefly some theoretical background on the stability analysis of travelling wave solutions. For more detailed discuss we refer to the review article [51] and the references therein. 5.1.1 Stability Criteria for Travelling Waves We consider partial differential equations of the form Ut = A(∂x)U, x ∈ R (5.1.1) where A(∂x) is a nonlinear differential operator. Travelling waves are solutions of the form U(x, t) = Q(x− ct). We intro- duce a moving frame ξ = x− ct, then the travelling wave solution Q(ξ) is a stationary solution of problem Ut = A(∂ξ)U + c∂ξU, ξ ∈ R. (5.1.2) That is, U = Q(ξ) satisfies 0 = A(∂ξ)U + c∂ξU, ξ ∈ R. (5.1.3) We linearize equation (5.1.2) about the steady state Q(ξ) and denote the resulting equation by Ut = LU := D(ξ, ∂ξ)U + c∂ξU, ξ ∈ R (5.1.4) where D(ξ, ∂ξ) is linear and therefore L is a linear operator. It is easy to see that LQξ = 0, (5.1.5) so zero is always in the spectrum of L. This is an immediate consequence of the fact that any translation of a travelling wave is also a travelling wave. 74 Chapter 5. Stability of the Travelling Wave Solutions Therefore, one cannot expect to get asymptotic stability of the waves in a strict sense, i.e., if the initial data of the form U(ξ, 0) = Q(ξ) + ²U0(ξ) (5.1.6) then in general one cannot expect that ‖U(ξ, t)−Q(ξ))‖∞ → 0 as t→∞. (5.1.7) Rather, one should expect that ‖U(ξ, t)−Q(ξ + γ))‖∞ → 0 as t→∞ (5.1.8) for some suitably chosen phase shift γ. An absolute instability occurs if perturbations grow in time at every fixed point in the domain. Convective instabilities are characterized by the fact that, even though the overall norm of the perturbation grows in time, per- turbations decay locally at every fixed point in the unbounded domain; in other words, the growing perturbation is transported, or convected towards infinity. When a convective instability occurs, we refer to the situation if every unstable mode travels to either left or right but not simultaneously to the left and right as transient instability. Transient instabilities can be stabilized by replacing the ‖ · ‖∞ norm in (5.1.8) by ‖ · ‖η defined as ‖u‖2η = ∫ ∞ −∞ |eηxu(x)|2dx (5.1.9) for appropriate value of η ([51, 52]). For example, if η < 0 and the unstable modes travel towards +∞, then term eηx decays exponentially as x → +∞ and may reduce the growth of function u(x) or even cause it to decay. We call patterns convectively stable (transient unstable) if they are un- stable with ‖ · ‖∞ but stable with respect to a suitably chosen exponentially weighted norm ‖ · ‖η. Otherwise, if the wave is not affected by any weighted norm, we say that it is remnantly unstable. Whether a travelling wave is stable or not depends essentially on the spectrum of the linear operator L, i.e., the eigenvalue problem λU = LU. (5.1.10) Roughly speaking, a travelling wave is (absolute) stable if the spectrum (ex- cept 0) of L lies in the open left half plane and zero is a simple eigenvalue. 75 Chapter 5. Stability of the Travelling Wave Solutions Remnant instabilities are captured by the so-called absolute spectrum: the wave experiences a remnant instability if and only if the absolute spectrum has points in the right half plane. Thus a wave is convectively stable (tran- sient unstable) if and only if its absolute but not its essential spectrum is contained in the open left half plane. 5.1.2 Spectrum on Unbounded Domains The spectrum of differential operator L is the union of essential spectrum Σess and point spectrum Σpt: Σ = Σess ∪ Σpt. (5.1.11) We rewrite eigenvalue problem (5.1.10) as a first order ODE system: Vξ = A(ξ, λ)V. (5.1.12) Suppose A(ξ, λ)→ A±(λ) as x→ ±∞ where A±(λ) do not depend on ξ. We define the dispersion relation: d(λ, ν) := det(A±(λ)− Iν). (5.1.13) The essential spectrum of L is then given by Σess(L) = {λ| d(λ, ik) = 0, k ∈ R}. (5.1.14) The point spectrum of L can be calculated by the Evans function defined by D(λ) := det[v1(λ), ..., vn(λ)] (5.1.15) where v1, ...vn are the eigenvectors corresponding to the stable eigenvalues of A+(λ) and A−(λ), i.e, the eigenvectors such that the corresponding eigenval- ues have negative real part for A+(λ) and those such that the corresponding eigenvalues have positive real part for A−(λ). Those λ that satisfy D(λ) = 0 are in the point spectrum of operator L. We refer to [2, 48] and a review article [52] for more details about the Evans function. We remark that there is no point spectrum to the right of the essential spectrum if A+(λ) = A−(λ) since all vi will then be linearly independent. The concept of absolute spectrum also plays an important role in deter- mining the stability. We label eigenvalues of A±(λ) according to their real part, and repeated with their multiplicity, Reν1(λ) ≥ ... ≥ Reνn(λ). (5.1.16) 76 Chapter 5. Stability of the Travelling Wave Solutions We define the absolute spectrum as Σabs = {λ|Reνm(λ) = Reνm+1(λ)} (5.1.17) where m is the asymptotic morse index of the matrices A±(λ) as Reλ→∞. We remark that the absolute spectrum is not a real spectrum in that it is not defined as the set of complex numbers for which a certain operator is not invertible [52]. Regarding absolute spectrum, if we ignore point spectrum, then L is transiently unstable if, and only if, its absolute but not its essential spectrum spectrum is contained in the open left half plane. This instability can be stabilized by an exponentially weighted norm. 5.1.3 Spectrum on Large Bounded Domains In many applications, it appears impossible, or at lest very difficult, to in- vestigate the stability of travelling waves by analytical means. Instead, the problem is usually truncated onto a large but bounded domain where spec- trum can be approximated numerically. The essential spectrum will no longer exist on the bounded domain. In- stead, it is replaced by an accumulation of discrete eigenvalues. These eigen- values are in general different from the continuous spectrum if separated boundary conditions are used. In this case, no matter how large the domain size is, the discrete spectrum will never be close to the spectrum of the op- erator on the real line. It turns out that the discrete eigenvalues converge to the absolute spectrum Σabs as the domain size L tends to infinity if separated boundary conditions are used. In the case of periodic boundary conditions, the set of discrete eigenvalues do converge to the essential spectrum Σess as L tends to infinity. If a travelling wave is Σess unstable but Σabs stable, it is called transient instability and it could be stabilized by exponentially weighted norm. Such a travelling wave is stable on any bounded domain even though it is unsta- ble on the real line. The reader might well ask why a transient instability can be stabilized on a bounded domain? It is because the transient unsta- ble eigenmodes transport perturbations toward either +∞ or −∞ and this process can continue forever on infinite intervals. On bounded intervals, the perturbations are damped at the finite distance boundary conditions. 77 Chapter 5. Stability of the Travelling Wave Solutions 5.1.4 A Simple Example In this section, we illustrate the computation of essential spectrum and ab- solute spectrum using a simple convection-diffusion equation. We consider equation ut = uxx + cux. (5.1.18) The associated eigenvalue problem is given by λu = uxx + cux. (5.1.19) We rewrite this equation as a first order ODE system:( u1 u2 ) x = A(λ) ( u1 u2 ) = ( 0 1 λ −c )( u1 u2 ) . (5.1.20) The dispersion relation d(λ, ν) is given by d(λ, ν) = det(A(λ)− Iν) = det (−ν 1 λ −c− ν ) = ν2 + cν − λ. (5.1.21) Now replace ν by ik and set d(λ, ik) = 0: λ = −k2 + ick, k ∈ R. (5.1.22) This equation gives the essential spectrum and it is described by the parabola c2x+ y2 = 0. (5.1.23) To find the absolute spectrum, we solve for the spatial eigenvalues ν from the dispersion relation which gives ν1,2 = −1 2 c± √ 1 4 c2 + λ. (5.1.24) Then Reν1 = Reν2 requires that λ ≤ −1 4 c2. (5.1.25) This gives the absolute spectrum Σabs = (−∞,−14c2]. Since A(λ) does not depend on x there is no point spectrum. 78 Chapter 5. Stability of the Travelling Wave Solutions 5.2 Linear Stability In this section, we consider numerically the linear stability of the travelling wave solutions derived before. We focus on the more general coupled motion problem which involves both mean curvature motion and surface diffusion. The analysis for the equal order problem is similar and its numerical results will be shown in the next section. Recall that the PDAE formulation for the coupled motion is given by X1t · ~n = κ X2t · ~n = −κss X3t · ~n = −κss (5.2.1) X iσ ·X iσσ = 0 i = 1, 2, 3. Let X i = (ui(σ) + ct, vi(σ)) represent the parameterized travelling wave solution shown in section 2.4. We introduce the moving frame ũ = u − ct. Then system (5.2.1) becomes X1t · ~n = κ− C · ~n X2t · ~n = −κss − C · ~n X3t · ~n = −κss − C · ~n (5.2.2) X iσ ·X iσσ = 0 i = 1, 2, 3 where C = (c, 0)T and c is the travelling wave speed. Obviously, X i0 = (ui(σ), vi(σ)) is a steady state solution of above system. To analyze the linear stability of the travelling wave solution, we consider a small perturbation to the steady state X i0, i.e., X1 = X10 + ²X 1 1 X2 = X20 + ²X 2 1 X3 = X30 + ²X 3 1 and linearize system (5.2.2) around X i0 to get a linear system: X1t · ~n0 = (κ′ − C · ~n) X2t · ~n0 = (−κ′ss − C · ~n) X3t · ~n0 = (−κ′ss − C · ~n) (5.2.3) 0 = (X i0σ ·X iσσ) + (X iσ ·X i0σσ) i = 1, 2, 3, 79 Chapter 5. Stability of the Travelling Wave Solutions where κ′, κ′ss represent the linearized expression for κ and κss respectively. Here we also dropped the subscripts that indicate the O(²) term of X i. The associated eigenvalue problem is given by λX1 · ~n0 = (κ′ − C · ~n) λX2 · ~n0 = (−κ′ss − C · ~n) λX3 · ~n0 = (−κ′ss − C · ~n) (5.2.4) 0 = (X i0σ ·X iσσ) + (X iσ ·X i0σσ) i = 1, 2, 3. The last algebraic equation appearing in (5.2.4) is another reason why standard theory does not apply. When discretized, this leads to a generalized eigenvalue problem. The discretized junction conditions are easily handled in this generalized eigenvalue problem. As a bounded problem truncated from the unbounded domain the follow- ing artificial conditions are implemented at the domain boundary (at σ = 1): uiσ(1) = K0 i = 1, 2, 3 v1σ(1) = 0 v2(1) = v3(1) (5.2.5) κ2(1) = κ3(1). We shall give some explanation of these conditions. Here K0 is a constant determining the grid spacing. The value ofK0 does not appear in the stability analysis since it vanishes after linearization. The last two conditions make the top two surface diffusion curves periodic. This is reasonable since both curves follow the same evolution law and approach zero asymptotically. We remark that these conditions are different from those we used for numerical simulations in Chapter 4. For the equal order grain boundary motion problem, we can do the same procedure to linearize the system. However, for the boundary conditions, it is impossible to implement the standard periodic condition since the asymptotic values at σ = 1 are different for any two of the curves. We use following condition as a substitute for the periodic condition: v1(1) = v2(1)− 1. (5.2.6) The constant value 1 is the vertical distance between the first and the second curve at σ = 1. Therefore, the boundary conditions for the equal order 80 Chapter 5. Stability of the Travelling Wave Solutions problem are given by uiσ(1) = K0 i = 1, 2, 3 v3σ(1) = 0 (5.2.7) v1(1) = v2(1)− 1. We remark that the second condition in both (5.2.5) and (5.2.7) are serv- ing as separated boundary conditions for the associated curves. All artificial boundary conditions and triple junction conditions are lin- earized around the steady state solution. The resulting linear systems are discretized on a staggered grid which give two generalized eigenvalue prob- lems for the mixed order problem and the equal order problem, respectively. The stability of the travelling wave solution on the truncated bounded domain requires that all discrete eigenvalues λ satisfy Reλ < 0 or λ = 0. This PDAE formulation does not fit into the framework that we discussed in section 5.1 because of the presence of the triple junction and the algebraic equations. However, the numerical results give us confidence that both trav- elling waves are convectively stable and the results can be interpreted in the analytic framework discussed above. 5.3 Numerical Results on Stability In this section, we show the numerically computed spectrum for both the mixed order quarter loop problem and the equal order grain boundary motion problem and interpret the associated stability of the travelling wave solutions. 5.3.1 Stability for the Grain Boundary Motion The resulting generalized eigenvalue problem for the equal order problem is solved numerically using the MATLAB function eig and the distribution of the discrete eigenvalues is shown in Figure 5.1. The accumulation of the eigenvalues shows two different patterns in Fig- ure 5.1. First, an elliptic shaped curve and second, a line segment on the negative real axis. As is typical, only those eigenvalues with small magnitude can well resolve the spectrum of the unbounded problem. In our simulations, 81 Chapter 5. Stability of the Travelling Wave Solutions −1600 −1400 −1200 −1000 −800 −600 −400 −200 0 −10 −8 −6 −4 −2 0 2 4 6 8 10 x y Eigenvalues for the Simple Problem numerical result x+(y/c)2=0 −5 −4.5 −4 −3.5 −3 −2.5 −2 −1.5 −1 −0.5 0 −1 −0.5 0 0.5 1 1.5 x y Eigenvalues for the Simple Problem numerical result x+(y/c)2=0 point(−c2/4,0) Figure 5.1: Eigenvalue distribution for the problem of grain boundary mo- tion. Left: distribution of eigenvalues; Right: Close-up near zero. The dots are numerical eigenvalues. The solid line and the star point are predictions. The solid parabola is given by x+ (y/c)2 = 0 and the star point is given by (−c2/4, 0). Note that both the middle line segment and the circle in the left picture are composed of numerical eigenvalues. no significant change for the part of the spectrum close to origin was found by increasing the size of the computational domain or decreasing grid spacing. It seems that the circular part near the origin is a good approximation to the essential spectrum of curvature motion equation (see system 5.2.4), Xt · ~n0 = κ− C · ~n. The essential spectrum of this equation can be calculated analytically by considering following linear equation: vt = vσσ − cvσ (5.3.1) which is obtained by taking σ →∞ and noting the fact that ~n0 → ( 0 1 ) , κ→ vσσ, as σ →∞. The essential spectrum of equation (5.3.1) is then given by Σess = {λ|λ = −k2 − ick, k ∈ R}, (5.3.2) 82 Chapter 5. Stability of the Travelling Wave Solutions or equivalently, the parabola x + (y/c)2 = 0 in the xy−plane. This curve is plotted as the solid line in Figure 5.1. As we increase the domain size, the ellipse opens wider and gets even closer to the solid parabola. Our feeling is that this phenomena can be explained by the periodic conditions that we used for the top two curvature motion curves. When periodic boundary conditions are used the numerical eigenvalues approximate the essential spectrum. This is proved true for equations without junction, and it also seems to be true for our problem when triple junctions are present. The second accumulation of points converges to a line segment on the left real axis. This line turns out to be an approximation of the absolute spectrum for equation (5.3.1) which is given analytically by Σabs = {λ ≤ −c 2 4 , λ ∈ R}. (5.3.3) The point (−c2/4, 0) is marked by the star point in Figure 5.1. The right end of the line segment comes closer to the star point as we increase the domain size and refine the grid. The reason that we believe it approximates the absolute spectrum is because of the separated boundary condition used in (5.2.7) for the third curvature motion curve. This again coincides with the conclusion that separated boundary conditions lead to absolute spectrum. Any point in the spectrum of the bounded operator is a point spectrum. Besides those lying in one of the two continuous patterns and the zero eigen- value we observe no other point spectrum in the diagram. It is reasonable considering the fact that the resulting linear system of (5.2.4) has asymptotic constant coefficient (does not depend on σ) as σ →∞. Therefore, the Evans function is not zero for any λ which means there is no point spectrum. The zero eigenvalue seems to be of multiplicity two. One of them corresponds to a horizontal shift and the other corresponds to a vertical shift. Two eigenvectors are plotted in Figure 5.2. We choose eigenvalues that are close to the origin since only eigenvectors corresponding to small magnitude eigenvalues will well resolve the eigenfunctions on the unbounded domain. These eigenvectors include perturbations to the travelling wave for both x coordinate and y coordinate. For better illustration of their behaviour, the eigenvectors ν are added to the original travelling wave solution X0 and the new vectors X0 + ν are plotted. Travelling wave solutions with different wave speed are also investigated. All have similar behaviour to those shown above. 83 Chapter 5. Stability of the Travelling Wave Solutions −10 −5 0 5 10 15 20 25 −1.5 −1 −0.5 0 0.5 1 1.5 eigenvalue=−4.4594 −15 −10 −5 0 5 10 15 20 −1 −0.5 0 0.5 1 1.5 2 eigenvalue=−4.365+0.99984i Figure 5.2: Typical eigenvectors are plotted after being added to the trav- elling wave solution. Left: when eigenvalue is real number; Right: when eigenvalue is complex number. We are interested in the stability of the travelling wave on the unbounded domain. It seems that the numerical results partly reproduced the essential spectrum and the absolute spectrum of the original problem. If this is true, we can say the travelling wave is at least convectively stable since its absolute spectrum is contained in the left half plane. There is no rigorous theory to support the conclusion that our method reproduces the essential spectrum or absolute spectrum. However, the numerical results give evidence that the problem is stable on any bounded domain and therefore the travelling wave is convectively stable. In other words, the travelling wave is pointwise stable and any unstable eigenmode, if exist, is transported to either −∞ or ∞ and it can be stabilized by an appropriately chosen exponential weight. 5.3.2 Stability for the Coupled Surface and Grain Boundary Motion In this section, we consider the stability of the travelling wave solution for the coupled problem derived by Kanel et al. in [30] (see also in section 2.4.2). The distribution of eigenvalues for the coupled surface and grain boundary motion is shown in Figure 5.3. The behaviour of the eigenvalues for this coupled problem is similar to the equal order case. 84 Chapter 5. Stability of the Travelling Wave Solutions −5 −4 −3 −2 −1 0 x 104 −15 −10 −5 0 5 10 15 x y Eigenvalues for Coupled Problem numerical result x+(y/c)4=0 −20 −15 −10 −5 0 −4 −3 −2 −1 0 1 2 3 4 x y Eigenvalues for Coupled Problem numerical result x+(y/c)4=0 point(−c2/4,0) Figure 5.3: Eigenvalue distribution for the coupled surface and grain bound- ary motion. Left: distribution of eigenvalues; Right: Close-up near zero. The dots are numerical eigenvalues. The solid line and the star point are predictions. The solid parabola is given by x+(y/c)4 = 0 and the star point is given by (−c2/4, 0). The eigenvalues show two different patterns. One of them is more like a parabola but it closes up at the far left end. The other is a relatively short line segment on the negative real axis (almost invisible in the left picture). We can use the same reasoning as above to explain them. Since periodic boundary conditions are used for the top two surface diffusion curves the numerical parabola now approximates the essential spectrum of vt = −vσσσσ − cvσ (5.3.4) which is obtained by taking σ →∞ for the linearized surface diffusion motion in system (5.2.4). The analytical essential spectrum is given by Σess = {λ|λ = −k4 − ick, k ∈ R} (5.3.5) or equivalently the curve x+(y/c)4 = 0 in the xy−plane. This curve is plotted as the solid line in Figure 5.3. As we increase the domain size the numerical parabola (dotted line) opens wider and resolves the analytical parabola (solid line) better. The line segment on the real axis approximates the absolute spectrum of equation (5.3.1) with c being the travelling wave speed of the coupled 85 Chapter 5. Stability of the Travelling Wave Solutions problem. Note that it does not reproduce the absolute spectrum of equation (5.3.4) since the separated boundary conditions are applied to the curvature motion curve, the second order equation not the fourth order equation. The point (−c2/4, 0) is marked by the star point in Figure 5.3. The right end of the line segment approaches to this point as we increase the domain size and refine the grid. Different from the equal order problem, the coupled problem is of mixed order. Periodic boundary conditions are applied to the fourth order equa- tions, while separated boundary conditions are used for the second order equation. Therefore, numerical results exhibit the essential spectrum of the fourth order equation and the absolute spectrum of the second order equa- tion. There is no point spectrum in the diagram except the zero eigenvalue and those falling into one of the continuous patterns. Again the zero eigenvalue is of multiplicity two. One corresponds to a horizontal shift and the other corresponds to a vertical shift. Two eigenvectors corresponding to small eigenvalues are added to the travelling wave solution and plotted in Figure 5.4. −30 −20 −10 0 10 20 30 −1.4 −1.2 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 eigenvalue = −2.0718 −30 −20 −10 0 10 20 30 −1.2 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 eigenvalue = −3.7163+2.3759i Figure 5.4: Typical eigenvectors are plotted after being added to the trav- elling wave solution. Left: when eigenvalue is real number; Right: when eigenvalue is complex number. The conclusion on the stability of this travelling wave solution is the same as the equal order problem. Numerical results show stability on any bounded 86 Chapter 5. Stability of the Travelling Wave Solutions domain and therefore we believe the travelling wave is convectively stable. In other words, it is stable for every fixed point in the real line, even though the perturbation when measured by standard ‖ · ‖∞ norm could become unbounded. No existing theory covers the spectral analysis for the PDAE case with junction. Our numerical results show that many features the theory predicts for simpler problems appear here, i.e., when an unbounded domain is trun- cated to a bounded domain, the spectrum of the bounded problem converges to the essential spectrum of the unbounded problem when periodic boundary conditions are used; and it converges to the absolute spectrum when sepa- rated boundary conditions are used. However, rigorous analysis is still an open problem. To conclude this section, we show a 2D stability test for the quarter loop problem. In Figure 5.5, a large perturbation to a travelling wave solution is given. The travelling wave structure redevelops after some time. This and other runs we conducted with similar behaviour indicate the 2D stability of these travelling wave solutions and corroborate the conclusions above. −5 0 5 10 15 20 −0.4 −0.2 0 0.2 Initial Status −5 0 5 10 15 20 −0.4 −0.2 0 0.2 Result at time 1 −5 0 5 10 15 20 −0.4 −0.2 0 0.2 Result at time 5 −5 0 5 10 15 20 25 30 35 40 −0.4 −0.2 0 0.2 Result at time 20 Figure 5.5: Results for initial data that is a perturbation of a travelling wave. 87 Chapter 5. Stability of the Travelling Wave Solutions 5.4 3D Stability In this section we consider the linear stability of the travelling waves with 3D perturbations. We focus only on the more complicated mixed order problem. Let X = (x(α, β), y(α, β), z(α, β)) represent a surface in 3D. One has the following formula for mean curvature: H = G11B22 − 2G12B12 +G22B11 2(G11G22 −G212) where G11 = ∂X ∂α · ∂X ∂α ,G22 = ∂X ∂β · ∂X ∂β ,G12 = ∂X ∂α · ∂X ∂β B11 = ∂2X ∂α2 · ~n,B22 = ∂ 2X ∂β2 · ~n,B12 = ∂ 2X ∂α∂β · ~n and ~n denotes the unit normal given by ~n = ∂X ∂α × ∂X ∂β ‖∂X ∂α × ∂X ∂β ‖ . This expression can also be rewritten as H = (Xα ·Xα)Xββ − 2(Xα ·Xβ)Xαβ + (Xβ ·Xβ)Xαα 2(Xα ·Xα)(Xβ ·Xβ)− 2(Xα ·Xβ)2 · ~n. Recall that the motion by surface diffusion has normal velocity being equal to the surface Laplacian of mean curvature and the surface Laplacian operator is defined by ∆s = ∇s · ∇s where ∇s = ∇− n∂n. (5.4.1) The tangential gradient ∇s applied to a function H defined on a surface X = (x(α, β), y(α, β), z(α, β)) is given by ∇sH = g11∂H ∂α Xα + g 12∂H ∂β Xα + g 21∂H ∂α Xβ + g 22∂H ∂β Xβ where gij are the components of the inverse matrix of (gij) = ( Xα ·Xα Xα ·Xβ Xα ·Xβ Xβ ·Xβ ) . 88 Chapter 5. Stability of the Travelling Wave Solutions The surface Laplacian operator is then given by ∆sH = 1√ g ( ∂ ∂α (g11 √ gHα + g 12√gHβ) + ∂ ∂β (g21 √ gHα + g 22√gHβ)). (5.4.2) where g = det(gij). We consider a surface parameterized as X = (x(α, β), y(α, β), β), where parameter β is the z coordinate. The coupled surface and grain boundary motion in 3D is then described by X1t · ~n = H X2t · ~n = −∆sH X3t · ~n = −∆sH (5.4.3) X iα ·X iαα = 0 i = 1, 2, 3, where H represents the mean curvature. Since the third component β is fixed, we have Xt · ~n = xt · nx + yt · ny where nx, ny represent the components of the unit normal. Most of the junction and boundary conditions can be extended in a straightforward way to the 3D case. We show only the junction condition that reflects the balance of mass flux. It is now interpreted as ∇sH2 · (X2β × ~n) = ∇sH3 · (X3β × ~n) (5.4.4) where the superscripts represent curve indices. Let X0 = (x0(α, t), y0(α, t)) represent the travelling wave solution in 2D. A travelling wave solution for the 3D problem can be expressed in the form of X0 = (x0, y0, z) = (x0(α, t), y0(α, t), β) by extending along the z direction (see Figure 5.6). We remark that the travelling wave speed in 3D is only half of that in 2D after such an extension. This comes from the fact that X0α = (x0α, y0α, 0), X0αα = (x0αα, y0αα, 0), X0β = (0, 0, 1), X0ββ = (0, 0, 0), X0αβ = (0, 0, 0), ~n = (nx, ny, 0). 89 Chapter 5. Stability of the Travelling Wave Solutions x y z X0 = (x0, y0) X0 = (x0, y0, z) Figure 5.6: Extension of the 2D travelling wave solution to 3D. Therefore, the mean curvature of the surface is given by H = X0αα 2(X0α ·X0α) · ~n = x0αα · nx + y0αα · ny 2(x20α + y 2 0α) , which is only half of the curvature of the 2D curve (x0, y0) given by κ = x0αα · nx + y0αα · ny (x20α + y 2 0α) . For simplicity, the wave speed is still denoted by c. We shall now repeat the procedure for the 2D linear stability analysis. We consider a small perturbation to the steady state solution in the moving frame x̃ = x− ct: x̃ = x0 + ²x̃1, ỹ = y0 + ²ỹ1, and assume x̃1 = e iωzx1(α), ỹ1 = e iωzy1(α), (5.4.5) where ω is wave number. 90 Chapter 5. Stability of the Travelling Wave Solutions The PDAE system (5.4.3) is then linearized around (x0, y0, z) and the resulting linear system for x1, y1 is discretized to get a generalized eigenvalue problem with wave number ω in the coefficient matrix. We again use periodic boundary conditions for the top surface that evolves by surface diffusion and separated boundary conditions for the bottom sur- face that evolves by mean curvature motion. The discrete eigenvalue problem is solved numerically and the distribution of eigenvalues with wave number ω = 1 is shown in Figure 5.7. −2.5 −2 −1.5 −1 −0.5 0 x 104 −6 −4 −2 0 2 4 6 x y Eigenvalues for Coupled Problem in 3D numerical result x+(y/c)4/2+ω2(y/c)2+ω4/2=0, ω=1 −10 −8 −6 −4 −2 0 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 x y Eigenvalues for Coupled Problem in 3D numerical result x+(y/c)4/2+ω2(y/c)2+ω4/2=0, ω=1 point (−c2/2−ω2/2,0), ω=1 Figure 5.7: Eigenvalue distribution for the coupled surface and grain bound- ary motion in 3D with wave number ω = 1. Left: distribution of eigen- values; Right: Close-up near zero. The dots are numerical eigenvalues. The solid line and the star point are predictions. The solid parabola is given by x + (y/c)4/2 + ω2(y/c)2 + ω4/2 = 0 and the star point is given by (−c2/2− ω2/2, 0). Similar to the 2D case, the eigenvalues form two continuous patterns: one parabola and one small line segment. The parabola approximates the essential spectrum of equation ut = −1 2 uσσσσ + ω 2uσσ − cuσ − 1 2 ω4u (5.4.6) where ω is the wave number introduced in equation (5.4.5). This equation is again obtained by taking σ → ∞ for the linearized equation of the motion 91 Chapter 5. Stability of the Travelling Wave Solutions by surface diffusion. The essential spectrum of above equation is given by Σess = {λ|λ = −1 2 k4 − ω2k2 − 1 2 ω4 − ick, k ∈ R} (5.4.7) or equivalently the curve x+ (y/c)4/2+ ω2(y/c)2+ ω4/2 = 0 which is drawn in Figure 5.7. The line segment on the negative real axis approximates the absolute spectrum of equation ut = 1 2 uσσ − cuσ − 1 2 ω2u (5.4.8) which is given by Σabs = {λ ≤ −1 2 c2 − 1 2 ω2, λ ∈ R}. (5.4.9) The point (−c2/2− ω2/2, 0) is marked by the star point in Figure 5.7. The approximation becomes more accurate as we increase the domain size and refine the grid. The conclusion on the stability of this travelling wave is similar to the 2D case. Since all numerical eigenvalues are contained in the open left half plane it is stable on any bounded domain. Note that unless ω = 0 there is no zero eigenvalue any more. All eigenvalues move further to the left as we increase the wave number ω. Therefore, the numerical results show that the travelling wave solution is convectively stable with 3D perturbations. Remark on 3D computations: For real 3D simulations, one may use the following equations to describe the motion of a surface X = (x(α, β), y(α, β), z(α, β)): X1t · ~n = H X2t · ~n = −∆sH (5.4.10) X3t · ~n = −∆sH The mean curvatureH and the surface Laplacian of mean curvature are given by H = (Xα ·Xα)Xββ − 2(Xα ·Xβ)Xαβ + (Xβ ·Xβ)Xαα 2(Xα ·Xα)(Xβ ·Xβ)− 2(Xα ·Xβ)2 · ~n, 92 Chapter 5. Stability of the Travelling Wave Solutions and ∆sH = 1√ g ( ∂ ∂α (g11 √ gHα + g 12√gHβ) + ∂ ∂β (g21 √ gHα + g 22√gHβ)), where g = det(gij) and (gij) = ( Xα ·Xα Xα ·Xβ Xα ·Xβ Xβ ·Xβ ) . In previous linear stability analysis, we need only one constraint since the parameter β is fixed as the z coordinate. However, in order to produce a good mesh for numerical computation, one usually should allow the motion of the grid nodes in both directions. Therefore, one would need two constraints to maintain a suitable parametrization during the evolution. However, it does not seem clear how to impose appropriate constraints for 3D computations. 5.5 Conclusions Numerical results show that both the equal order problem and the coupled problem are stable when truncated onto a bounded domain. Therefore, both problems seem to be convectively stable. We also investigated the numer- ical stability for the coupled problem with 3D perturbations and the same conclusion applies. No existing spectral theory covers the case of our PDAE system with junc- tions. Numerical results make us suspect that related conclusions will also apply for this more complicated situation: first, the essential spectrum of the original problem is determined by the essential spectrum of the limit equation at infinity; second, when an unbounded domain is truncated to a bounded domain, the spectrum of the bounded problem converges to the essential spectrum of the unbounded problem when periodic boundary conditions are applied, and it converges to the absolute spectrum when separated bound- ary conditions are applied. It is also interesting to notice that our numerical results show mixed spectra due to the mixed type of far field conditions at the artificial domain boundary. 93 Chapter 6 Extension to a Class of High Order Problems In this chapter, we extend our discussion to a class of high order three- phase boundary motion problems, including system formulations, numerical simulations and linear well-posedness analysis. In the general class of high order problems we consider here, three phase boundaries evolve under certain motion laws with specified normal velocities. They meet at a triple junction with appropriate conditions applied. The discussion in Chapter 2 is extended to describe the new problems. A system of partial differential algebraic equations(PDAE) is proposed. The approach for the linear well-posedness analysis in Chapter 3 also applies in this chapter. All problems in our class are shown to be well-posed if all three boundaries evolve under the same evolution law. For problems involving more than one type of motion we show the well-posedness case by case for some examples. 6.1 Mathematical Formulation We consider a class of three-phase boundary motion problems in two dimen- sions as shown in Figure 6.1. There are three phases in this model, phases I, II and III. Their interfaces Γ1,Γ2 and Γ3 are undergoing given evolution laws with parabolic normal velocities given by curvature or its even order derivatives with respect to arc length. They meet at a triple junction with prescribed angles θ1, θ2 and θ3. Depending on the type of motions for the interfaces, extra junction conditions may needed. Following the idea of Chapter 2, we represent Γ1,Γ2 and Γ3 byX 1, X2 and X3 respectively and consider the following generalized three-phase problem 94 Chapter 6. Extension to a Class of High Order Problems I III Γ1 Γ2 Γ3 θ1 θ3 θ2 II Figure 6.1: A three-phase geometry. described by a PDAE formulation, X1t · ~n = (−1)m1κ(2m1) X2t · ~n = (−1)m2κ(2m2) X3t · ~n = (−1)m3κ(2m3) (6.1.1) X iσ ·X iσσ = 0 i = 1, 2, 3, where κ(2mi), i = 1, 2, 3, represent the derivative of curvature κ of order 2mi with respect to arc length s and mi are non-negative integers satisfying, 0 ≤ m1 ≤ m2 ≤ m3. (6.1.2) The evolution laws are fulfilled by the first three equations for all three in- terfaces. The last algebraic equation has the same function as before to maintain a uniform parametrization during the evolution. To propose appropriate junction conditions, we first consider the equal order case when m1 = m2 = m3 = m. The simplest equal order case is the problem of grain boundary motion with m = 0 and the junction conditions have been discussed in previous chapters as well as in [9]. The six junction conditions are given by X1(0) = X2(0) = X3(0) X1s ·X2s = cos θ3 X1s ·X3s = cos θ2. 95 Chapter 6. Extension to a Class of High Order Problems When m = 1, the system (6.1.1) is reduced to a fourth order problem which has been discussed in [25]. The continuity of the chemical potentials and the balance of mass flux give two more conditions (see [25] for more details): γ1κ 1 + γ2κ 2 + γ3κ 3 = 0 γ1κ 1 s = γ2κ 2 s = γ3κ 3 s where γi represents the surface energy of the interface Γi. The superscripts represent the indices of the interfaces. Motivated by the conditions above, we consider the following junction conditions for the equal order problems when m ≥ 1: X1(0) = X2(0) = X3(0) X1s ·X2s = cos θ3 X1s ·X3s = cos θ2 γ1κ (2i−2) 1 + γ2κ (2i−2) 2 + γ3κ (2i−2) 3 = 0 i = 1, · · · ,m γ1κ (2i−1) 1 = γ2κ (2i−1) 2 = γ3κ (2i−1) 3 i = 1, · · · ,m where the superscripts in the brackets represent the order of derivatives with respect to arc length. Note that the curve indices of κ are represented by the subscripts instead of the usual superscripts. Each equal order problem has 3m+ 6 junction conditions in total. The junction conditions for the mixed order problems are more compli- cated. The coupled surface and grain boundary motion studied in previous chapters is one of the simplest cases. The discussion there suggests the fol- lowing junction conditions for the mixed higher order problems: X1(0) = X2(0) = X3(0) X1s ·X2s = cos θ3 X1s ·X3s = cos θ2 γ1κ (2i−2) 1 + γ2κ (2i−2) 2 + γ3κ (2i−2) 3 = 0 i = 1, · · · ,m1 γ1κ (2i−1) 1 = γ2κ (2i−1) 2 = γ3κ (2i−1) 3 i = 1, · · · ,m1 γ2κ (2i−2) 2 + γ3κ (2i−2) 3 = 0 i = m1 + 1, · · · ,m2 γ2κ (2i−1) 2 = γ3κ (2i−1) 3 i = m1 + 1, · · · ,m2 γ3κ (2i−1) 3 = 0 i = m2 + 1, · · · ,m3. When m1 = 0,m2 = m3 = 1, these conditions are reduced to the junction conditions for the mixed order quarter loop problem with the assumption that γ2 = γ3. 96 Chapter 6. Extension to a Class of High Order Problems We remark that these high order problems are more mathematically in- teresting rather than physically interesting. Some of them may not have physical meanings. 6.2 Linear Well-posedness In this section, we use the same approach discussed in Chapter 3 to show the well-posedness of the generalized three-phase boundary motion problems described above. For convenience, we take γ1 = γ2 = γ3 in this section and the rest of this chapter. This implies equal angle conditions and similarity of the phases in the higher order conditions. We first linearize the problems around fixed straight line solutions at the junction. We then consider the resulting linear system for all equal order problems and some examples of the mixed order case, respectively. 6.2.1 Linearization of the System We again consider a perturbation expansion around the tangential direction at the triple junction for each curve, i.e., X1 = d1σ + ²X̄ 1 +O(²2) X2 = d2σ + ²X̄ 2 +O(²2) X3 = d3σ + ²X̄ 3 +O(²2) where all notations have the same meaning as before. The resulting linear system for these new problems is given by X1t · d⊥1 = (−1)m1X1(2m1+2) · d⊥1 X2t · d⊥2 = (−1)m2X2(2m2+2) · d⊥2 X3t · d⊥3 = (−1)m3X3(2m3+2) · d⊥3 0 = di ·X iσσ i = 1, 2, 3 (6.2.1) where (2mi+2) represent derivative of order 2mi+2 with respect to the scaled arc length σ and the superscripts outside the brackets are curve indices. To show the linearization of the junction conditions we consider the more 97 Chapter 6. Extension to a Class of High Order Problems general mixed order case and the linearized conditions are given by X1(0) = X2(0) = X3(0) d1 ·X2σ + d2 ·X1σ − (d1 · d2)(d1 ·X1σ + d2 ·X2σ) = 0 d1 ·X3σ + d3 ·X1σ − (d1 · d3)(d1 ·X1σ + d3 ·X3σ) = 0 X1(2m) · d⊥1 +X2(2m) · d⊥2 +X3(2m) · d⊥3 = 0, m = 1, · · · ,m1 (6.2.2) X1(2m+1) · d⊥1 = X2(2m+1) · d⊥2 = X3(2m+1) · d⊥3 , m = 1, · · · ,m1 X2(2m) · d⊥2 +X3(2m) · d⊥3 = 0, m = m1 + 1, · · · ,m2 X2(2m+1) · d⊥2 = X3(2m+1) · d⊥3 , m = m1 + 1, · · · ,m2 X3(2m+1) · d⊥3 = 0, m = m2 + 1, · · · ,m3. The linearized junction conditions for the equal order case can be obtained as part of (6.2.2). Without loss of generality, we specify one of the tangential directions, say d1 = (0,−1)T . Further we assume di,1, di,2 6= 0 for i = 2, 3. Then system (6.2.1) has solutions in the form u1 = ∑m1+1 j=1 C1je λ1jσ v1 = C1,m1+2 ui = ∑mi+1 j=1 Cije λijσ + Ci,mi+2, i = 2, 3 vi = −ki( ∑mi+1 j=1 Cije λijσ) + 1 ki Ci,mi+2, i = 2, 3 (6.2.3) where ki = di1/di2 and λij are roots of equation λ 2mi+2 = (−1)mis satisfying pi 2 < arg(λij) < 3pi 2 . Here s stands for the Laplace transform variable and it is chosen such that | arg(s)| < pi 2 . (6.2.4) Thus both s and the roots of s have positive real part. When these solutions are substituted into the linearized boundary con- ditions (6.2.2), a square matrix M of size (m1 +m2 +m3) + 6 is obtained. Whether the problem is well-posed or not depends on the determinant of matrix M . We shall discuss this determinant for all equal order problems and some mixed order problems respectively. 98 Chapter 6. Extension to a Class of High Order Problems 6.2.2 Equal Order Problems We first consider the case when m1 = m2 = m3. We denote the value of mi by m (the order of each equation is 2m + 2). For simplicity, we assume the angles between the curves at the triple junction are θi = 2pi/3. Then one has d1 = 0 −1 d2 = − √ 3 2 1 2 d3 = √ 3 2 1 2 . The coefficient matrix M is given by 1 · · · 1 0 −1 · · · −1 −1 0 · · · 0 0 0 · · · 0 1 −√3 · · · −√3 √ 3 3 0 · · · 0 0 1 · · · 1 0 0 · · · 0 0 −1 · · · −1 −1 0 · · · 0 1 0 · · · 0 0 √3 · · · √3 − √ 3 3 λ1 · · · λm+1 0 2λ1 · · · 2λm+1 0 0 · · · 0 0 λ1 · · · λm+1 0 0 · · · 0 0 2λ1 · · · 2λm+1 0 λ1 2 · · · λm+12 0 −2λ12 · · · −2λm+12 0 −2λ12 · · · −2λm+12 0 λ1 3 · · · λm+13 0 2λ13 · · · 2λm+13 0 0 · · · 0 0 λ1 3 · · · λm+13 0 0 · · · 0 0 2λ13 · · · 2λm+13 0 ... ... ... ... ... ... ... ... ... ... ... λ1 2m · · · λm+12m 0 −2λ12m · · · −2λm+12m 0 −2λ12m · · · −2λm+12m 0 λ1 2m+1 · · · λm+12m+1 0 2λ12m+1 · · · 2λm+12m+1 0 0 · · · 0 0 λ1 2m+1 · · · λm+12m+1 0 0 · · · 0 0 2λ12m+1 · · · 2λm+12m+1 0 . Here we use the same set of λi for all three curves. Note that the fifth row and the sixth row have been multiplied through by −2√3/3 and 2√3/3, respec- tively. For the convenience of the determinant computation, we rearrange 99 Chapter 6. Extension to a Class of High Order Problems the matrix above. Appropriate row operations yield 2 √ 3 3 2 3 0 0 −2 0 0 −2 0 0 −2 0 0 1 0 −1 1 0 −1 1 0 · · · −1 1 0 0 0 −1 1 0 −1 1 0 −1 1 0 −1 0 A1 A2 · · · Am+1 0 λ21A1 λ 2 2A2 · · · λ2m+1Am+1 ... ... ... ... 0 λ2m1 A1 λ 2m 2 A2 · · · λ2mm+1Am+1 (6.2.5) where Ai are 3× 3 block matrices in the form of Ai = 1 −2 −2 λi 2λi 0 λi 0 2λi . (6.2.6) The determinant of matrix M is given recursively by M0 = −8 √ 3λ21 Mm = 12λ 2 m+1 m∏ i=1 (λ2m+1 − λ2i )3 ·Mm−1, where Mi is the determinant of M when m = i. Now recall that λj are non- repeated roots of (−1)ms that have negative real part and whose arguments are in (pi/2, 3pi/2). Therefore, λi 6= 0 and λi 6= ±λj for i 6= j. This guarantees that the determinant of M is not zero for any m. Thus all equal order problems are well-posed when θi = 2pi/3. Physically, γ1 = γ2 = γ3 implies θi = 2pi/3. However, a similar derivation can be applied to the case of arbitrary angles and the determinant is given by M0 = (sin(θ2 + θ3)− sin(θ2)− sin(θ3))λ21 sin θ2 sin θ3 cos θ2 cos θ3 Mm = 3λ2m+1 cos θ2 cos θ3 m∏ i=1 (λ2m+1 − λ2i )3 ·Mm−1. The determinant is not zero provided 0 < θi < pi. We can further see that when θ1 = pi and θ2, θ3 > 0, above problem is also well-posed. Now recall 100 Chapter 6. Extension to a Class of High Order Problems that we have assumed di,1, di,2 6= 0 for i = 2, 3. This excludes the case when θi = pi/2 or pi. If only one angle, say θ2 = pi/2, one can rotate the system such that Γ2 is pointing to the direction (0,−1)T . Hence above discussion still applies and the problem is well-posed if 0 < θ1, θ3 < pi . Now we consider the case when θ2 = θ3 = pi/2. In this case, the solution for system (6.2.1) becomes u1 = ∑m+1 j=1 C1je λjσ v1 = C1,m+2 ui = Ci1 i = 2, 3 vi = ∑m+2 j=2 Cije λjσ + Ci,m+2 i = 2, 3. (6.2.7) The associated determinant is now given by M0 = 2λ 2 1 Mm = −3λ2m+1 m∏ i=1 (λ2m+1 − λ2i )3 ·Mm−1. Therefore, the problem is still well-posed. The case of θi = pi can be discussed similarly. In summary, all equal order problems are well-posed provided 0 < θi ≤ pi. 6.2.3 Mixed Order Problems The same procedure can be applied to the mixed order problems. However, no common determinant formula is found. We can still verify case by case that all tested problems are well-posed. For example, ifm1 = 1,m2 = m3 = 2 and θi = 2pi/3, the coefficient matrix M is given by 101 Chapter 6. Extension to a Class of High Order Problems 2 √ 3 3 2/3 0 0 0 −2 −2 −2 0 0 0 0 1 0 −1 −1 1 1 1 0 0 0 0 0 −1 1 1 0 0 0 −1 −1 −1 0 0 0 1 1 −2 −2 −2 −2 −2 −2 0 0 0 λ1 λ2 2ω1 2ω2 2ω3 0 0 0 0 0 0 λ1 λ2 0 0 0 2ζ1 2ζ2 2ζ3 0 0 0 λ21 λ 2 2 −2ω21 −2ω22 −2ω23 −2ζ21 −2ζ22 −2ζ23 0 0 0 λ31 λ 3 2 2ω 3 1 2ω 3 2 2ω 3 3 0 0 0 0 0 0 λ31 λ 3 2 0 0 0 2ζ 3 1 2ζ 3 2 2ζ 3 3 0 0 0 0 0 −2ω41 −2ω42 −2ω43 −2ζ41 −2ζ42 −2ζ43 0 0 0 0 0 −2ω51 −2ω52 −2ω53 2ζ51 2ζ52 2ζ53 (6.2.8) where the roots λij are denoted by λi, ωi, ζi respectively for the three curves. The determinant of the matrix above after substituting in the values for λi, ωi, ζi is given by |M | ≈ −2660.4 i (1.4 + 2.0 s1/12 + 5.7 s2/12 + 12.0 s3/12 + 2.8 s4/12) s39/12. Recall that both s and the roots of s have positive real part, therefore the sum in the brackets has positive real part. Thus the determinant of M is not zero and the corresponding problem is well-posed with at most algebraic growth in time. Similarly, when m1 = m2 = 1,m3 = 3, the determinant of the resulting matrix is given by |M | ≈ (4369.6 + 18918.6 s1/8 + 8739.3 s2/8 + 2854.6 s3/8 + 1809.9 s4/8)s27/8. This determinant is not zero for any s satisfying the requirements and the corresponding problem is well-posed. Some of the mixed order problems are more complicated. For example, 102 Chapter 6. Extension to a Class of High Order Problems when m1 = 1,m2 = 2,m3 = 3, the determinant of M is given by |M | ≈ (1567.5− 7568.4 s2/24 − 11936.6 s3/24 − 10703.4 s4/24 − 16384.0 s5/24 − 20488.5 s6/24 − 33059.1 s7/24 − 51680.4 s8/24 − 49546.0 s9/24 − 42813.5 s10/24 − 39554.5 s11/24 − 25840.2 s12/24 − 11585.2 s13/24 − 7568.4 s14/24 − 6992.3 s15/24 − 3134.9 s16/24)s29/8. This determinant can be zero for some s close to zero noticing the first positive term. However, it is bounded away from zero if s satisfies Re(s) > 1 and condition (6.2.4). Condition (6.2.4) guarantees that any power of s has positive real part provided the power is less than one. Thus this problem is still well-posed but with possible exponential growth in the solution. Similar results are obtained for other mixed order examples we tested. They are all well-posed, possibly with exponential growth. 6.3 Numerical Simulation In this section we consider several numerical simulations for the high order problems for both the equal order case and the mixed order case. The problem of the grain boundary motion discussed in section 4.4.3 gives the simplest example of equal order case. Here we consider a similar problem but with higher order derivatives. In this problem, all curves are enclosed in a unit circle (see Figure 6.3) and the normal velocity of each curve is equal to κssss. The PDAE system that describes this problem is given by X it · ~n = κssss i = 1, 2, 3 X iσ ·X iσσ = 0 i = 1, 2, 3. (6.3.1) This is a sixth order problem with m = 2 and the twelve junction conditions are X1(0) = X2(0) = X3(0) X1s ·X2s = cos θ3 X1s ·X3s = cos θ2 κ1 + κ2 + κ3 = 0 (6.3.2) κ1s = κ 2 s = κ 3 s κ1ss + κ 2 ss + κ 3 ss = 0 κ1sss = κ 2 sss = κ 3 sss 103 Chapter 6. Extension to a Class of High Order Problems where θ2 = θ3 = 2pi/3. At the domain boundary, the following conditions are applied to each curve: |X i(1)| = 1 X is · T (X i(1)) = 0 κiss = 0 κisss = 0. (6.3.3) The first two conditions indicate that the end of each curve is always attached perpendicularly to the circle. The term T (X(1)) represents the tangential direction of the unit circle at point X(1). X 1 2 X 1 1 X 1 0 X 1 −1 (κ1 −1 ) X 3 2 X 3 1 X 3 0 X 3 −1 (κ3 −1 ) X 2 2 X 2 0 X 2 −1 (κ2 −1 ) Junction X 2 1 Figure 6.2: Sketch of the ghost values at the triple junction for the discretiza- tion of the sixth order problem with m = 2. Two ghost points (X0 and X−1) and one ghost value (κ−1) are needed for each curve. This system is discretized using standard finite difference schemes on a staggered grid. Note that we need two ghost points and one ghost value of curvature (five ghost values in total) for each curve (see Figure 6.2). The numerical results are shown in Figure 6.3. 104 Chapter 6. Extension to a Class of High Order Problems −1 −0.5 0 0.5 1 −0.5 0 0.5 x y time = 0 −1 −0.5 0 0.5 1 −0.5 0 0.5 x y time = 0.00031 −1 −0.5 0 0.5 1 −0.5 0 0.5 x y time = 0.002 −1 −0.5 0 0.5 1 −0.5 0 0.5 x y time = 0.01 Figure 6.3: Simulation of equal order problem with m = 2 (sixth order). Angles between curves at the triple junction are 2pi/3. The second numerical example is a mixed order problem with m1 = 0,m2 = m3 = 2 which is given by X1t · ~n = κ X it · ~n = κssss i = 2, 3 (6.3.4) X iσ ·X iσσ = 0 i = 1, 2, 3. 105 Chapter 6. Extension to a Class of High Order Problems This system needs ten junction conditions which are given by X1(0) = X2(0) = X3(0) X1s ·X2s = cos θ3 X1s ·X3s = cos θ2 κ2 + κ3 = 0 (6.3.5) κ2s = κ 3 s κ2ss + κ 3 ss = 0 κ2sss = κ 3 sss where θ2 = θ3 = (pi/2 + arcsin 1 4 ). Inspired by the coupled surface and grain boundary motion which admits travelling wave solutions, we expect a travelling wave solution for this new problem. Therefore, we start the simulation with the travelling wave solution for the coupled problem given in Chapter 2 and hope that the travelling wave solution for the new problem, if exist, will stay close. To keep the curves flat, the conditions at the far field domain boundary are given by X i(1, t) = given value i = 1, 2, 3 κis = 0 i = 2, 3 κiss = 0 i = 2, 3 κisss = 0 i = 2, 3. (6.3.6) The discretization of the junction and boundary conditions again needs some ghost points and ghost values which can derive analogously with pre- vious discussions. The numerical results are shown in Figure 6.4. The shape of the curves changes at the beginning but it stabilizes after a long enough time and moves to the right at a constant speed. The speed of the triple junction versus time is shown in Figure 6.5. It clearly shows the constant speed which is around 1.78. This result confirms our expectation of the ex- istence of travelling wave solutions for mixed high order problems. However, we have not been able to find them analytically. We remark that all the motions discussed in this chapter that are higher than order four are area preserving. Therefore, no travelling wave solution is possible if mi ≥ 1 for all three curves. This point can be shown by the next example. 106 Chapter 6. Extension to a Class of High Order Problems −20 −10 0 10 20 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 Result at time 0 −10 0 10 20 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 Result at time 5 20 30 40 50 60 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 Result at time 25 50 60 70 80 90 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 Result at time 40 Figure 6.4: Simulation of mixed order problem with m1 = 0,m2 = m3 = 2. 5 10 15 20 25 1.55 1.6 1.65 1.7 1.75 1.8 1.85 1.9 1.95 Speed of triple junction time sp ee d Figure 6.5: Plot of the junction speed versus time for the mixed order problem with m1 = 0,m2 = m3 = 2. The junction speed converges approximately to 1.78. 107 Chapter 6. Extension to a Class of High Order Problems We consider another mixed order problem with m1 = 1,m2 = m3 = 2. This problem can be described by the following PDAE system: X1t · ~n = −κss X it · ~n = κssss i = 2, 3 (6.3.7) X iσ ·X iσσ = 0 i = 1, 2, 3. The eleven junction conditions are given by X1(0) = X2(0) = X3(0) X1s ·X2s = cos θ3 X1s ·X3s = cos θ2 κ1 + κ2 + κ3 = 0 (6.3.8) κ1s = κ 2 s = κ 3 s κ2ss + κ 3 ss = 0 κ2sss = κ 3 sss where θ2 = θ3 = (pi/2 + arcsin 1 4 ). We again start the simulation with the travelling wave solution of the coupled surface and grain boundary motion. The conditions at the domain boundary are given by X i(1, t) = given value i = 1, 2, 3 κis = 0 i = 1, 2, 3 κiss = 0 i = 2, 3 κisss = 0 i = 2, 3. (6.3.9) The discretization of this system is similar to previous examples. The nu- merical results are shown in Figure 6.6. We are not expecting any travelling wave solutions for this problem because of the reason stated above. 108 Chapter 6. Extension to a Class of High Order Problems −20 −10 0 10 20 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 −20 −10 0 10 20 −1.2 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 −20 −10 0 10 20 −1.5 −1 −0.5 0 0.5 −20 −10 0 10 20 −1.5 −1 −0.5 0 0.5 Figure 6.6: Simulation of mixed order problem with m1 = 1,m2 = m3 = 2. 6.4 Conclusions In this chapter, we considered a class of high order three-phase boundary motion problems. The PDAE formulation is extended easily to describe these new problems. The extension of other formulations or methods, such as the parabolic formulation, the finite element methods ([5, 6, 57, 58]) and the MBO methods ([22, 37, 36, 38, 50]) is difficult, or impossible. All equal order problems are shown to be well-posed using the same technique discussed in Chapter 3. Some examples of mixed order problems are also shown to be well-posed and some of them may admit travelling wave solutions. 109 Chapter 7 Thesis Summary and Future Work 7.1 Thesis Summary In this thesis, we studied the coupled surface and grain boundary motion in the context of a quarter loop geometry and other related problems. Differ- ent formulations were presented and compared in Chapter 2. The cartesian formulation is the simplest but cannot handle the non-single-valued case. The parabolic formulation is given in a parametric form and it is convenient for analysis because of its parabolic character in all components. However, its numerical performance is not as good as the PDAE formulation. The PDAE formulation introduced in this thesis uses an algebraic equation to preserve the arc length up to scaling and outstanding numerical performance was achieved as shown in Chapter 4. It was applied to several problems of physical interest, giving more accurate results than other methods did. It maintains uniform resolution and allows for convenient global regridding. This could be necessary as curves in a network expand or shrink quickly. The PDAE formulation can also be used to efficiently simulate the motion of a curve that involves only mean curvature motion or surface diffusion as shown in Chapter 4. In Chapter 3, the well-posedness of the parabolic formulation and the PDAE formulation were analyzed in a reduced, linear setting. Both formu- lations are well-posed with appropriate angle conditions at the triple junc- tion. The parabolic formulation requires extra artificial tangential condi- tions. It is shown that these conditions are arbitrary and will affect only the parametrization of the solution. The PDAE formulation includes algebraic equations in the system and thus may introduce extra difficulty for analy- sis. However, the approach we used can be applied to both formulations and our analysis gives the conditions that match those that more complicated nonlinear analysis gives in the cases studied in [9, 25]. 110 Chapter 7. Thesis Summary and Future Work The stability of travelling waves to 2D and 3D perturbations were also discussed in this thesis. The basic approach is to linearize the system around the travelling waves and approximate the spectrum of the resulting linear system numerically. The numerical results gave evidence that both problems are convectively stable. However, the numerical computation has to be done on a finite domain and there is still a theoretical gap between the numerical eigenvalues on a finite domain and the spectrum on an infinite domain when triple junctions and algebraic constraints are present. We suspect that rele- vant theories can also be applied to our problem, but it is still an interesting open question. We also introduced a class of high order three-phase boundary motion problems in Chapter 6. The PDAE formulation was extended naturally to the whole class of problems which can then be approximated by finite difference methods. This again shows the superiority of the PDAE formulation. The problems were classified into two categories and their well-posedness were considered separately. For the equal order case, i.e., when all three curves move according to the same evolution law, all the problems were shown to be well-posed. When more than one motion are involved in the system, we investigated their well-poseness case by case. All the cases that we considered are well-posed. 7.2 Future Work The well-posedness analysis described in Chapter 3 is based on the lineariza- tion of the system around straight line solutions. Rough analysis of the reason to investigate only at the junction is given in this thesis but there is still room to make the investigation more rigorous. In Chapter 5, we have shown some interesting numerical results on the stability of the travelling waves. These results gave some predictions re- garding the spectral theory for problems with triple junctions. The existing spectral theory may be extended to such problems. This is a much more challenging problem. For the higher order three-phase problems discussed in Chapter 6, we have shown the possible existence of travelling wave solutions for some mixed order problems. It would be interesting to show the existence of such solutions analytically. Finally, regarding the numerical methods for curve evolution problems, 111 Chapter 7. Thesis Summary and Future Work the threshold dynamics algorithms (MBO algorithms) discussed in [22, 37, 36, 38, 50] seem to be interesting. These methods exhibit excellent stability properties and can also be used for multiphase problems involving topolog- ical changes. They have been used for mean curvature motions with triple junctions ([50]) and surface diffusion problems with triple junctions ([22]). However, both cases involve only one type of motion and it is still a difficult question on how to implement these techniques for mixed order problems. Even for the equal order problems there are unanswered questions. For exam- ple, for the motion by surface diffusion discussed in [22], it is not clear what conditions are implemented at the triple junction. It would be interesting to investigate these problems. 112 Bibliography [1] I. Ahmed, H.-S. Chu, W. B. Ewe, E.-P. Li, and Z. Chen. Modeling and simulation of nano-interconnects for nanophotonics. Electronics Pack- aging Technology Conference, 2007. EPTC 2007. 9th, pages 436–440, Dec. 2007. [2] J. Alexander, R. Gardner, and C. Jones. A topological invariant arising in the stability of travelling waves. Journal für die Reine und Ange- wandte Mathematik, 410:167–212, 1990. [3] U. M. Ascher and L. R. Petzold. Computer Methods for Ordinary Dif- ferential Equations and Differential-Algebraic Equations. SIAM, 1998. [4] E. Bansch, P. Morin, and R. H. Nochetto. A finite element method for surface diffusion: the parametric case. Journal of Computational Physics, 203(1):321–343, Feb. 2005. [5] J. W. Barrett, H. Garcke, and R. Nürnberg. On the variational ap- proximation of combined second and fourth order geometric evolution equations. SIAM Journal on Scientific Computing, 29(3):1006–1041, 2007. [6] J. W. Barrett, H. Garcke, and R. Nürnberg. A parametric finite ele- ment method for fourth order geometric evolution equations. Journal of Computational Physics, 222(1):441–467, 2007. [7] W.-J. Beyn and J. Lorenz. Stability of traveling waves: dichotomies and eigenvalue conditions on finite intervals. Numerical Functional Analysis and Optimization, 20(3):201–244, 1999. [8] A. Brokman, R. Kris, W. W. Mullins, and A. J. Vilenkin. Analysis of boundary motion in thin films. Scripta Metallurgica et Materialia, 32(9):1341–1346, May 1995. 113 Bibliography [9] L. Bronsard and F. Reitich. On three-phase boundary motion and the singular limit of a vector-valued Ginzburg-Landau equation. Archive for Rational Mechanics and Analysis, 124(4):355–379, Dec. 1993. [10] L. Bronsard and B. T. R. Wetton. A numerical method for tracking curve networks moving with curvature motion. Journal of Computa- tional Physics, 120(1):66–87, Aug. 1995. [11] M. Burger. Surface diffusion including adatoms. Communications in Mathematical Sciences, 4(1):1–51, 2006. [12] J. Cahn and J. Taylor. Overview 113: Surface motion by surface diffu- sion. Acta Metallurgica et Materialia, 42:1045–1063, 1994. [13] P. P. Camanho, C. G. Davila, and M. F. de Moura. Numerical simulation of mixed-mode progressive delamination in composite materials. Journal of Composite Materials, 37(16):1415–1438, 2003. [14] S. J. Chapman, S. D. Howison, and J. R. Ockendon. Macroscopic models for superconductivity. SIAM Review, 34(4):529–560, 1992. [15] D. L. Chopp and J. Sethian. Motion by intrinsic Laplacian of curvature. Interfaces and Free Boundaries, 1:1–18, 1999. [16] U. Czubayko, V. G. Sursaeva, G. Gottstein, and L. S. Shvindlerman. Influence of triple junctions on grain boundary motion. Acta Materialia, 46(16):5863–5871, Oct. 1998. [17] K. Deckelnick, G. Dziuk, and C. M. Elliott. Computation of geometric partial differential equations and mean curvature flow. Acta Numerica, 14:139–232, 2005. [18] Q. Du. Numerical approximations of the Ginzburg–Landau models for superconductivity. Journal of Mathematical Physics, 46(9):095109, 2005. [19] Q. Du, M. D. Gunzburger, and J. S. Peterson. Analysis and approxima- tion of the Ginzburg-Landau model of superconductivity. SIAM Review, 34(1):54–81, 1992. [20] C. Dunn, F. Daniels, and M. Bolton. Transactions of the American Institute of Mining, Metallurgical and Petroleum Engineers, 185:708, 1949. 114 Bibliography [21] G. Dziuk, E. Kuwert, and R. Schatzle. Evolution of elastic curves in Rn: Existence and computation. SIAM Journal on Mathematical Analysis, 33(5):1228–1245, 2002. [22] S. Esedoglu, S. Ruuth, and R. Tsai. Threshold dynamics for high order geometric motions. Interfaces and Free Boundaries, to appear. [23] M. Furtkamp, G. Gottstein, D. A. Molodov, V. N. Semenov, and L. S. Shvindlerman. Grain boundary migration in Fe-3.5% Si bicrystals with [001] tilt boundaries. Acta Materialia, 46(12):4103–4110, July 1998. [24] A. V. Galina, V. E. Fradkov, and L. S. Shvindlerman. Influence of mobility of triple grain junctions on boundary migration. Physics Metals Metallogr, 63(6):165, 1987. [25] H. Garcke and A. Novick-Cohen. A singular limit for a system of degen- erate Cahn-Hilliard equations. Advanced Differential Equations, 5:401– 434, 2000. [26] G. Gottstein and L. Shvindlerman. Grain boundary migration in metals: thermodynamics, kinetics, applications. CRC Press, 1999. [27] W. Hergert, A. Ernst, and M. Dane. Computational Materials Science: From Basic Principles to Material Properties. Springer, 2004. [28] T. Hou, J. Lowengrub, and M. Shelley. Boundary integral methods for multicomponent fluids and multiphase materials. Journal of Computa- tional Physics, 169:302–362, 2001. [29] J. Kanel, A. Novick-Cohen, and A. Vilenkin. A traveling wave solu- tion for coupled surface and grain boundary motion. Acta Materialia, 51(7):1981–1989, Apr. 2003. [30] J. Kanel, A. Novick-Cohen, and A. Vilenkin. Coupled surface and grain boundary motion: a travelling wave solution. Nonlinear Analy- sis, 59(8):1267–1292, Dec. 2004. [31] J. Kanel, A. Novick-Cohen, and A. Vilenkin. Coupled surface and grain boundary motion: nonclassical traveling wave solution. Advanced Dif- ferential Equations, 9:299–327, 2004. 115 Bibliography [32] J. Kanel, A. Novick-Cohen, and A. Vilenkin. A numerical study of grain boundary motion in bicrystals. Acta Materialia, 53(2):227–235, Jan. 2005. [33] B. B. Karki, R. M. Wentzcovitch, S. de Gironcoli, and S. Baroni. First- principles determination of elastic anisotropy and wave velocities of MgO at lower mantle conditions. Science, 286(5445):1705–1707, 1999. [34] M. Khenner, A. Averbuch, M. Israeli, and M. Nathan. Numerical sim- ulation of grain-boundary grooving by level set method. Journal of Computational Physics, 170(2):764–784, July 2001. [35] U. F. Mayer. Numerical solutions for the surface diffusion flow in three space dimensions. Computational and Applied Mathematics, 20:361–379, 2001. [36] B. Merriman, J. K. Bence, and S. J. Osher. Diffusion generated motion by mean curvature. in Computational Crystal Growers Workshop, J. Taylor, ed., American Mathematical Society, Providence, Rhode Island, pages 73–83, 1992. [37] B. Merriman, J. K. Bence, and S. J. Osher. Motion of multiple junctions: a level set approach. Journal of Computational Physics, 112(2):334–363, 1994. [38] B. Merriman and S. J. Ruuth. Diffusion generated motion of curves on surfaces. Journal of Computational Physics, 225(2):2267–2282, 2007. [39] K. Mikula and D. Sevcovic. Computational and qualitative aspects of evolution of curves driven by curvature and external force. Computing and Visualization in Science, 6(4):211–225, 2004. [40] D. Min and H. Wong. A model of migrating grain-boundary grooves with application to two mobility-measurement methods. Acta Materialia, 50(20):5155–5169, Dec. 2002. [41] W. Mullins. Two-dimensional motion of idealized grain boundaries. Journal of Applied Physics, 27:900–904, 1956. [42] W. Mullins. Theory of thermal grooving. Journal of Applied Physics, 28:333–339, 1957. 116 Bibliography [43] W. Mullins. The effect of thermal grooving on grain boundary motion. Acta Metallurgica, 6:414–427, 1958. [44] S. Osher and R. Fedkiw. Level Set Methods and Dynamic Implicit Sur- faces, volume 153 of Applied Mathematical Sciences. Springer, 2003. [45] S. Osher and N. Paragios. Geometric Level Set Methods in Imaging, Vision, and Graphics. Springer, 2003. [46] S. Osher and J. A. Sethian. Fronts propagating with curvature- dependent speed: algorithms based on hamilton-jacobi formulations. Journal of Computational Physics, 79(1):12–49, 1988. [47] Z. Pan and B. T. Wetton. A numerical method for coupled surface and grain boundary motion. European Journal of Applied Mathematics, 19:311–327, 2008. [48] R. L. Pego and M. I. Weinstein. Eigenvalues, and instabilities of solitary waves. Philosophical Transactions: Physical Sciences and Engineering, 340(1656):47–94, 1992. [49] S. Ruuth. Efficient Algorithms for Diffusion-Generated Motion by Mean Curvature. PhD thesis, University of British Columbia, 1996. [50] S. J. Ruuth. Efficient algorithms for diffusion-generated motion by mean curvature. Journal of Computational Physics, 144(2):603–625, 1998. [51] B. Sandstede. Stability of travelling waves. Handbook of Dynamical Systems II (B Fiedler, ed.). North-Holland, pages 983–1055, 2002. [52] B. Sandstede and A. Scheel. Absolute and convective instabilities of waves on unbounded and large bounded domains. Physica D: Nonlinear Phenomena, 145(3-4):233–277, Nov. 2000. [53] J. Sethian. Level set methods and fast marching methods : Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vi- sion, and Materials Science, volume 3 of Cambridge Monographs on Applied and Computational Mathematics. Cambridge University Press, 1999. 117 Bibliography [54] P. Smereka. Semi-implicit level set methods for curvature and surface diffusion motion. Journal of Scientific Computing, 19(1):439–456, Dec. 2003. [55] V. Solonnikov. Boundary value problems in physics. Proceedings of the Steklov Institute of Mathematics, LXXXIII, 1965. [56] J. M. Stockie and B. T. R. Wetton. Stability analysis for the immersed fiber problem. SIAM Journal on Applied Mathematics, 55(6):1577–1591, 1995. [57] B. Sun and Z. Suo. A finite element method for simulating interface motion–II. large shape change due to surface diffusion. Acta Materialia, 45(12):4953–4962, Dec. 1997. [58] B. Sun, Z. Suo, and W. Yang. A finite element method for simulat- ing interface motion–I. migration of phase and grain boundaries. Acta Materialia, 45(5):1907–1915, May 1997. [59] R. Sun and C. Bauer. Measurement of grain boundary mobilities through magnification of capillary forces. Acta Metallurgica, 18:635–638, 1970. [60] R. Sun and C. Bauer. Tilt boundary migration in NaCl bicrystals. Acta Metallurgica, 18:639–647, 1970. [61] Z. Suo. Motions of microscopic surfaces in materials. Advances in Ap- plied Mechanics, 33:193–294, 1997. [62] A. J. Vilenkin, R. Kris, and A. Brokman. Breakup and grain growth in thin-film array. Journal of Applied Physics, 81(1):238–245, 1997. [63] T. Wang, N. Moll, K. Cho, and J. D. Joannopoulos. Deliberately de- signed materials for optoelectronics applications. Physical Review Let- ters, 82(16):3304–3307, Apr 1999. [64] I. Webman, J. Jortner, and M. H. Cohen. Numerical simulation of elec- trical conductivity in microscopically inhomogeneous materials. Physical Review B, 11(8):2885–2892, Apr 1975. [65] H. Zhang and H. Wong. Coupled grooving and migration of inclined grain boundaries: Regime I. Acta Materialia, 50(8):1983–1994, May 2002. 118 Bibliography [66] H. Zhang and H. Wong. Coupled grooving and migration of inclined grain boundaries: Regime II. Acta Materialia, 50(8):1995–2012, May 2002. [67] W. Zhang and I. Gladwell. Thermal grain boundary grooving with anisotropic surface free energy in three dimensions. Journal of Crys- tal Growth, 277(1-4):608–622, Apr. 2005. [68] W. Zhang and J. H. Schneibel. Numerical simulation of grain- boundary grooving by surface diffusion. Computational Materials Sci- ence, 3(3):347–358, Jan. 1995. [69] X. Zhang, J.-S. Chen, and S. Osher. A multiple level set method for modeling grain boundary evolution of polycrystalline materials. Subm- mited. [70] H.-K. Zhao, T. Chan, B. Merriman, and S. Osher. A variational level set approach to multiphase motion. Journal of Computational Physics, 127(1):179–195, Aug. 1996. 119
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Simulation and analysis of coupled surface and grain...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Simulation and analysis of coupled surface and grain boundary motion Pan, Zhenguo 2008
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Simulation and analysis of coupled surface and grain boundary motion |
Creator |
Pan, Zhenguo |
Publisher | University of British Columbia |
Date Issued | 2008 |
Description | At the microscopic level, many materials are made of smaller and randomly oriented grains. These grains are separated by grain boundaries which tend to decrease the electrical and thermal conductivity of the material. The motion of grain boundaries is an important phenomenon controlling the grain growth in materials processing and synthesis. Mathematical modeling and simulation is a powerful tool for studying the motion of grain boundaries. The research reported in this thesis is focused on the numerical simulation and analysis of a coupled surface and grain boundary motion which models the evolution of grain boundary and the diffusion of the free surface during the process of grain growth. The “quarter loop” geometry provides a convenient model for the study of this coupled motion. Two types of normal curve velocities are involved in this model: motion by mean curvature and motion by surface diffusion. They are coupled together at a triple junction. A front tracking method is used to simulate the migration. To describe the problem, different formulations are presented and discussed. A new formulation that comprises partial differential equations and algebraic equations is proposed. It preserves arc length parametrization up to scaling and exhibits good numerical performance. This formulation is shown to be well-posed in a reduced, linear setting. Numerical simulations are implemented and compared for all formulations. The new formulation is also applied to some other related problems. We investigate numerically the linear stability of the travelling wave solutions for the quarter loop problem and a simple grain boundary motion problem for both curves in two dimensions and surfaces in three dimensions. The numerical results give evidence that they are convectively stable. A class of high order three-phase boundary motion problems are also studied. We consider a region where three phase boundaries meet at a triple junction and evolve with specified normal velocities. A system of partial differential algebraic equations (PDAE) is proposed to describe this class of problems by extending the discussion for the coupled surface and grain boundary motion. The linear well-posedness of the system is analyzed and numerical simulations are performed. |
Extent | 2043312 bytes |
Subject |
Numerical simulation Grain boundary motion Surface diffusion Linear stability |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2008-10-30 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0066768 |
URI | http://hdl.handle.net/2429/2733 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2009-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2009_spring_pan_zhenguo.pdf [ 1.95MB ]
- Metadata
- JSON: 24-1.0066768.json
- JSON-LD: 24-1.0066768-ld.json
- RDF/XML (Pretty): 24-1.0066768-rdf.xml
- RDF/JSON: 24-1.0066768-rdf.json
- Turtle: 24-1.0066768-turtle.txt
- N-Triples: 24-1.0066768-rdf-ntriples.txt
- Original Record: 24-1.0066768-source.json
- Full Text
- 24-1.0066768-fulltext.txt
- Citation
- 24-1.0066768.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0066768/manifest