AIAA 2013-4174

Joint Propulsion Conferences

July 14 - 17, 2013, San Jose, CA
49th AIAA/ASME/SAE/ASEE Joint Propulsion Conference

Numerical Simulation of the Zefiro 9 Performance

Using a New Dynamic SRM Ballistic Simulator

R. Bertacin1, F. Ponti2, E. Corti3 and D. Fedele4

University of Bologna, Forlì, FC, 47121, Italy

A. Annovazzi5
Avio Space Propulsion, Colleferro, Rome, 00034, Italy

This work describes the application of a new three-dimensional ballistic model, named
Downloaded by UNIV OF CALIFORNIA LOS ANGELES on May 1, 2014 | [Link] | DOI: 10.2514/6.2013-4174

ROBOOST (ROcket BOOst Simulation Tool), developed at the Laboratory Propulsion and
Mechanics of the University of Bologna (Department of Industrial Engineering), to the Vega
- Zefiro 9 Solid Rocket Motor, manufactured by the Avio company in Colleferro (Rome).
The code uses an original graphical approach and a point-by-point description of the
propellant burning surface regression. The main purpose of the newly developed model is to
investigate non homogeneous behaviors of the surface regression rate or non-isotropic
characteristics of the grain. A zero-dimensional unsteady thermo-dynamic model, coupled
with a mono-dimensional quasi-steady one, computes the internal fluid-dynamics of the
combustion chamber and contributions by igniter, nozzle erosion and thermal protections
ablation are also considered. Comparisons with reference curves and experimental data, in
terms of volume and surface regression and mean pressure time evolution, have been
performed and final results are presented and discussed.

Nomenclature
2
A = cross section area, m
A* = nozzle throat area, m2
cp = specific heat at constant pressure, J/kg K
cv = specific heat at constant volume, J/kg K
c* = propellant characteristic velocity, m/s
MW = molar weight, kg/kmol
ṁ = mass flow rate, kg/s
p = pressure, Pa
Qloss = heat power losses, W
R = ideal gas constant, J/kg K
r = regression rate, m/s
S = surface, m2
T = temperature, K
t = time, s
V = volume, m3
 = specific heats ratio, adim.
 = density, kg/m3

1
PhD Candidate, Department DIN, via Fontanelle 40, [Link]@[Link], Student Member AIAA.
2
Associate Professor, Department DIN, via Fontanelle 40, [Link]@[Link].
3
Assistant Professor, Department DIN, via Fontanelle 40, enrico.corti2@[Link].
4
PhD Candidate, Department DIN, via Fontanelle 40, [Link]@[Link], Student Member AIAA.
5
Senior Engineer, AVIO Space Propulsion Design Department, via Ariana km 5.2, [Link]@[Link].

1
American Institute of Aeronautics and Astronautics

Copyright © 2013 by the American Institute of Aeronautics and Astronautics, Inc. All rights reserved.
 Subscripts
abl = TP material ablation
b = propellant burn rate
I = node number
ign = igniter
n = nozzle
p = solid propellant or propellant combustion gas
TP = motor thermal protections
0 = combustion chamber gas

I. Introduction
Downloaded by UNIV OF CALIFORNIA LOS ANGELES on May 1, 2014 | [Link] | DOI: 10.2514/6.2013-4174

