Chinese Journal of Aeronautics (2015) xxx, xxx-xxx

Chinese Society of Aeronautics and Astronautics & Beihang University

Chinese Journal of Aeronautics

JOURNAL OF

AERONAUTICS

cja@buaa.edu.cn www.sciencedirect.com

Uncertainty analysis and design optimization of hybrid rocket motor powered vehicle for suborbital flight

Zhu Hao a, Tian Hui a *, Cai Guobiao a, Bao Weimin b

a School of Astronautics, Beihang University, Beijing 100191, China b China Aerospace Science and Technology Corporation, Beijing 100048, China

Received 2 July 2014; revised 21 October 2014; accepted 30 December 2014

KEYWORDS

Design optimization; Hybrid rocket motor; Kriging model; Uncertainty analysis; Uncertainty-based design optimization

Abstract In this paper, we propose an uncertainty analysis and design optimization method and its applications on a hybrid rocket motor (HRM) powered vehicle. The multidisciplinary design model of the rocket system is established and the design uncertainties are quantified. The sensitivity analysis of the uncertainties shows that the uncertainty generated from the error of fuel regression rate model has the most significant effect on the system performances. Then the differences between deterministic design optimization (DDO) and uncertainty-based design optimization (UDO) are discussed. Two newly formed uncertainty analysis methods, including the Kriging-based Monte Carlo simulation (KMCS) and Kriging-based Taylor series approximation (KTSA), are carried out using a global approximation Kriging modeling method. Based on the system design model and the results of design uncertainty analysis, the design optimization of an HRM powered vehicle for suborbital flight is implemented using three design optimization methods: DDO, KMCS and KTSA. The comparisons indicate that the two UDO methods can enhance the design reliability and robustness. The researches and methods proposed in this paper can provide a better way for the general design of HRM powered vehicles.

© 2015 The Authors. Production and hosting by Elsevier Ltd. on behalf of CSAA & BUAA. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

* Corresponding author. Tel.: +86 10 82338117. E-mail addresses: zhuhao@buaa.edu.cn (H. Zhu), tianhui@buaa. edu.cn (H. Tian).

Peer review under responsibility of Editorial Committee of CJA.

http://dx.doi.org/10.1016/j.cja.2015.04.015

1000-9361 © 2015 The Authors. Production and hosting by Elsevier Ltd. on behalf of CSAA & BUAA.

This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

1. Introduction

With the increasing demands for green, nontoxic and cheap propulsion technologies, hybrid rocket motors (HRMs) show great potential as they are less complex and cheaper than liquid rocket motors (LRMs), and more easily throttled and restarted than solid rocket motors (SRMs).1-3 It makes sense to develop sub-orbit vehicles with HRMs which have such advantages as safety, cheapness and non-toxicity, since the near space of 30-100 km altitude is becoming increasingly important in scientific research and military applications in recent years. Therefore, there are many academic studies and projects about sub-orbit vehicles with HRMs recently.4-7

It is necessary and important to apply design optimization methods in the aerospace vehicle design process in order to improve the design level and efficiency. In the traditional design optimization methods, the input parameters are considered as deterministic values to simplify the modeling process. However, it may be inconsistent with the objective reality. Therefore, the studies on the uncertainties in the aerospace vehicle design process have important theoretical and practical values to improve the overall design level.

Compared with the traditional SRMs or LRMs, HRMs have both a liquid oxidizer feeding system and a solid fuel combustion chamber, so the system design model of HRMs has more input variables and model parameters. Moreover, since the combustion mechanism of HRMs is not fully researched at present, there are more uncertainties in the design process of HRMs. The uncertainties probably result in the fact that the optimal design results under deterministic design optimization (DDO) are infeasible or unreliable in the following manufacturing process. Nevertheless, the current studies on design optimization of HRMs or its applications typically focus on the DDO method,4'5 so it is necessary to study the uncertainties and develop uncertainty-based design optimization (UDO) methods to enhance the design reliability and robustness. Therefore, an approach to the uncertainty analysis and design optimization of HRM powered vehicles is proposed in this paper, based on our former work about the conceptual design of HRM powered rockets.8,9

