# Robust Design Optimization of an Aerospace Vehicle Prolusion System

Authors: Muhammad Aamir Raza; Wang Liang
Journal: Mathematical Problems in Engineering
 (2011)

---

## Abstract

This paper proposes a robust design optimization methodology under design uncertainties of an aerospace vehicle propulsion system. The approach consists of 3D geometric design coupled with complex internal ballistics, hybrid optimization, worst-case deviation, and efficient statistical approach. The uncertainties are propagated through worst-case deviation using first-order orthogonal design matrices. The robustness assessment is measured using the framework of mean-variance and percentile difference approach. A parametric sensitivity analysis is carried out to analyze the effects of design variables variation on performance parameters. A hybrid simulated annealing and pattern search approach is used as an optimizer. The results show the objective function of optimizing the mean performance and minimizing the variation of performance parameters in terms of thrust ratio and total impulse could be achieved while adhering to the system constraints.


---

## Body

## 1. Introduction

Uncertainties are unavoidable in all real world scientific and engineering problems and products. The products are often encountered with variations in various stages of design, manufacturing, service, and/or degradation during storage or operation. In design phase, these variations present challenges to the design process and have a direct effect on the product quality, performance, and reliability. The uncertainties arise during design phases are either aleatory or epistemic. The aleatory or random uncertainties mostly account for initial operating conditions such as temperature, pressure, velocity, loading conditions, and manufacturing tolerances. The aleatory uncertainties can be represented by well-established and precisely known probability distribution techniques. On the other hand epistemic uncertainties account for the model error due to lack of design knowledge or information, numerical error, incomplete understanding of the system, or unpredictable behaviour during operation of some parameters having fix values [1–4]. If unaccounted during design process, these uncertainties can cause loss in quality, deterioration of performance, and product reliability. Therefore, uncertainty representation, aggregation, and propagation are the challenges in design optimization process [1, 2].The robust design optimization (RDO) strategy accounts for the effects of variation and uncertainties by simultaneously optimizing the objective function and minimizing performance parameter variations [5]. The RDO does not expunge the uncertainties rather it makes the product performance insensitive to these variations. The main difficulty in optimization under uncertainty is how to deal with an uncertainty space that is huge and frequently leads to very large-scale optimization models [4]. To guarantee solution quality, robustness assessment is often used and is considered as a key part of robust optimization [6]. In this paper, we computed variance using first-order Taylor approximation and percentile difference approach (PDM).The architecture of aerospace propulsion system consists of a single chamber solid propellant dual thrust solid rocket motor (DTSRM) having two levels of thrust, namely, boost phase thrust (Fb) and sustain phase thrust (Fs). Fb is required to accelerate the vehicle from zero velocity to a certain stabilized velocity and make the vehicle to reach a certain altitude with high Mach number in quite a short period of time. Then it sustains the vehicle at a constant velocity for a longer time with low level of thrust, Fs. The performance of such propulsion system is determined by internal ballistics through either motor grain burnback analysis with CAD modeling or employing analytical expressions [7, 8]. Literature on single thrust solid rocket motor design using traditional optimization methods without uncertainties is inundated. Sforzini used pattern search (PS) technique for automated design of solid rockets utilizing a special internal ballistics model [9]. Brooks carried out optimization by generating parametric design data for evaluating several characteristics of the star grain geometry [10]. Clegern adopted an expert system knowledge-based approach for computer-aided motor design [11]. Anderson used genetic algorithm (GA) for motor optimization [12]. Recently, Kamran and Guozhu proposed a hyperheuristic approach to minimize the motor mass [13]. Literature on dual grain geometry optimization and design however is scarce. Hu et al. presented a design study on high thrust ratio approach using single chamber dual-thrust solid rocket motor [14]. Dunn and Coats suggested use of dual propellant to achieve dual thrust [15].


## 2. Robust Design Optimization Methodology

### 2.1. Robust Design Optimization Methodology

The uncertainty based robust design optimization (UBRDO) is used in the preliminary design phase to find a robust solution that is insensitive to resulting changes and variations in subsequent design phases without expunging the source of uncertainty. This goal is achieved by simultaneously “optimizing the mean performance” and “reducing the performance variation”, subject to the robustness of constraints [16]. The key components of UBRDO methodology used in present research are as follows:(i)
objective robustness: maintaining robustness in the objective function,(ii)
feasibility robustness: maintaining robustness in the constraint,(iii)
robustness assessment: estimating variance of the performance parameters,(iv)
optimization approach: to perform the optimization,(v)
sensitivity analysis: parametric analysis of performance response.The framework of the proposed UBRDO methodology is presented in Figure1.Figure 1
Framework of proposed UBRDO.The main steps in executing the robust strategy are as follows:Step 1.
Define the design problem in terms of design variables (X).Step 2.
Initialization of population size by Latin hypercube sampling (LHS), function evaluation and stopping criterion.Step 3.
Decision for next generation subpopulation and its function evaluation.Step 4.
Hybrid optimization through simulated annealing and pattern search to obtain optimal solution.Step 5.
Execute DTSRM analysis code.Step 6.
Uncertainty propagation by perturbing the upper and lower bounds for the sampled design variables through worst-case scenario.Step 7.
Execute grain regression at each design point to calculate the effect of variation of objective function.Step 8.
Run optimization algorithm to obtain robust solution.Step 9.
Robustness assessment of objective function and constraint.Step 10.
Sensitivity analysis of performance parameters.


### 2.2. Robust Design Optimization Formulation