I N the design phase of a new Solid Rocket Motor (SRM) the most crucial point is represented by the prevision of
its future performance in terms of pressure1 and, hence, thrust profiles2. The main factor that leads the combustion
pressure is the time evolution of the propellant volume or, in other terms, of its burning surface and regression rate
amount. In the theoretic situation, the propellant burn rate is assumed uniformly distributed over the burning surface
and the performance calculation results quite simple. In the real case, instead, a lot of factors occur in the burning
environment and pressure and flow conditions result not uniform within the chamber3. Their main effect is on the
propellant burn rate which can locally vary from the expected theoretical amount and this alterations are normally
reported in literature with different designations, such as HUMP factor, BARF (Burning Anomaly Rate Function),
etc4,5.
The pressure drop along the motor axis, for example, causes the grain to burn faster at head-end due to the higher
local pressure, while the high speed of gas flow improves the erosive effect toward the aft-end. Moreover, acoustic
and flow field oscillations, especially in the ignition transient, could alterate the local grain performances and
introduce further instabilities in the combustion process. Hence, the capability to reconstruct instant by instant the
geometry of the grain volume during the entire combustion duration is fundamental to predict the internal ballistics
of the SRM.
A variety of techniques have been proposed to numerically represent the geometric information that defines a three-
dimensional solid and its surface. One possible solution is the analytical method6, as implemented in the Solid
Propellant rocket motor Performance computer program (SPP). This is a popular SRM simulation software with
three available approaches (two-dimensional, axisymmetric and three-dimensional), all of them based on the
subdivision of the generic section or the entire surface in pieces. Using a boolean geometry each piece is described
by combination or intersection of primitive curves (segments or arches) or solids (spheres, planes, cylinders, etc.),
whose parametric equations are evolved during surface regression.
Another way to describe the evolution of
propellant surface is the Phase-based Analytical
Method, in which a numerical layering technique
is needed and applied to each 2-D cross section7.
This approach divides the burning surface
evolution process into different consecutive
phases, within which the cross section perimeter
and port area can be described mathematically.
Each phase describes time periods during which
the set of parametric equations is the same and
ends when the set is modified. In order to define
chamber volume and burning surface area, the
section port area and perimeter are multiplied by
the grain segment length. Finally, one of the most
recent techniques is the so called Level-Set Figure 1. Layout of the developed code. The goal of the right
Method, in which the generic grain cross section part is the geometrical simulation of the grain consumption,
is described by an interface, i.e. the intersection while the left module reproduces the internal ballistics of the
8
of a level set function with a plane . The burning chamber.
perimeter represents the zero level set of the level

2
American Institute of Aeronautics and Astronautics
 set function, and its evolution describes the time-based motion of the interface. Moreover, the resulting initial value
partial differential equation (PDE) for the evolution of the level set function is similar to an Hamilton-Jacobi
equation and it is solved using an entropy-satisfying schemes borrowed from the solutions of hyperbolic
conservation laws. The main advantage related to this method is an easy evaluation of surface curvatures and
normal, and a natural evolution of topology.
Recently, a new approach has been proposed 9,10 and it is based on a so called Minimum Distance Function (MDF)
which, starting from a typical Computer-Aided-Design (CAD) geometry file, is able to reconstruct the surface
evolution by calculating the zero-level contour of the time-dependent MDF on each motor cross-section, in a similar
way to the Level-Set-Method.
All these described methods allow a good tracking of the surface regression during simulations and could also
implement non-uniform burn rate distributions, which however must be smoothed and quite extensive, i.e. not
spatially confined, in order to be efficiently implemented in the layering technique. Moreover, the section-by-section
description of the propellant volume used by the analytical methods and by the level-set theory, make these burn rate
Downloaded by UNIV OF CALIFORNIA LOS ANGELES on May 1, 2014 | [Link] | DOI: 10.2514/6.2013-4174

distributions to have some form of prevalent orientation, for example in the axial direction rather than in the
azimuthal one.
The work described in this paper try to get over these limitations and proposes the last development phase of a new
approach11 to reproduce the combustion process and the internal ballistics of large-scale SRMs. The objective is to
create a simulation tool, adaptable to different motor configurations and able to assure anyway an accurate three-
dimensional representation of the burning surface regression and a reliable one-dimensional estimation of the
chamber internal fluid dynamics.
Moreover, a key requisite is the ability to implement heterogeneous distributions of the grain burn back all over the
surface and eventually to bring macroscopic discontinuities (i.e. cavities or air entrapments) into the propellant
volume. The method used to generate the initial geometry and simulating its evolution during the combustion burn
back is the discrete representation of the grain surface and its controlled motion point-by-point.
Compared to the level-set approach the proposed method could introduce some inaccuracies in terms of curvatures
and normal evolutions but, on the other side, allows a more “physical” description of the surface regression process
and versatility in the insertion of the burn rate heterogeneities,
which in this case could also be spatially constricted.

II. Code Overview

The developed code has been programmed entirely in the
Matlab environment and represents the final evolution of some
previous versions based on an two-dimensional approach,
applied respectively on an axial- and on a cross-section of the
motor.
The software ROBOOST has a modular structure, as shown in
Fig. 1, and additional modules will be implemented in future
developments. Two main subsystems could be identified: the
grain burn back part and the internal ballistics one. Most of the
code algorithm is organized to work automatically on multi-
core machines to better exploit the computer performances and
reduce the computational time of simulations.

A. Definition of the Grain Geometry

