International Journal of Aerospace Engineering

Volume 2016, Article ID 3765796, 8 pages

http://dx.doi.org/10.1155/2016/3765796

## Monte Carlo Uncertainty Quantification Using Quasi-1D SRM Ballistic Model

^{1}Department of Aerospace Science and Technology, SPLab, Politecnico di Milano, 20156 Milan, Italy^{2}Space Propulsion Design Department, AVIO S.p.A., 00034 Colleferro, Italy

Received 5 January 2016; Accepted 23 March 2016

Academic Editor: William W. Liou

Copyright © 2016 Davide Viganò et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

Compactness, reliability, readiness, and construction simplicity of solid rocket motors make them very appealing for commercial launcher missions and embarked systems. Solid propulsion grants high thrust-to-weight ratio, high volumetric specific impulse, and a Technology Readiness Level of 9. However, solid rocket systems are missing any throttling capability at run-time, since pressure-time evolution is defined at the design phase. This lack of mission flexibility makes their missions sensitive to deviations of performance from nominal behavior. For this reason, the reliability of predictions and reproducibility of performances represent a primary goal in this field. This paper presents an analysis of SRM performance uncertainties throughout the implementation of a quasi-1D numerical model of motor internal ballistics based on Shapiro’s equations. The code is coupled with a Monte Carlo algorithm to evaluate statistics and propagation of some peculiar uncertainties from design data to rocker performance parameters. The model has been set for the reproduction of a small-scale rocket motor, discussing a set of parametric investigations on uncertainty propagation across the ballistic model.

#### 1. Introduction

In the current time frame, private space companies, supplying manned or unmanned flight services, are replacing governmental involvement in such missions. Different options are explored for what can be considered the rise of space access commercialization, spanning from the optimization of current launch options (namely, solid and liquid propulsion) to the development of advanced concepts such as hybrid propulsion, air-launched alternatives, or reusable stages. In any case, the development is targeting the reduction of costs. In this context, solid propulsion has an active role for boosting phases, initial stages, or embarked systems. In a paper describing the roadmap for solid propulsion published in 2010 the authors underlined appealing features such as the cost-effectiveness, the reliability of such technology, the capability of high thrust-to-weight ratio and the high propellant density; they also addressed well known limits, namely, reduced specific impulse and scarce flexibility **[1].

In general a solid rocket motor is not throttleable, limiting the capability of in-flight corrections. Variable thrust solutions given by the actuation of a pintle-nozzle are still matter of experimentation and are not implemented in large-scale systems. The thrust profile is decided during the design process and is strictly related to the pressure-time history of the combustion chamber. Once ignited, solid rocket motors proceed till the exhaustion of all the propellant stored in the combustion chamber. Shutdown can be achieved using destructive techniques or injection of flame suppressors. Reignition is not possible anyway **[2].

Deviations from the nominal behavior directly influence the mission profile and may be caused by multiple reasons. The commercially used composite solid propellant is one of the primary elements to be considered. This material is a complex mixture of oxidizer salts, metallic powders, and a binder which are compounded, cast, and cured to form the grain. The initial shape of the propellant charge evolves during combustion and affects pressure-time history **[3]. Internal ballistics consists of the interaction of several details, which are not limited to the nominal propellant properties. Real performance may be altered by casting process, raw material lot variability, or environmental factors. This kind of behavior imposes strict requirements on acceptance criteria and reliability. For example, unexpected variations of the propellant properties due to casting effects have been highlighted by different experimental and modeling works and were found to be responsible for the hump behavior in pressure traces of small-scale rocket motors **[4–7]. Another example is represented by rocket burnout. Ideally speaking, the pressure level should drop as the burning surface reaches the propellant liner. The actual effect can also include either a pressure spike, commonly referred to as Friedman’s Curl, followed by pressure drop, or a longer burnout transient caused by grain misalignment **[3, 8]. In the case of strap-on boosting units, differential thrust or burning time might be critical at stage detachment.

The analysis of internal rocket ballistics is requested for prediction and design of rocket system performance. Different degrees of complexity can be adopted, grouped into four categories: simple, engineering, full-up, and research models **[9]. In the first group we find approaches based on equilibrium thermodynamics and zero-dimensional geometries and empirical models for burning rate and loss quantification. The second category mediates some simplifications with detailed description of relevant aspects, targeting a practical application. Examples of such codes are represented by the works published by Greatrix’s team. These solvers do not include full solution of combustion process but can model different aspects of internal ballistics such as star-grained geometry or transient burning **[10, 11]. Full-up models deal with complete physics and, in general, have commercial nature. Out of this category we mention the Rocstar code, a multiphysics multiscale computing framework used for solid rocket simulation **[12]. Research codes are more focused on physics detailed investigation than on product development.

The present paper focuses on the development of POLIRocket, an engineering model for the simulation of rocket internal ballistics and performance prediction. The code implements a quasi-1D model, including evolution of grains with complex geometry, erosive burning phenomena, compressible gas dynamics, and nonuniform propellant ballistics. In the aforementioned context a tool based on Monte Carlo method was adopted for the analysis of performance uncertainties. This tool is able to treat statistically all the major sources of uncertainties for a solid rocket motor and, on this basis, their propagation towards performances.

#### 2. Uncertainty in Rocket Motor Model