The generalized form of the optimization problem without considering robust strategy can be given by(2.1)min⁡xf(x)s.t.LB≤gi(x)≤UBlb≤x≤ub,
where f(x) is the objective function, x is the vector of design variables, and gi(x) is the ith constraint. When considering the robust optimization strategy under uncertainty the objective function “f(x)” in (2.1) is replaced by mean and standard deviation. It can be formulated as follows:
(2.2)min⁡df(μ,σ)=(wμf(d)+vσf(d))s.t.LB+(kσ(gi(d,z)))≤E((gi(d,z)))≤UB-(kσ(gi(d,z)))∀i,lbi+(kσ(xi))≤(di)≤ubi-(kσ(xi))fori=1,2,…,nddv,lbi≤di≤ubifori=1,2,…,nrdv,
where, μf is the mean value of the objective function; σf the standard deviation of the objective function; d the vector of deterministic design variables x; d the mean values of the uncertain design variables x; nrdv the number of the random design variables; nddv the numbers of the deterministic design variables; z the vector of nondesign input random vectors; w the weighting coefficients objectives μf; v the weighting coefficients objectives σf; gi(d,z) the ith constraint; E(gi(d,z)) the expectation of design mean; σ(gi(d,z)) the standard deviation of the ith constraint; LB and UB the vectors of lower and upper bounds of constraints gi’s; lb and ub the vectors of lower and upper bounds of the design variables; σ(xi) the vector of standard deviations of the random variables; k the adjusting constant.The focus of this study is to evolve insensitive optimized design in the presence of aleatory and epistemic uncertainties in the design parameters of dual thrust solid rocket motor. In the current work, the 3D grain configuration geometry was modeled as an inherent uncertain variable described with a normal probability distribution. The propellant burning rate was modeled as epistemic uncertain variable, since its uncertainty originates due to the lack of knowledge in a physical model, and they were represented as an interval with specified bounds. The space-filling designs strategies such as Latin hypercube sampling (LHS) are helpful when there is little or no information about the primary effects of data on responses, as in the case of epistemic uncertainties. The aim of this sampling technique is to spread the points as evenly as possible around the operating space. These designs literally fill out then-dimensional space with points that are in some way regularly spaced [17]. A dynamic penalty function is embedded to handle the violations in weighted sum of costs. A symbolic problem statement can be expressed as follows:
(2.3)min⁡f(x)=f(x)+h(kcit)∑i=1mmax⁡{0,gi(x)},
where f(x) is the objective function, h(k) is a dynamically modified penalty value, and kcit is the current iteration number of the algorithm. The function gi(x) is a relative violated function of the constraint [18].


### 2.3. MeasuringRobustness Assessment

In RDO,robustness assessmentis a measure of solution quality of performance parameters due to presence of uncertainty in design variables. In present formulation of objective function modeling, we have considered two robustness assessment measures:(i)
variance by mean and standard deviation,(ii)
percentile difference.The variance is estimated by a standard first-order Taylor series approximation due to its simplicity. Using this approximation, the variance of performance functionY=g(X) at the mean values, μX of X, is given by
(2.4)σγ2≅∑i=1n(∂g∂Xi|μX)2σXi2+2∑i<j∑∂g∂Xi|μX∂g∂Xj|μXσXiXj2,
where σXi2 is the variance of Xi and σXiXj2 is the covariance of Xi and Xj.IfX is mutually independent, the variance σγ2 of performance function Y is given by:(2.5)σγ2≅∑i=1n(∂g∂Xi|μX)2σXi2.
The upper bounds of standard deviation of uncertainty in Y are given by(2.6)σγ≅∑i=1n(∂g∂Xi|μX)σXi.
The estimation of mean neglecting the second-order sensitivities is given by(2.7)μγ=f(μx1,μx2,…,μxn).
The worst-case deviation Δx1w,Δx2w,…,Δxnw corresponds to the uncertainties in x1,x2,…,xn, respectively; then the worst-case estimation ΔYw in Y is given by the absolute sum of(2.8)ΔYw=∑i=1n|∂g∂Xi||Δxiw|.
The second method employs percentile difference approximation and is given as(2.9)Δyα1α2=yα2-yα1,
where yα1 and yα2 are the two values of performance parameters (Y) given as
(2.10)Prob{Y≤yαi}=αi(i=1,2).Hereα1 is the cumulative distribution function (CDF) of Y at the left tail and its value is taken as 0.05. α2 is the CDF of Y right tail and its value is taken as 0.95. yαi is the value of Y that corresponds to CDF αi and such a value is called a percentile value. The variance by percentile difference depicts broader picture relative to the standard deviation. In addition to variance, it also provides the skewness of distribution and the probability level at which design robustness could be achieved. The goal is to reduce the percentile difference for robustness. Normally, minimizing the mean of objective function can shift the location of the distribution towards left, while minimizing the percentile performance difference is conducive to shrinking the range of the distribution [19].


### 2.4. Uncertainty Propagation Modeling

We integrate a worst-case variation and statistical approach in our UBRDO approach to propagate the effect of uncertainties. A simulation-based method of First-Order Orthogonal Design Matrices (FOODMs) is proposed to propagate uncertainty using the mean and worst-case estimation of variance. FOODMs of +1’s and −1’s are used whose rows and columns are orthogonal. The last column is the actual variable settings, while the first column (all ones) enables us to measure the mean effect in the linear equation:(2.11)Y=Xβ+ε.
All control factors are assigned to FOODMs as shown in Table 1. The two levels called “high” and “low”, denoted by +1 and −1, are worst-case variations (uncertainties) in design variable that is, Xi+ΔXi and Xi-ΔXi, respectively, (i=1,2,…,n). The high level “+1” corresponds to upper limit of perturbed variable and correspondingly low level “−1” corresponds to lower limit of perturbed variable. The response yi is simulation output corresponding to each test run with perturbed design variables. The output mean is calculated with nominal values of design variables at current design point. The “worst-case variation” of output is estimated using simulation runs as
(2.12)ΔYw=max⁡(|yi-y̅|),
where yi is the current perturbed output and y̅ is the nominal value of the design point.Table 1
Layout of FOODM.
Control factors (design variables)TestX1X2X3X4X5X6X7X8Xny1+1+1+1+1+1+1+1+1·y12+1−1+1−1+1−1+1−1·y23+1+1−1−1+1+1−1−1·y34+1−1−1+1+1−1−1+1··5+1+1+1+1−1−1−1−1··6+1−1+1−1−1+1−1+1··7+1+1−1−1−1−1+1+1··8+1−1−1+1−1+1+1−1·ynn·········yn