In the first subsystem the SRM configuration is acquired and
converted in the geometric domain, used afterwards for the
surface regression simulation. Because a wide range of Figure 2. Example of adopted triangle mesh.
Computer-Aided-Design softwares are commercially available The original CAD project (upper) is converted
and they offer three-dimensional drafting capabilities, the motor in a triangle surface mesh. Yellow elements are
configuration is generated outside the code and then uploaded virtually exposed to the combustion process
in the simulator using the common stereo-lithography file and are moved during simulation, the gray ones
format (STL). are inhibited and represent the casing profile.

3
American Institute of Aeronautics and Astronautics
 Borrowing a typical technique from Compuetr Graphics, the STL motor geometry is converted in a more thick
surface discrete polygonal. In particular, a triangle mesh is preserved in order to obtain an unambiguous definition of
the normal vector for each surface element (also defined as face). For simple configurations the mesh generation
procedure could be implemented directly within the Matlab software or as an alternative, for complex grain shapes it
could be used an external mesher, e.g. ANSYS, Comsol, etc., uploading later the vertices and triangulations arrays
into the code.
Normally, an unstructured mesh is used and the number of triangles vary in relation to the grain complexity,
curvatures and, obviously, the desired grade of accuracy. In fact, an high accuracy requires more elements for better
approximate the surface, but on the other side will increase the simulation computational time.
The next step, is the distinction of inhibited and exposed triangles. The first ones are virtually in contact with the
internal thermal protections of the motor casing and therefore are maintained fixed during simulations, while
through the controlled motion of the second ones the regression of the burning surface is reconstructed. The
presence of the inhibited surface is not required for the surface burnback process, but is necessary for the calculation
Downloaded by UNIV OF CALIFORNIA LOS ANGELES on May 1, 2014 | [Link] | DOI: 10.2514/6.2013-4174

of the grain volume and its consumption, through a signed sum of volumes of single tetrahedrons. Moreover, in the
initialization phase it is necessary to define the casing two-dimensional profile which will be used later for the
imposition of the wall condition to the mesh vertices.

B. Surface Evolution
In order to respect the so called Huygens’ Principle, applied to the
thermal wave propagation inside the propellant due to the
combustion process, the motion of the burning surface takes place
by shifting the mesh vertices along their respective normal
directions. Vertices that lie on the boundary between the exposed
and the inhibited surface are forced to move along the casing
profile in order to prevent anomalous detachments of the inhibited
elements.
Compared to the level-set method or other section-based
techniques, on a three-dimensional polygonal mesh the calculation
of the normal-to-vertex is quite critical and requires a combination
Figure 3. Example of mesh motion. The of the normals of the triangle elements connected to a single node.
initial surface (green) evolve in a new one Different solutions12 are reported in literature, all dependent on
(blue) depending on single vertices the particular field of application and purpose of the study. Most
displacement, carried out along their of them use either the internal angle of connected triangles or their
respective normal vectors (red arrows). area as weight function,
and have been all tested
and evaluated. At present, the code implements a particular linear
combination of adjacent normals-to-face in which the angle between two
vectors is used as weight variable in the calculation. This formulation
prevents redundancies when two or more adjacent triangles are coplanar,
but at the same time gives the same weight value to each normal direction
(in a similar way to an arithmetic mean).
Once defined the normal direction, the vertex displacement is calculated
according to the instantaneous local amount of the propellant burn rate and
the simulation time step, as shown in Fig. 3. Obviously the surface
regression rate is assumed constant for the duration of the advancement
time step. The adopted meshing approach allows to define point-by-point
the surface motion, applying alterations to the burn rate mean value both in
axial and in azimuthal direction, based on the local estimated instantaneous Figure 4. Remeshing procedures
pressure or rheological characteristics of the propellant. outline. The goal of this group of
During simulations, in order to preserve the coherence of the mesh, e.g. the subroutines is of cheching and
observance of the Euler’s Formula, a set of remeshing procedures has been preserving the mesh quality and
implemented at each time iteration to manage the mesh motion and quality. resolution.
As represented in Fig. 4, different checks are scheduled and the proposed

4
American Institute of Aeronautics and Astronautics
 sequence of operations represents the solution which assures a better evolution of the geometric domain. Because