The main problem in UDO is the low efficiency of the uncertainty analysis when the system design model is complicated. The approximate model technology is one of the most popular methods to solve this problem. Kriging model is a widely used approximate model for its advantages such as unbiased estimator at the training sample point, desirably strong nonlinear approximating ability, flexible parameter selection of the model and accurate global approximation abil-ity.10,11 An approach that applying the Kriging model to two uncertainty analysis methods, i.e., Monte Carlo simulation (MCS) and Taylor series approximation (TSA), is proposed in this paper. Both newly formed methods are applied to the design optimization of the HRM powered vehicle for suborbital flight and the design results with high reliability and robustness are obtained.

2. System design model

The HRM powered sub-orbit vehicle is a ballistic rocket with an aerodynamic stable shape. The system design process involves many disciplines including structure, propulsion,

aerodynamic, launching dynamics and trajectory. Each discipline is analyzed to find out possible mathematical relationships between design variables and performance parameters, such as the rocket lift-off mass MR or the rocket body length LR, and develop a feasible multidisciplinary design model of the rocket system.

2.1. Rocket structure design

The structure of the HRM powered rocket consists of head (containing payloads), fins, HRM and the linking structures,8 as shown in Fig. 1. The rocket lift-off mass MR can be obtained by

Mr = mm + ms + m„;

where Mpay is the payload mass. The HRM mass Mm is deduced by HRM design. The rocket structure mass Ms consists of head mass, fin mass and linking structures mass. It is related to rocket diameter DR and defined as

Ms = 75Dr - 7.5 (2)

The rocket body length LR can be obtained by

Lr = Lm + Lh (3)

where Lh is the rocket head length and it is 1 m in this paper. The HRM length Lm is also deduced by HRM design.

2.2. HRM design

HRM is the main part of the rocket. Its mass and dimension almost determine the mass, dimension and trajectory of the rocket. A wheel port grain is selected in the HRM with a pro-pellant combination of 98% hydrogen peroxide (HP) and hydroxyl-terminated polybutadiene (HTPB). A nitrogen gas pressure feed subsystem is used. The oxidizer mass flow rate is controlled to keep constant by an ideal venturi section. The HRM is designed as shown in Ref.8 The propellant, including the solid fuel and liquid oxidizer, constitutes the

Fig. 1 HRM powered rocket structure.

main mass and dimension of the HRM, therefore it determines MR and LR indirectly.

2.3. Launch simulation

A 35 m length ramp is used in the launch subsystem. The dynamic equation of the ramp launch process can be obtained based on Newton's Second Law as follows:

V = F/M - g sin uo - kg cos uo (4)

where V is the rocket velocity, F the thrust, M the rocket mass that changes as the HRM combusts, g gravity constant and u0 the launch elevation angle which needs to be optimized. The friction coefficient k is 0.05. The aerodynamics drag is ignored because the rocket velocity is slow and the launch phase time is short.

2.4. Trajectory simulation

A two degree of freedom (DOF) mass point trajectory equation is used with no control law as follows:

Table 2 Design variables.

Design variable Lower Upper Distribution Relative

limit limit limit deviation

Grain outer 0.2 0.5 Normal 0.005

diameter Dp (m)

Initial web 0.02 0.05 Normal 0.01

thickness e1 (m)

Number of wheel 3 10 Not

port n uncertain factor

Initial thrust Fi 10 20 Uniform 0.05

Initial chamber 1 3 Uniform 0.05

pressure pc

Initial oxidizer to 2 5 Uniform 0.05

fuel ratio a

Nozzle expansion 3 10 Normal 0.015

ratio e

Launch elevation 65 85 Normal 0.005

angle U0 (°)

X — (F - CDqSM)X/V/M Y — (F - CDqSM)Y/V/M - g

where X and Y are the horizontal and vertical distances of the position from the launch point, q is the dynamic pressure and CD represents the drag coefficient. Rocket body section area is used as the reference area SM.

2.5. Dynamic computation

A U.S. 76 standard atmosphere model with no atmosphere motion is used here. The drag coefficient CD of ''Titan II'' rocket12 is used as shown in Table 1, given a 20% increase as the sub-orbit rocket in this paper has fins.

Based on the analysis and deduction above, the multidisci-plinary design model of the rocket system is established. There are 43 input parameters, including 8 design variables (shown in Table 2) and 35 model parameters (shown in Table 3), in the mathematical model through which the performance response parameters, such as MR and LR, are computed.

3. Design uncertainty analysis