### 2.5. Hybrid Optimization Methodology

Simulated Annealing (SA) is a metaheuristics algorithm, analogous to the annealing processes of metals to generate succeeding solutions by means of local search procedure. Following a predefined selection criterion, some of them are accepted while others will be rejected. Metropolis et al. applied the same idea to simulate atoms in equilibrium. The “Metropolis Algorithm” has also been applied to solving combinatorial optimization problems [20]. SA algorithm uses a probabilistically determined sequence on each iteration to decide whether a new point is selected or not [21–24].Pattern search (PS) is a direct search-derivative free optimization algorithm that does not require gradient of the problem to be optimized. Therefore PS can be used on functions that are not continuous or differentiable. Hooke and Jeeves were the first to propose PS as a direct line search optimization method [25].A pattern search algorithm computes a sequence of points that get closer to the optimal point. At each step, the algorithm searches a set of points, called a mesh, around the current point, that is, the point computed at the previous step of the algorithm. The algorithm forms the mesh by adding the current point to a scalar multiple of a fixed set of vectors called a pattern. At each step, the algorithm polls the points in the current mesh by computing their objective function values. If the algorithm finds a point in the mesh that improves the objective function at the current point, the poll is then called successful and that point becomes the current point at the next step; otherwise the iteration continues.The hybrid approach, SAPS, is based on exploring and exploiting the local search to multiple individuals in population in between generations of an SA. The hybridization of SA and direct search algorithms has shown successful in other difficult optimization problems [26–28]. In SAPS, a two-stage hybrid method SAPS (Simulated Annealing with Pattern Search) is investigated. The first stage is explorative, employing a traditional SA to identify promising areas of the search space. The best solution found by the simulated annealing (SA) is then refined using a pattern search (PS) method during the second exploitative stage. The aforementioned optimization methods are incorporated to find the optimal solution that cannot be considered as robust one [29]. Therefore here we refer to optimal solution (XSAPS*) as nonrobust one. The hybrid SAPS optimization framework is shown in Figure 2.Figure 2
Framework of hybrid optimization SAPS.


## 2.1. Robust Design Optimization Methodology

The uncertainty based robust design optimization (UBRDO) is used in the preliminary design phase to find a robust solution that is insensitive to resulting changes and variations in subsequent design phases without expunging the source of uncertainty. This goal is achieved by simultaneously “optimizing the mean performance” and “reducing the performance variation”, subject to the robustness of constraints [16]. The key components of UBRDO methodology used in present research are as follows:(i)
objective robustness: maintaining robustness in the objective function,(ii)
feasibility robustness: maintaining robustness in the constraint,(iii)
robustness assessment: estimating variance of the performance parameters,(iv)
optimization approach: to perform the optimization,(v)
sensitivity analysis: parametric analysis of performance response.The framework of the proposed UBRDO methodology is presented in Figure1.Figure 1
Framework of proposed UBRDO.The main steps in executing the robust strategy are as follows:Step 1.
Define the design problem in terms of design variables (X).Step 2.
Initialization of population size by Latin hypercube sampling (LHS), function evaluation and stopping criterion.Step 3.
Decision for next generation subpopulation and its function evaluation.Step 4.
Hybrid optimization through simulated annealing and pattern search to obtain optimal solution.Step 5.
Execute DTSRM analysis code.Step 6.
Uncertainty propagation by perturbing the upper and lower bounds for the sampled design variables through worst-case scenario.Step 7.
Execute grain regression at each design point to calculate the effect of variation of objective function.Step 8.
Run optimization algorithm to obtain robust solution.Step 9.
Robustness assessment of objective function and constraint.Step 10.
Sensitivity analysis of performance parameters.


## 2.2. Robust Design Optimization Formulation

The generalized form of the optimization problem without considering robust strategy can be given by(2.1)min⁡xf(x)s.t.LB≤gi(x)≤UBlb≤x≤ub,
where f(x) is the objective function, x is the vector of design variables, and gi(x) is the ith constraint. When considering the robust optimization strategy under uncertainty the objective function “f(x)” in (2.1) is replaced by mean and standard deviation. It can be formulated as follows:
(2.2)min⁡df(μ,σ)=(wμf(d)+vσf(d))s.t.LB+(kσ(gi(d,z)))≤E((gi(d,z)))≤UB-(kσ(gi(d,z)))∀i,lbi+(kσ(xi))≤(di)≤ubi-(kσ(xi))fori=1,2,…,nddv,lbi≤di≤ubifori=1,2,…,nrdv,
where, μf is the mean value of the objective function; σf the standard deviation of the objective function; d the vector of deterministic design variables x; d the mean values of the uncertain design variables x; nrdv the number of the random design variables; nddv the numbers of the deterministic design variables; z the vector of nondesign input random vectors; w the weighting coefficients objectives μf; v the weighting coefficients objectives σf; gi(d,z) the ith constraint; E(gi(d,z)) the expectation of design mean; σ(gi(d,z)) the standard deviation of the ith constraint; LB and UB the vectors of lower and upper bounds of constraints gi’s; lb and ub the vectors of lower and upper bounds of the design variables; σ(xi) the vector of standard deviations of the random variables; k the adjusting constant.The focus of this study is to evolve insensitive optimized design in the presence of aleatory and epistemic uncertainties in the design parameters of dual thrust solid rocket motor. In the current work, the 3D grain configuration geometry was modeled as an inherent uncertain variable described with a normal probability distribution. The propellant burning rate was modeled as epistemic uncertain variable, since its uncertainty originates due to the lack of knowledge in a physical model, and they were represented as an interval with specified bounds. The space-filling designs strategies such as Latin hypercube sampling (LHS) are helpful when there is little or no information about the primary effects of data on responses, as in the case of epistemic uncertainties. The aim of this sampling technique is to spread the points as evenly as possible around the operating space. These designs literally fill out then-dimensional space with points that are in some way regularly spaced [17]. A dynamic penalty function is embedded to handle the violations in weighted sum of costs. A symbolic problem statement can be expressed as follows:
(2.3)min⁡f(x)=f(x)+h(kcit)∑i=1mmax⁡{0,gi(x)},
where f(x) is the objective function, h(k) is a dynamically modified penalty value, and kcit is the current iteration number of the algorithm. The function gi(x) is a relative violated function of the constraint [18].