inhibited elements attached to the burning portion of surface, i.e. with at least one vertex in common with an
exposed triangle, are indirectly deformed and compressed by the displacements of the moving nodes, overlapping
phenomenon can happen. Even if not harmful at the beginning for the grain volume calculation, in the last part of
the simulation, increasing the number of overlapped triangles, this problem can produce considerable errors in
addition to an useless increase of mesh array dimensions. For this reason the first subroutines solves this problem
removing overlapped triangulations. Then, the second, third and fourth check analyze the domain resolution, i.e. the
length of single edges or triangle areas, and preserve it in the user-defined range by collapsing or splitting
anomalous elements.
The next subroutine, instead, could be enabled or not and execute a sort of smoothing over the burning surface 13. In
particular this procedure, identifies the mesh features, analyze the sharp corners evolution and corrects possible
surface spikes. However, unlike other approaches, such as level-set methods, when a sharp corner is internal-burning
it preserves its feature and, at present, it is not automatically smoothed to a semicircular arc. This introduce a certain
Downloaded by UNIV OF CALIFORNIA LOS ANGELES on May 1, 2014 | [Link] | DOI: 10.2514/6.2013-4174

error in the calculation of instantaneous grain volume if compared with the real physical regression process, but this
inaccuracy becomes zero if the evaluation is performed with others analytical solvers, such as SPP, that use a
primitives intersection technique. Anyway, future improvements will be focused on the solution of this source of
geometrical inaccuracy.
The next check detects possible residual unconnected vertices and verify the two-manifold coherence of the entire
surface mesh. Finally, the last function imposes the wall condition to all moving vertices and specifies as new
inhibited elements all exposed triangulations that have reached the casing profile with all their connected vertices.
An additive stability check is put in parallel to the surface regression subroutines in order to resize, when needed, the
simulation time step, depending on the resolution of the mesh, in order to prevent discontinuities of the moving
surface, such as self-intersections between triangle elements.

C. Internal Ballistics Model

The internal ballistics module represents the most recent part developed within the software and, for this reason, it
will be the primary argument of future improvements. It implements a solution proposed by Avio and normally used
in its experimental data analysis. In detail, the ballistic algorithm consists of the combination of a zero-dimensional
unsteady model and a one-dimensional quasi-steady one.

Figure 5. Graphical representation of the 0D unsteady ballistic model.

Mass contributions due to burning grain, igniter and thermal protections
ablation are considered. The erosion of nozzle throat section is also
implemented.

The zero-dimensional algorithm computes the time evolution of thermo-dynamic parameters inside the combustion
chamber using an averaged approach and is able to simulate also the ignition and tail-off transients. Moreover, it is
based on the following assumptions:
 The gas in the flow cannel respects to the perfect gas law;

5
American Institute of Aeronautics and Astronautics
  The internal surfaces of the combustion chamber, except the propellant one, are not adiabatic;
 The chamber gas thermodynamic parameters (T0, p0, 0) are only time dependent;
 The physics (, MW) of the gas within the chamber and those generated by igniter and propellant can vary;
 The heat production due to the combustion process is assumed occurring close to the propellant surface and
outside the control volume;
 The flow within the motor chamber is quite subsonic;
 The kinetic energy of the chamber flow and the velocities of the gas leaving the burning surface and at the
nozzle inlet can be considered negligible compared to the thermal energy released by the propellant;
 The fluid is considered inviscid.
Downloaded by UNIV OF CALIFORNIA LOS ANGELES on May 1, 2014 | [Link] | DOI: 10.2514/6.2013-4174

Through the combination and rearrangement of the continuity and energy balance equations, the time dependencies
of chamber pressure and temperature can be obtained as follows:

. . . . cV dV0
m ign c Pign Tign  m p c Pp T p  m TP c PTP Tabl  m n c Pn Tn  p0  Qloss
dp 0 R dt
 (1)
dt cV
V0
R
dT0 1  dp0 dV p . . . . 
 V0  p0 0  0  m ign  m p  mTP  m n  (2)
dt  0V0 R  dt dt  0  

In particular, as introduced previously, the instantaneous time-derivative of the chamber volume is assumed equal to
the variation of the propellant volume in a single simulation step, with the hypothesis that no deformations are
applied at the casing structure and its initial volume is maintained constant. This value is obtained directly from the
regression module without the need of computing the instantaneous burning surface and mean burn rate amount, and
represents also the mass of gas introduced in the chamber from the propellant surface.
At present, the integration of these ordinary differential equations (ODE) is performed using the Euler’s Algorithm,
and the stability of the solution is assured by dividing the main time step in multiple sub-integrations, but future
improvements will replace it with a faster fourth-order Runge Kutta method.
Contributions due to the igniter, burning propellant, thermal protections and nozzle are defined separately within
specific subroutines. In detail, the igniter performance is implemented through the interpolation of available
experimental data, in terms of mass flow rate and gas properties. The physics of gases produced by the grain
combustion is computed during simulations from parametric thermo-chemical maps generated through the NASA
C.E.A. code, while the exhaust mass flow rate through the nozzle
is obtained from the equation:

