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.
Efficiency of the characteristic velocity and that of thrust coefficient are figures of merit of combustion chamber and nozzle. These values are user-defined constants in the present implementation. The code can be easily extended to account for further sources of losses such as two-dimensional nozzle exhaust, two-phase flow, or chemical kinetics using compact semiempirical formulations [23].
Code parameters have been set to reproduce the ballistic behavior of a BARIA SSRM (small-scale rocket motor), used for the experimental evaluation of propellants by Avio S.p.A. (Colleferro, Italy). A scheme of the rocket is reported in Figure 2. Internal grain is a simple central-perforated grain without lateral inhibition, designed to have a quasi-neutral burning surface and, consequently, quasi-constant pressure trace. Graphite convergent nozzle with different diameters is used to change the mean operating pressure. Igniter and pressure gauge are placed at the motor head-end. The propellant cartridge is 290 mm long. Inner diameter is 100 mm and propellant web thickness is 30 mm.
Validation tests were performed by comparing internal pressure trace from model results with experimental data obtained from BARIA rocket motor, finding a general good agreement for different pressure levels. Thrust is disregarded since data are not recorded during tests. One example is reported in Figure 3. Typical simulation parameters are listed in the next section.
4. Simulations for Uncertainty Propagation
An MC method was implemented for the analysis of uncertainty propagation towards performances, starting from construction parameters. The solver for BARIA internal ballistics was used. The aforementioned motor model was modified introducing a divergent section with area ratio of 4 to gain some sensitivity on thrust and specific impulse. Presented simulations are performed assuming a nozzle diameter of 25.25 mm, which corresponds to an expected mean operating pressure of about 45–50 bar. Gaussian distributions were assumed for Vieille’s parameters and , and efficiency, hump effect, grain length (), and offset of axis perforation (). Nominal values and relevant standard deviations of input populations, derived from comparison with available data, are reported in Table 1. Erosive combustion is not considered in the present work.
Performance parameters monitored in this work consisted in burning time, gravimetric specific impulse, total impulse, mean thrust, MEOP (Maximum Expected Operating Pressure), and mean pressure. Pressure traces were analyzed using a thickness-over-time postprocessing method [8]. Some of them (mean thrust, gravimetric specific impulse, and total impulse) are important for system-level considerations while other data gain particular interest in mission profile evaluation (burning time, mean pressure, and Maximum Expected Operative Pressure). The axial geometry was discretized with 500 evenly spaced grid points. Sensitivity study, considering up to 1000 intervals, is reported for the burning time in Figure 4.
MC method accuracy depends on the chosen population size. Generally, a bigger population will give more accurate statistics results; on the other hand the computational effort will grow accordingly. A compromise is usually requested. For this reason an analysis using sample mean value and estimated standard deviation was conducted to choose a proper population size. Burning time was arbitrarily chosen as monitoring variable. The convergence of the numerical method is obtained when both mean value and confidence interval become insensitive from the size of the input population. If this happens, the standard deviation of the results is generated by the sole propagation of the original uncertainty. Representative results are reported in Table 2.
We select 200 data points since the algorithm reaches a satisfactory convergence for our purposes. In the investigated range the variation of the mean value changes about 0.1% while the variation on the confidence interval is less than 3%. The choice of 200 points represents a compromise between accuracy and computational time, even though better results could be obtained for a larger population. With these parameters, the Matlab® run of the MC procedure took about 8/9 hours on a desktop PC.
Simulations were performed varying one parameter per each simulation and keeping a nominal value for the other input parameters. Populations of 200 samples were generated using Gaussian distributions for propellant ballistics ( and ), characteristic velocity , nozzle efficiency , hump level, propellant length (which turns into a variation on loaded mass), and propellant axis offset. Standard deviations used for the generation are reported in Table 1. Simulations were postprocessed to extract, from output populations, the relative standard deviation with respect to the mean value, in percent. Data are reported in Table 3. The last column reports results for a simulation where all input parameters contribute to uncertainty propagation.
Results allow understanding how model input parameters influence accuracy of performance predictions. The uncertainty of coefficients and does not propagate into gravimetric specific impulse and total impulse statistics; rather, it has substantial influence over the variability of other performance data. The preexponential coefficient amplifies its initial relative standard deviation, while the pressure exponent slightly decreases it. The uncertainty on characteristic velocity causes variability in all mapped parameters, less evident in burning time and more pronounced for thrust and pressure data. Thrust coefficient refers to nozzle expansion, so performances connected to the combustion chamber are not affected. It is important to notice that major uncertainties on impulses are caused by efficiencies. The reproducibility of the hump effect causes only a significant alteration in MEOP statistics while it does not procure other correlated effects; this is consistent since the hump is a local effect. For the tested range, variability of grain geometry does not have significant effect on gravimetric and volumetric specific impulses. The offset of grain axis causes a change of the pressure tail-off, influencing burning time and, as a consequence, total impulse computation. When tail-off becomes longer, an increased portion of the pressure curve is below the threshold value used by the TOT method to identify the burnout instant and is not considered in the relevant integral. Grain length is connected to propellant mass load as well as to geometric shape. The ballistic model amplifies the uncertainty of this parameter, propagating mainly on thrust and pressure.
The uncertainty resulting from single sources can be combined and compared with the standard deviation from the global MC simulation. According to TS method, the global standard deviation results from the formula , under the hypothesis of complete independency between uncertainty sources [16]. If not, a correlation term should be included. Table 4 reports the comparison between global standard deviation of parameters obtained from TS and MC. It appears that uncertainty contributions to most of the performance parameters are independent of each other, according to the present model. Only a minor correlation is visible for impulses. Conversely, dependence is observed in MEOP statistics and the correlation between error sources shrinks the actual uncertainty, with respect to the TS prediction.
5. Parametric Estimation of Main Input Uncertainties
The evaluation of incremented uncertainty level for some fundamental rocket parameters derived from experimentation is reported in the present paragraph. The level of accuracy for model input data is fundamental for the determination of the prediction confidence bounds of performances. The scalability of such interval with respect to input quality is of interest when experimental campaigns are planned or literature data are available. In the following paragraph the standard deviation for the population of three input parameters is varied, with respect to a reference value reported in Table 1. All sets of simulations were run with only one source of uncertainty per time. The other parameters were assumed nominal.
The first group of simulations was relevant to the ballistic parameters and . Resulting statistics of rocket performance are listed in Tables 5 and 6 for input populations having 50%, 100%, and 200% the respective reference standard deviation. In both cases, the input variability is reduced across the model. As already observed, there is a substantial independence of impulses from ballistic coefficients. Conversely, there is a more than linear correlation between input and output uncertainty of the other parameters for the preexponential coefficient and a less than linear effect when variability of pressure exponent is considered.
A similar analysis was conducted on the efficiency of the characteristic velocity. Resulting uncertainty on performance parameters is reported in Table 7. The width of the input confidence interval decreases for burning time and for impulses, while it is amplified for pressure and thrust data. Similar results were obtained in the previous section. Variation of output standard deviations is in line with input ones. This observation is valid for all tested parameters in this model.
6. Conclusion
The work presented a model for the description of rocket internal ballistics, based on Shapiro’s quasi-1D formulation. The tool is a robust and simple code capable of dealing with a wide variety of phenomena that occur inside a solid rocket motor. Current implementation is based on Matlab and can easily be extended to include new features. Thanks to its capability of running fast simulations, the code is suitable for use in a MC algorithm. The numerical model made during this work proved out to be suitable for the approximation of the BARIA motor, a small-scale rocket motor used at industrial level. In the present framework, the code was used to perform an investigation on uncertainty propagation of input parameters, based on MC method.
Uncertainty on propellant ballistic coefficients were demonstrated to influence burning time, pressure, and thrust, while no effect was observed on impulses. Efficiency of characteristic velocity propagated across all mapped quantities, while thrust coefficient variability affected only quantities correlated to nozzle expansion. The only relevant effect of hump reproducibility was observed on the MEOP statistics. Variability of central perforation propagated to burning time and total impulse. The grain length was connected to propellant loaded mass and burning surface, thus influencing pressure, burning rate, thrust, and total impulse.
Future development of the POLIRocket code consists of model refinement, without losing the simplicity of the original framework. Thanks to Shapiro’s equations, the inclusion of wall friction, heat exchange, chemical reactions, and phase changes is possible. The statistical framework developed so far can be applied in industrial environment by comparing predictions to resulting production uncertainty based on real firing data. Moreover, the indications obtained by the code about uncertainty propagation can be used to understand which effort for better knowledge of input data should be developed to improve global prediction accuracy. In this respect, parallel execution of Monte Carlo method will enable faster evaluations with larger populations.
Nomenclature
| Area | |
| Preexponential coefficient of Vieille law | |
| Coefficient after hump correction | |
| Time coefficient | |
| Characteristic velocity | |
| Hump correction factor | |
| Hump peak coefficient | |
| Specific impulse | |
| Total impulse | |
| Volumetric specific impulse | |
| Specific heat ratio | |
| Grain length | |
| Mach number | |
| Mass flow rate | |
| Size of population | |
| Exponential coefficient of Vieille law | |
| Pressure | |
| Mean pressure | |
| Radial distance of consumed propellant | |
| Grain radius | |
| Burning rate | |
| Estimate of the standard deviation | |
| Temperature | |
| Time | |
| Burn time | |
| Speed | |
| Cartesian coordinate | |
| Sample mean for statistical analysis | |
| Coefficients for hump correction | |
| Efficiencies | |
| Population mean value | |
| Standard deviation | |
| MC: | Monte Carlo |
| MEOP: | Maximum expected operating pressure |
| SRM: | Solid rocket motor |
| TOT: | Thickness over time |
| TS: | Taylor series |
| Thrust coefficient. |
Competing Interests
The authors declare that they have no competing interests.