## 2.3. MeasuringRobustness Assessment

In RDO,robustness assessmentis a measure of solution quality of performance parameters due to presence of uncertainty in design variables. In present formulation of objective function modeling, we have considered two robustness assessment measures:(i)
variance by mean and standard deviation,(ii)
percentile difference.The variance is estimated by a standard first-order Taylor series approximation due to its simplicity. Using this approximation, the variance of performance functionY=g(X) at the mean values, μX of X, is given by
(2.4)σγ2≅∑i=1n(∂g∂Xi|μX)2σXi2+2∑i<j∑∂g∂Xi|μX∂g∂Xj|μXσXiXj2,
where σXi2 is the variance of Xi and σXiXj2 is the covariance of Xi and Xj.IfX is mutually independent, the variance σγ2 of performance function Y is given by:(2.5)σγ2≅∑i=1n(∂g∂Xi|μX)2σXi2.
The upper bounds of standard deviation of uncertainty in Y are given by(2.6)σγ≅∑i=1n(∂g∂Xi|μX)σXi.
The estimation of mean neglecting the second-order sensitivities is given by(2.7)μγ=f(μx1,μx2,…,μxn).
The worst-case deviation Δx1w,Δx2w,…,Δxnw corresponds to the uncertainties in x1,x2,…,xn, respectively; then the worst-case estimation ΔYw in Y is given by the absolute sum of(2.8)ΔYw=∑i=1n|∂g∂Xi||Δxiw|.
The second method employs percentile difference approximation and is given as(2.9)Δyα1α2=yα2-yα1,
where yα1 and yα2 are the two values of performance parameters (Y) given as
(2.10)Prob{Y≤yαi}=αi(i=1,2).Hereα1 is the cumulative distribution function (CDF) of Y at the left tail and its value is taken as 0.05. α2 is the CDF of Y right tail and its value is taken as 0.95. yαi is the value of Y that corresponds to CDF αi and such a value is called a percentile value. The variance by percentile difference depicts broader picture relative to the standard deviation. In addition to variance, it also provides the skewness of distribution and the probability level at which design robustness could be achieved. The goal is to reduce the percentile difference for robustness. Normally, minimizing the mean of objective function can shift the location of the distribution towards left, while minimizing the percentile performance difference is conducive to shrinking the range of the distribution [19].


## 2.4. Uncertainty Propagation Modeling

We integrate a worst-case variation and statistical approach in our UBRDO approach to propagate the effect of uncertainties. A simulation-based method of First-Order Orthogonal Design Matrices (FOODMs) is proposed to propagate uncertainty using the mean and worst-case estimation of variance. FOODMs of +1’s and −1’s are used whose rows and columns are orthogonal. The last column is the actual variable settings, while the first column (all ones) enables us to measure the mean effect in the linear equation:(2.11)Y=Xβ+ε.
All control factors are assigned to FOODMs as shown in Table 1. The two levels called “high” and “low”, denoted by +1 and −1, are worst-case variations (uncertainties) in design variable that is, Xi+ΔXi and Xi-ΔXi, respectively, (i=1,2,…,n). The high level “+1” corresponds to upper limit of perturbed variable and correspondingly low level “−1” corresponds to lower limit of perturbed variable. The response yi is simulation output corresponding to each test run with perturbed design variables. The output mean is calculated with nominal values of design variables at current design point. The “worst-case variation” of output is estimated using simulation runs as
(2.12)ΔYw=max⁡(|yi-y̅|),
where yi is the current perturbed output and y̅ is the nominal value of the design point.Table 1
Layout of FOODM.
Control factors (design variables)TestX1X2X3X4X5X6X7X8Xny1+1+1+1+1+1+1+1+1·y12+1−1+1−1+1−1+1−1·y23+1+1−1−1+1+1−1−1·y34+1−1−1+1+1−1−1+1··5+1+1+1+1−1−1−1−1··6+1−1+1−1−1+1−1+1··7+1+1−1−1−1−1+1+1··8+1−1−1+1−1+1+1−1·ynn·········yn


## 2.5. Hybrid Optimization Methodology

Simulated Annealing (SA) is a metaheuristics algorithm, analogous to the annealing processes of metals to generate succeeding solutions by means of local search procedure. Following a predefined selection criterion, some of them are accepted while others will be rejected. Metropolis et al. applied the same idea to simulate atoms in equilibrium. The “Metropolis Algorithm” has also been applied to solving combinatorial optimization problems [20]. SA algorithm uses a probabilistically determined sequence on each iteration to decide whether a new point is selected or not [21–24].Pattern search (PS) is a direct search-derivative free optimization algorithm that does not require gradient of the problem to be optimized. Therefore PS can be used on functions that are not continuous or differentiable. Hooke and Jeeves were the first to propose PS as a direct line search optimization method [25].A pattern search algorithm computes a sequence of points that get closer to the optimal point. At each step, the algorithm searches a set of points, called a mesh, around the current point, that is, the point computed at the previous step of the algorithm. The algorithm forms the mesh by adding the current point to a scalar multiple of a fixed set of vectors called a pattern. At each step, the algorithm polls the points in the current mesh by computing their objective function values. If the algorithm finds a point in the mesh that improves the objective function at the current point, the poll is then called successful and that point becomes the current point at the next step; otherwise the iteration continues.The hybrid approach, SAPS, is based on exploring and exploiting the local search to multiple individuals in population in between generations of an SA. The hybridization of SA and direct search algorithms has shown successful in other difficult optimization problems [26–28]. In SAPS, a two-stage hybrid method SAPS (Simulated Annealing with Pattern Search) is investigated. The first stage is explorative, employing a traditional SA to identify promising areas of the search space. The best solution found by the simulated annealing (SA) is then refined using a pattern search (PS) method during the second exploitative stage. The aforementioned optimization methods are incorporated to find the optimal solution that cannot be considered as robust one [29]. Therefore here we refer to optimal solution (XSAPS*) as nonrobust one. The hybrid SAPS optimization framework is shown in Figure 2.Figure 2
Framework of hybrid optimization SAPS.