p0 An
m n  (3)
c

Finally, the contribution due to the ablation of Thermal Protections

(TP), also called liner, is obtained from

 TP  TP  rabl  STP

m (4) Figure 6. Portion of the axial flow as
described in the 1D quasi-steady
where the density and ablation rate law of the particular thermal ballistic model. The combustion chamber
protection material is extracted from experimental data or is modeled as a sequence of sections with
calculated by Avio through the CMA commercial software14, while different port areas.

6
American Institute of Aeronautics and Astronautics
 the evolution of the exposed PT surface is once again computed directly from the three-dimensional regression
module.
When a SRM has a low length-to-diameter ratio the only zero-dimensional approach is sufficient to describe its
internal ballistics. Otherwise, with high length-to-diameter ratios the pressure drop along the motor length becomes
not negligible and must be considered for a correct calculation of the local propellant burn rate, together with other
contributions such as the erosive combustion.
Therefore, the one-dimensional quasi-steady model
has been added in order to estimate the longitudinal
distribution of the thermodynamic parameters
whose mean values are defined in the zero-
dimensional calculation.
Applying the common finite-difference approach,
the internal motor bore, defined by the triangle
Downloaded by UNIV OF CALIFORNIA LOS ANGELES on May 1, 2014 | [Link] | DOI: 10.2514/6.2013-4174

mesh, is divided into consecutive segments and the

calculation of Mach number and other physical
variables is performed within each one. Moreover,
because the radial component of the flow field is
neglected, each section between two consecutive
segments is converted into an equivalent circular
one and the motor chamber can be compared to a
cylindrical section-variable pipe.
Using a set of equations defined in Shapiro15 for an
Figure 7. Subdivision of the computational domain in
isentropic flow, the variation of four dependent
the 1D ballistic model. In each node thermodynamic
variables (Mach, pressure, density and temperature)
parameters are computed as linear combination of port
between two consecutive chamber sections is
area variation, mass adduction and thermal exchanges.
computed as a linear combination of effects due to
three independent known variables, i.e. the port area
variation, mass adduction and thermal exchanges. Obviously, this assumption represents a computational
simplification because in the real case the three contributions happen at the same time.
The solution of the channel flow field is performed through iterations, by imposing a first attempt pressure value on
the first node and solving the others in cascade. Convergence is achieved when the weighted average pressure along
the motor bore, i.e. in each channel segment i,

V  p i i
p 1D  i (5)
V i
i

equals the solution of the zero-dimensional model. The so calculated pressure profile is then implemented in the
regression module to determine the burn rate distribution for the
next time integration, using the classical Vieille’s Law, while the
flow field velocity is used to define the erosive burning
contribution through the Lawrence and Beddini16 equation (an
evolution of the Lenoir-Robillard model).

III. The Zefiro 9 SRM

Acquiring the simulation domain from an external CAD software,
the developed code proves to be very flexible and is able to
implement a great range of SRM geometries, in terms of motor
dimensions and internal configurations. After a first validation Figure 8. The Zefiro 9 SRM. This
phase performed on the small-scale motor configuration, also propulsion system constitutes the third-stage
called Baria, used by Avio for the ballistic characterization of the equips the new italian launcher Vega,
manufactured by Avio.

7
American Institute of Aeronautics and Astronautics
 solid propellant, the code has been applied to the real Zefiro 9 SRM, shown n Fig. 8.
This motor, entirely manufactured by Avio, represents the third stage of the new european launcher Vega, whose
maiden flight has been performed successfully in February 2012 from the Guyana Space Center (French Guyana). It
is about 3 meters long, with a maximum diameter of about 1.9 meters and a total propellant mass of about 10 tons.
Its main characteristics are an external carbon-epoxy filament-wound casing, a consumable igniter, an
electromechanical thrust vector control system and a very high density Al-HTPB-AP propellant mixture.
Its internal bore has a quite original configuration with a very variable cross section, circular in the fore and central
part and finocyl-shaped in the rear part, near the nozzle inlet. The motor solid geometry has been generated using the
Catia CAD software and its conversion to triangle surface mesh has been completed within the mesh generator of
the ANSYS software package.
In detail, a simulation domain of about 52.000 triangle elements is normally used and their size is not uniform but
vary depending on the local surface curvature and on the type of element, i.e. it is exposed or inhibited. A maximum
resolution of 3÷5 mm over the burning surface is so allowed while a minimum of about 50÷60 mm is applied over
Downloaded by UNIV OF CALIFORNIA LOS ANGELES on May 1, 2014 | [Link] | DOI: 10.2514/6.2013-4174

