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 levelset 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
twodimensions 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 FaceO 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 FirstOrder 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 pressuretime pro le for a solid rocket motor . . . . . . . . . . . . . . . 6
2.1 Star Grain De nition [27] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2 FaceO 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 7pointed 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 vepointed star using mesh of 10 10 with 9th order
polynomial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
5.7 Normalized residuals for vepointed 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 1D 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 coneshaped 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 1D ow throughout both the combustion cham
ber and the nozzle. For 2 and 3D 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 RatioMach 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
indepth 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 pressuretime 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 TaylorCulick
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 quasiviscous. 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 TaylorCulick
pro le was examined to determine the timedependent 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 TaylorCulick 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 TaylorCulick 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 TaylorCulick
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 prelaunch 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 crosssection 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 thrusttime 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 6DOF) 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 crosssection. 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 3D 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. 3D nitevolume gasparticle solver for internal aerodynamics [43]
2.2 FaceO setting Method
Another popular method for surface propagation is called the FaceO 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: FaceO 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 UrbanaChampaign 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 faceo 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 markerpoint
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 2D, if a uniform velocity eld is
assumed, the level set will simply convect in the direction of the zaxis with no shape changes.
However, if a nonuniform 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
11 0 0 11101
22 1 0 02101
22 1 1 02101
01 0 0 11001
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 2D
or 3D 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 HamiltonJacobi 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 2D rstorder 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 7pointed 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 7pointed 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 7pointed 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 7pointed 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 HamiltonJacobi 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 1D polynomials and satisfy the di erential equation in a pointwise fashion.
This makes nite di erencing schemes simple to implement and still allow for highorder
implementation. However, implementing highorder 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 nonlocal manner, satisfying the equation across the entire domain which make
them suitable for highorder. However, they are implicit in time. The desired formulation
would be to have local highorder 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 HamiltonJacobi 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 xderivative of using data from either the left or the right (1 and
2) and q represents the yderivative 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 xdirection 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)
i1 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 ydirection 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 zdirection 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 2D level set is essentially being convected
in the zdirection, 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 CourantFriedrichsLewy (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 pseudotime
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 psuedotime 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 vepointed 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 vepointed 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 zdirection,
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 2D, 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 2D 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 vepointed 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 rootmeansquare 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 tailo 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 tailo 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 crosssection. 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) Nonzero 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 (nonuniform) wavespeed
should be implemented to make inclusion of second order propellant burning e ects such
as multiphase ow and erosive burning possible. It would also be desirable to test the
threedimensional 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., \3D Grain Design and Ballistic Analysis Using the SPP97
Code," AIAA paper 19973340, 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 SP8051, 1971.
[5] Foster, W. and Jenkins, R., \Analysis of Advanced Solid Rocket Motor Ignition Phe
nomena," NASA CR199427, 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., \TwoDimensional Axisymmetric Analysis of SRM
Ignition Transient," AIAA paper 19932311, 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 TimeDependent 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 TaylorCulick Pro le,"
AIAA paper 20075797, 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
20107079, 2010.
[14] Maicke, B. A. and Majdalani, J., \On the Compressible HartMcClure Mean Flow
Motion in Simulated Rocket Motors," AIAA paper 20107077, 2010.
55
[15] Akiki, M. and Majdalani, J., \QuasiAnalytical Approximation of the Compressible
Flow in a Planar Rocket Con guration," AIAA paper 20107080, 2010.
[16] Chedevergne, F., Casalis, G., and Majdalani, J., \DNS Investigation of the Taylor
Culick Flow Stability," AIAA paper 20075796, 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 19960761, 1996.
[20] Coats, D. and Dunn, S., \Improved Motor Stability Predictions for 3D Grains Using
the SPP Code," AIAA paper 19973251, 1997.
[21] French, J., \Tangential Mode Instability of SRMs with Even and Odd Numbers of
Slots," AIAA paper 20023612, 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 20024341, 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 1999408, 1999.
[27] Jenkins, R. and Hart eld, R., \Hybrid Particle SwarmPattern Search Optimizer for
Aerospace Applications," AIAA paper 20107078, 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
20115798, 2011.
[29] Woltosz, W., The Application of Numerical Optimization Techniques to SolidPropellant
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 ThreeDimensional 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 IgnitiontoBurn Out Analysis of
SRM Internal Ballistic and Performances," AIAA paper 20085141, 2008.
[34] Cavallini, E., Favini, B., Giacinto, M. D., and Serraglia, F., \SRM Internal Ballistic
Numerical Simulation by SPINBALL Model," AIAA paper 20095512, 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 20024298, 2002.
[37] Sforzini, R., \An Automated Approach to Design of Solid Rockets Utilizing a Special
Internal Ballistics Model," AIAA paper 801135, 1980.
[38] Sforzini, R., \Design and Performance Analysis of SolidPropellant Rocket Motors Using
a Simpli ed Computer Program," NASA CR129025.
[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 19983702, 1998.
[41] French, J. and Coats, D., \Automated 3D Solid Rocket Combustion Stability Analysis,"
AIAA paper 19992797, 1999.
[42] Coats, D., Dunn, S., and Berker, D., \Improvements to the Solid Performance Program,"
AIAA paper 20034504, 2003.
[43] French, J. and Tullos, J., \Solid Rocket Motor Grid Generation and CFD for Internal
Ballistcs Analysis," AIAA paper 20054162, 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 UrbanaChampaign, 2008.
57
[46] Dick, W., Fiedler, R., and Heath, M., \Building Rocstar: Simulation Science for Solid
Propellant Rocket Motors," AIAA paper 20064590, 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
20053991.
[48] Fiedler, R., Wasistho, B., and Brandyberry, M., \Full 3D Simulation of Turbulent Flow
in the RSRM," AIAA paper 20064587, 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 PQM1," AIAA
paper 20064592, 2006.
[50] Osher, S. and Sethian, J., \Fronts Propagating with Curvature Dependent Speed: Al
gorithms Based on HamiltonJacobi 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
HiddenSurface 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 RungeKutte discontinuous Galerkin method for
convectiondominated 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 HamiltonJacobi 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
HamiltonJacobi 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
TimeExplicit Discontinuous Galerkin Method," Mathematical Modeling of Natural Phe
nomenon, Vol. 10, No. 10, 2010, pp. 1{42.
[64] Persson, P. and Peraire, J., \SubCell Shock Capturing for Discontinuous Galerkin
Methods," AIAA paper 2006112, 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 (1D) nonoverlapping 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
k1
D
k
D
k+1
x
r
k1
=x
l
k
x
r
k
=x
l
k1
Figure B.1: Stencil for 1D example
From Equation B.4, the local state is de ned using a higherorder 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.
Reexamining 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 illconditioned 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 higherorder polynomials. Using instead something
like an orthonormal basis (threeterm recursion relation) will lead to a more wellconditioned
matrix.
The threeterm 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 nonheirarchical, 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 GaussLobatto 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, reexamining
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 FirstOrder Finite Element to Discontinuous Galerkin:
A Classic Fluid Dynamics Example
The twodimensional 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 reenters 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 zaxis, the
vortex has dissipated (nonphysical) 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