## 3. Design Methodology

### 3.1. Grain Geometry

In this paper, single chamber dual grain geometry is considered for motor test case. The 3D finocyl grain consists of fin geometry in the aft followed by regular hollow tubular cross-section that runs through the entire length till the forward end of the motor. The two grain geometries are smoothly interconnected by a transition taper concave fin section (see Figure3). Figure 4 shows the geometrical parameters of fin geometry used in this analysis. The ballistic design analysis of this 3D grain configuration requires parametric evaluation of burning surface areas including geometrical variables, such as hollow tubular grain inner diameter (dtu), fin radius (f), inner radius (ri), minimum web (w), maximum web (W), fin height (l), number of fins (N), and transition section taper angle (ϕ).Figure 3
Schematic of dual thrust solid rocket motor.Figure 4
Finocyl grain geometry parameters.The burning surface area of the 3D grain geometry is calculated by grain regression analysis through computer program and is given by(3.1)Abk=Vk+1-Vkwk+1-wk,
where k is the web step, V the volume of propellant, and w the web at respective position of geometry and propellant mass is calculated simply by
(3.2)mp=Vkρp.


### 3.2. Performance Prediction

Internal ballistic calculation is the core of any propulsion system to evaluate its performance. The performance of propulsion system utilizing solid rocket motor is calculated using a zero-dimension ballistic model. The dual levels of thrust are achieved by change in burning areas of grain geometry during burning. The chamber pressurePc is calculated by equating mass generated in the chamber to the mass ejected through nozzle throat using (3.3) [30]: (3.3)Pc=(ρpac*AbAt)1/(1-n),
where At denotes the throat area, Ab the burning surface area, a the propellant burning rate coefficient, n the pressure sensitivity index, c* the characteristic velocity, and ρp the propellant density. Thrust is determined by using (3.4) as follows: (3.4)F=PcAtCF,F=AtCF(ρpac*AbAt)1/(1-n).
Since At is assumed to be constant during burning, therefore dual thrust is obtained by change in Ab of the propellant grain. The thrust coefficient (CF) is given by (3.5): (3.5)CF=2γ2γ-1(2γ+1)(γ+1)/(γ-1)[1-(PePc)(γ-1)/γ]+Pe-PambPcε,
where γ is ratio of the specific heat capacity, Pe is nozzle exit pressure, Pamb is ambient pressure and ε is nozzle area expansion ratio. The total impulse is given by
(3.6)IT=CFc*mp.


## 3.1. Grain Geometry

In this paper, single chamber dual grain geometry is considered for motor test case. The 3D finocyl grain consists of fin geometry in the aft followed by regular hollow tubular cross-section that runs through the entire length till the forward end of the motor. The two grain geometries are smoothly interconnected by a transition taper concave fin section (see Figure3). Figure 4 shows the geometrical parameters of fin geometry used in this analysis. The ballistic design analysis of this 3D grain configuration requires parametric evaluation of burning surface areas including geometrical variables, such as hollow tubular grain inner diameter (dtu), fin radius (f), inner radius (ri), minimum web (w), maximum web (W), fin height (l), number of fins (N), and transition section taper angle (ϕ).Figure 3
Schematic of dual thrust solid rocket motor.Figure 4
Finocyl grain geometry parameters.The burning surface area of the 3D grain geometry is calculated by grain regression analysis through computer program and is given by(3.1)Abk=Vk+1-Vkwk+1-wk,
where k is the web step, V the volume of propellant, and w the web at respective position of geometry and propellant mass is calculated simply by
(3.2)mp=Vkρp.


## 3.2. Performance Prediction

Internal ballistic calculation is the core of any propulsion system to evaluate its performance. The performance of propulsion system utilizing solid rocket motor is calculated using a zero-dimension ballistic model. The dual levels of thrust are achieved by change in burning areas of grain geometry during burning. The chamber pressurePc is calculated by equating mass generated in the chamber to the mass ejected through nozzle throat using (3.3) [30]: (3.3)Pc=(ρpac*AbAt)1/(1-n),
where At denotes the throat area, Ab the burning surface area, a the propellant burning rate coefficient, n the pressure sensitivity index, c* the characteristic velocity, and ρp the propellant density. Thrust is determined by using (3.4) as follows: (3.4)F=PcAtCF,F=AtCF(ρpac*AbAt)1/(1-n).
Since At is assumed to be constant during burning, therefore dual thrust is obtained by change in Ab of the propellant grain. The thrust coefficient (CF) is given by (3.5): (3.5)CF=2γ2γ-1(2γ+1)(γ+1)/(γ-1)[1-(PePc)(γ-1)/γ]+Pe-PambPcε,
where γ is ratio of the specific heat capacity, Pe is nozzle exit pressure, Pamb is ambient pressure and ε is nozzle area expansion ratio. The total impulse is given by
(3.6)IT=CFc*mp.


## 4. Dual Thrust Motor Test Case

### 4.1. Design Objectives and Constraints