the inhibited one.

Figure 9. Conversion of the motor CAD geometry to the triangle mesh. The original motor
configuration has been created using the Catia code, while the surface mesh has been generated
within the ANSYS software package.

The Z9 burnout time is of about 120 s but, at present, the duration of a single simulation is in the order of some
hours on a Intel Core i7 CPU machine. This time results from the used spatial resolution and a simulation time step
in the order of some cents of second must be used in order to assure the stability of the calculation , i.e. observing a
sort of Courant-Friedrichs-Lewy (CFL) condition, and the accuracy of results.
Because experimental data about the burning surface time evolution normally are not available, in this case the
comparison of simulated results is achieved assuming as reference curves the outputs of a particular version of the
SPP code, normally used in the Avio Solid Propulsion Division. In detail, this software couples the performance of
the SPP algorithm with another code called GEOMGDE, developed by Avio, and represents the main tool adopted
by the company for the performance estimation in the design of new propulsion systems.
The goal of the described work is to have an absolute difference between the new simulator results and the reference
SPP-GEOMGDE curves below 3%, a common threshold used for the validation of new codes. Fig. 11 and Fig. 12
show the results of the comparison in terms of burning surface and grain volume derivative.
Both the curves are referred to the burned propellant web, in order to filter the ballistic model contribution and
analyze only the performance of the grain regression module. The error profiles are calculated as the instantaneous
difference between the two curves and highlight how the developed three-dimensional software is able to maintain
globally a gap with the reference data in the range of ±2%. In particular, the attention must be mainly focused on the
second comparison because the instantaneous volume derivative is the parameter that directly affect the next internal
ballistic calculation.

8
American Institute of Aeronautics and Astronautics
Downloaded by UNIV OF CALIFORNIA LOS ANGELES on May 1, 2014 | [Link] | DOI: 10.2514/6.2013-4174

Figure 10. The Zefiro 9 SRM grain regression. Some different simulation frames are
proposed. The yellow portion represents the evolving burning surface, while the gray one is the
initial casing profile.

9
American Institute of Aeronautics and Astronautics
Downloaded by UNIV OF CALIFORNIA LOS ANGELES on May 1, 2014 | [Link] | DOI: 10.2514/6.2013-4174

Figure 11. Time evolution of grain burning surface. The simulated curve is compared with a
reference one obtained by Avio using its SPP/GEOMGDE simulation code.

Figure 12. Time evolution of grain volume derivative. The simulated curve is compared
with a reference one obtained by Avio using its SPP/GEOMGDE simulation code.

The isolated spikes in the first part of the simulation are due to the remeshing procedures and correspond to the time
instant when the finocyl parts reach the casing wall. In this situation, many triangle elements, which were previously
exposed, become simultaneously inhibited and overlap, and the execution of the checking subroutines introduces
some oscillations in the volume computation. Future improvements to these procedures will be focused on reducing
their alteration effect in this simulation phase, in order to obtain a more regular profile.
A further comparison is now proposed and the evolution of chamber pressure is analyzed. In this case, the
experimental profile of the QM3 firing test is assumed as reference curve, while the output of the zero-dimensional
ballistic model is analyzed for the evaluation.
Contributions due to igniter, nozzle erosion and ablation of thermal protections are considered in the simulation. In
detail, the igniter performance is implemented by interpolation of available experimental data and the same method
is used for the nozzle throat consumption. Because the particular configuration of the Zefiro 9 motor, the thermal
protection coating of the carbon casing results exposed to the chamber hot gas since the beginning of the combustion
and its contribution is not negligible. As mentioned above, the exposed TP surface is calculated directly from the
three-dimensional simulation domain and, using specific semi-empirical formulations developed by Avio, it defines
both the ablation mass flow rate and the thermal flux absorbed. Fig. 13 shows the evolution of TP surface calculated
during the simulation.
An empirical web-dependent HUMP function, obtained from experimental firing tests, has been introduced in the
propellant burn rate calculation, in order to achieve a more realistic pressure profile. As shown in Fig. 14 the
simulated curve trend is appreciable and very close to the experimental one, despite the initial oscillation
corresponding to the peaks observed in the volume derivative profile. A mean percentage error in the range of ±2%
can be observed.