3.1. Design uncertainties

In design phase, specifically the simulation-based computational design, there are three sources contributing to the total uncertainty of computational simulation, namely model input uncertainty, model uncertainty, and model error.13 The model

Table 1 Drag coefficient of USA "Titan II rocket''

0 6 Ma < 0.8 0.29

0.8 6 Ma < 1.068 Ma - 0.51

Ma p 1.068 0.091 + 0.5 Ma-1

input uncertainties constitute the main part of the total uncertainties, so only they are considered in this paper. The model input uncertainties are the uncertainty factors in the design variables and model parameter. Most of them are emerged from the engineering realization process. Except for the number of wheel port n, all the other 42 input parameters are uncertain factors in the hybrid rocket system design model. The most accurate methodology to classify and quantify these uncertain factors is probability theory, which is applied in this study. The uncertainties in the three types of the input parameters, including physical or chemistry characteristic parameters such as fuel density, machining parameters such as grain outer diameter and parameters obtained by testing(such as injector pressure drop coefficient), are classified to be normal distribution. The uncertainties in the other parameters, which are either generated from the cognitive incompletion of the physical word such as the regression rate, or determined by these parameters directly or indirectly such as the initial thrust, are classified to be uniform distribution. According to design criteria and engineering experience14-16, the relative limit deviations of the design uncertainties are confirmed as shown in Tables 2 and 3. The standard deviations of uncertainties with normal distribution can be evaluated in terms of "6r" principle.

At each step of the optimization process, the uncertainties of design variables can be quantified as

; = -x'A

xU — x'A

where A is the relative limit deviation; XL and XU are the upper limit deviation and the lower limit deviation of the design variable when its value is x .

3.2. Sensitivity analysis

It makes the uncertainty analysis complicated and time-consuming with a lot of input uncertainties, so a sensitivity analysis method17 is used to filter out the insignificant model

Table 3 Model parameters.

Model parameter Symbol Mean Limit deviation Distribution

Regression rate precision coefficient x1 1 0.02 Uniform

Weld coefficient x2 0.8 0.03 Uniform

Pressure oscillation coefficient x3 1.2 0.03 Uniform

Volume fraction of remained oxidizer x4 0.05 0.05 Uniform

Volume fraction of initial air cushion x5 0.05 0.05 Uniform

Combustion efficiency x6 0.96 0.015 Uniform

Nozzle efficiency x7 0.93 0.015 Uniform

Initial tank temperature (K) x8 293.15 0.0682 Uniform

Injector pressure drop coefficient x9 0.2 0.0325 Normal

Tube pressure loss coefficient x10 0.2 0.08 Normal

Initial gas bottle pressure (MPa) x11 30 0.03 Normal

Fuel density (kg/m3) x12 1218 0.015 Normal

Oxidizer density (kg/m3) x13 1440 0.015 Normal

Chamber heat insulation layer density (kg/m3) x14 1000 0.015 Normal

Chamber and nozzle shell density (kg/m3) x15 7750 0.015 Normal

Chamber head density (kg/m3) x16 7750 0.015 Normal

Tank shell density (kg/m3) x17 2850 0.015 Normal

Gas bottle shell density (kg/m3) x18 1750 0.015 Normal

Semi-minor axis length ratio of chamber head x19 3 0.01 Normal

Semi-minor axis length ratio of oxidizer tank head x20 2 0.01 Normal

Semi-minor axis length ratio of gas bottle head x21 2 0.01 Normal

Chamber heat insulation layer thickness (m) x22 0.003 0.025 Normal

Nozzle heat insulation layer thickness (m) x23 0.015 0.025 Normal

Injector panel thickness (m) x24 0.005 0.02 Normal

Minimum machining thickness of chamber and nozzle shell material (m) x25 0.0015 0.066 Normal

Minimum machining thickness of oxidizer tank shell material (m) x26 0.0015 0.066 Normal

Minimum machining thickness of gas bottle shell material (m) x27 0.005 0.025 Normal

Nozzle half expansion angle (°) x28 15 0.012 Normal

Nozzle half convergence angle (°) x29 45 0.012 Normal

Yield limit of chamber and nozzle shell material (MPa) x30 1080 0.01 Normal

Yield limit of chamber head material (MPa) x31 1080 0.01 Normal