There can be different objective functions depending upon the mission requirements. In case of solid rocket motor, the designers have always probed for high total impulse, minimum motor mass, and high reliability, and so forth. In our case, the objective function is to minimize the mean performance of average boost phase to sustain phase thrust ratio (FRave) and total impulse (It) of the aerospace vehicle propulsion system and their variations.The context of objective function formulation under uncertainty, the variations in grain geometric parameters (design variables), and propellant burning rate (nondesign parameter) are considered as aleatory and epistemic uncertainties respectively. The variations in design variables are realized as lower and upper bounds whereas nondesign parameter has a fix value with an attributed tolerance by expert elicitation. The design and nondesign variablesX and Z are given by
(4.1)X=f(dtu,l,f,ri,w,W,η,N,ϕ),Z=f(rg).
The formulation of objective function of (4.1) according to (2.2) will be given as (4.2)d*=argmind⁡(wE(FRave)+(1-w)σ(FRave)),d*=argmind⁡(wE(It)+(1-w)σ(It))fori=1,2,…,8,lb+kσ(x)≤(d)≤ub-kσ(x),μz*=argminμz⁡(wE(FRave)+(1-w)σ(FRave)),μz*=argminμz⁡(wE(It)+(1-w)σ(It))fori=1,Zl≤μz≤Zu,
subject to the following constraints:
(4.3)Cj(X)≥0(j=1,2,…,7),Cj(X)≤0(j=8,…,10),
where C is given by (4.4) as follows:
(4.4)C1:FRave≥5.5,C2:It≥930kN-sec,C3:Pb_min⁡≥14MPa,C4:Ps_min⁡≥2MPa,C5:tb≥10.5sec,C6:mp≥400kg,C7:rg_min⁡≥4mm/sec,C8:Pb_max⁡≤16MPa,C9:Ps_max⁡≤3MPa,C10:rg_max⁡≤8.5mm/sec,
and bound to
(4.5){Lowerbound=min⁡(Xi)Upperbound=max⁡(Xi)}(i=1,2,…,8).The upper bound (UB) and lower bound (LB) of design variables are listed in Table2. The propellant properties and nozzle parameters for the ballistics analysis are listed in Table 3.Table 2
LB and UB of design variables.
ParametersLBUBUnitsTubular grain inner diameter (dtu)105130mmFin radius (f)618mmInner radius (ri)5070mmMin. web (w)1535mmMax. web (W)75100mmFin height (l)90120mmNumber of fins (N)412—Taper fin angle (ϕ)1825degreeTable 3
Propellant and nozzle parameters.
ParametersValueUnitsPropellant density(ρp)1745kg/m3Grain burning rate(rg)8.5@7 MPa ± 1mm/secCharacteristic velocity(c*)1560m/secNozzle throat diameter(dt)90mmNozzle area expansion ratio (ε)9—


## 4.1. Design Objectives and Constraints

There can be different objective functions depending upon the mission requirements. In case of solid rocket motor, the designers have always probed for high total impulse, minimum motor mass, and high reliability, and so forth. In our case, the objective function is to minimize the mean performance of average boost phase to sustain phase thrust ratio (FRave) and total impulse (It) of the aerospace vehicle propulsion system and their variations.The context of objective function formulation under uncertainty, the variations in grain geometric parameters (design variables), and propellant burning rate (nondesign parameter) are considered as aleatory and epistemic uncertainties respectively. The variations in design variables are realized as lower and upper bounds whereas nondesign parameter has a fix value with an attributed tolerance by expert elicitation. The design and nondesign variablesX and Z are given by
(4.1)X=f(dtu,l,f,ri,w,W,η,N,ϕ),Z=f(rg).
The formulation of objective function of (4.1) according to (2.2) will be given as (4.2)d*=argmind⁡(wE(FRave)+(1-w)σ(FRave)),d*=argmind⁡(wE(It)+(1-w)σ(It))fori=1,2,…,8,lb+kσ(x)≤(d)≤ub-kσ(x),μz*=argminμz⁡(wE(FRave)+(1-w)σ(FRave)),μz*=argminμz⁡(wE(It)+(1-w)σ(It))fori=1,Zl≤μz≤Zu,
subject to the following constraints:
(4.3)Cj(X)≥0(j=1,2,…,7),Cj(X)≤0(j=8,…,10),
where C is given by (4.4) as follows:
(4.4)C1:FRave≥5.5,C2:It≥930kN-sec,C3:Pb_min⁡≥14MPa,C4:Ps_min⁡≥2MPa,C5:tb≥10.5sec,C6:mp≥400kg,C7:rg_min⁡≥4mm/sec,C8:Pb_max⁡≤16MPa,C9:Ps_max⁡≤3MPa,C10:rg_max⁡≤8.5mm/sec,
and bound to
(4.5){Lowerbound=min⁡(Xi)Upperbound=max⁡(Xi)}(i=1,2,…,8).The upper bound (UB) and lower bound (LB) of design variables are listed in Table2. The propellant properties and nozzle parameters for the ballistics analysis are listed in Table 3.Table 2
LB and UB of design variables.
ParametersLBUBUnitsTubular grain inner diameter (dtu)105130mmFin radius (f)618mmInner radius (ri)5070mmMin. web (w)1535mmMax. web (W)75100mmFin height (l)90120mmNumber of fins (N)412—Taper fin angle (ϕ)1825degreeTable 3
Propellant and nozzle parameters.
ParametersValueUnitsPropellant density(ρp)1745kg/m3Grain burning rate(rg)8.5@7 MPa ± 1mm/secCharacteristic velocity(c*)1560m/secNozzle throat diameter(dt)90mmNozzle area expansion ratio (ε)9—


## 5. Results and Discussion

### 5.1. Comparison of Results