10
American Institute of Aeronautics and Astronautics
Downloaded by UNIV OF CALIFORNIA LOS ANGELES on May 1, 2014 | [Link] | DOI: 10.2514/6.2013-4174

Figure 13. Time evolution of Z9 Thermal Protections (TP). The plotted curve represents the
surface of the casing TP coating exposed to the chamber hot gases.

Figure 14. Combustion chamber mean pressure. The simulated pressure profile, computed
using the 0D unsteady ballistic model, is compared with the experimental one. An empirical
HUMP function has been introduced in the propellant burn rate calculation.

About the one-dimensional ballistic model, because the Zefiro 9 SRM has a low length-to-diameter ratio, no evident
pressure drop along the motor axis is expected and the only zero-dimensional analysis could be sufficient.
Fig. 15 shows the pressure distribution within the motor channel as obtained from the simulation. In the quasi-steady
phase of the combustion profile the gap between the head-end and bottom-end sections remains lower than 1 bar and
it could be neglected. Each marker identifies the middle-section of a single chamber segment in which physical
quantities are solved through the ballistic module. In Fig. 15 are also reported the mass flow rate contributions
referred to each domain cell, due to the propellant combustion process and the TP coating ablation.
An evenly spaced segmentation is currently implemented within the code, but a specific analysis has begun to
optimize the spatial ballistic domain and to achieve a more smoothed trend of profiles between the cylindrical and
the slotted part of the channel.

11
American Institute of Aeronautics and Astronautics
 Finally, Fig. 16 shows the used segmentation of the internal motor bore. As described above, for each cross-section
an equivalent circular profile has been defined, in order to compare the combustion chamber to a simplified circular
section-variable channel.
In particular, for each obtained cell the amount of exposed TP coating is calculated as the difference between the
casing internal surface of the segment and the sum of area portions referred to included inhibited triangulations.
Obviously, when an equivalent circular section reaches the motor case profile, i.e. no exposed triangulations are
intersected, it is hold on the profile for the rest of the simulation.
Downloaded by UNIV OF CALIFORNIA LOS ANGELES on May 1, 2014 | [Link] | DOI: 10.2514/6.2013-4174

Figure 15. Simulated pressure, gas velocity and mass flow rate distributions along the Zefiro
9 length in quasi-steady conditions. As expected, the pressure drop between the head-end and
the aft-end is lower than 1 bar and could be neglected.

Figure 16. Axial segmentation of internal bore. The blue lines represent the circular cross-sections
used for the chamber segmentation. An evenly spaced distribution is currently implemented.

12
American Institute of Aeronautics and Astronautics
 IV. Conclusion
A new three-dimensional ballistic code, called ROBOOST (Rocket BOOst Simulation Tool), has been developed.
Using a triangle moving mesh to define the grain volume, the code is able to simulate the regression of propellant
burning surface of a wide range of SRM configurations. Surface displacement
is applied to each mesh vertex and, for this reason, any heterogeneous distribution of propellant burn rate can be
easily applied. The SRM internal ballistics is simulated using the coupling of a zero-dimensional unsteady model
and a one-dimensional quasi-steady one. This solution allows to describe the combustion gas properties along the
chamber length not only in the steady-state phase, but also in the transient phases. Contributions due to igniter,
nozzle erosion and thermal protections ablation are also implemented.
In this paper the developed code has been applied to the Zefiro 9 SRM, manufactured by Avio, and simulated results
Downloaded by UNIV OF CALIFORNIA LOS ANGELES on May 1, 2014 | [Link] | DOI: 10.2514/6.2013-4174

have been compared with experimental and reference curves. Globally, an error in the range of ±2% has been
obtained in terms of surface and volume time evolutions, and a very appreciable trend has been observed for the
pressure profile.
Future developments will be focused on the reduction of these errors especially by improving the simulation
performances of the regression module and, in particular, of the remeshing procedures. Once reached the desired
accuracy, the next step will be represented by the analysis of the experimental HUMP function through the
implementation in the code of specific burn rate three-dimensional distributions, computed both from experimental
tests and from fluid-dynamic simulations of propellant casting processes.
Moreover, the code will be applied to other SRM configurations, also multi-segment, always manufactured by Avio.