Yield limit of oxidizer tank shell material (MPa) x32 490 0.01 Normal

Yield limit of gas bottle shell material (MPa) x33 1760 0.01 Normal

Ramp length (m) x34 35 0.0001 Normal

Drag precision coefficient x35 1 0.1 Normal

parameter uncertainties in the UDO approach of hybrid rocket design. According to the probability distributions, the model parameter uncertainties are sampled 1000 times respectively with a Latin hypercube sampling (LHS) method.18 Then the main performance parameters are computed by the system design model. The sensitivities of the rocket's main performance parameters, including MR, the top altitude of trajectory Ymax, the maximum rocket axes acceleration Nxmax, and the rocket length to diameter ratio L/D are obtained by the sensitivity analysis method. The top ten model parameter uncertainties that have the most significant effects on MR, Ymax, Nxmax, and L/D are shown respectively from Figs. 2-5. The length of the bar in the figures represents the sensitivity of the performance parameters to model parameter uncertainties. The positive value of bar means that the performance parameter increases when the model parameter uncertainty increases and vice versa.

When the oxidizer mass flow rate is constant, the regression rate determines the burning time, which directly determines the oxidizer mass Mo; Mo has an important influence on MR and Lr as discussed in Section 2.2. Therefore, the fuel regression rate precision coefficient (x1) has the greatest impact on MR and L/D indirectly as shown in Figs. 2 and

x25 ZZJ

1 , \ 1

■0.10 -0.05 0 0.05

Fig. 2 Sensitivity of MR to model uncertain factors.

3. The oxidizer density determines the oxidizer mass which influences the rocket structure mass. The yield limit of oxi-dizer tank shell material determines the tank wall thickness which influences the tank mass Mt consequently. In addition, Nxmax and Ymax are mainly determined by the rocket velocity, attitude and altitude at the burnout point which is influenced by rocket structure mass according to Tsiolkovski formula. Therefore oxidizer density (x13) and the yield limit of oxidizer tank shell material (x32) have

x4 =□

xlO ■1

■0.20 -0.15 -0.10 -0.05 0 0.05

Fig. 3 Sensitivity of L/D to model uncertain factors.

x25 —

x23 x30

x32 xl3 -

-0.02 -0.01 0 0.01 0.02

Fig. 4 Sensitivity of Nxmax to model uncertain factors.

x23 x30

x32 xl3 -

-0.02 -0.01 0 0.01 0.02

Fig. 5 Sensitivity of Ymax to model uncertain factors.

the greatest impact on Nxmax and Ymax indirectly as shown in Figs. 4 and 5. At the current level of science and technology, the uncertainty of material density and strength can be quantified accurately, while the fuel regression rate model of HRM is not investigated clearly enough. Most uncertain data of regression rate is mainly obtained through certain test of a special HRM and it is not suitable for all HRMs. Therefore uncertain factors in the regression rate are the main source of the rocket performance uncertainties using an HRM.

Through the above sensitivity analysis, three model parameter uncertainties, including the fuel regression rate precision coefficient (x1), oxidizer density (x13) and the yield limit of oxidizer tank shell material (x32) (shown in Table 3 with boldfaced words), together with 7 design variables (except for n) are selected as the uncertain factors in the following UDO approach.

4. Kriging-based UDO

The Kriging-based uncertainty analysis methods are proposed in this section. The general UDO process and two uncertainty analysis method including MCS and TSA are discussed, then a global approximation Kriging modeling approach is presented. Two newly formed uncertainty analysis methods including the Kriging-based Monte Carlo simulation (KMCS) and Kriging-based Taylor series approximation (KTSA) are proposed sequentially.

4.1. UDO

The tradition optimization can be expressed in a mathematics model as