The solid rocket motors are designed to provide required pressure-time (P-t) and thrust-time (F-t) profiles. The P-t and F-t profiles for DTSRM test case are shown in Figure 5 for the robust results. The peak pressure during boost phase is within the maximum limit imposed by the constraint. Tables 4 and 5 show the comparison of robust and nonrobust results achieved by UBRDO and SAPS hybrid optimization approaches, respectively, for grain design variables and required motor performance parameters. Table 6 shows the comparison of robustness assessment for robust and nonrobust optimization approach. The two measures of variance modeled in Section 2 reveal that the robustness is achieved with minimum dispersion and adhering to the targeted mean by the robust approach. The variance of output mean calculated at current design points with robust solution exhibited very less variance compared to nonrobust approach of hybrid SAPS. The performance parameters achieved by the UBRDO strategy are well within the range. All these values have been achieved by adhering to and obeying the bounds of design variables, propellant properties and nozzle parameters given in Tables 2 and 3, and constraints imposed by (4.4). The nonrobust result however shows violation for the maximum boost pressure constraint limit.Table 4
Comparison of design variables.
Design variablesRobustOptimalUnitsTubular grain inner diameter (dtu)116120mmFin radius (f)911mmInner radius (ri)5760mmMin. web (w)2325mmMax. web (W)8186mmFin height (l)105110mmNumber of fins (N)77—Taper fin angle (ϕ)1922degTable 5
Comparison of performance parameters.
ParametersRobustOptimalFRave5.806.11It (kN-sec)936.75965.55mp (kg)404.10420.10Pb_max (MPa)14.5515.94Table 6
Robustness assessment by variance.
MethodFRavePb_max (MPa)It (kN-sec)mp (kg)Percentile differenceoptimalrobustoptimalrobustoptimalrobustoptimalrobusty0.055.935.7815.3514.41947.17938.22405.17404.35y0.956.515.9016.4814.81978.34948.69433.32407.33Δy0.950.050.580.121.130.4031.1710.4728.152.98First-order Taylor seriesσ0.160.050.410.1112.753.5513.411.45μ (targeted)5.514.50930400μ (achieved)6.115.8015.9414.55965.55936.75420.10404.10Figure 5
F-t and P-t profiles of DTSRM.Beside mean-variance framework, robustness is also measured using percentile difference plots of performance parameters. It provides the extent to which the probability level of the design robustness is achieved. The comparison of data in Table6 confirms the robustness in all the performance parameters evaluated through robust approach. The results of percentile difference approach are plotted in Figures 6(a)–6(d). On examining the figure, shrinking and swelling of performance parameter data under uncertainty for robust and nonrobust optimization, respectively, can be seen clearly. The goal of shrinking the data dispersion and shifting the data distribution towards left has been achieved in all performance parameters through UBRDO.Percentile difference plots of performance parameters.
(a)
Average thrust ratio(b)
Total impulse(c)
Maximum boost pressure(d)
Propellant massFurthermore, to investigate the effects of uncertainties of design variables on performance parameters deviation, Monte Carlo simulation is performed by selecting a sample of 500 runs on random basis for a worst-case deviation ofΔ=±8%. The scatter plots are shown in Figures 7(a)–7(d) for FRave, It, Pb_max, and mp. The scatter plots reflect the anarchy of nonrobust solution afflicted due to the presence of uncertainty. The data points of robust solution have smoothly settled around the mean value confirming the robustness and insensitivity of the solution.Scatter plots of performance parameters.
(a)
Average thrust ratio(b)
Total impulse(c)
Maximum boost pressure(d)
Propellant mass


### 5.2. Parametric Sensitivity Analysis

Performance parameters sensitivity analysis is used to compute the response of variation of performance parameters with respect to variation of design variable. The design variables(dtu,l,f,ri,w,W,η,N,ϕ) related to grain geometries and nondesign parameters of propellant burning rate of star and cylindrical grains (rg_max,rg_min) are considered to analyze the sensitivity of performance parameters FRave, It, mp, Pb_max, and Pb_ave. The sensitivity of these design parameters on performance has been analyzed using the values given in Tables 2 and 3 and the constraints impose on performance parameters given by (4.4). The behavior of sensitivity is depicted in Figures 8(a)–8(d).Sensitivity plots of performance parameters.
(a)
Number of fins (N)(b)
Major web (W)(c)
Inner cylindrical grain diameter (dtu)(d)
Fin radius (f)By examining the sensitivity plots, theFRave decreases towards upper bounds of N, W, and f sharply and slowly with dtu. The It and mp decrease sharply with increase in value of N, f, and dtu and increase with W. The Pb_max increases sharply with increasing W and dtu, decreases gradually with increases in N and remains unaffected by f.From performance view point, special focus should be paid on the selection ofN and W as both have influence on FRave and It. The Pb_max should be taken care of as it violates constraint limit towards upper limits of W and dtu.Sensitivity of Propellant Burning Rate
The propellant burning rate parameter is considered as nondesign parameter. The affect of burning rate data uncertainty is analyzed by setting a tolerance value of ±1 which is considered as worst-case deviation. The sensitivity affects burning rate on pressure system of DTSRM are shown in Figure9 in detail. Towards higher values, it violates the constraint limit of Pb_max, whereas at lower limits, it violates the constraint limit of Ps_min. As higher pressure leads towards motor bursting and at lower pressure combustion extinguishing phenomenon can occur, both affects are undesirable for stable working of motor.Figure 9
Sensitivity analysis of burning rate.


## 5.1. Comparison of Results