Acknowledgments
Authors would like to thanks the Avio s.p.a. Space Division, and Ing. Adriano Annovazzi as first, for
supporting us in this research project and supplying us the necessary data.

References
1
Kallmeyer, T.E., Sayer, L.H., “Differences between Actual and Predicted Pressure-time Histories of Solid Rocket Motors”,
AIAA 82-1094, 18th AIAA/ASME JPC, 1982, Cleveland, Ohio.
2
Thepenier, J., Ribéreau, D., and Giraud, E., “Application of Advanced Computational Softwares in Propellant Grain
Analysis: A Major Contribution of Future SRM Development for Space Application”, International Astronautical Federation,
IAF 98-S.2.06, October 1997, Turin, Italy.
3
Hessler, R.O., Glick, R.L., “Consistent Definitions for Burning-Rate Measurements in Solid Rocket Motors”, Combustion,
Explosion and Shock Waves, vol. 36, no. 1, 2000.
4
Friedlander, M.P., Jordan, F.W., “Radial Variation of Burning Rate in Center Perforated Grains”, AIAA 84-1442, 20th
AIAA/SAE/ASME JPC, 1984, Cincinnati, Ohio.
5
Maggi, F., De Luca L.T.,m Bandera, A., Subith, V.S., Annovazzi, A., “Burn-Rate Measurement on Small-Scale Rocket
Motors”, Defence Science Journal, vol. 56, no. 3, July 2006, pp. 353-367.
6
Cai, Q., Bao, F., Liu, Y., Hu, H., and Wei, H., “Solid Rocket Grain Configuration Variational Design Method Based on
Geometric Constraint Sketch”, 2nd WRI World Congress on Software Engineering, December 19-20, 2010, Wuhan, China.
7
NASA SP-8076, “Solid Propellant Grain Design and Internal Ballistics”, March 1972.
8
Cavallini, E., Bianchi, D., Favini, B., Di Giacinto, M., “Propellant Effects on SRM Upper Stage Internal Ballistics and
Performance with Nozzle Erosion Characterization”, AIAA 2012-3887, 48th AIAA/ASME/SAE/ASEE JPC, Atlanta, Georgia.
9
Willcox, M. A., Brewster M. Q., Tang K. C., Stewart, D. S., and Kuznetsov, I., “Solid Rocket Motor Internal Ballistics
Simulation Using Three-Dimensional Grain Burnback”, Journal of Propulsion and Power, Vol. 23, No. 3, May-June 2007, pp.
575-584.
10
Willcox, M.A., Brewster, M.Q., Tang, K.C., Stewart, D.S., “Solid Propellant Grain Design and Burnback Simulation Using
a Minimum Distance Function”, Journal of Propulsion and Power, vol. 23, no. 2, 2007, pp. 465-475.
11
Bertacin, R., Ponti, F., Annovazzi, A., “A New Three- Dimensional Ballistic Model for Solid Rocket Motor Non-
Homogeneous Combustion”, AIAA 2012-3974, 48th AIAA/ASME/SAE/ASEE JPC, 2012, Atlanta, Georgia.

13
American Institute of Aeronautics and Astronautics
 12
Guoy, D., Wilmarth, T., Alexander, P., Jiao, X., Campbell, et al., “Parallel Mesh Adaptation for Highly Evolving
Geometries with Application to Solid Propellant Rockets”, 16th International Meshing Roundtable, Session 5B, 2008, pp. 515-
534
13
Morigi, S., “Geometric Surface Evolution with Tangential Contribution”, Journal of Computational and Applied
Mathematics, Vol. 233, Issue 5, January 2010, pp. 1277-1287.
14
Bianchi, D., “Modeling of Ablation Phenomena in Space Applications”, Ph.D. thesis, 2006/07, University of Rome “La
Sapienza”.
15
Shapiro, A.H., “The Dynamics and Thermodynamics of Compressible Fluid Flow”, 1957, The Ronald Press Company,
New York.
16
Lawrence, W., Matthews, D., Deverall, L., “The Experimental and Theoretical Comparison of the Erosive Burning
Characteristics of composite Propellants”, AIAA 68-531, 1968.
Downloaded by UNIV OF CALIFORNIA LOS ANGELES on May 1, 2014 | [Link] | DOI: 10.2514/6.2013-4174

14
American Institute of Aeronautics and Astronautics