Application of the Level Set Method to Solid Rocket Motor Simulation by Kevin M. Albarado A thesis submitted to the Graduate Faculty of Auburn University in partial ful llment of the requirements for the Degree of Master of Science Auburn, Alabama August 4, 2012 Keywords: rocket propulsion, solid rocket, level set, discontinous Galerkin Copyright 2012 by Kevin M. Albarado Approved by Roy Hart eld, Chair,Woltosz Professor of Aerospace Engineering John Burkhalter, Professor Emeritus of Aerospace Engineering Andrew Shelton, Assistant Professor of Aerospace Engineering George Flowers, Professor of Mechanical Engineering, Dean of Graduate School Abstract The body of work encompassed in this thesis merges two advanced concepts for develop- ing ow solutions with level sets to develop an accurate and e cient method for simulating solid rocket motor grain regression. The level-set method has been implemented using the discontinous Galerkin numerical method (DGM) to represent the burning surface area and chamber volume as a function of time. The combination of DGM with geometrically arbi- trary solid motors presents a unique and novel environment for solid rocket motor analysis, giving the user more freedom for design and a more accurate result. This thesis provides a complete background and introduction to both solid rocket motor research completed and currently underway as well as an introduction to computational uid dynamics and the level set method. Development of the LSM/DGM approach to the grain regression problem in two-dimensions is presented with complementary comparisons to analytical approaches. This work represents the rst known implementation of the discontinous Galerkin method with a level set to solve the grain regression problem. ii Acknowledgments I would rst like to thank my committee for the opportunity and encouragement to pursue this work. I would like thank Dr. Hart eld rst and foremost for his knowledge and expertise in the rocket propulsion eld. He provided the foundation from which this work developed. Dr. Shelton provided a unique perspective on how to approach this work, and also helped tremendously in implementation of the approach developed in this thesis. Without his input, the end product of this thesis would not have been realized. Finally, I thank Dr. Burkhalter for advice on life decisions and continuous support and encouragement. Without an open door policy from all three members of my committee, this work would not have been completed. I would also like to thank my family and friends for their continued support and encouragement as I make the transition from academic life to the \real world". iii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Performance Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 Ignition Phenomena . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.3 Chamber Aerodynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.3.1 Analytical Approaches . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.3.2 Full Scale CFD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.4 Motor Stability Considerations . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.5 Propellant Solid Mechanics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.6 Automated Grain Design and Performance Control . . . . . . . . . . . . . . 8 2 Existing Internal Ballistics Models . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.1 Analytical Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2 Face-O setting Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 3 Eulerian Approaches to Surface Tracking . . . . . . . . . . . . . . . . . . . . . . 13 4 First Order Grain Regression Program . . . . . . . . . . . . . . . . . . . . . . . 19 4.1 Development and Implementation . . . . . . . . . . . . . . . . . . . . . . . . 19 4.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 5 Higher Order Grain Regression Program . . . . . . . . . . . . . . . . . . . . . . 28 5.1 Development and Implementation . . . . . . . . . . . . . . . . . . . . . . . . 28 5.1.1 Discontinuous Galerkin Method . . . . . . . . . . . . . . . . . . . . . 28 5.1.2 Stabilizing the Di erential Equation . . . . . . . . . . . . . . . . . . . 31 iv 5.1.3 Motor Case Implementation . . . . . . . . . . . . . . . . . . . . . . . 41 5.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 6 Capabilities and Further Development . . . . . . . . . . . . . . . . . . . . . . . 52 6.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 6.2 Recommendations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 A Analytical method for a star grain . . . . . . . . . . . . . . . . . . . . . . . . . 60 B Derivation for Discontinuous Galerkin Approach . . . . . . . . . . . . . . . . . . 63 C Comparison of First-Order Finite Element to Discontinuous Galerkin:A Classic Fluid Dynamics Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 D DG/LSM Fortran Source . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 v List of Figures 1.1 Generic SRM Schematic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Typical pressure-time pro le for a solid rocket motor . . . . . . . . . . . . . . . 6 2.1 Star Grain De nition [27] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2 Face-O setting Method: Advective Reconstruction vs. Wavefront Reconstruction 12 3.1 A Closed Propagating Curve . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 3.2 Convex and Concave curvature with Lagrangian propagation (left) and actual propagation (right) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.3 Volume of Fluid Representation . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.4 A Discretized Curve and its Complementary Level Set . . . . . . . . . . . . . . 17 3.5 Winding Number Diagram . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 4.1 Smearing functions for arbitrary near the zeroeth level set . . . . . . . . . . . 22 4.2 A Neutral Burning Seven Pointed Star . . . . . . . . . . . . . . . . . . . . . . . 23 4.3 Neutral Burning Star Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 4.4 Smeared Heaviside function for 7-pointed star at ignition . . . . . . . . . . . . . 24 4.5 Motor Surface near the Case Boundary . . . . . . . . . . . . . . . . . . . . . . . 25 vi 4.6 Grid Sensitivity Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 4.7 O center Circular Perforated Grain Regression . . . . . . . . . . . . . . . . . . 27 4.8 Pressure Time Pro le for an O -Center Cylindrical Perforated Grain . . . . . . 27 5.1 Stencil Used to Develop DG Derivatives . . . . . . . . . . . . . . . . . . . . . . 31 5.2 Level set for a circle of radius = 0.25 . . . . . . . . . . . . . . . . . . . . . . . . 32 5.3 Level set propagation for a circle of radius = 0.25 . . . . . . . . . . . . . . . . . 33 5.4 Level set propagation for a circle of radius = 0.25 with source term . . . . . . . 34 5.5 Five Pointed Star Level Set . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 5.6 Normalized residual for ve-pointed star using mesh of 10 10 with 9th order polynomial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 5.7 Normalized residuals for ve-pointed star for varying dissipation constants . . . 37 5.8 Level set propagation for a ve pointed star with dissipation constant set to 0.006 38 5.9 Level set propagation for a ve pointed star with dissipation constant set to 0.009 39 5.10 Level set propagation for a ve pointed star at t = 0:2 . . . . . . . . . . . . . . 40 5.11 Motor Surface near the Case Boundary . . . . . . . . . . . . . . . . . . . . . . . 41 5.12 Heaviside function for motor case . . . . . . . . . . . . . . . . . . . . . . . . . . 42 5.13 Five pointed star burning into the case wall . . . . . . . . . . . . . . . . . . . . 43 5.14 Comparison of polynomial order with error . . . . . . . . . . . . . . . . . . . . . 45 vii 5.15 Burn area versus burn distance comparison between analytical method and DG schemes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 5.16 Comparison of polynomial order with error compiled using analytical area . . . 46 5.17 Five pointed star with grain cross with grain regression lines . . . . . . . . . . . 47 5.18 Comparison of ve pointed star simulated with DG versus analytical solution . 48 5.19 Propellant remaining for ve pointed star at various burn distances for a circular motor case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 5.20 Propellant remaining for ve pointed star at various burn distances for a square motor case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 5.21 Comparison of pressure pro les for star grain in a circular case and a square case 51 6.1 Example utilization of sparse matrices . . . . . . . . . . . . . . . . . . . . . . . 54 A.1 Star Grain De nition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 B.1 Stencil for 1-D example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 B.2 Modal versus Nodal Graphical Representation . . . . . . . . . . . . . . . . . . . 66 B.3 Lobatto quadrature for varying degrees of freedom . . . . . . . . . . . . . . . . 69 C.1 Euler Solution to Isentropic Vortex using Finite Element . . . . . . . . . . . . . 73 C.2 Euler Solution to Isentropic Vortex using Discontinuous Galerkin . . . . . . . . 74 viii Chapter 1 Introduction The roots of solid rocketry can be dated back to the 13th century, when the Chinese rst used gunpowder as a solid fuel propellant in warfare. The Arabs, Indians, and Mongolese shortly followed the Chinese in developing some of the earliest rocket weaponry. The earliest rockets were simple bamboo rods lled with gunpowder and carried incendiary material and shrapnel. While history is lled with examples of rocketry use in warfare, solid rocket motor design would go relatively unchanged for hundreds of years. During the American Revolution, William Congreve developed rockets for use by the British military. The advancement made by Congreve in these rockets were the cone-shaped cavity in the end of the propellant. This cavity provided a higher burn area resulting in higher pressure and mass ow rate, in turn producing much higher thrust levels. This was a signi cant advancement leading engineers to begin experimenting with grain modi cation for performance enhancement. In the late 1800?s, solid rockets saw improvement in propellants as Paul Vieille and Alfred Nobel began experimenting with explosive components such as nitroglycerine. In the early 1900?s, most of the solid rocket motors produced were used in rockets for military use due to the long shelf life of the propellant. Improvements in solid rocket motors were restricted to mainly propellant development, led by John Parsons, Frank Malina, and Theodore von Karman until the 1960?s. Parsons, Malina, and von Karman made signi cant advances in development of composite solid fuels that enhanced the stability of the motors along with advancing performance. During the 1960?s, the idea of a Space Shuttle had developed, and the design of the two largest solid rockets in history began. Due to the unprecedented size of these motors, known as the Solid Rocket Boosters (SRB), numerous advancements in solid rocket science were required. During initial development of the SRBs, a simpli ed 1 treatment of the grain regression problem was given to the analysis. The programs used in this analysis, while accurate, were somewhat tailored speci cally to the exact problem at hand. These same tools built during the development of the Shuttle booster program are still in use today in industrial settings, such as Solid Performance Program [1]. SPP still remains the industry standard solid rocket motor analysis tool. A solid rocket motor is simple enough in design and implementation (which is what makes it advantageous over more complicated rocket motors) as it is only composed of four main components: propellant, igniter, case, and nozzle (see Figure 1.1). Igniter Propellant Motor Case Nozzle Figure 1.1: Generic SRM Schematic Despite its simplicity, many physical phenomena occur during a solid motor ring. From ignition to burnout, engineers and scientists have made entire careers out of studying just a handful of the myriad of topics present in solid rocket motor research such as: per- formance estimation, ignition phenomena, grain regression simulation techniques, chamber aerodynamics, motor stability considerations, propellant solid mechanics, propellant manu- facturing, and automated grain design amongst many more. While this body of work focuses solely on grain regression simulation, a brief discussion of the remaining topics is necessary to fully understand the underlying physics in the unique environment of a solid rocket motor. 2 1.1 Performance Estimation In order to analyze a solid rocket motor, some estimation of performance must be available. For a rocket motor, the performance parameters of interest are chamber pressure and thrust. A simple derivation for chamber pressure and thrust can be developed using lumped parameter analysis. The general thrust equation for any chemical propulsion system is a manipulation of the momentum equation assuming steady, uniform ow at the intake and exit planes of the motor nozzle: T = _m(ue ui) +Ae (pe pa) (1.1) For a rocket motor, there is no incoming jet of momentum, so this equation simpli es to just T = _mue +Ae (pe pa) (1.2) Analyzing Equation 1.2, it should be noted that the unknowns in this equation are _m, ue, and pe, and as such are the variables to estimate with a rocket analysis program. These unknowns are simply functions of burn area and propellant properties. The continuity equation requires that the mass production rate of the propellant grain at an time must equal the mass egress rate of the nozzle. This is represented mathematically as Abrb b = poAtc (1.3) Equation 1.3 assumes uniform, steady, and 1-D ow throughout both the combustion cham- ber and the nozzle. For 2- and 3-D ows, mass weighted averages apply. The burn rate, rb, is empirically de ned to be a function of propellant properties rb = apno (1.4) 3 Substituting Equation 1.4 into 1.3 and rearranging yields an equation for the total chamber pressure. po = A b Ata bc 1 1 n (1.5) In this equation, it is assumed that the propellant burn rate constant, propellant density, characteristic velocity, and burn rate exponent (a, b, c , and n respectively) are knowns as well as the throat area, At. This implies that chamber pressure is function of burn area only for a given propellant choice. With total chamber pressure calculated, the mass ow rate is calculated as _m = poAtc (1.6) Using isentropic nozzle performance equations, the exit pressure can be calculated by solving the Area Ratio-Mach number relation for Mach number, and using isentropic relations to nd the pressure ratio. Ae At = 1 Me 2 + 1 1 + 12 M2e +1 2( 1) (1.7) po pe = 1 + 12 M2e 1 (1.8) Equations 1.2, 1.5, 1.6, 1.7, and 1.8 provide the necessary conditions to calculate thrust for a given burn area in a solid rocket motor. The primary work for this thesis is to determine accurately the burn area as a function of time for any generic solid motor. For a more in-depth derivation of the performance estimation see references [2, 3]. 1.2 Ignition Phenomena The ignition of a solid rocket motor has as much of an impact on overall performance as grain design and manufacturing. To start a solid rocket motor, an igniter introduces (usually at the head end) high temperature and high pressure gases or ames. These ames spread 4 throughout the motor gradually heating up the propellant surface. When the propellant reaches the ignition temperature, the propellant burns, releasing large amounts of thermal energy into the chamber. Typically at this point a phenomenon known as pressure overpeak occurs (see Figure 1.2). The whole process takes place in only a few hundred milliseconds for the largest motors such as the Space Shuttle booster and is much shorter in smaller motors, but is signi cant in computing the total impulse for the motor. This process has been studied extensively both experimentally and analytically. Refer- ence [4] lists numerous characteristics of igniters including e ects of temperature and pres- sure on ignition time as well as types of igniters and performance characteristics. Foster and Jenkins developed a computer model and performed validation experiments to determine whether a single port igniter or multiport igniter would be more e ective at igniting the head end star on the Space Shuttle Redesigned Solid Rocket Motor [5]. Their work was a collaborative e ort between Auburn University and NASA?s Marshall Space Flight Center. Cho and Baek [6], Bai et al. [7], and Johnston [8] all developed numerical simulations of the internal aerodynamics during ignition in order to analyze ignition transients in an axisym- metric motor. However, Cho and Baek took this work a step further by including radiation e ects to determine quantitatively what e ect radiation has on ignition. They found that radiation plays an important role in the heat ux to the propellant surface, and without radiation a longer ignition delay is required. 1.3 Chamber Aerodynamics Chamber aerodynamics can play a large roll in motor performance. While the local velocity gradients and temperature gradients are not large compared to conditions in the nozzle, they play a roll on local pressure as well as the local burn rate. If the burn rate gradient along the surface is large enough, some areas of the motor will burn more rapidly than others, leading to a divergence from predicted burn area. Aerodynamic predictions have been performed in both analytical approaches and full scale CFD. 5 Time Pr es sur e Pressure Overpeak Steady State Operation Tailoff Ignition Transient Figure 1.2: Typical pressure-time pro le for a solid rocket motor 1.3.1 Analytical Approaches Numerous assumptions can reasonably be made for solid rocket motor chamber aero- dynamics. The set of assumptions made speci cally can be referred to as the Taylor-Culick Pro le. Taylor and Culick independently determined the conditions and assumptions for solid rocket motors in the 1950?s and 60?s [9, 10]. The assumptions are that the ow solution is incompressible, rotational, axisymmetric, and quasi-viscous. Majdalani et al. [11, 12, 13, 14, 15] has produced a series of publications detailing the analytical Taylor- Culick pro le and its many applications. In [11], basic development of the Taylor-Culick pro le was examined to determine the time-dependent velocity eld for a ow- eld with sidewall injection. Majdalani and Saad examined how an arbitrary headwall injection af- fected the ow eld solution (this particular e ect would be observed in a hybrid rocket motor) [12]. Majdalani, Xu, Lin, and Liao used the Homotopy Analysis Method to exam- ine the Taylor-Culick pro le with regressing and injecting sidewalls [13]. This development proved to be vital in development of an analytical aerodynamic model that closely resembles a solid rocket motor. 6 1.3.2 Full Scale CFD The use of Computational Fluid Dynamics (CFD) has been more sparse than its an- alytical counterpart, but has been instrumental in solving unsteady aerodynamic problems in solid rocket motors. Chedevergne, Casalis, and Majdalani [16] used Direct Numerical Simulation (DNS) to investigate the validity of the Taylor-Culick pro le. Apte and Yang used Large Eddy Simulation (LES) to investigate the e ects of turbulence in a rocket motor with respect to chemical combustion rates [17]. 1.4 Motor Stability Considerations Pressure uctuations within solid motors can lead to undesired and sometimes catas- trophic e ects on the motors performance. This phenomena has been studied extensively as both a means to predict and a means to prevent. Vortex shedding has been shown to be a direct cause to these pressure uctuations [18, 19]. The developers of SPP [1] have worked extensively to include modules in their program that would allow for motor stability analysis. Coats and Dunn [20] linked SSP (an o shoot analytical approach utilizing the Taylor-Culick pro le) to SPP to automate stability predictions. This work was further developed by French to investigate tangential mode instabilities in grains with even and odd numbers of slots [21]. 1.5 Propellant Solid Mechanics The study of propellant solid mechanics involves examining the e ect of stress and strain on the propellant block from both resting during pre-launch and accelerating during launch. Slump is the term used to describe the displacement of the interior propellant as opposed to the propellant bonded to the case wall from resting or accelerating during launch. This e ect has been known to alter motor performance. Fiedler et al. investigated this e ect directly via simulation by coupling a solid mechanics model to a solid motor grain regression model [22]. Renganathan, Rao and Jana investigated e ects of slump on segmented motors under 7 storage conditions using FEA and compared this analysis with experimental results [23]. Lajczok [24] performed an analysis on propellant properties during ignition. This study determined an e ective propellant modulus to accurately predict propellant deformation under the immense pressure gradient prevalant shortly after ignition. 1.6 Automated Grain Design and Performance Control Automated grain cross-section design and performance control has recently become a hot topic in the SRM eld with the advances in computing technology. The "Achilles heel" of solid rocket motors has always been a lack of variability in the thrust-time pro le. One primary reason for choosing a liquid rocket over a solid rocket motor is the added capability of a liquid to throttle. Nunn and Cha n developed a method for throttling a solid by way of a variable sized nozzle throat and gas injector in the combustion zone forward of the nozzle but aft of the propellant grain [25]. Of course, the issue of throttling a solid motor only becomes relevant when o design performance is required, and has been achieved in the past using pintles. If a motor is designed properly (i.e. designed for a speci c mission), throttling can be avoided. Anderson rst showed how this could be possible by linking a fully encompassing missile analysis tool (i.e. complete with aerodynamics, propulsion, mass properties, guidance and control, and 6-DOF) to a genetic algorithm [26]. This tool proved e ective for designing entire missile interceptor systems with as much detail as was desired quickly and e ectively. The idea for automated design using a genetic algorithm inspired Jenkins [27] and Albarado [28] to develop optimization programs to design uniform cross section solid motors for speci c missions. Their work entailed development of a particle swarm/pattern search hybrid optimization scheme similar to that of Woltosz [29]. 8 Chapter 2 Existing Internal Ballistics Models Many of the common techniques for modeling internal ballistics (grain regression) are analytical in nature. As Barrere demonstrates[30], the grain cross section for a uniform solid rocket motor can be decomposed into a handful of parameters that geometrically de ne the motor. Using this prede ned set of parameters, analytical equations can be developed to describe the internal ballistics of the motor. SPP expanded on this idea to produce the rst method for analyzing motors that are nonuniform along the grain length (prominently the Shuttle boosters). This method was modi ed heavily in order to accurately model geometrically complex motors, and the true accuracy of the program was never fully realized. Only as recently as 2005 was a new method developed for performing the grain regression problem [31]. Willcox found that a concept known as a Signed Distance Function (SDF) could locate a surface within a domain, and determine where a surface will be at some prescribed burn distance later. Cavallini showed how an extension of the SDF called the Level Set Method (LSM) could be utilized to include erosive burning and internal aerodynamics with their programs Solid Propellant rocket motor Internal Ballistics (SPINBALL) and the Grain REGression model (GREG) [32, 33, 34]. To date, these two programs collectively represent one of the most accurate rocket motor analysis tools in literature and o ers the computationally most e cient method for solving the grain regression problem. The present work aims to improve on their development of the grain regression problem by improving the propagation routine developed in GREG using a higher order integration scheme known as the Discontinous Galerkin Method. 9 2.1 Analytical Methods Classically, solid rocket motor regression is performed via analytical techniques exploit- ing the geometric method for de nition. The limiting factor for almost all solid rocket motor grain simulations exploiting this idea is their dependence on a speci c method for de ning the grain cross-section. There is no standard for parametrically de ning rocket motors (such as wing span and chord length for aircraft wings), and thus these analytical techniques, while accurate, fail to be robust. Nevertheless, these analytical approaches will provide so- lutions for comparison to the LSM approach, and warrant further discussion. Barrere [30] is credited with developing the rst systematic analytical techniques for star and wagon wheel con gurations. The star grain, as de ned by Barrere, is shown in the schematic below in Figure 2.1. A complete description of the Barrere star grain analytical method is given in Appendix A. Figure 2.1: Star Grain De nition [27] 10 Using the same parameters, Hart eld et al. [35] review the analytical solution to the long spoke wagon wheel originally found in [30] and describe the analytical solution for the short spoke wagon wheel. Reference [36] gives an analytical description for how to solve grain regression in a slotted grain solid rocket motor. Sforzini [37, 38] and Ricciardi [39] developed expressions for performing grain regression in star grains and truncated star grains (slotted tube grains). Coats and Dunn [1] expanded on the analytical methods approach and developed ana- lytical equations for simple 3-D shapes that would be found in a typical solid rocket motor. These shapes are cone, cylinder, prism, sphere, and torus. From these shapes, designs such as nocyls, tapered stars, and forward stars can be designed and analyzed using SPP. This code has been expanded upon further to include the following advancements: 1. Motor stability and combustion stability predictions [20, 40, 41] 2. Advanced motor de nition options including segmented motors, dual propellant mo- tors, and case insulation [42] 3. Surface mesh generation [42] 4. Ignition transient calculations [42] 5. 3-D nite-volume gas-particle solver for internal aerodynamics [43] 2.2 Face-O setting Method Another popular method for surface propagation is called the Face-O setting Method (FOM). This method, formalized by Jiao [44] is of Lagrangian formulation, meaning that the surface or curve of interest is explicitly propagated in a discrete manner. In order to perform this method, each panel or facet is propagated outward. Then the vertices connect- ing the faces are reconstructed using the method described by Jiao [44] with appropriate logic to handle the intersections of faces. To improve mesh quality, the vertices are then 11 curr ent int er face new int er face wavefr ont adve ct ive Figure 2.2: Face-O setting Method: Advective Reconstruction vs. Wavefront Reconstruc- tion redistributed along the new surface. The vertex reconstruction can be completed in one of two manners: advective and wavefront motion. An illustration of each is shown below in Figure 2.2. Advective reconstruction would be more suitable for concave surfaces while wavefront reconstruction would be more suitable for convex surfaces. FOM was developed further by researchers at the University of Illinois at Urbana-Champaign into a program called Rocstar [45, 46, 47, 48, 49]. Rocstar was developed as an analysis tool to perform everything from basic SRM performance analysis to o -design analysis that includes e ects such as muliphase core ow, propellant solid mechanics (slumping) and erosive burning. The Rocstar software has been used to model propellant solid mechanics problems in the Space Shuttle boosters [46], turbulent ow e ects in the RSRM [48], and propellant slumping in the Titan IV [49]. 12 Chapter 3 Eulerian Approaches to Surface Tracking In order to ensure that an SRM code can handle arbitrary motor geometries, i.e. the analysis method has no dependence on geometric de nition, a discrete surface tracking method must be employed. There are two schools of thought to surface tracking that can be employed for the grain regression problem: Lagrangian and Eulerian. A Lagrangian formula- tion, like the face-o setting method used in Rocstar, is fraught with numerical and stability issues. Consider the closed curve given in Figure 3.1a. s s s s n n n n n (a) Closed Curve (b) Marker Particles Figure 3.1: A Closed Propagating Curve Using marker particles along the boundary (see Figure 3.1b), the Lagrangian formulation of surface propagation would dictate that each particle move normal to the curvature of the surface. For discretized treatments, this creates numerical issues in areas where particles begin to cluster together or pass by each other. This formulation also creates accuracy issues in areas where particles begin to spread out. Figure 3.2 demonstrates these numerical issues in convex and concave surfaces, showing interpretation of the particles versus what should physically occur. In order to alleviate these problems, subroutines must be developed 13 Original Surface Propagated Surface Original Surface Propagated Surface Convex Curvature Concave Curvature Figure 3.2: Convex and Concave curvature with Lagrangian propagation (left) and actual propagation (right) to determine which surfaces are convex and which are concave. Appropriate logic to round out convex curvature and create sharp points for concave curvature must also be employed. In cases where the front velocity is a function of the curvature, concave surfaces will form sharp edges, leading to singularities in the front velocity. When this occurs, a shock in the solution forms and smoothing must be applied in order to continue with the solution. The methodology used in [44, 45] has accounted for the issues associated with the marker-point method; however, more elegant solutions to this problem are available. Another popular surface tracking method is \volume of uid" (VOF) method which is an Eulerian formulation and uses a stationary mesh and tracks the motion of the interior region. Each cell is assigned a value based on the percentage of the element volume that is occupied by the surface or uid of interest. This value is 0 in elements where the uid is not present, between 0 and 1 where the uid or surface cuts through a cell, and 1 where the entire cell is occupied (see Figure 3.3a). VOF handles the surface implicitly; however, because the only information stored for each cell is the percentage volume occupied, information regarding the curvature of the surface edge is not easily ascertainable. Figure 3.3b shows the interpretation by the VOF method of the curve given in Figure 3.3a. From the Figure, it is readily apparent how the initial curvature is lost in the VOF interpretation. The Level Set Method (LSM) is an implicity method much like the VOF method to studying the evolution of boundary between two regions. The method was formalized by 14 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.3 0.0 0.0 0.0 0.0 0.2 0.4 0.3 0.1 0.2 0.4 0.7 1.0 1.0 1.0 0.8 Solid Fluid (a) Numerical Representation 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.3 0.0 0.0 0.0 0.0 0.2 0.4 0.3 0.1 0.2 0.4 0.7 1.0 1.0 1.0 0.8 Solid Fluid (b) Advection Front Figure 3.3: Volume of Fluid Representation Osher and Sethian in 1979 [50, 51]. The primary di erence between LSM and VOF is the information stored within the domain. In VOF, the percentage of the cell occupied by a uid or region of interest is stored. For a level set, the minimum distance to a closed surface within some domain is stored along with a sign indicating interior or exterior to the surface. By storing the distance do the surface rather than just the volume of the cell occupied by a region, more information regarding the evolution of the initial surface can be inferred. The level set method has been proven useful in numerous areas of engineering including graphics, computational uid dynamics, material sciences, and other elds. Rather than explicitly tracking surface propagation, the signed distance function, simply referred to as the level set, is propagated which in turn implicitly propagates the surface. Propagation of the level set is handled using a transport equation of the form @u @t + @u @x = 0 (3.1) where is the characteristic wave speed for the posed problem. This approach o ers a convenient and inexpensive alternative to the Lagrangian formulations which can become 15 cumbersome with the logic and subroutines necessary to handle convex and concave curva- ture. The level set equation has some properties that can be taken advantage of to develop a stable scheme for integration. For one, the magnitude of the gradient of the signed distance function is equal to one everywhere in a general sense. In 2-D, if a uniform velocity eld is assumed, the level set will simply convect in the direction of the z-axis with no shape changes. However, if a non-uniform velocity eld is desired, there are relatively straightforward and inexpensive methods for implementing them (see [52, 53]). The LSM approach propagates a function representing the distance to a curve in dis- cretized space as opposed to propagating a discretized curve. It is this characteristic which gives the level set method an advantage over Lagrangian approaches. LSM handles convex and concave surface implicitly, eliminating the need for logic and subroutines to handle each case. The level set of a curve describes the surface as a function of distance from the original curve. At each time step, the zeroeth level set becomes the curve itself. To initialize the level set for a curve, begin with a discretized curve and overlay a discretized mesh on top of it. The level set is then initialized by nding the signed minimum distance (SDF) from each point in discretized space to the curve. Figure 3.4 represents a discretized curve and its complementary level set function at initialization. Determining the minimum distance from each grid point to the initial surface is simple enough, using the distance formula: DFi;j = min q (xi;j xk)2 + (yi;j yk)2 8 x; y (3.2) Of course, using Equation 3.2 will always yield a positive value, but for grid points located inside of the surface, the desired outcome is a negative value. In order to make this deter- mination, the winding number of the point with respect to the curve is calculated [54]. The winding number is calculated by summing the angles subtended by each of edges as seen by the point of interest. The summation is always either 0 or n . If the sum is n then 16 +k +i +j (a) Curve in Discretized Space 00 0 1 20001 -1-1 0 0 1-1-101 -2-2 -1 0 0-2-101 -2-2 -1 -1 0-2-101 0-1 0 0 1-1001 10 1 1 20012 21 2 2 31123 11 1 2 31112 (b) Signed Minimum Distance Figure 3.4: A Discretized Curve and its Complementary Level Set the point lies within the curve, and if the sum is 0, the point lies outside the boundaries of the surface. A graphic representation of this is shown in Figure 3.5, where A has a winding number of 1, and B has a winding number of 0. This method for determining the signed minimum distance function is performed the same regardless of whether the problem is 2-D or 3-D in nature. The signed minimum distance is what allows LSM to handle convex and concave curves implicitly, in that the gradient of the level set will determine whether surfaces develop into sharp points (concave) or round valleys (convex). ?? ?? ???=360 Sur rounded ???=0 Not Sur rounded A B Figure 3.5: Winding Number Diagram 17 According to Osher and Sethian [50, 51], the Level Set equation can be written as 8 >< >: t +rb (~x;t) ~r = 0 (~x;t = 0) = 0; 8~x2 (3.3) Di erentiating Equation 3.3 with respect to time yields @ @t + X xi@xi@t = 0 (3.4) @ @t + ~V ~r = 0 (3.5) Here it is important to observe that the front propagation is related only to the normal component of the velocity eld because the gradient, ~r , at each point in ~x is de ned normal to the level set. In fact, as Cavallini [32] shows, when writing the velocity vector in terms of its normal and tangential components, and substituting into Equation 3.5, Equation 3.6 can be arrived at, which is equivalent to the level set equation as described by Osher and Sethian above. @ @t +Vn ~r = 0 (3.6) The level set equation, Equation 3.6 is a partial di erential equation that falls into a class of equations known as Hamilton-Jacobi equations, where t +V ~r = 0 ! t +H( x; y; z) = 0 (3.7) where H is the Hamiltonian de ned for Equation 3.6 as H(x;y;z; x; y; z;t) = V ~r (3.8) 18 Chapter 4 First Order Grain Regression Program 4.1 Development and Implementation The rst order two dimensional burn regression program was developed as a test bed for the level set method. This program was developed using a similar setup to that found in references [32, 33]. From a mathematical sense, the approach taken to solid rocket motor simulation is the same approach taken in the studies of uid dynamics except that the equations of motion are di erent. Equation 3.3 must to be rewritten into a scheme that can be easily integrated. For the rst order code, it was written of the form n+1j;k = nj;k tnrbjnj;k max max D x nj;k;0 2;min D+x nj;k;0 2 + max max D y nj;k;0 2;min D+y nj;k;0 2 12 (4.1) where D x nj;k = n j;k n j 1;k x ; D +x n j;k = nj+1;k nj;k x D y nj;k = n j;k n j;k 1 y ; D +y n j;k = nj;k+1 nj;k y This scheme is rst order in space and time. To enhance the order for the time integration, back and error compensation was utilized as described by Dupont and Liu[55]. To perform back and forth error compensation, the initial level set is propagated forward in time one time step through normal operation. The result of this initial propagation is then propagated backwards in time, yielding a new level set at the same time step as the initial level set. The 19 di erence in the original level set and the new level set is then used to correct the original level set. The next time step is nally calculated using this corrected level set. Mathematically this appears as: ^ n+1 = f ( n;+rnb ) (4.2) ^ n = f ^ n+1; rnb (4.3) n = n + n ^ n 2 (4.4) n+1 = f ;+rnb (4.5) While this method works for propagating the propellant surface, what is of importance to rocket motor simulation is the surface area and chamber volume evaluation at each time step, as this determines chamber pressure and propellant weight. The volume and area can be calculated simply as: V = Z h( (x))dx (4.6) A = Z ( (x)) ~r dx (4.7) where h is the Heaviside step function h( ) = 8 >< >: 0 0 1 > 0 (4.8) and ( ) is the Dirac delta ( ) = dh( )d (4.9) 20 Because Equation 4.8 is a step function, Equation 4.9 cannot be evaluated easily in a numerical scheme. Furthermore, these evaluations, while sound mathematically, do not translate well to numerical schemes as theorized. To handle this problem, some smearing was applied to h and consequently as follows. h( ) = 8 >>> >< >>>> : 0 < 1 2 + 2 + 1 2 sin j j 1 > (4.10) ( ) = 8 >< >: 1 2 + 1 2 cos j j 0 > (4.11) In Equations 4.10 and 4.11, is a tuning parameter de ned by some grid spacing char- acteristic length. Figure 4.1 demonstrates the smearing function for the Heaviside function and Dirac delta. Contributions for both surface area and chamber volume are made only by those cells contained within the boundary of the case. 21 Figure 4.1: Smearing functions for arbitrary near the zeroeth level set 4.2 Results The 2-D rst-order accurate program (henceforth referred to as AUBurnSim2D) was developed to demonstrate the e ectiveness of LSM at accurately capturing the burn re- gression physics. The following results serve as a justi cation for moving forward with the development of a higher order scheme. Neutral Burning Star According to the Barrere model [30], it is possible to analytically design a motor which exhibits a neutral burning pro le for phase I burning. Using this approach, a neutral burning motor was designed and the analytical result was comparied with the result from AUBurn- sim2D for validation. Motor performance was simpli ed to lumped parameter analysis and local burn rates, erosive burning, and chamber aerodynamics were ignored in order to test 22 geometric competence of the level set approach. The motor tested was a 7-pointed star, shown in Figure 4.2. Figure 4.2b was generated using a mesh size of 200 200. (a) Motor Surface and Case (b) Signed Minimum Distance Figure 4.2: A Neutral Burning Seven Pointed Star (a) Burn Area vs Burn Distance (b) Pressure vs Time Figure 4.3: Neutral Burning Star Results The results of this simulation are shown in Figure 4.3. From the Figure, it is apparent that LSM performs well at predicting burn area (and therefore pressure) over time. However, there are a few key problems to note about Figure 4.3. While LSM was able to accurately 23 trend the motor (neutral burn at the beginning followed by a progressive burn; burn time is correct), there was an underprediction in total pressure throughout burn of approximately 15 psia. Secondly, near burnout a dip in pressure is observed when using LSM, which disagrees with the analytical solution. Both of these e ects can be attributed to the smearing required when calculating area and chamber volume. Shown in Figure 4.4 is a contour of the smeared Heaviside function for the 7-pointed star at ignition. From Equation 4.10, the function of the Heaviside equation is to step from a value of 0 to a value of 1 at the boundary, but in a smooth manner. This causes an apparent bandwidth at the boundary of approximately 2 (dark region dividing the lighter sections in Figure 4.4). Figure 4.4: Smeared Heaviside function for 7-pointed star at ignition A closer look at this bandwidth when interacting with the motor case reveals the source of the issue with the dip in pressure near burnout. Figure 4.5 is a schematical representation of what is occurring at the case wall. The actual motor surface is still completely within the boundary, and as such should be completely contributing to the surface area. However, because the surface must be represented numerically as this bandwidth, we see that part of the numerical surface is outside the case boundary. As such, areas of the numerical that should be contributing to the burn area have been neglected. The obvious x to this issue is 24 to slim down the bandwidth as much as possible with a ner mesh, which will in turn drive up memory useage and wall clock time. Case Actua l Sur face Numerical Sur fac e Negle ct ed Sur fac es C ont ribut ing Sur fac e Figure 4.5: Motor Surface near the Case Boundary A grid sensitivity study was performed using this neutral 7-pointed star to show how a ner mesh will in fact reduce error but will drive up memory usage and wall clock time. The mesh sizes studied varied between 50 50 and 200 200. The results of this grid study are shown below in Figure 4.6. The grid study showed that the RMS error decreased with ner and ner meshing. In fact, it should expected that a ner mesh will in fact produce a more accurate result. But also as was expected, the total runtime increased with a ner mesh size. From Figure 4.6, the lowest RMS error achievable is quickly realized by increasing mesh size, but the lowest RMS achievable was still around only 10%. It can be concluded from this study that attempts to resolve the grid further are not warranted as the gain in RMS error will be unattractively low with increasing memory useage and wall clock time. 25 0 50 100 150 200 250 300 35010 11 12 13 14 15 16 17 18 19 20 Memory Useage (kB) RMS error % (a) Memory Useage vs RMS Error 0 50 100 150 200 250 300 3500 1 2 3 4 5 6 Memory Useage (kB) Wall Clock Time (s) (b) Memory Useage vs Wall Clock Figure 4.6: Grid Sensitivity Results O center Right Circular Perforated Grain This second case using AUBurnSim2D is an o center right CP grain to demonstrate the robustness of LSM at handling more complex grains. The initial surface is simple enough, but the burning of the surface poses a problem for analytical methods when a small portion of the surface burns into the motor case before the rest of the surface. This case represents a common problem when manufacturing these motors: perfectly centering the bore in the grain. The case presented here represents an exaggerated version of this problem. Shown in Figure 4.7 is the time evolution of the o center CP grain. From the Figure, the surface initially burns into the motor case around at time of t = 5:0 seconds. A CP grain should have a progressive burn pro le until the surface burns into motor case, at which point the burn becomes regressive. From Figure 4.8, the result is exactly that. At time t = 5:0 seconds, the pressure changes from progressive to regressive. 26 (a) t = 0 s (b) t = 5 s (c) t = 8.5 s (d) t = 12 s Figure 4.7: O center Circular Perforated Grain Regression Figure 4.8: Pressure Time Pro le for an O -Center Cylindrical Perforated Grain 27 Chapter 5 Higher Order Grain Regression Program 5.1 Development and Implementation The following sections detail the methods used to develop the higher order grain regres- sion program. A discussion of the discontinuous Galerkin method is provided along with a detailed discussion of how to implement DG with Hamilton-Jacobi equation, as will be needed for the level set method. The next section discusses some of the modi cations that were needed in order to stabilize the di erential equation for integration. The nal sec- tion of development and implementation discusses the method in which the motor case was incorporated into the level set area calculation. 5.1.1 Discontinuous Galerkin Method The discontinuous Galerkin method (DGM) combines avors of nite di erence, nite element, and nite volume methods. Finite di erence schemes approximate the solution locally with 1-D polynomials and satisfy the di erential equation in a pointwise fashion. This makes nite di erencing schemes simple to implement and still allow for high-order implementation. However, implementing high-order schemes cannot be implemented with geometric variance. Finite volume schemes also approximate the solution locally, using a cell average as opposed to polynomials. Finite volume schemes are robust and fast just as the nite di erencing schemes due to local solution approximation. But, nite volume schemes cannot be modi ed to be higher order in general. Lastly, nite element methods de ne the solution in a non-local manner, satisfying the equation across the entire domain which make them suitable for high-order. However, they are implicit in time. The desired formulation would be to have local high-order approximation that can work with complex geometries. 28 Hesthaven [56] provides much of the technical background required to get started with DG. A comprehensive synopsis on development of DG can be found in Appendix B. Much of the discussion found in Appendix B emenates from a variety of sources, primarily [56] but also [57, 58, 59, 60]. Appendix C gives an example comparison of DG to a conventional nite volume method for solving the classical isentropic vortex problem. The basic transport equation @u @t + @f(u) @x = 0 (5.1) in a DG sense becomes Z Dk bki @u k h @t dx Z Dk @bki @kf k hdx+ I @Dk bkif uk; h ;uk;+h ;nk; x ds = 0 (5.2) In matrix vector form, this can be rewritten as JkMdeu k dt S T efk +Laef a k +L bef b k = 0 (5.3) Since we want to determine the time derivative term for each cell, the nal form of the algebraic equation becomes deuk dt =J kMST efk JkM Laef ak +Lbef bk (5.4) Equation 5.4 is useful when the exact ux and numerical ux is straightforward to calculate, such as the case with the Euler equations. For states such as the level set equation, determining the gradient is more di cult and relies on the auxiliary variable approach like the rst order LSM. However, from Equation 5.4, the tools needed to calculate the derivative of the level set equation are already available. Looking at Equation 5.5 Z b@u@tdx+ Z b@F@xdx = 0 (5.5) 29 In order to use with the level set, DG must be successfully implemented using a Hamiltonian. Cheng and Shu [61] discuss how to incorporate the basic Hamiltonian into DG schemes. Yan and Osher [62] describe how this Hamiltonian term can be replaced with a modi ed Hamiltonian term to incorporate arti cial viscosity. Z bdudtdx+ Z b ^H(ux)dx = 0 (5.6) where ^H(p1;p2;q1;q2) = H p 1 +p2 2 ; q1 +q2 2 12 (p1 p2) 12 (q1 q2) (5.7) This is just one of many ways to represent the modi ed Hamiltonian term. The Gudonov scheme in Equation 4.1 is also of Hamilton-Jacobi form where ^H takes on the form of ^H(p1;p2;q1;q2) = tnrbjnj;k max max (p1;0)2;min (p2;0)2 + max max (q1;0)2;min (q2;0)2 1 2 (5.8) where p represents the x-derivative of using data from either the left or the right (1 and 2) and q represents the y-derivative of . Development of the derivatives in the x- and y- directions can be developed using the same operators used in Equation 5.4. In integral form, the derivatives in the x-direction would take the form Z Dk bki @u k h @x left = Z Dk @bki @xu k hdx+ I bkif uk; h ;uk+1; h ;nk; x ds (5.9) Z Dk bki @u k h @x right = Z Dk @bki @xu k hdx+ I bkif uk+1;+h ;uk;+h ;nk; x ds (5.10) and in matrix vector form are expressed as 30 p1 = JkSTui +JkLxrightui+1left JkLxleftuileft (5.11) p2 = JkSTui +JkLxrightuiright JkLxleftui 1right (5.12) i-1 i+1i Right Data Lef t C ell Left Data Right Cell Left Data Current C ell Right Data Current C ell Figure 5.1: Stencil Used to Develop DG Derivatives Using Figure 5.1 as a reference, Equation 5.11 can be interpreted as the derivative is equivalent to taking the derivative of the basis function for the cell plus lifted data from the left edge of the cell to the right minus lifted data from the left edge of the current cell. Similarly, Equation 5.12 can be interpreted as being equivalent to taking the derivative of the basis function for the cell plus lifted data from the right edge of the current cell minus lifted data from the right edge of the cell to the left. Derivatives in the y-direction are performed in exactly the same manner by replacing the sti ness matrices and lifting matrices appropriately and using data from upper and lower surfaces of cells above and below the current cell. The idea of using data from the left or data from the right is analogous to upwinding used for conservation laws. 5.1.2 Stabilizing the Di erential Equation To properly couple discontinuous Galerkin with the level set equation, additional mea- sures must be taken to ensure that the numerical scheme is stable and accurate. Using a 31 circle as an example, shown in Figure 5.2, sharp gradients occur in the level set when points in the domain are equidistant to multiple locations on the surface (such as the center of the circle). With any high order scheme, sharp gradients will cause numerical instabilities if not treated properly. In the case of the level set, a sharp gradient at the center of the circle will continue to grow sharper, corrupting the solution. (a) Isometric View (b) Side View Figure 5.2: Level set for a circle of radius = 0.25 Another issue that it is di cult to avoid is a case where the gradient is zero. In the example of the circle, this will happen if the center point, where the sharp gradient occurs, is located within a cell, rather than at the boundary of a cell. In that case, the sharp point will be smoothed out, e ectively changing the state at that point from a sharp gradient to a gradient of zero. Using the basic time evolution scheme described in Chapter 5.1.1 for discontinuous Galerkin schemes, zero gradients (such as shown in the center of the circle shown in Figure 5.3, will develop incorrectly. From the Figure, the zero gradient at the center e ectively propagates outward toward the domain edges, steadily attening the level set. Remember, the level set should simply convect in the z-direction for some constant wave speed. Thus, the result in Figure 5.3b does not make physical sense. Osher and Fedkiw [52] provide a number of corrections for this, one being a method known as reinitialization. Reinitialization is a process performed periodically to \reset" 32 (a) Initial Surface (b) Propagated with basic ODE Figure 5.3: Level set propagation for a circle of radius = 0.25 the gradient of the level set to its original shape (or close to it). The principal behind reinitialization is that only behavior at or near the zeroeth level set need to be conserved, and that behavior elsewhere can be reset without issue. The basic equation for reinitialization is the standard di erential equation with constraint added, mathematically represented as: h(x) = sign( ) (jr j 1) (5.13) This constraint essentially takes advantage of the fact that the magnitude of the gradient of the level set should be equal to 1 everywhere. The downside to reinitialization is that it is costly due to the fact that it is an iterative process that is separate from the time evolution. However, it stands to reason that if the constraint were added to the original di erential equation to be solved, smaller corrections to the level set would occur simultaneously with time evolution, and reduce the cost associated with reinitialization. The equation of motion thus becomes @ @t +jr j= sign( ) (jr j 1) (5.14) The same circle that was propagated without this correction in Figure 5.3 is given below in Figure 5.4. From the Figure, the level set was correctly propagated from state 5.4a to 5.4b. 33 (a) Initial Surface (b) Propagated Surface Figure 5.4: Level set propagation for a circle of radius = 0.25 with source term At this point, Equation 5.14 represents a robust equation for propagating simple, well behaved level sets, such as that for the circle. However, for more complex grains, such as a ve pointed star shown below in Figure 5.5, more work must be done on the di erential equation. Speci cally, convex regions of the star (areas where the zeroeth contour will remain sharp) will cause numerical issues as time evolves. (a) Top Down View (b) Isometric View Figure 5.5: Five Pointed Star Level Set 34 To demonstrate the numerical instability, the normalized residual was tracked during time evolution. It is worth noting here that the normalized residual at every time step should be 1. This is due to the special property of the level set which states that the gradient must be equal to one everywhere at all times. Since the 2-D level set is essentially being convected in the z-direction, the residual should be the same at every time step. The normalized residual for discontinuous Galerkin is calculated as Rn = qP i P krhs 2 i;kwk R1 (5.15) By tracking this quantity over time, a sense of stability can be measured. To determine the proper time step, the Courant-Friedrichs-Lewy (CFL) condition was satis ed as follows CFL = t d 2 x (5.16) where is the maximum wavespeed for the problem at hand. The traditional CFL condition essentially states that information at one point should not travel farther in a single time step than the adjacent cell. For DG, this de nition is expanded to include the total degrees of freedom for a cell (d2) [56]. The desired simulation time for this was set to a pseudo-time of 0.5, meaning that the level set values at simulation end should be exactly 0 + 0:5. The residual for this simulation is shown below in Figure 5.6. From the Figure, the simulation grew unstable very quickly, since the simulation only reached a psuedo-time of 0.03 out of 1.0. To relieve this e ect, an explicit arti cial viscos- ity term was added to the di erential equation as discussed by [63] and [64], resulting in Equation 5.17. @ @t +jr j= sign( ) (jr j 1) + xr (r ) (5.17) 35 0 0.005 0.01 0.015 0.02 0.025 0.03 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 Psuedo?Time Residual Figure 5.6: Normalized residual for ve-pointed star using mesh of 10 10 with 9th order polynomial In Equation 5.17, x is the di usion constant which is directly controlled by the user. Arti cial viscosity, implicit or explicit, reduces the gradients in the state, e ectively smooth- ing the state. Whether this is physical or not, it helps keep the solution stable. For the level set equation, it is desirable to keep the sharp ridges in the state as well maintained as possible; however sharp gradients tend to cause numerical instability. For this reason, the di usion term was only used when jr j was excessively large (jr j >> 1) or nearly 0. When either of these conditions were met, the tuning parameter was set to the desired value, otherwise the di usion term was not added to Equation 5.17. In order to tune the dissipation constant and the following steps were taken. a136 Set the CFL number to a su ciently low number to drive the timestep down. A typical value is around 0.2. a136 Begin with a low dissipation constant, and track the residual. set to 0 will give some insight into how much dissipation will be required. 36 a136 Slowly raise the dissipation constant until the residual remains approximately 1 or slightly less than 1. a136 Drive the CFL number back up as high as possible while remaining stable The results of this process for the ve pointed star in Figure 5.5 are shown below in Figure 5.7. 0 0.1 0.2 0.3 0.4 0.5 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 Pseudo?Time Normal ized Resid ual ? x = 0.000 ? x = 0.006 ? x = 0.009 Figure 5.7: Normalized residuals for ve-pointed star for varying dissipation constants As a comparison demonstrating why it is necessary to keep the residual close to 1, the evolution of x = 0:006 was compared to x = 0:009 in Figures 5.8 and 5.9 respectively. From Figure 5.8b, small wrinkles start to form which directly lead to the rise in the residual observed in Figure 5.7. This is where tracking the residual and ne tuning the dissipation constant becomes important. Although the simulation in Figure 5.8 continued to run without excessive rises in the residual, the simulation had already been corrupted, leading to a poor solution. From Figure 5.9,the level set simply convects in the z-direction, maintaining its original shape. While a few wrinkles do develop, there is enough dissipation to keep them from growing and compromising the solution. Figure 5.10 shows the zeroeth level set for each case at an approximate time of 0.2 after start. Clearly from Figure 5.10a, the zeroeth level set, the area that is most important in the domain, was corrupted without 37 (a) t=0.022 (b) t=0.110 (c) t=0.200 (d) t=0.285 (e) t=0.373 (f) t=0.460 Figure 5.8: Level set propagation for a ve pointed star with dissipation constant set to 0.006 38 (a) t=0.022 (b) t=0.110 (c) t=0.200 (d) t=0.285 (e) t=0.373 (f) t=0.460 Figure 5.9: Level set propagation for a ve pointed star with dissipation constant set to 0.009 39 su cient arti cial dissipation. When su cient arti cial dissipation is provided, the solution remains well behaved, and the zeroeth level set evolves correctly. ?2 ?1 0 1 2 ?2 ?1.5 ?1 ?0.5 0 0.5 1 1.5 2 (a) = 0:006 ?2 ?1 0 1 2 ?2 ?1.5 ?1 ?0.5 0 0.5 1 1.5 2 (b) = 0:009 Figure 5.10: Level set propagation for a ve pointed star at t = 0:2 40 5.1.3 Motor Case Implementation While it is important to evolve the level set correctly over time, it is equally important to accurately determine burn area (in 2-D, perimeter) at each time step for a rocket motor simulation program. There are issues associated with determining area accurately, particu- larly when it comes to determining area as the surface interferes with the motor case wall as was discussed in Section 4.2. Revisiting an explicitly de ned case, shown in Figure 5.11, the bandwidth associated with smoothing the Dirac delta function, used to determine where the surface is, causes an unfavorable drop in calculated area near the boundary if grid spacing is relatively large. Case Actua l Sur face Numerical Sur fac e Negle ct ed Sur fac es C ont ribut ing Sur fac e Figure 5.11: Motor Surface near the Case Boundary This e ect could be further exaggerated in DG grid where nodes are based on Lobatto quadrature as shown in Figure B.3 of Appendix B, and can be sparse in some areas and more densely populate other areas of a cell. Recall the equation for burn area: A = Z ( (x)) ~r dx (5.18) 41 This equation can be modi ed to account for an arbitrarily de ned case simply by de ning a a step function de ning locations inside and outside of the case. This new step function is incorporated into the burn area calculation as follows. A = Z hcase ( case) ( (x)) ~r dx (5.19) In this instance, hcase is de ned as hcase = 8 >>> >< >>> >: 1 case > 0 1 2 1 + cos " "< case < 0 0 case < 0 (5.20) where a negative signed distance function de nes the region exterior to the case, and the positive signed distance function de nes the region interior to the case. By implicitly de ning the case boundary, two advantages are gained: 1) the smeared surface is allowed to exceed the explicit boundary and continue to contribute to the area calculation, and 2) arbitrary case geometries can be implemented. An example of a 2-D cylindrical case and square case are given in Figure 5.12. (a) Cylindrical Case (b) Square Case Figure 5.12: Heaviside function for motor case 42 To demonstrate, the same ve-pointed star given above was simulated with a circular case of radius equal 1 as well as with a square case with side length of 2. Figure 5.13 demonstrates the star as it burns into the two di erent case walls. In Figure 5.13, the raised surface represents the remaining propellant in the case. (a) Cylindrical Case (b) Square Case Figure 5.13: Five pointed star burning into the case wall 43 5.2 Results The discussion in Chapter 3 and Section 5.1.1 has led to the development of a robust and e cient approach to solving the grain regression problem. Two example cases were tested with this modi ed DG approach in order to con rm both robustness and accuracy. The rst case is a circular grain cross section in a circular case, and was tested with varying number of cells using varying degree polynomials. A comparison was made between these higher order results and the rst order program developed in Chapter 4. This comparison was made for a circular grain cross section of radius 0.7 with a outer case diameter of 2. The number of cells were varied from 9 to 625 in higher order schemes and from 81 to 10,000 in the rst order scheme. A true comparison of the schemes looks not at the number of cells but rather the degrees of freedom per cell. Since rst order cells contain only one degree of freedom, the number of cells is equivalent to the degrees of freedom. However, for higher order cells, the degrees of freedom can be calculated using Equation 5.21, and the total work required scales with degrees of freedom. DOF = (d+ 1)2 Ncells (5.21) The higher order schemes varied from degree 3 polynomial up to degree 9 polynomial, giving total degrees of freedom ranging from 144 up to 10,000, which is more on par for the degrees of freedom for the rst order scheme. The results for the simulation are shown below in Figure 5.14. It should be noted that the error was a root-mean-square based error term between the simulation and the analytical solution for the entire level set at time t. In other words, Figure 5.14 gives a measure of a particular schemes ability to propagate the level set correctly. This comparison will only serve as a validation for using DG over conventional rst order schemes, and does not speak to the accuracy of the scheme?s to correctly determine area as a function of time. 44 10 0 10 1 10 2 10 3 10 4 10 ?4 10 ?3 10 ?2 10 ?1 10 0 10 1 Number of Cells Error P =1 P =4 P =6 P =10 (a) Error vs Number of Cells 10 1 10 2 10 3 10 4 10 ?4 10 ?3 10 ?2 10 ?1 10 0 10 1 Degrees of Freedom Error P =1 P =4 P =6 P =10 (b) Error vs Degrees of Freedom Figure 5.14: Comparison of polynomial order with error From Figures 5.14a and 5.14b, it is clearly cheaper to use higher order integration schemes over conventional rst order methods from the standpoint of number of cells and total degrees of freedom. For higher order polynomials, the reqiured number of cells is lower to achieve the same error at lower orders. The same trend is observed for error versus degrees of freedom. The second judgement for accuracy of the scheme is ability to determine area as a function of time. Figure 5.15 gives various burn area pro les for di erent con gurations of number of cells and polynomial orders. From Figure 5.15, two of the con guration have no issues propagating the level set; however, some issues arise when calculating area as the surface approaches the case boundary. In every case, a tail-o begins to form near the boundary around a distance of 0.3 due to case radius of 1 and inner radius of 0.7. Where this tail-o begins depends on the number of cells as well as the polynomial order. However, combining both high order and high number cells is not necessarily cheaper anymore. In Figure 5.16, error was compiled by comparing the analytical solution for burn area to the higher order schemes. From the Figure, the total degrees of freedom apparently dictate the error in area. This is directly related to the smoothed Heaviside function bandwidth problem previously discussed. The bandwidth is 45 0 0.050.1 0.150.2 0.250.3 0.35 4 4.5 5 5.5 6 6.5 Burn Distance Burn Area Analytical Solution I = 9, P = 4 I = 81, P = 6 I = 225, P = 10 Figure 5.15: Burn area versus burn distance comparison between analytical method and DG schemes proportional to characteristic length and inversely proportional to order polynomial. This is a situation where solution adaptivity could be useful. Modifying the mesh to use smaller cells or higher order polynomials near the case, and large cells/lower order polynomials away far from the case could alleviate this issue without penalty in computational cost. 10 0 10 1 10 2 10 3 0 0.02 0.04 0.06 0.08 0.1 0.12 Ncells Error P=4 P=6 P=10 (a) Error vs Number of Cells 10 2 10 3 10 4 10 5 0 0.02 0.04 0.06 0.08 0.1 0.12 Degrees of Freedom Error P=4 P=6 P=10 (b) Error vs Degrees of Freedom Figure 5.16: Comparison of polynomial order with error compiled using analytical area The second test case is a more customary example of a star grain that would be found in typical rockets. It is a ve pointed star with a maximum inner radius of 7 inches, a 46 minimum inner radius of 5 inches, a case diameter of 20 inches, and a grain length of 24 inches. The fuel type was assumed to be PBAN/AP/AL. Shown in Figure 5.17 is the grain cross section (red), the motor case (blue), and the successive grain regression lines (black) e ectively representing the level set function. Figure 5.17: Five pointed star with grain cross with grain regression lines This motor was tested using 81 total cells with a 9th order polynomial. Su cient domain space was established to ensure that the zeroeth level set could not occupy areas both within the motor case and outside the global domain simultaneously. This helps to establish con dence that the calculation for area does not con ict with the domain boundary and corrupt the calculation for area. Figure 5.18 shows the burn area as a function of burn distance and the chamber pressure as a function of time. From Figure 5.18a, there is an underprediction in burn area of about 30 square inches throughout the burn distance. Likewise there is an underprediction in pressure versus time as well. This underprediction in area should be expected according to [52] due to the rst order nature of the smeared heaviside and dirac delta functions. These two equations are rst order accurate no matter 47 the integration method. Therefore, there is a current limitation as to how accurate an area calculation can be using the level set method. 0 1 2 3 4 5 200 400 600 800 1000 1200 1400 1600 Burn Distance (in) Burn Area (in 2 ) Analytical Predicted (a) Burn area versus burn distance 0 2 4 6 8 10 0 50 100 150 200 250 300 350 400 450 500 Time (s) Pressure (psi) Analytical Predicted (b) Pressure versus time Figure 5.18: Comparison of ve pointed star simulated with DG versus analytical solution Despite the limitations on the area calculation, there are some positives that can be noted on these results. One key observation to make is the ability of this approach to correctly account for the circular case. For this test case, the analytical solution should (and does) end at a burn distance of 3 inches where phase 2 burning ends and phase 3 burning begins. In most analytical schemes, phase 3 burning is neglected. The DG scheme should also burn out or at least begin tailing o at the same burn distance of 3 inches. And indeed, from Figure 5.18a the tailo begins to occur at 3 inches. Also, it would be expected that the nal burnout would occur at a burn distance of 5 inches because the minimum inner radius of the grain started 5 inches from the motor case. Figure 5.19 shows the remaining propellant in the motor case at various burn distances. At a burn distance of 3 inches (Figure 5.19d), part of the propellant has reached the case wall and will no longer contribute to the burn area, as evidenced previously in Figure 5.18a. To demonstrate the generality of this SRM code, the same star grain was tested in a square case with a side length of 20 inches. The grain evolution for this test case can be 48 (a) Burn Distance = 0 in. (b) Burn Distance = 1 in. (c) Burn Distance = 2 in. (d) Burn Distance = 3 in. (e) Burn Distance = 4 in. (f) Burn Distance = 5 in. Figure 5.19: Propellant remaining for ve pointed star at various burn distances for a circular motor case 49 (a) Burn Distance = 0 in. (b) Burn Distance = 1 in. (c) Burn Distance = 2 in. (d) Burn Distance = 3 in. (e) Burn Distance = 4 in. (f) Burn Distance = 5 in. Figure 5.20: Propellant remaining for ve pointed star at various burn distances for a square motor case 50 seen in Figure 5.20. From the Figure, it is apparent that geometric generality is conserved. For the circular case, at a burn distance of 3 inches, every star point reached the case wall. However, for a square case at the same burn distance, only one star point reached the case wall. In fact, after a burn distance of 5 inches, a relatively large amount of propellant still remains as was expected. Knowing the propellant regression, it should be expected that the burn area as a function of burn distance should be identical to that of the circular case for the rst 3 inches. Shown in Figure 5.21 is a plot comparing the burn area and pressure pro les for the circular case to the pro les for the square case. 0 1 2 3 4 5 200 400 600 800 1000 1200 1400 1600 Burn Distance (in) Burn Area (in 2 ) Circular Case Square Case (a) Burn area versus burn distance 0 2 4 6 8 10 0 50 100 150 200 250 300 350 400 450 500 Time (s) Pressure (psi) Circular Case Square Case (b) Pressure versus time Figure 5.21: Comparison of pressure pro les for star grain in a circular case and a square case As was expected, the pro les are identical until the surface reaches the motor case. At the motor case, the burn area continues to increase until more star points reach the case wall thus causing the burn area to decrease. This result is useful from the standpoint that this approach can handle completely arbitrary propellant shapes, both interior or exterior. 51 Chapter 6 Capabilities and Further Development 6.1 Conclusions This thesis developed a high delity motor simulation approach that is robust and arbitrary, as opposed to analytical techniques that require derivation of the solution for each type of grain cross-section. The methodology used combined previous work in the elds of surface tracking and propagation with higher order integration schemes to develop a fast and computationally e cient approach to grain surface evolution. Problems inherent to propagation of the signed distance function were identi ed, and successfully eliminated. The level set method had already proven to be an arbitrary approach to the grain regression problem, and the addition of the discontinuous Galerkin method to the level set approach developed into the robust and computationally e cient method discussed. The end product was a solid foundation for an approach to SRM simulation that has the potential to be the cheapest and most robust method available. The following remarks are concluded from the work performed in this thesis: a136 The discontinuous Galerkin approach to surface evolution proved to be more e cient than traditional rst order schemes that have previously been coupled with the level set method. a136 The marriage between DG and LSM did require some modi cation of the level set di erential equation. Inclusion of the slope correction source term and the explicit arti cial source term did introduce complexity to the approach, but proved to be vital in properly evolving the level set over time. The source term should always be 52 incorporated into nite volume type schemes. Explicit arti cial viscosity is required for high order methods due to oscillations. a136 While it was shown that surface tracking is more accurate and e cient with both decreasing number of cells and higher order polynomial, the limiting factor on error in area calculation proved to be total degrees of freedom. a136 An approach to incorporating arbitrary motor case boundaries was developed and proved to be more robust and as accurate as previous work involving discretized meth- ods. a136 The method for calculating burn area is only rst order accurate, and as such was limited in accuracy to total degrees of freedom and problem complexity. This issue is should be addressed regardless of scheme used for integration. 6.2 Recommendations While a solid framework for an accurate solid rocket motor grain regression tool was developed, several key topics should be pursued for further improvement. A more accurate method for calculating area must be developed in order to move away from the rst order limitation currently implemented. However, a higher order scheme for developing the area calculation will increase computational e ciency in both work and wallclock. It was not discussed in this work, but a decrease in wallclock time can be achieved for higher order schemes by taking advantage of sparse matrix multiplication. In its current implementation, any manipulation of the sti ness matrix (for example) requires N2 calculations. But the sti ness matrix is sparse (see Figure 6.1), and would allow sparse multiplication of the large matrices that are used at higher order. If sparse muliplication can be taken advantage of, the time required to complete a simulation can be further reduced. This savings in time can be utilized in one of two ways. The rst utilization of time savings can trivially be taken at face value as time savings. The 53 0 5 10 15 20 25 0 5 10 15 20 25 non?zero = 110 (a) Non-zero entries for sti ness ma- trix of order 5 First Order Higher Order 0 1 Normal ized Wallclo ck Available Wallclock (b) Sample wallclock savings Figure 6.1: Example utilization of sparse matrices second utilization of time savings could be used to develop a more accurate calculation of area, include second order e ects, or any other avor of improvements to the simulation, while still remaining as cheap as rst order schemes. Solution adaptivity could also be taken advantage of to provide more accuracy. Moving forward, the current implementation of DG/LSM should be extended from two dimensions to three dimensions. Development of variable (non-uniform) wavespeed should be implemented to make inclusion of second order propellant burning e ects such as multi-phase ow and erosive burning possible. It would also be desirable to test the three-dimensional code against real world motors as opposed to simpler analytical solutions. Some nal improvements to the program would be to include the following: ignition and burnout transients, chamber aerodynamics, and propellant solid mechanics. A successful achievement of these recommendations would produce an end product believed to be the most robust and accurate simulation tool available. 54 Bibliography [1] Dunn, S. and Coats, D., \3-D Grain Design and Ballistic Analysis Using the SPP97 Code," AIAA paper 1997-3340, 1997. [2] Mattingly, J., Elements of Propulsion: Gas Turbines and Rockets, AIAA Education Series, 2006. [3] Sutton, G. and Biblarz, O., Rocket Propulsion Elements, Wiley & Sons, 2001. [4] \Solid Rocket Motor Igniters," NASA SP-8051, 1971. [5] Foster, W. and Jenkins, R., \Analysis of Advanced Solid Rocket Motor Ignition Phe- nomena," NASA CR-199427, 1995. [6] Cho, I. H. and Baek, S. W., \Numerical Simulation fo Axisymmetric Solid Rocket Motor Ignition Transient with Radiation E ect," Journal of Propulsion and Power, Vol. 16, No. 4, 1999, pp. 725{728. [7] Bai, S., Han, S., and Pardue, B., \Two-Dimensional Axisymmetric Analysis of SRM Ignition Transient," AIAA paper 1993-2311, 1993. [8] Johnston, W., \Solid Rocket Motor Internal Flow During Ignition," Journal of Propul- sion and Power, Vol. 11, No. 3, 1995, pp. 489{496. [9] Taylor, G., \Fluid Flow in Regions Bounded by Porous Surfaces," Proceedings of the Royal Society, London, Series A, Vol. 234, No. 1199, 1956, pp. 456{475. [10] Culick, F., \Rotational Axisymmetric Mean Flow and Damping of Acoustic Waves in a Solid Propellant Rocket," AIAA Journal, Vol. 4, No. 8, 1966, pp. 1462{1464. [11] Majdalani, J. and Moorhem, W. V., \Improved Time-Dependent Flow eld Solution for Solid Rocket Motors," AIAA Journal, Vol. 36, No. 2, 1998. [12] Majdalani, J. and Saad, T., \Energy Steepened States of the Taylor-Culick Pro le," AIAA paper 2007-5797, 2007. [13] Majdalani, J., Xu, H., Lin, Z.-L., and Liao, S.-J., \Exact HAM Solutions for the Viscous Rotational Flow eld in Channels with Regressing and Injecting Sidewalls," AIAA paper 2010-7079, 2010. [14] Maicke, B. A. and Majdalani, J., \On the Compressible Hart-McClure Mean Flow Motion in Simulated Rocket Motors," AIAA paper 2010-7077, 2010. 55 [15] Akiki, M. and Majdalani, J., \Quasi-Analytical Approximation of the Compressible Flow in a Planar Rocket Con guration," AIAA paper 2010-7080, 2010. [16] Chedevergne, F., Casalis, G., and Majdalani, J., \DNS Investigation of the Taylor- Culick Flow Stability," AIAA paper 2007-5796, 2007. [17] Apte, S. and Yang, V., \Unsteady Flow Evolution and Combustion Dynamics of Homo- geneous Solid Propellant in a Rocket Motor," Combustion and Flame, Elsevier Science Inc., 2002. [18] Dotson, K., Koshigoe, S., and Pace, K., \Vortex Shedding in a Large Solid Rocket Motor Without Inhibitors at the Segment Interfaces," Journal of Propulsion and Power, Vol. 13, No. 2, 1997. [19] Lupoglazo , N. and Vuillot, F., \Parietal vortex shedding as a cause of instability for a long solid propellant motors { Numerical simulations and comparisons with ring tests," AIAA paper 1996-0761, 1996. [20] Coats, D. and Dunn, S., \Improved Motor Stability Predictions for 3-D Grains Using the SPP Code," AIAA paper 1997-3251, 1997. [21] French, J., \Tangential Mode Instability of SRMs with Even and Odd Numbers of Slots," AIAA paper 2002-3612, 2002. [22] Fiedler, R., Jiao, X., Haselbacher, A., Geubelle, P., Guoy, D., and Brandyberry, M., \Simulations of Slumping Propellant and Flexing Inhibitors in Solid Rocket Motors," AIAA paper 2002-4341, 2002. [23] Renganathan, K., Rao, B., and Jana, M., \Slump Estimation of Cylindrical Segment Grains of a Typical Rocket Motor under Vertical Storage Condition," Trends in Applied Research, Vol. 1, No. 1, 2006, pp. 97{104. [24] Lajczok, M., \E ective Propellant Modulus Approach For Solid Rocket Motor Ignition Structural Analysis," Computers and Structures, Vol. 56, No. 1, 1995, pp. 101{104. [25] Nunn, R. and Cha n, L. C., \Method and Means for Controlling the Thrust in a Solid Propellant Rocket Motor," 1971. [26] Anderson, M., Burkhalter, J., and Jenkins, R., \Design of an Air to Air Interceptor Using Genetic Algorithms," AIAA paper 1999-408, 1999. [27] Jenkins, R. and Hart eld, R., \Hybrid Particle Swarm-Pattern Search Optimizer for Aerospace Applications," AIAA paper 2010-7078, 2010. [28] Albarado, K., Hart eld, R., Hurston, B., and Jenkins, R., \Solid Rocket Motor Per- formance Matching Using Pattern Search/Particle Swarm Optimization," AIAA paper 2011-5798, 2011. [29] Woltosz, W., The Application of Numerical Optimization Techniques to Solid-Propellant Rocket Motor Design, Master?s thesis, Auburn University, 1977. 56 [30] Barrere, M., Jaumotte, A., de Veubeke, B. F., and Vandenkerckhove, J., Rocket Propul- sion, Elsevier Publishing Company, 1960. [31] Willcox, M., Brewster, M. Q., Tang, K., Stewart, D. S., and Kuznetsov, I., \Solid Rocket Motor Internal Ballistics Simulation Using Three-Dimensional Grain Burnback," Journal of Propulsion and Power, Vol. 23, No. 3, 2007. [32] Cavallini, E., Modeling and Numerical Simulation of Solid Rocket Motors Internal Bal- listics, Ph.D. thesis, Sapienze Universita di Roma, Rome, Italy, 2009. [33] Favini, B., Cavallini, E., and Giacinto, M. D., \An Ignition-to-Burn Out Analysis of SRM Internal Ballistic and Performances," AIAA paper 2008-5141, 2008. [34] Cavallini, E., Favini, B., Giacinto, M. D., and Serraglia, F., \SRM Internal Ballistic Numerical Simulation by SPINBALL Model," AIAA paper 2009-5512, 2009. [35] Hart eld, R., Jenkins, R., Burkhalter, J., and Foster, W., \A Review of Analytical Methods for Predicting Grain Regression in Tactical Solid Rocket Motors," Journal of Spacecraft and Rockets, Vol. 41, No. 4, 2004, pp. 689{693. [36] Hart eld, R., Burkhalter, J., Jenkins, R., Anderson, M., and Witt, J., \Analytical Development of a Slotted Grain Solid Rocket Motor," AIAA paper 2002-4298, 2002. [37] Sforzini, R., \An Automated Approach to Design of Solid Rockets Utilizing a Special Internal Ballistics Model," AIAA paper 80-1135, 1980. [38] Sforzini, R., \Design and Performance Analysis of Solid-Propellant Rocket Motors Using a Simpli ed Computer Program," NASA CR-129025. [39] Ricciardi, A., \Generalized Geometric Analysis of Right Circular Cylindrical Star Per- forated and Tapered Grains," AIAA Journal of Propulsion and Power, Vol. 8, No. 1, 992, pp. 51{58. [40] French, J., \Three Dimensional Combustion Stability Modeling for Solid Rocket Mo- tors," AIAA paper 1998-3702, 1998. [41] French, J. and Coats, D., \Automated 3-D Solid Rocket Combustion Stability Analysis," AIAA paper 1999-2797, 1999. [42] Coats, D., Dunn, S., and Berker, D., \Improvements to the Solid Performance Program," AIAA paper 2003-4504, 2003. [43] French, J. and Tullos, J., \Solid Rocket Motor Grid Generation and CFD for Internal Ballistcs Analysis," AIAA paper 2005-4162, 2005. [44] Jiao, X., \Face o setting: A uni ed approach for explicit moving interfaces," Journal of Computational Physics, 2006. [45] Fiedler, R., Rocstar 3 User?s Guide, University of Illinois at Urbana-Champaign, 2008. 57 [46] Dick, W., Fiedler, R., and Heath, M., \Building Rocstar: Simulation Science for Solid Propellant Rocket Motors," AIAA paper 2006-4590, 2006. [47] Jiao, X., Zheng, G., Lawlor, O., Alexander, P., Campbell, M., Heath, M., and Fiedler, R., \An Integration Framework for Simulations of Solid Rocket Motors," AIAA paper 2005-3991. [48] Fiedler, R., Wasistho, B., and Brandyberry, M., \Full 3-D Simulation of Turbulent Flow in the RSRM," AIAA paper 2006-4587, 2006. [49] Fiedler, R., Namazifard, A., Campbell, M., Xu, F., Aravas, N., , and Sofronis, P., \Detailed Simulations of Propellant Slumping in the Titan IV SRMU PQM-1," AIAA paper 2006-4592, 2006. [50] Osher, S. and Sethian, J., \Fronts Propagating with Curvature Dependent Speed: Al- gorithms Based on Hamilton-Jacobi Formulations," Journal of Computational Physics, Vol. 79, No. 2, 1988, pp. 12{ 49. [51] Sethian, J., \Curvature and the Evolution of Fronts," Communication of Mathematical Physics, Vol. 101, No. 4, 1985. [52] Osher, S. and Fedkiw, R., Level Set Methods and Dynamic Implicit Surfaces, Texts in Applied Mathematics, Springer, 2003. [53] Adalsteinsson, D. and Sethian, J., \The Fast Construction of Extension Velocities in Level Set Methods," Journal of Computational Physics, Vol. 148, pp. 2{22. [54] Sutherland, I. E., Sproull, R. F., and Schumacker, R. A., \A Characterization of Ten Hidden-Surface Algorithms," ACM Comput. Surv., Vol. 6, No. 1, 1974, pp. 1{55. [55] Dupont, T. and Liu, Y., \Back and forth error compensation and correction methods for removing errors induced by uneven gradients of the level set function," Journal of Computational Physics, Vol. 190, 2003, pp. 311{324. [56] Hesthaven, J. and Warburton, T., Nodal Discontinuous Galerkin Methods, Texts in Applied Mathematics, Springer, 2008. [57] Cockburn, B. and Shu, C., \The Runge-Kutte discontinuous Galerkin method for convection-dominated problems," Journal of Scienti c Computing, Vol. 16, No. 3, 2001, pp. 173{261. [58] Karniadakis, G. and Sherwin, S., Spectral/hp Element Methods for CFD, Oxford Uni- versity Press, 1999. [59] Gautschi, W., Orthogonal Polynomials: Computation and Approximation, Oxford Uni- versity Press, 2004. [60] Golub, G. and van Loan, C., Matrix Computations, John Hopkins University Press, 1996. 58 [61] Cheng, Y. and Shu, C.-W., \A discontinuous Galerkin nite element method for directly solving the Hamilton-Jacobi equations," Journal of Computational Physics, Vol. 223, 2007, pp. 398{415. [62] Yan, J. and Osher, S., \A local discontinuous Galerkin method for directly solving Hamilton-Jacobi equations," Journal of Computational Physics, Vol. 230, 2011, pp. 232{ 244. [63] Kl ockner, A., Warburton, T., and Hesthaven, J., \Viscous Shock Capturing in a Time-Explicit Discontinuous Galerkin Method," Mathematical Modeling of Natural Phe- nomenon, Vol. 10, No. 10, 2010, pp. 1{42. [64] Persson, P. and Peraire, J., \Sub-Cell Shock Capturing for Discontinuous Galerkin Methods," AIAA paper 2006-112, 2006. 59 Appendix A Analytical method for a star grain Beginning with Figure A.1, the star grain is de ned with the following parameters: a136 Rp { maximum inner radius with no llet a136 Ri { minimum inner radius a136 N { number of star points a136 f { llet radius a136 " { fraction of the angle allowed for a single star point taken up by the structure of the star point (angular fraction) Figure A.1: Star Grain De nition An analytical expression for the grain regression for Phase I burning can be determined from these 5 parameters (from the Figure, Phase I burning ends when the star point burns into the lower boundary. From Figure A.1, the height to the origin of the llet radius is H1 = Rp sin " N (A.1) 60 The star point half angle is de ned as 2 = tan 1 H1 tan "N H1 Ri tan "N ! (A.2) The three burn perimeters (S1, S2, and S3) can be calculated as follows. S1 = H1sin 2 (y +f) cot 2 (A.3) S2 = (y +f) 2 2 + " N (A.4) S3 = (Rp +y +f) N " N (A.5) The total burn perimeter is then calculated as S = 2N (S1 +S2 +S3) (A.6) The port area is calculated as Ap1 = 12H1 Rp cos " N +H1 tan 2 12S21 tan 2 (A.7) Ap2 = 12(y +f)2 2 2 + " N (A.8) Ap3 = 12 (Rp +y +f)2 N " N (A.9) Ap = 2N (Ap1 +Ap2 +Ap3) (A.10) As previously mentioned, this set of equations de nes phase I burn. Phase II burning begins with the following condition. y +f = H1cos 2 (A.11) Phase II burn perimeter can be calculated with the following set of equations. = 2 2 + " N (A.12) = tan 1 p (y +f)2 H21 H1 ! 2 (A.13) S2 = (Y +f)( ) (A.14) 61 S3 = (Rp +y +f) N " N (A.15) S = 2N (S2 +S3) (A.16) It should be noted that the rst burn perimeter, S1, is not present in phase II because the star point has burned to the boundary, making S1 equal to zero. Phase II port area is calculated as follows. Ap1 = 12H1 Rp cos " N + q (y +f)2 H21 (A.17) Ap2 = 12(y +f)2( ) (A.18) Ap3 = 12 (Rp +y +f)2 N " N (A.19) A = 2N (Ap1 +Ap2 +Ap3) (A.20) The end of phase II burning occurs when the burn distance is equal to the web thickness which can be calculated as Web = Ro (Rp +f) (A.21) 62 Appendix B Derivation for Discontinuous Galerkin Approach Consider the linear scalar transport equation @u @t + @f(x) @x = 0; x2[L;R] = (B.1) f(u) = au (B.2) Initial conditions are de ned as u(x;0) = u0(x) (B.3) and boundary conditions are given as an in ow boundary u(L;t) = g(t) if a 0 u(R;t) = g(t) if a 0 To continue, the domain, , is approximated in K (1-D) non-overlapping elements, demon- strated as a stencil in Figure B.1, and the global state is approximated in a piecewise con- tinuous. x2Dk :: u(x;t)?uh(x;t) = KM k=1 ukh(x;t) (B.4) x l 1 =L x r K =R D k-1 D k D k+1 x r k-1 =x l k x r k =x l k-1 Figure B.1: Stencil for 1-D example From Equation B.4, the local state is de ned using a higher-order local basis function: ukh(x;t) = NX j=0 eukj(t)bkj(x) (B.5) Here, each eukj is the modal coe cient for the expansion, and bkj is the basis function a function space PN. The local residual for the approximate solution is calculated as Rkh(x;t) = @u k h @t + @fkh @x (B.6) 63 The residual is required to vanish in a Galerkin sense, so the residual becomes Z Dk bkiRkhdx = Z Dk bki @u k k @t dx+ Z Dk bki @f k h @x dx = 0; i = 0;:::;N (B.7) Thus the solutions on each element are not required to be continuous at element inter- faces (hence the name discontinuous Galerkin). Using the divergence theorem to integrate Equation B.7 results in Z Dk bki @u k h @t dx+ Z Dk @ @x bk if k h dx Z Dk @bki @xf k i dx = 0 (B.8) Z Dk bki @u k h @t dx+ I @Dk bkifkhnkxds Z Dk @bki @xf k i dx = 0 (B.9) Because the state was approximated to be piecewise continuous, special care must be given to the second term in Equation B.9. As is performed in nite element methods, a numerical ux function will replace the ux at the interface between adjacent cells. x2@Dk :: fkhnkx!f uk; h ;uk;+h ;nk; x (B.10) De ning the interface ux in this manner allows for the following: a136 Boundary conditions set using uk;+h a136 Numerical ux egressing from the right edge is the negative of the ux egressing from the left edge of adjacent cell !f uk; h ;uk;+h ;nk; x = f uk;+h ;uk; h ;nk;+x a136 Numerical ux of two identical states is simply the ux !f (u;u;nx) f(u)nx This leaves the nal equation as Z Dk bki @u k h @t dx+ I @Dk bkif uk; h ;uk;+h ;nk; x ds Z Dk @bki @xf k i dx = 0 (B.11) A large part of DG is determining exactly what to use as the bases function, bki (x). These will typically be some sort of polynomial: linear connections, cubics, Taylor series expansion, etc. For a more comprehensive listing of polynomials and the theory behind them, see [59]. Consider a single cell k in Dk. Using a standard interval of x ! 2 ( 1;1) yields access to a class of functions known as recursion polynomials, where polynomial of degree N relies on the solution to polynomial degree (N 1). The state space for this standard interval is of the form uh( )2PN (B.12) The polynomial function space is restricted to the given domain and vanishes outside of this domain. We now need to de ne the polynomial formula to t our space. Expanding the 64 state, we arrive at two possible choices: modal vs. nodal. uh( ) = NX m=0 eum m( ) = NX n=0 bunln( ) (B.13) In Equation B.13, the rst representation is modal (bi! i) and the second representation is nodal (bi!li). Each representation is mathematically equivalent, but there are advantages and disadvantages to both. Graphically, these two representations look entirely di erent (see Figure B.2 for an example of P5). From Figure B.2a, there are a total of six characteristic curves (a constant line plus one for each of the 5 polynomials) that combine in a modal sense to form the state for the cell. Figure B.2b gives the same P5 represented with 6 nodes (circles). In the nodal sense, when the basis is queried at a node, a curve is generated yielding a value of 1 at that node, and 0 at every other node. Re-examining the modal representation, uh( ) = NX m=0 eum m( ) (B.14) it is apparent that we can determine the modal coe cients using an L2 projection. Z 1 1 i( )u( )d = Z 1 1 i( )uh( )d = Z 1 1 ( ) NX m=0 eum m( ) ! d = NX m=0 Z 1 1 i( ) m( )d eum = NX m=0 Mimeum (B.15) where Mim is the mass matrix. Because a solution will require inversion of the mass ma- trix, an ill-conditioned matrix will lead to numerical inaccuracy or possibly even instability. Therefore, it is important to intelligently choose the basis function to represent the state. Simply choosing any polynomial basis, say monomials, will lead to poorly conditioned ma- trix, leading to losses in accuracy with higher-order polynomials. Using instead something like an orthonormal basis (three-term recursion relation) will lead to a more well-conditioned matrix. The three-term recurrence relation is expressed as p m+1 m+1( ) = ( m) m( ) p m 1 m 1( ) (B.16) 65 ?1 ?0.5 0 0.5 1 ?2 ?1 0 1 2 ? 0 ?1 ?0.5 0 0.5 1 ?2 ?1 0 1 2 ? 1 ?1 ?0.5 0 0.5 1 ?2 ?1 0 1 2 ? 2 ?1 ?0.5 0 0.5 1 ?2 ?1 0 1 2 ? 3 ?1 ?0.5 0 0.5 1 ?2 ?1 0 1 2 ? 4 ?1 ?0.5 0 0.5 1 ?2 ?1 0 1 2 ? 5 (a) Modal Basis ?1 ?0.5 0 0.5 1 ?1 ?0.5 0 0.5 1 1.5 | 0 ?1 ?0.5 0 0.5 1 ?1 ?0.5 0 0.5 1 1.5 | 1 ?1 ?0.5 0 0.5 1 ?1 ?0.5 0 0.5 1 1.5 | 2 ?1 ?0.5 0 0.5 1 ?1 ?0.5 0 0.5 1 1.5 | 3 ?1 ?0.5 0 0.5 1 ?1 ?0.5 0 0.5 1 1.5 | 4 ?1 ?0.5 0 0.5 1 ?1 ?0.5 0 0.5 1 1.5 | 5 (b) Nodal Basis Figure B.2: Modal versus Nodal Graphical Representation 66 where m = 0;:::;N 1, and 1 := 0; 0 := 1p 0 In these recurrence relations, and provide special subclasses of recurrence polynomials such as Chebychev and Legendre. For Legendre polynomials these coe cients are m = 0; m = 2 if m = 0 1=(4 m 2) if m 1 (B.17) Looking at the nodal representation of Equation B.13, uh( ) = NX n=0 bunln( ) (B.18) we can use Lagrange interpolating polynomials through a given set of nodal locations, i. ln( ) = NY i=0i6=n i n i (B.19) Evaluating the Legendre polynomials at the same nodal locations that we set for these Lagrange polynomials yields the Vandermonde matrix, which, as will be shown, forms the foundation for the operators required for use in DG. The use of Lagrange polynomials and Legendre polynomials o er some powerful properties. Lagrange interpolating polynomials are non-heirarchical, meaning that each polynomial is of degree N. This leads to the property demonstrated in Figure B.2b. On the other hand, Legendre polynomials are heirarchical, where each polynomial is of 1 degree higher than the previous polynomial degree, shown in Figure B.2a. Lagrange polynomials are convenient for boundary data extraction while Legendre polynomials provide an accurate way to perform integration. Vij := j( i) (B.20) The Vandermonde matrix also o ers a convenient way to switch back and forth between modal and nodal representation. However, one important observation to make from Equa- tion B.21 is that the Vandermonde matrix must well conditioned enough to be invertible without penalty. Because we have chosen to use orthonormal basis, the conditioning of the Vandermonde matrix is entirely dependent on nodal locations. This leads into a discussion on quadrature. ~uh := uh(~ ) = Veu = bu (B.21) Quadrature is a weighted sum of function values at predetermined points within a standard domain. Numerous quadrature rules exist for approximating the exact integral of 67 a function. In general, this approximation would look like Z 1 1 f( )d ? NpX j=1 !jf( j) (B.22) where ! is the quadrature weights, and is the quadrature nodal locations. For Gaussian quadrature, nodal locations will fall in f 2P2Np 1 with 1 < 1!Np < 1. This quadrature rule, while well conditioned, exlcudes the endpoints of the standard interval, so we choose instead to use Gauss-Lobatto quadrature, nodal locations will fall in f2P2Np 3 with 1 = 1 and Np+1 = 1. Mathematically, this is expressed as Equation B.23. A sample of Gauss- Lobatto quadrature cells are shown in Figure B.3 for varying degree polynomials. Z 1 1 f( )d = !0f( 0) + NpX i=1 !if( i) +!Np+1f( Np+1) (B.23) Looking back now at the original problem to solve, @u @t + @f(u) @x = 0 (B.24) and rewriting in descretized form, we nd Z Dk bki @u k h @t dx Z Dk @bki @kf k hdx+ I @Dk bkif uk; h ;uk;+h ;nk; x ds = 0 (B.25) In this form, it is apparent that discontinuous Galerkin is a local weighted residual statement for each element based on the standard element. In order for this form to be useful, a transformation from the standard element to the physical element is required. This can be achieved using the following set of equations. x2Dk :: x =Mk( ) (B.26) Mk( ) = xka + 1 + 2 xkb xka (B.27) x2Dk :: dx =Jkd ; Jk = x k b x k a 2 (B.28) where Jk is the Jacobian of the transformation. Using the coordinate transformation, the state decomposition can be rewritten as ukh(x) = uh Mk( ) = NX j=0 eukj(t)bj( ) (B.29) 68 (a) d = 2 (b) d = 3 (c) d = 4 (d) d = 5 (e) d = 6 (f) d = 7 Figure B.3: Lobatto quadrature for varying degrees of freedom 69 It is also assumed that the ux can be decomposed in the same manner. Now, re-examining Equation B.25, the time derivative term can be rewritten as Z Dk bki @u k n @t dx = Z Dk bki (x) @@t NX j=0 eukj(t)bkj(x) ! dx = NX j=0 Z Dk bki (x)bkj(x)dx deuk j dt = NX j=0 Z 1 1 bi( )bj( )Jkd deuk j dt =Jk NX j=0 Mijdeu k j dt (B.30) The mass matrix, Mij, is an operator that can determined before time evolution begins rather than at each iteration. Following a similar process as the time derivative, the spatial derivative term can be rewritten as Z Dk dbki dx = NX j=0 ST ij efkj (B.31) where ST is the sti ness matrix and is de ned as Sij := Z 1 1 bi( )dbj( )d d (B.32) The numerical ux term in Equation B.25 must be handled signi cantly di erent. Because the setup thus far has allowed for discontinuities at the boundaries of cells, some numerical ux function must be used in order to interpret jumps at the boundaries. Writing this ux as an expansion f (s) = f uk; h ;uk;+h ;nk; x s = NX j=0 ef skjbj(s) (B.33) It can be shown that using this expansion, the cell boundary term can be written as I @Dk bkif uk; h ;uk;+h ;nk; x ds = NX m=0 (La)im ef akm NX n=0 (Lb)im ef bkn (B.34) where La and Lb are lifting matrices. The lifting matrices are called so because they lift lower dimensional boundary information to higher dimensional element data. So in matrix vector form, Equation B.25 can be written as JkMdeu k dt S T efk +Laef a k +L bef b k = 0 (B.35) 70 The matrix vector form (element operator form) given above is the useful form of DG. For more information regarding formulation of the operator matrices and time evolution using DG, refer to [57, 56, 58, 59, 60]. 71 Appendix C Comparison of First-Order Finite Element to Discontinuous Galerkin: A Classic Fluid Dynamics Example The two-dimensional isentropic vortex problem is often posed in comparing numerical methods because it is well behaved and smooth and the exact solution is known. The isentropic vortex is created by introducing perturbations to the velocity and temperature in a uniform ow eld. The exact solution to the convection of a vortex is such that the vortex strength will be maintained but will move with a bulk velocity. The given freestream conditions for this problem are as follows: u1 = M1cos (C.1) v1 = M1sin (C.2) p1 = 1 (C.3) T1 = 1 (C.4) For the problem given here, the freestream Mach number was set to 0.5 and the speci c heat ratio was set to 1.4. This results in a freestream density of 1 and an acoustic speed of 1. The perturbations to develop the initial vortex are proportional (in a way) to the distance from the vortex origin (set at the domain origin for simplicity). Perturbations are then calculated as follows: u = yf (C.5) v = + xf (C.6) T = 12 f2 (C.7) where ( x; y) = (x xo;y yo) (C.8) f = b2 exp a(1 r2) (C.9) r2 = x2 + y2 (C.10) 72 The vortex strength parameters were set as a = 0:5, b = 5. The computational mesh set on the domain (x;y)2( 5;5) ( 5;5). Using nite element method, the mesh size was set to 50 50, while for DG, the mesh size requirement was only 4 4. The vortex was set to convect, and the boundary conditions were set as periodic so that as the vortex leaves the domain on say the left boundary, it re-enters the domain on the right boundary. With the domain and Mach number set as speci ed, and a time span of 20 seconds, the nal location of the center of the vortex should be exactly the initial location. Results for the nite element method are shown in Figure C.1. ?6 ?4 ?2 0 2 4 6 ?6 ?4 ?2 0 2 4 6 0.4 0.5 0.6 0.7 0.8 0.9 1 (a) Initial Density ?6 ?4 ?2 0 2 4 6 ?6 ?4 ?2 0 2 4 6 0.5 0.6 0.7 0.8 0.9 1 (b) Final Density (c) Initial Velocity Field (d) Final Velocity Field Figure C.1: Euler Solution to Isentropic Vortex using Finite Element From Figures C.1a and C.1b, it is apparent that the rst order scheme for time evolution was unable to fully resolve the vortex after 20 seconds. Looking at the scale on the z-axis, the vortex has dissipated (non-physical) and there appear to be slight wrinkles in the solution for the density. From the velocity elds (Figures C.1c and C.1d), it appears as though some 73 of the velocity has also dissipated. Looking ahead at the results for DG, it is apparent that the features of the ow were maintained accurately with little to no dissipation (Figure C.2). It should be noted that in both Figures C.2c and C.2d, the solution appears di erent from those in Figure C.1. The artifacts present are a function of plotting only. The important note to take away from the results is the similarities between the initial and nal states. ?5 0 5 ?5 0 5 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 (a) Initial Density ?5 0 5 ?5 0 5 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 (b) Final Density (c) Initial Velocity Field (d) Final Velocity Field Figure C.2: Euler Solution to Isentropic Vortex using Discontinuous Galerkin 74 Appendix D DG/LSM Fortran Source Main Evolution Source program main INTEGER : : NI , NJ, d ,Np, Np2 , steps , kk , i , j , k , itmax , kchek , islow , output , n stage INTEGER,ALLOCATABLE,DIMENSION( : ) : : East , West , North , South INTEGER,ALLOCATABLE,DIMENSION( : ) : : Ekx ,Wkx, Nky , Sky INTEGER,ALLOCATABLE,DIMENSION( : , : ) : : Kx,KY REAL : : epsilon , staring , ending , start time , end time , pi , burnstop , res sum , res1 , lambda REAL : : CFL, dx , dy , ex , ey ,Hmax, Hmin , Jx , Jy , Jk , dstep , y , xmax , xmin , ymax , ymin , l r e f ,gamma, alpha (3 ,3) REAL,ALLOCATABLE,DIMENSION( : ) : : ww,www, p1 , p2 , q1 , q2 , delta , umode , Hmode, ypos , area , rhs , res , cmode REAL,ALLOCATABLE,DIMENSION( : , : ) : : Sx , Sy ,V,Bnm REAL,ALLOCATABLE,DIMENSION( : , : , : ) : : Lx , Ly , phx , phy ,Ham, Hea , caseheavi REAL,ALLOCATABLE,DIMENSION( : , : , : ) : : u , uold CHARACTERa4220 : : DIRNAME, fname if ( islow . eq . 1 ) call system ( ? c l e a n s e ? ) ! clean output folder ! read in input settings open(unit=11, file=? i n p u t s e t t i n g s . txt ? ) read(11 ,a42) NI , NJ read(11 ,a42) d read(11 ,a42) e p s i l o n read(11 ,a42)CFL read(11 ,a42) xmin , ymin read(11 ,a42)xmax , ymax read(11 ,a42) ex , ey read(11 ,a42)Hmax, Hmin read(11 ,a42) l r e f ,gamma read(11 ,a42) burnstop , itmax read(11 ,a42) islow , output close (11) ! determine necessary domain settings dx = (xmax xmin )/NI dy = (ymax ymin )/NJ l r e f=l r e f a42min(dx , dy )/da42a422.0 Np = d+1; Np2 = 4a42Np 2; Jx = 0.5a42dy ; Jy = 0.5a42dx ; Jk = .25a42dxa42dy ; lambda=1.0+gamma dstep = CFL/( lambdaa42da42a422./ min (dx , dy ) ) e p s i l o n=e p s i l o n a42min (dx , dy )/( d+1) steps = i n t ( min( dxa42NI , dya42NJ)/ dstep ) pi=acos ( 1.0); ! a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42 ! ALLOCATE MEMORY FOR MATRICES 75 ALLOCATE( ypos ( steps ) ) ; ypos =0; ALLOCATE( area ( steps ) ) ; area =0; ALLOCATE( r e s ( steps ) ) ; r e s =0; ALLOCATE( Sx (Npa42a422 ,Npa42a422 ) ) ; Sx=0; ALLOCATE( Sy (Npa42a422 ,Npa42a422 ) ) ; Sy=0; ALLOCATE(Kx(Np , 2 ) ) ; Kx=0; ALLOCATE(Ky(Np , 2 ) ) ; Ky=0; ALLOCATE(Lx(Npa42a422 ,Np , 2 ) ) ; Lx=0; ALLOCATE(Ly(Npa42a422 ,Np , 2 ) ) ; Ly=0; ALLOCATE(Bnm(Np2a42a422 ,Npa42a422 ) ) ;Bnm=0; ALLOCATE(ww(Np2a42a422 ) ) ;ww=0; ALLOCATE(www(Npa42a422 ) ) ;www=0; ALLOCATE(u(Npa42a422 ,NI , NJ ) ) ; u=0; ALLOCATE( uold (Npa42a422 ,NI , NJ ) ) ; uold=u ; ALLOCATE( caseheavi (Npa42a422 ,NI , NJ ) ) ; caseheavi =0; ALLOCATE(umode(Np2a42a422 ) ) ; umode=0; ALLOCATE(Hmode(Np2a42a422 ) ) ;Hmode=0; ALLOCATE( cmode (Np2a42a422 ) ) ; cmode=0; ALLOCATE( de l ta (Np2a42a422 ) ) ; de l ta =0; ALLOCATE( East (NI ) ) ; East =0; ALLOCATE(West(NI ) ) ; West=0; ALLOCATE( North (NJ ) ) ; North=0; ALLOCATE( South (NJ ) ) ; South=0; ALLOCATE(Ekx(NI ) ) ; Ekx=1; ALLOCATE(Wkx(NI ) ) ;Wkx=2; ALLOCATE(Nky(NJ ) ) ; Nky=1; ALLOCATE( Sky (NJ ) ) ; Sky=2; ALLOCATE( p1 (Npa42a422 ) ) ; p1=0; ALLOCATE( p2 (Npa42a422 ) ) ; p2=0; ALLOCATE( q1 (Npa42a422 ) ) ; q1=0; ALLOCATE( q2 (Npa42a422 ) ) ; q2=0; ALLOCATE( phx (Npa42a422 ,NI , NJ ) ) ; phx=0; ALLOCATE( phy (Npa42a422 ,NI , NJ ) ) ; phy=0; ALLOCATE(Ham(Npa42a422 ,NI , NJ ) ) ;Ham=0; ALLOCATE(Hea(Npa42a422 ,NI , NJ ) ) ; Hea=0; ALLOCATE( rhs (Npa42a422 ) ) ; rhs =0; ! a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42a42 if (Np. l t . 1 0 ) write(DIRNAME, 9 1 ) ?DG ? ,Np if (Np. ge . 1 0 ) write(DIRNAME, 9 2 ) ?DG ? ,Np call chdir (DIRNAME) ! read in stiffness matrix open(unit=22, file=? S t i f f n e s s . dat ? ) do i =1,Npa42a422 read( 2 2 ,a42) ( Sx ( i , j ) , j =1,Npa42a422) enddo Sx=Jx/Jka42Sx do i =1,Npa42a422 read( 2 2 ,a42) ( Sy ( i , j ) , j =1,Npa42a422) enddo close (22) Sy=Jy/Jka42Sy ! read in Lifting Matrices open(unit=23, file=? LiftingX . dat ? ) do k=1,2 do i =1,Npa42a422 read( 2 3 ,a42) ( Lx( i , j , k ) , j =1,Np) enddo enddo close (23) Lx=Jx/Jka42Lx 76 open(unit=24, file=? LiftingY . dat ? ) do k=1,2 do i =1,Npa42a422 read( 2 4 ,a42) ( Ly( i , j , k ) , j =1,Np) enddo enddo close (24) Ly=Jy/Jka42Ly ! read in tall vandermonde matrix open(unit=25, file=? TallVandermonde . dat ? ) do i =1,Np2a42a422 read( 2 5 ,a42) (Bnm( i , j ) , j =1,Npa42a422) enddo close (25) ! read in weights for tall vandermonde matrix open(unit=26, file=? TallWeights . dat ? ) do i =1,Np2a42a422 read(26 ,a42)ww( i ) enddo close (26) ! read in weights for regular lobatto points open(unit=26, file=? LobWeights . dat ? ) do i =1,Npa42a422 read(26 ,a42)www( i ) enddo close (26) ! read in initial state matrix call chdir ( ? . . ? ) open(unit=27, file=? i n i t i a l s t a t e . dat ? ) do i =1,NI do j =1,NJ read( 2 7 ,a42) ( u(k , i , j ) , k=1,Npa42a422) enddo enddo close (27) open(unit=28, file=? c a s e h e a v i . dat ? ) do i =1,NI do j =1,NJ read( 2 8 ,a42) ( caseheavi (k , i , j ) , k=1,Npa42a422) enddo enddo close (28) ! Setup north , south , east , and west boundaries and nodes do i =1,Np Kx( i ,1)=1+( i 1)a42Np; Kx( i ,2)=( i 1)a42Np+Np; Ky( i ,1)= i ; Ky( i ,2)=Npa42a422 (Np i ) enddo do i =1,NI East ( i )= i+1 West( i )=i 1 enddo do j =1,NJ North ( j )= j+1 77 South ( j )=j 1 enddo East (NI)=NI ; West(1)=1 North (NJ)=NJ South (1)=1 Ekx(NI)=2; Wkx(1)=1; Nky(NJ)=2; Sky (1)=1; uold=u ; ! set up alpha for staging RK TVD n stage =3; alpha (1 ,:)=(/ 1.0 , 0.75 , 1 . 0 / 3 . 0 /) alpha (2 ,:)=(/ 0.0 , 0.25 , 2 . 0 / 3 . 0 /) alpha (3 ,:)=(/ 1.0 , 0.25 , 2 . 0 / 3 . 0 /) call cpu time ( s t a r t t i m e ) ! wallclock y=0; kk=1; kchek =10; ! begin time evolution do while ( dstepa42kk . l e . burnstop ) res sum =0; ! begin stage integration do i j k =1, n stage ! perform single time evolution here do j =1,NJ do i =1,NI ! first derivatives p1= matmul(Sx , u ( : , i , j ))+matmul(Lx ( : , : , 2 ) , u(Kx( : , Ekx( i ) ) , East ( i ) , j ))& & matmul(Lx ( : , : , 1 ) , u(Kx( : , 1 ) , i , j ) ) ; p2= matmul(Sx , u ( : , i , j ))+matmul(Lx ( : , : , 2 ) , u(Kx( : , 2 ) , i , j ))& & matmul(Lx ( : , : , 1 ) , u(Kx( : ,Wkx( i ) ) , West( i ) , j ) ) ; q1= matmul(Sy , u ( : , i , j ))+matmul(Ly ( : , : , 2 ) , u(Ky( : , Nky( j ) ) , i , North ( j )))& & matmul(Ly ( : , : , 1 ) , u(Ky( : , 1 ) , i , j ) ) ; q2= matmul(Sy , u ( : , i , j ))+matmul(Ly ( : , : , 2 ) , u(Ky( : , 2 ) , i , j ))& & matmul(Ly ( : , : , 1 ) , u(Ky( : , Sky ( j ) ) , i , South ( j ) ) ) ; p1=p1 ; p2=p2 q1=q1 ; q2=q2 phx ( : , i , j )=0.5a42( p1+p2 ) ; phy ( : , i , j )=0.5a42( q1+q2 ) ; !form Hamiltonian Ham( : , i , j )=(phx ( : , i , j )a42a422.+phy ( : , i , j )a42a422 . )a42a420 . 5 ; Ham( : , i , j )=Ham( : , i , j )+gammaa42u ( : , i , j )/ s q r t (u ( : , i , j )a42a422.+ l r e f a42a422 . )a42(Ham( : , i , j ) 1.); Ham( : , i , j )=Ham( : , i , j ) 0.5a42(p1 p2 ) 0.5a42( q1 q2 ) ; ! calculate area if ( i j k . eq . 1 ) then umode=matmul(Bnm, u ( : , i , j ) ) ; Hmode=matmul(Bnm,Ham( : , i , j ) ) ; cmode=matmul(Bnm, caseheavi ( : , i , j ) ) ; d el ta =0; do k=1,Np2a42a422 if ( abs (umode( k ) ) . l e . e p s i l o n ) d el t a ( k)=0.5/ e p s i l o n +0.5/ e p s i l o n& &a42cos ( pi a42umode( k )/ e p s i l o n ) ; enddo area ( kk)=area ( kk)+Jka42sum( cmodea42Hmodea42 d el ta a42ww) endif 78 enddo enddo do j =1,NJ do i =1,NI ! second derivatives rhs= Ham( : , i , j ) if ( any ( abs ( phx ( : , i , j ) ) . gt .Hmax) )then p1=matmul(Lx ( : , : , 2 ) , 0 . 5 a42 ( phx (Kx( : , Ekx( i ) ) , East ( i ) , j )+phx (Kx( : , 2 ) , i , j )))& & matmul(Lx ( : , : , 1 ) , 0 . 5 a42 ( phx (Kx( : , 1 ) , i , j )+phx (Kx( : ,Wkx( i ) ) , West( i ) , j ) ) ) ; p1=(p1 matmul(Sx , phx ( : , i , j ) ) ) ; rhs=rhs+exa42p1 ; elseif ( any ( abs ( phx ( : , i , j ) ) . l t . Hmin) )then p1=matmul(Lx ( : , : , 2 ) , 0 . 5 a42 ( phx (Kx( : , Ekx( i ) ) , East ( i ) , j )+phx (Kx( : , 2 ) , i , j )))& & matmul(Lx ( : , : , 1 ) , 0 . 5 a42 ( phx (Kx( : , 1 ) , i , j )+phx (Kx( : ,Wkx( i ) ) , West( i ) , j ) ) ) ; p1=(p1 matmul(Sx , phx ( : , i , j ) ) ) ; rhs=rhs+exa42p1 ; endif if ( any ( abs ( phy ( : , i , j ) ) . gt .Hmax) )then q1=matmul(Ly ( : , : , 2 ) , 0 . 5 a42 ( phy (Ky( : , Nky( j ) ) , i , North ( j ))+phy (Ky( : , 2 ) , i , j )))& & matmul(Ly ( : , : , 1 ) , 0 . 5 a42 ( phy (Ky( : , 1 ) , i , j )+phy (Ky( : , Sky ( j ) ) , i , South ( j ) ) ) ) ; q1=(q1 matmul(Sy , phy ( : , i , j ) ) ) ; rhs=rhs+eya42q1 elseif ( any ( abs ( phy ( : , i , j ) ) . l t . Hmin) )then q1=matmul(Ly ( : , : , 2 ) , 0 . 5 a42 ( phy (Ky( : , Nky( j ) ) , i , North ( j ))+phy (Ky( : , 2 ) , i , j )))& & matmul(Ly ( : , : , 1 ) , 0 . 5 a42 ( phy (Ky( : , 1 ) , i , j )+phy (Ky( : , Sky ( j ) ) , i , South ( j ) ) ) ) ; q1=(q1 matmul(Sy , phy ( : , i , j ) ) ) ; rhs=rhs+eya42q1 endif ! calculate residual if ( i j k . eq . 3 ) res sum=s q r t ( res suma42a422+sum( rhsa42wwwa42rhs ) ) ; u ( : , i , j )=alpha (1 , i j k )a42 uold ( : , i , j )+alpha (2 , i j k )a42u ( : , i , j )+alpha (3 , i j k )a42 dstepa42rhs enddo enddo enddo call cpu time ( end time ) ! end single time evolution , perform checks , write output if necessary if ( islow . eq . 1 ) then if (mod( kk , output ) . eq . 0 )then call chdir ( ? output ? ) write( fname , 9 4 ) ? grain output ? , kk , ? . dat ? 94 format( a12 , i4 , a4 ) open(unit=41, file=fname ) write(41 ,a42) NI , NJ, d write(41 ,a42) xmin , xmax , ymin , ymax do j =1,NJ do i =1,NI write ( 4 1 ,a42) ( uold (k , i , j ) , k=1,Npa42a422) enddo enddo close (41) call chdir ( ? . . ? ) endif endif if ( kk . eq . 1 )then res1=res sum endif r e s ( kk)=res sum / res1 ; ypos ( kk)=kka42dstep ; if ( kk . eq . kchek )then write(a42 ,a42) kk , ypos ( kk ) , r e s ( kk ) 79 kchek=kchek +50; endif if ( r e s ( kk ) . gt . 2 )goto 51 kk=kk+1; uold=u ; enddo 51 write(a42 ,a42) ?Number of Steps : ? , kk 1 write(a42 ,a42) ? Burn Distance : ? , ypos ( kk 1) write(a42 ,a42) ? Burn Area : ? , area ( kk 1) write(a42 ,a42) ? Residual : ? , r e s ( kk 1) open(unit=81, file=? r e s i d u a l . dat ? ) write(81 ,a42) ? Variables = " t " ," r e s i d u a l " ," area " ? do i =1,kk 1 write(81 ,a42) ypos ( i ) , r e s ( i ) , area ( i ) enddo close (81) write(a42 ,93) ?Run Time = ? , end time start time , ? seconds ? ! write out final state matrix open(unit=31, file=? f i n a l s t a t e . dat ? ) do i =1,NI do j =1,NJ write ( 3 1 ,a42) ( u(k , i , j ) , k=1,Npa42a422) enddo enddo close (31) open(unit=31, file=? p o s t p r o c e s s s p e c s . dat ? ) write(31 ,a42) NI , NJ, d , xmin , xmax , ymin , ymax close (31) 91 FORMAT( a3 , i 1 ) 92 FORMAT( a3 , i 2 ) 93 FORMAT( a11 , 1 x , f14 . 8 , 1 x , a8 ) 99 FORMAT( e16 . 8 , 2 x , e16 . 8 , 2 x , e16 . 8 ) END Initialization Source program s h e l l ! shell program used for setting up and running level set propagation program INTEGER : : NI , NJ, d ,Np, Np2 , i , j , k , Ns REAL : : e p s i l o n c a s e , xmin , xmax , ymin , ymax , dump, xs , ys , dx , dy , ro , ri , lambda REAL : : beta , pi , f i l l e t , eps ang REAL,ALLOCATABLE,DIMENSION( : ) : : qx , qy , px , py , xi , yi REAL,ALLOCATABLE,DIMENSION( : , : ) : : xdomain , ydomain REAL,ALLOCATABLE,DIMENSION( : , : , : ) : : phi , heavi , c i CHARACTERa4220 : : gtype , ctype ! phi is for the grain ! heavi is for the case CHARACTERa4220 : : DIRNAME interface function s e t u p g r a i n ( xs , ys , dx , dy , l i s t x , l i s t y ,Np, Ns , gtype , xc , yc , c i )& 80 result ( d i s t ) real ,dimension(Npa42a422 ) : : l i s t x , l i s t y , d i s t real ,dimension(Ns , 2 , 4 ) : : c i real : : xs , ys , dx , dy , xc (Nsa426+1) , yc (Nsa426+1) integer : : Np, Ns charactera4220 : : gtype end function s e t u p g r a i n end interface interface function setup case ( xs , ys , dx , dy , l i s t x , l i s t y ,Np, ro , ctype , eps)& result ( heav ) real ,dimension(Npa42a422) : : l i s t x , l i s t y , heav real : : xs , ys , dx , dy , x , y , ro , eps charactera4220 : : ctype end function setup case end interface open(unit=11, file=? i n p u t s e t t i n g s . txt ? ) read(11 ,a42) NI , NJ read(11 ,a42) d read(11 ,a42)dump, e p s i l o n c a s e read(11 ,a42)dump read(11 ,a42) xmin , ymin read(11 ,a42)xmax , ymax close (11) open(unit=12, file=? design . txt ? ) read(12 ,a42) ro read(12 ,a42) rp read(12 ,a42) r i read(12 ,a42) lambda read(12 ,a42) Ns read(12 ,a42) f i l l e t read(12 ,a42) eps ang read(12 ,a42) gtype read(12 ,a42) ctype dx=(xmax xmin )/NI dy=(ymax ymin )/NJ Np=d+1 Np2=2a42Np 2 e p s i l o n c a s e=e p s i l o n c a s e a42min(dx , dy )/( d+1.0) pi=acos ( 1.0) !ALLOCATE MEMORY FOR STATE AND CASE MATRICES ALLOCATE( xdomain (NI+1,NJ+1)); xdomain=0; ALLOCATE( ydomain (NI+1,NJ+1)); ydomain=0; ALLOCATE( phi (Npa42a422 ,NI , NJ ) ) ; phi =0; ALLOCATE( heavi (Npa42a422 ,NI , NJ ) ) ; heavi =0; ALLOCATE( qx (Npa42a422 ) ) ; qx=0; ALLOCATE( qy (Npa42a422 ) ) ; qy=0; ALLOCATE( xi (Nsa426+1)); xi =0; ALLOCATE( yi (Nsa426+1)); yi =0; ALLOCATE( c i (Ns , 2 , 4 ) ) ; c i =0; !ALLOCATE(px(Np2a42a422));px=0; !ALLOCATE(py(Np2a42a422));py=0; if (Np. l t . 1 0 ) write(DIRNAME, 9 1 ) ?DG ? ,Np if (Np. ge . 1 0 ) write(DIRNAME, 9 2 ) ?DG ? ,Np call chdir (DIRNAME) open(unit=12, file=? LobattoPoints . dat ? ) do i =1,Npa42a422 read(12 ,a42) qx ( i ) , qy ( i ) enddo close (12) 81 !open( unit=13, file =?GaussPoints . dat ?) !do i=1,Np2a42a422 ! read(13 ,a42)px( i ) ,py( i ) !enddo ! close (13) call chdir ( ? . . ? ) ! set initial condition do j =1,NJ+1 do i =1,NI+1 xdomain ( i , j )=xmin+(i 1)a42dx ydomain ( i , j )=ymin+(j 1)a42dy enddo enddo ! setup initial discretized surface if ( gtype . eq . ? barrere ? )then call barrere ( rp , ri , f i l l e t , eps ang , Ns , xi , yi , c i ) elseif ( gtype . eq . ? s s t a r ? )then call hard star ( rp , ri , Ns , xi , yi ) elseif ( gtype . eq . ? h e l l f i r e ? ) then xi (1)= ro ; yi (1)=lambda ; endif ! setup initial state open(unit=21, file=? i n i t i a l s t a t e . dat ? ) open(unit=22, file=? c a s e h e a v i . dat ? ) do i =1,NI do j =1,NJ xs=xdomain ( i , j ) ys=ydomain ( i , j ) phi ( : , i , j )= s e t u p g r a i n ( xs , ys , dx , dy , qx , qy ,Np, Ns , gtype , xi , yi , c i ) ; heavi ( : , i , j )= setup case ( xs , ys , dx , dy , qx , qy ,Np, ro , ctype , e p s i l o n c a s e ) ; write ( 2 1 ,a42) ( phi (k , i , j ) , k=1,Npa42a422) write ( 2 2 ,a42) ( heavi (k , i , j ) , k=1,Npa42a422) enddo enddo close (21) close (22) 91 FORMAT( a3 , i 1 ) 92 FORMAT( a3 , i 2 ) ENDPROGRAM subroutine barrere ( rp , ri , f , eps , Ns , xi , yi , c i ) real : : rp , ri , f , eps , xi (6a42Ns+1) , yi (6a42Ns+1) ,H1, theta2 integer : : Ns , i , j , k , i i , j j , kk real : : x1 ( 4 , 2 ) , c1 ( 2 ) , c2 ( 2 ) , c3 ( 2 ) , c4 (2) real : : x ( 6 ) , y ( 6 ) , xr ( 6 ) , yr ( 6 ) , c i (Ns , 2 , 4 ) , ang , pi pi=acos ( 1.0) H1=rpa42 s i n ( pi a42eps /Ns) theta2=atan (H1a42tan ( pi a42eps /Ns )/(H1 r i a42tan ( pi a42eps /Ns ) ) ) x1 (1 ,1)=( rp+f )a42 cos ( pi /Ns ) ; x1(1 ,2)= ( rp+f )a42 s i n ( pi /Ns ) ; x1 (2 ,1)=( rp+f )a42 cos ( pi a42eps /Ns ) ; x1(2 ,2)= ( rp+f )a42 s i n ( pi a42eps /Ns ) ; x1 (3 ,1)= rpa42cos ( pi a42eps /Ns)+ f a42cos( pi/2+theta2 ) ; x1(3 ,2)= ( rpa42 s i n ( pi a42eps /Ns)+ f a42 s i n ( pi/2+theta2 ) ) ; x1 (4 ,1)= r i+f a42cos ( theta2 ) ; x1 (4 ,2)=0; c1 =(/0 ,0/) c2=(/rpa42cos ( pi a42eps /Ns), rpa42 s i n ( pi a42eps /Ns )/) c3=(/c2 (1) , c2 (2 )/ ) c4=(/c1 (1) , c1 (2 )/ ) k=1; do j =1,4 82 x ( k)=x1 ( j , 1 ) ; y ( k)=x1 ( j , 2 ) ; k=k+1 enddo do j =3,2, 1 x ( k)=x1 ( j , 1 ) y ( k)= x1 ( j , 2 ) k=k+1 enddo k=1 do i i =1,Ns ang = ( i i 1.0)a422.0a42 pi /Ns xr=xa42cos ( ang) ya42 s i n ( ang ) yr=xa42 s i n ( ang)+ya42cos ( ang ) c i ( i i ,1 ,1)= c1 (1)a42 cos ( ang) c1 (2)a42 s i n ( ang ) c i ( i i ,2 ,1)= c1 (1)a42 s i n ( ang)+c1 (2)a42 cos ( ang ) c i ( i i ,1 ,2)= c2 (1)a42 cos ( ang) c2 (2)a42 s i n ( ang ) c i ( i i ,2 ,2)= c2 (1)a42 s i n ( ang)+c2 (2)a42 cos ( ang ) c i ( i i ,1 ,3)= c3 (1)a42 cos ( ang) c3 (2)a42 s i n ( ang ) c i ( i i ,2 ,3)= c3 (1)a42 s i n ( ang)+c3 (2)a42 cos ( ang ) c i ( i i ,1 ,4)= c4 (1)a42 cos ( ang) c4 (2)a42 s i n ( ang ) c i ( i i ,2 ,4)= c4 (1)a42 s i n ( ang)+c4 (2)a42 cos ( ang ) do j j =1,6 xi ( k)=xr ( j j ) yi ( k)=yr ( j j ) k=k+1; enddo enddo xi ( k)=x (1) yi ( k)=y (1) return end subroutine subroutine hard star ( rp , ri , Ns , xi , yi ) real : : rp , ri , xi (6a42Ns+1) , yi (6a42Ns+1) , beta , dbeta integer : : Ns , i , j , k pi=acos ( 1.0) xi (1)= rp ; yi (1)=0 k=2; beta =0; dbeta=pi /Ns do i =1,Ns beta=beta+dbeta xi ( k)= r i a42cos ( beta ) ; yi ( k)= r i a42 s i n ( beta ) k=k+1 beta=beta+dbeta xi ( k)=rpa42cos ( beta ) ; yi ( k)=rpa42 s i n ( beta ) k=k+1 enddo return end subroutine function s e t u p g r a i n ( xs , ys , dx , dy , l i s t x , l i s t y ,Np, Ns , gtype , xc , yc , c i )& result ( d i s t ) real ,dimension(Npa42a422) : : l i s t x , l i s t y , d i s t real ,dimension(Ns , 2 , 4 ) : : c i real : : xs , ys , dx , dy , x , y , ro , ri , lambda , xc (Nsa426+1) , yc (Nsa426+1) integer : : Np, Ns integer : : i , j , k charactera4220 : : gtype interface function s t r a i g h t s t a r (x , y , Ns , xc , yc ) result ( phi ) 83 real : : x , y , phi , xc (300) , yc (300) integer : : Ns end function s t r a i g h t s t a r end interface do i =1,Npa42a422 x=xs +0.5a42( l i s t x ( i )+1.0)a42dx y=ys +0.5a42( l i s t y ( i )+1.0)a42dy if ( gtype . eq . ? s s t a r ? ) then d i s t ( i )= s t r a i g h t s t a r (x , y , Ns , xc , yc ) elseif ( gtype . eq . ? barrere ? ) then d i s t ( i )= d i s t t o b a r r e r e (x , y , Ns , xc , yc , c i ) elseif ( gtype . eq . ? h e l l f i r e ? ) then d i s t ( i )= d i s t t o h e l l f i r e (x , y , xc ( 1 ) , yc ( 1 ) ) endif enddo return end function function d i s t t o h e l l f i r e (x , y , ro , lambda ) result ( phi ) integer i , j , k real : : x , y , ro , lambda , r r=s q r t ( xa42a422.0+y a42a422.0) if ( r . gt . 0 . 5a42 ro )then phi=r 0.5a42ro 0.5a42lambda else phi=(1 r ) 0.5a42ro 0.5a42lambda endif end function function d i s t t o b a r r e r e (x , y , Ns , xc , yc , c i ) result ( phi ) integer : : Ns , i , j , k real : : x , y , phi , xc (Nsa426+1) , yc (Nsa426+1) , c i (Ns , 2 , 4 ) real : : p0 ( 2 ) , p1 ( 2 ) , p2 ( 2 ) , p ( 2 ) , c ( 2 ) , d(Nsa426) , s (Nsa426) , dmin , s i g p=(/x , y /) k=1 do i =1,Ns p0=(/ c i ( i , 1 , 1 ) , c i ( i , 2 , 1 ) / ) p1=(/xc ( k ) , yc ( k )/) p2=(/xc ( k+1) , yc ( k+1)/) call d i s t t o a r c (p , p0 , p1 , p2 , d( k ) , s ( k ) ) k=k+1 p0=(/ c i ( i , 1 , 2 ) , c i ( i , 2 , 2 ) / ) p1=(/xc ( k ) , yc ( k )/) p2=(/xc ( k+1) , yc ( k+1)/) call d i s t t o a r c (p , p0 , p1 , p2 , d( k ) , s ( k ) ) k=k+1 p1=(/xc ( k ) , yc ( k )/) p2=(/xc ( k+1) , yc ( k+1)/) call d i s t t o s e g (p , p1 , p2 , d( k ) , s ( k ) ) k=k+1 p1=(/xc ( k ) , yc ( k )/) p2=(/xc ( k+1) , yc ( k+1)/) call d i s t t o s e g (p , p1 , p2 , d( k ) , s ( k ) ) k=k+1 p0=(/ c i ( i , 1 , 3 ) , c i ( i , 2 , 3 ) / ) p1=(/xc ( k ) , yc ( k )/) 84 p2=(/xc ( k+1) , yc ( k+1)/) call d i s t t o a r c (p , p0 , p1 , p2 , d( k ) , s ( k ) ) k=k+1 p0=(/ c i ( i , 1 , 4 ) , c i ( i , 2 , 4 ) / ) p1=(/xc ( k ) , yc ( k )/) p2=(/xc ( k+1) , yc ( k+1)/) call d i s t t o a r c (p , p0 , p1 , p2 , d( k ) , s ( k ) ) k=k+1 enddo dmin=minval (d ) ; where(d /= dmin ) s =0.0 end where s i g=sum( s ) phi=dmina42sgn ( s i g ) end function function s t r a i g h t s t a r (x , y , Ns , xc , yc ) result ( phi ) integer : : Ns real : : x , y , phi , ro , ri , lambda , xc (300) , yc (300) integer : : i , j , k real : : x0 ( 2 ) , x1 ( 2 ) , x2 ( 2 ) , d(2a42Ns ) , s (2a42Ns ) , dmin ! interface ! subroutine dist to seg (p,p0 ,p1 ,d, s) ! real , intent (in ):: p(2) ,p0(2) ,p1(2) ! real , intent (out )::d, s ! end subroutine !end interface x0=(/x , y /) do i =1,2a42Ns x1=(/xc ( i ) , yc ( i )/) x2=(/xc ( i +1) , yc ( i +1)/) call d i s t t o s e g ( x0 , x1 , x2 , d( i ) , s ( i ) ) enddo dmin=minval (d ) ; where(d /= dmin ) s =0.0 end where s i g=sum( s ) phi=dmina42sgn ( s i g ) return end function ! function setup case ( xs , ys , dx , dy , l i s t x , l i s t y ,Np, ro , ctype , eps)& result ( heav ) real ,dimension(Npa42a422) : : l i s t x , l i s t y , heav real : : xs , ys , dx , dy , x , y , ro , eps integer : : Np integer : : i , j , k charactera4220 : : ctype interface function c i r c l e (x , y , ro , eps ) result ( heav ) real : : x , y , heav , ro , eps end function c i r c l e end interface interface function square (x , y , ro , eps ) result ( heav ) real : : x , y , heav , ro , eps end function square 85 end interface do i =1,Npa42a422 x=xs +0.5a42( l i s t x ( i )+1.0)a42dx y=ys +0.5a42( l i s t y ( i )+1.0)a42dy if ( ctype . eq . ? c i r c l e ? ) then heav ( i )= c i r c l e (x , y , ro , eps ) elseif ( ctype . eq . ? square ? )then heav ( i )=square (x , y , ro , eps ) else write(a42 ,a42) ? Check case type : ? , ctype endif enddo return end function function c i r c l e (x , y , ro , eps ) result ( heav ) real : : x , y , heav , ro , phi , eps , pi pi=acos ( 1.0) phi= (s q r t ( xa42a422.0+ya42a422.0) ro ) if ( phi . gt . 0 . 0 ) then heav =1.0 elseif ( phi . gt . eps )then heav =0.5a42(1+ cos ( phia42 pi / eps ) ) else heav =0.0; endif return end function function square (x , y , ro , eps ) result ( heav ) real : : x , y , ro , eps , x1 ( 5 ) , x2 ( 5 ) , p ( 2 ) , p1 ( 2 ) , p2 (2) integer : : i , j , k real : : d ( 4 ) , s ( 4 ) , pi , sig , phi , dmin pi=acos ( 1.0) x1=(/ 1, 1 , 1 , 1, 1/)a42ro x2=(/ 1, 1, 1 , 1 , 1/)a42ro p=(/x , y /) do i =1,4 p1=(/x1 ( i ) , x2 ( i )/) p2=(/x1 ( i +1) , x2 ( i +1)/) call d i s t t o s e g (p , p1 , p2 , d( i ) , s ( i ) ) enddo dmin=minval (d ) ; where(d /= dmin ) s =0.0 end where s i g=sum( s ) phi= dmina42sgn ( s i g ) if ( phi . gt . 0 . 0 ) then heav =1.0 elseif ( phi . gt . eps )then heav=(1+cos ( phia42 pi / eps ) ) else heav =0.0; endif return end function 86 ! subroutine d i s t t o a r c (p , p0 , p1 , p2 , d , s ) real , intent(in ) : : p ( 2 ) , p0 ( 2 ) , p1 ( 2 ) , p2 (2) real , intent(out ) : : d , s real ,dimension ( 2 ) : : u1 , u2 , v , w1 , w2 real : : pi , d1 , d2 , s1 , s2 , q1 , q2 , q pi=acos ( 1.0) u1=p1 p0 u2=p2 p0 v=p p0 q1=a n g l e o n c i r c l e ( u1 ) ; q2=a n g l e o n c i r c l e ( u2 ) ; q=a n g l e o n c i r c l e ( v ) ; if ( q2 . l t . q1 ) q2=2a42pi+q2 if ( q . l t . 0 . and . q . l t . q1 ) q=2.0a42 pi+q if ( q . gt . q1 . and . q . l t . q2 ) then d = s q r t ( v (1)a42a422.+ v (2)a42a422.) s q r t ( u1 (1)a42a422.+ u1 ( 2 )a42a422 . ) s = sign ( 1 . 0 , d) d = abs (d) else w1 = (/ u1 ( 2 ) , u1 (1)/)/ s q r t ( u1 (1)a42a422.+ u1 ( 2 )a42a422 . ) d1 = s q r t ( ( p(1) p1 (1))a42a422.+( p(2) p1 ( 2 ) )a42a422 . ) s1 = sign ( 1 . 0 , cross2d (p p1 , w1) ) w2 = (/ u2 ( 2 ) , u2 (1)/)/ s q r t ( u2 (1)a42a422.+ u2 ( 2 )a42a422 . ) d2 = s q r t ( ( p(1) p2 (1))a42a422.+( p(2) p2 ( 2 ) )a42a422 . ) s2 = sign ( 1 . 0 , cross2d (p p2 , w2) ) if ( d1 . l t . d2 ) then d = d1 s = s1 else d=d2 s=s2 endif endif end subroutine function a n g l e o n c i r c l e ( v ) result ( theta ) real : : v ( 2 ) , u ( 2 ) , theta , q u=(/1 ,0/) q = dot product (v , u) if ( q . ne . 0 ) then q = q /( s q r t ( v (1)a42a422.+ v ( 2 )a42a422 . )a42 s q r t (u(1)a42a422.+u ( 2 )a42a422 . ) ) endif theta = sgn ( cross2d (u , v ))a42 acos ( q ) end function subroutine d i s t t o s e g (p , p0 , p1 , d , s ) real , intent(in ) : : p ( 2 ) , p0 ( 2 ) , p1 (2) real , intent(out ) : : d , s real ,dimension(2) : : u , v0 , v1 , pb , vb real : : c1 , c2 , b u=p1 p0 v0=p p0 v1=p p1 c1=dot product ( v0 , u) c2=dot product (u , u) 87 if ( c1 . l e . 0 . 0 ) then d=s q r t ( v0 (1)a42a422.+ v0 ( 2 )a42a422 . ) s=sign ( 1 . 0 , cross2d ( v0 , u ) ) elseif ( c2 . l e . c1 ) then d=s q r t ( v1 (1)a42a422.0+ v1 ( 2 )a42a422 . 0 ) s=sign ( 1 . 0 , cross2d ( v1 , u ) ) else b=c1/c2 pb=p0+ba42u vb=p pb d=s q r t ( vb (1)a42a422.0+ vb ( 2 )a42a422 . 0 ) s=sign ( 1 . 0 , cross2d (vb , u ) ) endif end subroutine function sgn ( a ) result (b) real : : a , b b=sign ( 1 . 0 , a ) if (b . ge . 0 ) b=1.0 end function function cross2d (u , v ) result ( c ) real ,dimension ( 2 ) : : u , v real : : c c=u (1)a42v(2) u (2)a42v (1) end function 88