The solid rocket motors are designed to provide required pressure-time (P-t) and thrust-time (F-t) profiles. The P-t and F-t profiles for DTSRM test case are shown in Figure 5 for the robust results. The peak pressure during boost phase is within the maximum limit imposed by the constraint. Tables 4 and 5 show the comparison of robust and nonrobust results achieved by UBRDO and SAPS hybrid optimization approaches, respectively, for grain design variables and required motor performance parameters. Table 6 shows the comparison of robustness assessment for robust and nonrobust optimization approach. The two measures of variance modeled in Section 2 reveal that the robustness is achieved with minimum dispersion and adhering to the targeted mean by the robust approach. The variance of output mean calculated at current design points with robust solution exhibited very less variance compared to nonrobust approach of hybrid SAPS. The performance parameters achieved by the UBRDO strategy are well within the range. All these values have been achieved by adhering to and obeying the bounds of design variables, propellant properties and nozzle parameters given in Tables 2 and 3, and constraints imposed by (4.4). The nonrobust result however shows violation for the maximum boost pressure constraint limit.Table 4
Comparison of design variables.
Design variablesRobustOptimalUnitsTubular grain inner diameter (dtu)116120mmFin radius (f)911mmInner radius (ri)5760mmMin. web (w)2325mmMax. web (W)8186mmFin height (l)105110mmNumber of fins (N)77—Taper fin angle (ϕ)1922degTable 5
Comparison of performance parameters.
ParametersRobustOptimalFRave5.806.11It (kN-sec)936.75965.55mp (kg)404.10420.10Pb_max (MPa)14.5515.94Table 6
Robustness assessment by variance.
MethodFRavePb_max (MPa)It (kN-sec)mp (kg)Percentile differenceoptimalrobustoptimalrobustoptimalrobustoptimalrobusty0.055.935.7815.3514.41947.17938.22405.17404.35y0.956.515.9016.4814.81978.34948.69433.32407.33Δy0.950.050.580.121.130.4031.1710.4728.152.98First-order Taylor seriesσ0.160.050.410.1112.753.5513.411.45μ (targeted)5.514.50930400μ (achieved)6.115.8015.9414.55965.55936.75420.10404.10Figure 5
F-t and P-t profiles of DTSRM.Beside mean-variance framework, robustness is also measured using percentile difference plots of performance parameters. It provides the extent to which the probability level of the design robustness is achieved. The comparison of data in Table6 confirms the robustness in all the performance parameters evaluated through robust approach. The results of percentile difference approach are plotted in Figures 6(a)–6(d). On examining the figure, shrinking and swelling of performance parameter data under uncertainty for robust and nonrobust optimization, respectively, can be seen clearly. The goal of shrinking the data dispersion and shifting the data distribution towards left has been achieved in all performance parameters through UBRDO.Percentile difference plots of performance parameters.
(a)
Average thrust ratio(b)
Total impulse(c)
Maximum boost pressure(d)
Propellant massFurthermore, to investigate the effects of uncertainties of design variables on performance parameters deviation, Monte Carlo simulation is performed by selecting a sample of 500 runs on random basis for a worst-case deviation ofΔ=±8%. The scatter plots are shown in Figures 7(a)–7(d) for FRave, It, Pb_max, and mp. The scatter plots reflect the anarchy of nonrobust solution afflicted due to the presence of uncertainty. The data points of robust solution have smoothly settled around the mean value confirming the robustness and insensitivity of the solution.Scatter plots of performance parameters.
(a)
Average thrust ratio(b)
Total impulse(c)
Maximum boost pressure(d)
Propellant mass


## 5.2. Parametric Sensitivity Analysis

Performance parameters sensitivity analysis is used to compute the response of variation of performance parameters with respect to variation of design variable. The design variables(dtu,l,f,ri,w,W,η,N,ϕ) related to grain geometries and nondesign parameters of propellant burning rate of star and cylindrical grains (rg_max,rg_min) are considered to analyze the sensitivity of performance parameters FRave, It, mp, Pb_max, and Pb_ave. The sensitivity of these design parameters on performance has been analyzed using the values given in Tables 2 and 3 and the constraints impose on performance parameters given by (4.4). The behavior of sensitivity is depicted in Figures 8(a)–8(d).Sensitivity plots of performance parameters.
(a)
Number of fins (N)(b)
Major web (W)(c)
Inner cylindrical grain diameter (dtu)(d)
Fin radius (f)By examining the sensitivity plots, theFRave decreases towards upper bounds of N, W, and f sharply and slowly with dtu. The It and mp decrease sharply with increase in value of N, f, and dtu and increase with W. The Pb_max increases sharply with increasing W and dtu, decreases gradually with increases in N and remains unaffected by f.From performance view point, special focus should be paid on the selection ofN and W as both have influence on FRave and It. The Pb_max should be taken care of as it violates constraint limit towards upper limits of W and dtu.Sensitivity of Propellant Burning Rate
The propellant burning rate parameter is considered as nondesign parameter. The affect of burning rate data uncertainty is analyzed by setting a tolerance value of ±1 which is considered as worst-case deviation. The sensitivity affects burning rate on pressure system of DTSRM are shown in Figure9 in detail. Towards higher values, it violates the constraint limit of Pb_max, whereas at lower limits, it violates the constraint limit of Ps_min. As higher pressure leads towards motor bursting and at lower pressure combustion extinguishing phenomenon can occur, both affects are undesirable for stable working of motor.Figure 9
Sensitivity analysis of burning rate.


## 6. Conclusion

In this paper, we explored the robust design methodology based on an integrated approach of complex design scenario and optimization. The integration of 3D grain design, internal ballistics, hybrid optimization, and worst-case uncertainty formulation supplemented by efficient statistical methods has shown promising results in optimizing the grain geometry design variables and motor performance in the presence of both aleatory and epistemic uncertainties. The robustness assessment using efficient variance approach and simulation-based uncertainty modeling through worst-case estimation provides an insensitive robust design solution. The sensitivity analysis helps us in identifying those design variables that contributed significantly in ameliorating or deteriorating the performance parameters. The hybrid approach of Genetic Algorithm and Simulated Annealing has shown excellent performance about the efficacy as well as optimization view point. As expected, the values of performance parameters achieved by robust design are less than the ones achieved by optimal one but insensitive to variations. The important achievement that can be associated with proposed methodology is its ability to evaluate and optimize as well the dual grain design subject to performance constraints under a complex scenario of both aleatory and epistemic uncertainties. The proposed framework increased reliability and robustness of our design and is useful for propulsion systems that require both optimality and robustness and in the meanwhile it allows designers to make reliable decisions when there are uncertainties associated with design parameters.



---

*Source: 145692-2011-12-20.xml*