{find x

min f(x,P) 8

s.t. g(x, p) 6 0 ( )

XL 6 x 6 XU

where x is the design variables and p the model parameters, both of them are input parameters; xL and xU are the lower and upper boundaries of x, f(x,p) is the objective function and g(x,p) represents the constraint vector. In DDO, all of x;p, f(x,p) and g(x,p) are treated as deterministic parameters, and the optimal organization and the search strategy are based on the deterministic relationships as a result. The benefits are the simplification of the optimal process and the saving of computing time.

However, there are many uncertainties in the real world. In UDO, all design variables, model parameters and mathematical models are analyzed and the uncertain factors between them are separated. The uncertainties of the performance responses can be computed with different methods sequentially. The mathematical model of UDO can be expressed as

{find x

min F[if(x, p), rf(x, p)]

s.t. P[g(x, p) 6 0] p R ( )

xL 6 x 6 xU

where both x and p could be uncertain; if and Gf are the mean and standard deviation of the original optimization objective function f, and F is the reformulated optimization objective function with respect to if and Gf; P is the probability of the statement within the braces to be true, and R is the reliability vector specified for the constraint vector. The robustness of the system is achieved through minimizing if and Gf, and the reliability of constraints is accomplished through confidence level at which constraints are met with a higher probability.

Uncertainty analysis is a key procedure of UDO process. At this step, the uncertainty characteristics of the system responses propagated from the input uncertainties are quantitatively analyzed. There are many uncertainty analysis methods, including MCS, TSA, reliability analysis, etc.

MCS method is the most accurate solution for uncertainty problems based on probability theory. The mean value, standard deviation, distribution function and probability density of the responses are predicted statistically from the random sampling metric in MCS analysis. if and Gf are obtained as

J2f(x>)

f(xi)- if)

where f(xi) is the response function about the xt sample point. The prediction accuracy of MCS analysis is inversely proportional to the square root of the sampling number N, thus extensive sample number is needed to ensure the prediction accuracy. As a result, it is unacceptable for the optimization problems with a long-running time simulation program.

TSA method is one of the most efficient solutions for uncertainty problems with probability theory. Based on the firstorder Taylor series, if can be estimated as13

if — E(f(x)) « f(1x

where ix is the mean value of n-dimensional vector x. If the input variables are not related, Of can be estimated as13

where Oxi is the standard deviation of xi. The efficiency of TSA method is much higher than that of the MCS method, since it does not need repeated calculation. However, when f(x) is a nonlinear system whose first-order gradient cannot be obtained with analytical method, the application of TSA is restricted. How to obtain an accurate gradient value with low time cost is the key point to solve this problem.

4.2. Global approximation Kriging modeling

The Kriging model is an interpolation technique based on statistical theory. Its advantages, such as unbiased estimator at the training sample point, desirably strong nonlinear approximating ability, flexible parameter selection of the model and accurate global approximation ability, make it widely used in approximate models. It takes full account of the relevant characteristics of the variable space, containing the regression part and the nonparametric part

y(x) —f(x)+z(x)

where f(x) is deterministic function that is a global approximation of the design space represented by the polynomial of x; z(x) is a Gaussian stochastic process with zero mean and variance. Ref.10 showed the detailed developing process of Kriging model with N sampling points. When using Gauss function as the correlation function in Kriging model, the first-order derivative at point xi can be estimated by

dy(x) dxi

df(x) + dz(x)

The approximation accuracy of the approximate model is closely related to the quantity of sample points, thus a global approximation Kriging modeling method is developed by sequentially sampling the design space and update the

Kriging model in order to get higher approximation accuracy with fewer sample points. The process is shown in Fig. 6 and the main steps are shown as follows:

Step 1. Use LHS method to generate the initial training

points with a small scale.

Step 2. Establish the initial Kriging model.

Step 3. Use LHS method to generate the test points.

Step 4. Calculate the response values using the original

model and the Kriging model respectively, then compare

the results to get the approximation error of the Kriging

model.

Step 5. According to the results of Step 4, judge if the error satisfies the precision criterion. The adjusted multiple correlation coefficient R2a is used as the global precision criterion and the maximum relative error emax is used as the local precision criterion.19 If yes, end the iteration and output the present Kriging model; else, goto the next step. Step 6. Choose the test points at which the relative errors are larger than emax and put them into the training space to update the Kriging model in order to enhance its approximation accuracy at the design areas where the test points have large relative errors with the previous Kriging model. Step 7. Return to Step 2, continue the iteration until the Kriging model satisfies the desired precision.

There are many verifications about the approximation accuracy of Kriging model, so only the prediction precision of the first-order derivative is tested by a symmetric two-dimensional nonlinear function with multiple local extreme points20 shown as

y(x1; x2) — 2 + 4x1 + 4x2 - x? - x2 + 2sin(2x1) sin(2x2) x1; x2 2 [0.5,3.5]

The number of the initial training points is 10 to develop the Kriging model of the test function and the results are shown in Table 4. Along with the increase of the training point

Fig. 6 Procedures of global approximation Kriging modeling process.

Table 4 Iteration process of Kriging model for the test function.

Iteration Training points Test points Ra emax

0 10 20 0.6934 0.3302

1 27 20 0.9832 0.0555

2 36 20 0.9961 0.0384

3 37 20 0.9995 0.0189

number m, R2a increases to 1 gradually and emax decreases. When the number of training points reaches 37, Ra2 is close to 1 and emax nearly decreases to 0, indicating that Kriging model has reached a good prediction precision of the firstorder derivative as shown in Figs. 7 and 8.

4.3. Design optimization with KMCS

MCS is considered to be the most accurate method among the uncertainty analysis methods based on probability theory;

however its computation precision is directly related to sampling frequency, thus a great number of sampling points are needed when MCS is used. It is unacceptable for the simulation-based design optimization problems whose system design model is complicated and time-consuming. Therefore the KMCS method is developed in which the original complicated system design model f(x) used in MCS is surrogated by the approximate function f(x) established by the global approximation Kriging modeling method. The approximate mean value if and the standard deviation if of the system responses in KMCS can be computed using Eqs. (17) and (18). Fig. 9 shows the procedure of UDO with KMCS method.

X/(xi)

^(y(xi)- /

(a) First-order derivative of test function to (b) First-order derivative of Kriging model to *

Fig. 7 First-order derivative of test function and its Kriging model to x1.

(a) First-order derivative of test function to x2 (b) First-order derivative of Kriging model to x.

Fig. 8 First-order derivative of test function and its Kriging model to x2.

Fig. 9 Procedures of UDO with KMCS method.

4.4. Design optimization with KTSA

The Kriging model has a high prediction accuracy not only on the function value but also its first-order derivatives, thus the Kriging model is applied to obtaining the first-order gradient value of f(x) used in the TSA method and it can effectively solve the problems that generated when TSA is applied to multidimensional nonlinear systems. With a Kriging model accurate enough, the KTSA method can effectively expand the application field of TSA method. The mean value and the standard deviation of the system responses can be computed using Eqs. (19) and (20), where is computed by Eq. (15). The "6o" principle is used to compute constraint satisfaction probability as shown in Eqs. (21) and (22), where G is the reliability level, ig the mean of the constraint functions, og the standard deviation of the constraint functions, and k Sigma level. For the normal distribution function, when k = 6, constraint satisfaction probability is 99.9999998%; when k = 3, constraint satisfaction probability is about 99.9937%. k = 3 is chosen in this paper. Fig. 10 shows the procedure of UDO with KTSA method.

Ay = y(iJ

1g1 (X; p)-k'gi(X; p) P G1

(20) (21)

Fig. 10 Procedure of UDO with KTSA method.

1g2 (X; P)+ kog2 (X; P) 6 G2 (22)

5. Optimization, results and discussion

5.1. Design variables, target function and constraints

The general mission of sub-orbit rocket is to detect the high-altitude environment or supply the payload with a certain microgravity time. A HRM powered rocket with a mission to send a 50 kg payload to an altitude over 120 km is designed in this study. The payloads are often avionic equipment which is sensitive to acceleration, so the maximum rocket axial acceleration Nxmax is constrained. One of the characteristics of the rocket body with a series-wound structure HRM is its great length, so a value of 18 is chosen to constrain the rocket length to diameter ratio L/D, considering that a too high L/D is not good for the rocket structure strength. Both the design difficulty and the cost of a rocket vehicle are mainly determined by MR, so the target function is to minimize MR by satisfying the constraints above. All the design variables x and their boundaries xL, xU are shown in Table 2 and the three model uncertain parameters are shown in Table 3 with boldfaced words. The mathematics model of the DDO is shown in Eq. (23).

Table 5 Design results of HRM powered rocket.

Design variable Dp (m) e1 (m) n F¡ (kN) Pci (MPa) a¡ e UO (о) Run time(min)

DDO 0.267 0.0322 4 14.30 1.89 2.66 4.61 85.0 499

KMCS 0.277 0.0316 4 14.31 1.93 2.79 4.51 85.0 1887

KTSA 0.277 0.0417 4 13.65 1.81 2.33 4.56 85.0 281

Table 6 Statistical results of HRM powered rocket.

Target/constraints Method 1 r Maximum Minimum 1/r G

MR (kg) DDO 333.97 2.80 342.26 326.31 0.00837

KMCS 353.36 2.89 362.62 344.92 0.00819

KTSA 364.41 2.93 372.95 356.37 0.00803

Ymax (km) DDO 119.28 4.02 131.18 107.51 0.429

KMCS 127.04 4.23 139.47 114.85 0.954

KTSA 128.59 4.56 141.94 115.89 0.978

L/D DDO 17.97 0.14 18.40 17.58 0.589

KMCS 17.50 0.14 17.94 17.12 1

KTSA 17.57 0.13 17.92 17.22 1

Nxmax (g) DDO 9.23 0.25 9.77 8.67 1

KMCS 9.23 0.25 9.78 8.62 1

KTSA 8.23 0.23 8.78 7.78 1

(d) CDF of LID

Fig. 11 Cumulative distribution function (CDF) of the target and constraints.

find x

min Mr = f(x, p)

s.t. Ymax = g1 (x, p) p 120 km

L/D = g2(x,p) 6 18 ( )

Nxmax = g3(x,p) 6 10 g

xL 6 x 6 xU

In UDO approach, the robustness of the system is achieved through minimizing the mean value iM and standard deviation GMR of the target function while the reliability of constraints is accomplished through the confidence-level that the constraints are met with a higher probability, shown as Eq. (24).

find x

min F[iMr(x,p); GMR(xp)\

s.t. P\ [Ymax(x,p) p 120km] p 0.95

P2[L/D(x,p) 6 18] p 0.95 P3[Nxmax(x,p) 6 10g] p 0.95 xL 6 x 6 xU

5.2. Results and discussion

A modified differential evolution (MDE) algorithm21 is applied to implementing global optimization and improving the efficiency and quality of the optimization solution, then DDO, KMCS and KTSA methods are applied respectively to the design optimization of the HRM powered rocket for a suborbital flight. The design variables n and u0 are set to be constant and equal to the DDO optimal result to simplify the calculation. The comparisons of the results are shown in Table 5, and the statistical results considering all the uncertain input parameters are shown in Table 6 and Fig. 11. Compared to the DDO results, the two UDO methods achieve reliability requirements at a higher confidence level and r/i of MR reduces by 8.0% as shown in Table 6. The reason that reliability of Nxmax satisfying the constraint is 100% in DDO method is that optimal result is not at the boundary between feasible and unfeasible region as shown in Fig. 11. Although the mean value of the rocket lift-off mass is a little bigger than that of deterministic one, the reliability and robustness are enhanced obviously with two UDO method. Compared with KTSA, the results of KMCS are comparatively better while its efficiency is lower. The prediction precision of KTSA is not better than KMCS inherently since the former applies two approximate process, Kriging model and TSA method. As a result, the KTSA method is more suitable for the initial design optimization phase while the KMCS method is more applicable for the detailed design optimization phase. The results and comparisons prove that the uncertainty design optimization methods can also provide a better means for system conceptual design of aerospace vehicles.

6. Conclusions

(1) The multidisciplinary system design and analysis model of the HRM powered rocket is established and the input uncertain factors are quantified. The sensitivity analysis of the uncertain factors shows that among 42 uncertain factors the regression rate uncertainty has the most

significant effect on performances of the HRM powered rocket, thus it is necessary to accelerate the investigation on the combustion mechanism in HRMs.

(2) Two newly formed uncertainty analysis method including the KMCS and KTSA are carried out with a global approximation Kriging modeling method. The design optimization of the HRM powered rocket is carried out applying three methods including DDO, KMCS and KTSA. The results and comparisons show that the two UDO methods can provide design results with a higher reliability and robustness than DDO method and the KTSA method is more suitable for the initial design optimization phase while the KMCS method is more applicable for the detailed design optimization phase. The uncertainty design optimization methods can provide a better means for system concept design of aerospace vehicles.

(3) Due to the insufficiency of experiment or engineering data about the uncertain factors, there may be some inaccuracy on the related distributions and relative limit deviations. Our future work will focus on finding sufficient data to achieve accurate results. It is another important work to promote the engineering application of the UDO approach for HRM powered vehicles.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (No. 51305014) and China Postdoctoral Science Foundation (No. 2013M540842).

References

1. Altman D. Hybrid rocket development history. Proceedings of 27th AIAA/SAE/ASME/ASEE joint propulsion conference; 1991 Jun 24-26; Sacramento. Reston: AIAA; 1991.

2. Davydenko NA, Gollender RG, Gubertov AM, Mironov VV, Volkov NN. Hybrid rocket engines: the benefits and prospects. Aerosp Sci Technol 2007;11(1):55-60.

3. Li XT, Tian H, Cai GB. Numerical analysis of fuel regression rate distribution characteristics in hybrid rocket motors with different fuel types. Sci Chin Technol Sci 2013;56(7):1807-17.

4. Casalino L, Pastrone D. Optimal design of hybrid rocket motors for microgravity platform. J Propul Power 2008;24(3):491-8.

5. Kosugi Y, Oyama A, Fujii K, Kanazaki M. Multidisciplinary and multi-objective design exploration methodology for conceptual design of a hybrid rocket. Proceedings of Infotech@Aerospace; 2011 Mar 29-31; St. Louis, Missouri. Reston: AIAA; 2011.

6. Dyer J, Doran E, Dunn Z, Lohner K, Bayart C, Sadhwani A, et al. Design and development of a 100 km nitrous oxide/paraffin hybrid rocket vehicle. Proceedings of 43rd AIAA/ASME/SAE/ASEE joint propulsion conference and exhibit; 2007 Jul 8-11; Cincinnati. Reston: AIAA; 2007.

7. Kniffen R J, McKinney B, Estey P. Hybrid rocket development at the American Rocket Company. Proceedings of 26th AIAA| ASME/SAE/ASEE joint propulsion conference; 1990 Jul 16-18; Orlando. Reston: AIAA; 1990.

8. Cai GB, Zhu H, Rao DL, Tian H. Optimal design of hybrid rocket motor powered vehicle for suborbital flight. Aerosp Sci Technol 2013;25(1):114-24.

9. Rao DL, Cai GB, Zhu H, Tian H. Design and optimization of variable thrust hybrid rocket motors for sounding rockets. Sci Chin Technol Sci 2012;55(1):125-35.

10. Xuan Y, Xiang JH, Zhang WH, Zhang YL. Gradient-based Kriging approximate model and its application research to optimization design. Sci Chin Technol Sci 2009;52(4):1117-24.

11. Wang H, Wang SP, Mileta MT. Modified sequential Kriging optimization for multidisciplinary complex product simulation. Chin J Aeronaut 2010;23(5):616-22.

12. He LS. Ballistic missiles and launch vehicles design. Beijing: Beihang University Press; 2002.

13. Yao W, Chen XQ, Luo WC, van Tooren M, Guo J. Review of uncertainty-based multidisciplinary design optimization methods for aerospace vehicles. Prog Aerosp Sci 2011;47(6):450-79.

14. GJB 1026A-99, General specification for solid propellant rocket motors. Military standardization of People's Republic of China, 1999 Chinese.

15. GJB 2053A-2008, Specification for aluminum allot structural sheet for aerospace. Military standardization of People's Republic of China, 2008 Chinese.

16. Wen BC. Machine design handbook. Beijing: China Machine Press; 2010 Chinese.

17. Lucas JM. How to achieve a robust process using response surface methodology. J Qual Technol 1994;26(4):248-60.

18. Zhang RC, Wang ZJ. Design theory and data analysis of computer experiments. Chin J Appl Probab Stat 1994;10(4): 420-36 Chinese.

19. Xie YM. Research on robust optimization of sheet metal forming based on Kriging and grey relational analysisdissertation. Shanghai: Shanghai Jiao Tong University; 2007 [Chinese].

20. Rijpkema JJM, Etman LFP, Schoofs AJG. Use of design sensitivity information in response surface and Kriging metamodels. Optim Eng 2001;2(4):469-84.

21. Rao DL, Cai GB. Modified differential evolutionary algorithm for fast simulation optimization and its application. J Astronaut 2010;31(3):793-7 Chinese.

Zhu Hao is a postdoctoral researcher at School of Astronautics, Beihang University. His main research interests are design optimization method and general design of hybrid rocket motor powered vehicle.

Tian Hui is an associate professor at School of Astronautics, Beihang University. His main research interest is hybrid propulsion technologies.