The practical interest for uncertainty analysis in rocket motors is very wide, spanning from mission analysis for a single rocket to thrust imbalance evaluation for the design of passive and active control systems in case of multiple strap-on boosters (e.g., Ariane V). In general, we refer to rocket nominal performance, such as thrust, specific impulse, or MEOP (Maximum Expected Operating Pressure). In most of applications a rocket engine never works at its nominal parameters but close to them. Statistics helps in defining the bounds of confidence for predictions which should include both the nondeterministic component of a model derived by input parameter variability or environment and the structural uncertainty introduced by the model itself and its integration. The resulting quantification supports the process of risk assessment for complex missions in mutable scenarios. Interesting approaches to uncertainty estimation and propagation for models and parameters have been proposed by Oberkampf et al. and Roy and Oberkampf **[13, 14].

The propagation of uncertainty can be performed using both Taylor series (TS) and Monte Carlo (MC) methods. The former one is an analytical approach based on Taylor series expansion and requires the definition of sensitivity coefficients for each input parameter over the final result, starting from an analytical description of the problem. The second technique is a numerical statistic method for the analysis of complex models. Modern computers can run multiple instances of the same problem on the basis of probability density functions for input variables. MC methods do not need analytical differentiations and demonstrate flexibility in terms of magnitude, type of input uncertainty distributions, and nonlinearity of the model **[15]. In this work MC method was implemented by treating numerically the rocket internal ballistics as a “black box.” Construction parameters of the SRM were defined on the basis of Gaussian distributions of known standard deviation, obtaining the population of performance data as a result of different model runs. In order to reduce the number of iterations, latin hypersquares sampling was adopted. For input distributions sampled with points, this technique reduces the number of instances from combinations to calculations **[16]. The construction parameters (ballistic coefficients, characteristic velocity efficiency, nozzle efficiency, propellant hump, propellant axis offset, and propellant mass) are assumed aleatory variables and are characterized by known Gaussian distributions.

#### 3. Model for Internal Ballistics

The engineering model implemented in this work consists of a solver for internal ballistics simulation based on quasi-1D, quasi-steady, compressible, nonviscous flow equations, coupled with a zero-dimensional nozzle. This approach includes cross section propellant grain variation in both space and time. Local secondary behaviors, like vortexes and boundary layers, are not considered. Quasi-steady model can capture the evolution in time when the solid rocket motor operates under quasi-stationary condition, which happens in most cases apart from ignition and tail-off transients. Compressible fluid dynamics of the combustion chamber is described by Shapiro’s ordinary differential equations. The model is specific for the solution of a flow in a duct. A simple control volume between two sections at an infinitesimal distance is implemented deriving a series of basic physical equations in differential logarithmic form. The resulting ordinary differential equations describe the evolution of relevant gas properties in the frame of reference. The reader is encouraged to consult Shapiro’s book for details **[17]. The model used in this work is represented by (1) to (4), under the assumptions of quasi-one-dimensional, steady, nonviscous, continuous variable flow and perfect gas:In this reduced set of equations, flow Mach number, static temperature, and pressure are dependent variables. Injected mass flow rate and cross section are independent parameters which are updated by propellant combustion modeling. The injected mass flow is proportional to the local burning rate of propellant which is expressed through Vieille’s law . The simulation of nonconstant ballistic properties in the grain is also possible. In this respect, it is common practice to introduce a correction called hump effect. This parameter depends on and introduces a spatial variation of propellant preexponential ballistic coefficient. In the present code the function has the shape . The correction factor has symmetric parabolic fashion with equation and zero integral along the grain radius . Parabola coefficients are reported in (5) and are uniquely defined by the grain radius and the value of the correction factor in the propellant midpoint:The local static pressure is resolved by (4) using an iterative process, based on comparison between mass discharge from the nozzle and production by combustion. Initial guess is evaluated by means of a zero-dimensional model of the combustion chamber (whose properties are assumed constant along -axis, at this stage). Simulations have shown that motors with low L/D ratio do not require any iteration. Once propellant is fully burnt, motor tail-off is handled by a zero-dimensional and unsteady model to compute the combustion chamber emptying. Injected mass is treated as a monophasic mixture of known features, whose thermodynamic properties (temperature, molar mass, and specific heat ratios) are tabulated as a function of pressure using a thermochemical equilibrium code **[18]. The current version of the model solves Shapiro’s ODEs using a fourth-order Runge-Kutta solver with variable integration step **[19]. The current 1D computational domain is a simplified uniform grid, even though optimized grading might apply. Each cell contains a series of information about its geometry (position, port area, and perimeter), propellant (thickness burned both in radial and lateral directions), flow information (speed, pressure, and density), and a flag that defines if the cell represents a fixed duct, a burning propellant region, or, for special cases, a side-burning propellant. Time discretization is not uniform, due to variation of combustion velocity. Equation (6) relates and combustion velocity using the coefficient which is defined by the user:For each cell, information on propellant geometry is supplied by an external module which tabulates the evolution of the initial burning surface. Burning rate can vary along the axial direction **[20]. Local correction is performed by a model for erosive combustion using the Lenoir & Robillard semiempirical equation. This process is significant when high flow speed inside combustion chamber is locally registered **[3, 21, 22]. The model can handle general types of grain shapes and assumes that the regression of solid phase boundary is locally normal and uniform at the considered axial position. An example of evolution for a star-shaped section is reported in Figure 1. The procedure can handle also cases of grain axis misalignment or offset. The code can handle both side-inhibited grains or lateral combustion phenomena.