Kick Solid Rocket Motor Multidisciplinary Design Optimization Using Genetic Algorithm

In this paper, a multidisciplinary design optimization (MDO) approach of a solid propellant kick rocket motor is considered. A genetic algorithm optimization method has been used. The optimized kick solid rocket motor (KSRM) is capable of delivering a small satellite of 200 kg to a circular low earth orbit (LEO) of 600 km altitude. The KSRM should accelerate from the initial apogee velocity of 5000 m/s up to the orbital insertion velocity of 7560 m/s. The KSRM design variables and the orbital insertion trajectory profile variables were optimized simultaneously, whereas the mass characteristics of the payload deployment module were assigned. A depleted shutdown condition was considered, to avoid the necessity of a thrust termination device, resulting in a reduced total mass of the KSRM. The results show that the proposed optimization approach was able to find the convergence of the optimal solution with highly acceptable value for conceptual design phase. KeywoRdS: Kick solid rocket motor, Multidisciplinary design optimization, Genetic algorithm. INTRODUCTION Using small solid rocket motors (SRM) for space applications became very attractive, because of their advantages compared with liquid propellant driven rocket engines, especially for insertion of small payloads into a circular Low Earth Orbit (LEO). Simplicity, reliability, and easy to fabricate and operate, are some of the main characteristics of the SRMs. In this article, a specific design used for insertion of a payload into a specified orbit commonly known as kick solid rocket motors (KSRM) are analyzed and discussed. Previous works were performed on analysis of several options of upper stages for different launch vehicles (LV) (McGinnis and Joyner, 2005), and in design and optimization of upper stages for transfer from LEO to geostationary earth orbit (Motlagh and Novinzadeh, 2012). Casalino et al., (2011) carried out an optimization of a hybrid upper stage configuration with a special emphasis on the grain geometry. He Linshu and Murad (2005) developed an improved method for conceptual design of multistage solid rockets based on depleted shutdown condition, which means that the SRM has to burn all-contained propellant, resulting in a reduced gross mass. The genetic algorithm (GA) global optimization method is increasingly being used in optimization of aerospace and propulsion systems (Kamran et al., 2009; Tedford and Martins, 2010; Riddle et al., 2007; Rafique et al., 2009). Bayley and Hartfield (2007) used 1.School of Astronautics – Beihang University – Beijing – China Author for correspondence: Fredy Marcell Villanueva | Beihang University | XueYuan Road No.37, HaiDian District, New Main Building B1011 | Beijing 100191 – China | Email: marcell385@yahoo.com Received: 16/01/13 | Accepted: 30/04/13 294 Villanueva, F.M., Linshu, H. and Dajun, X. J. Aerosp. Technol. Manag., São José dos Campos, Vol.5, No 3, pp.293-304, Jul.-Sep., 2013 GA in multidisciplinary design optimization (MDO). Davis (2001), Goldberg (1989), and Fletcher (2000) provided a comprehensive description and step-by-step implementation of the GA in design and optimization of complex systems. The insertion process of a payload into the orbit of a typical solid propellant LV start after a ballistic flight phase, when the upper stage is in the apogee altitude, which correspond to the required circular orbital altitude. At this point, to reach  the required orbital velocity, a kick impulse is necessary, and for this purpose, a KSRM is required. Its propellant should be burnt completely (depleted shutdown condition), avoiding the use of a thrust termination device, resulting in a reduced total mass of the KSRM. The insertion accuracy is maintained by using an attitude control system grouped in a payload deployment module (PDM). Thus, the objective of this article is to describe the optimization approach of a KSRM for insertion of a small payload into a prescribed circular LEO orbit by using a well-performed GA optimization method as a global optimizer. OPTIMIZATION STRATEGY The MDO approach considered aims to find the KSRM optimal design and the corresponding trajectory to successfully accomplish the specified mission, which is an insertion of a small payload into a prescribed circular LEO orbit. To accomplish this task, a widely used GA optimization method is considered. The main advantage of this method relies on its independency of initial point to perform the optimization. The GA is considered to be a powerful heuristic global optimization tool in solving complex optimization problems that traditionally has been solved using the approximated analysis. An additional advantage of the GA is the ability to solve discrete and continuous variables. The operation sequence of the GA optimization method is shown in Fig. 1. Here, the procedure starts with the assignment of the design variables, the algorithm performs several operations as population initialization, selection, crossover and mutation until the optimal solution is reached, whereas all the considered constraints are satisfied. The drawback of the GA is its computational expense, because a large number of function evaluations are required before finding the optimal solution. The main characteristics of the GA are presented in Table 1. Table 1. Genetic algorithm characteristics. Variables Characteristics Generations 200 Population size 100 Stopping criteria function tolerance 10e-6 Population type double vector Selection stochastic uniform Crossover single point pc=0.8 Mutation uniform pm=0.25641 Reproduction elite count=2 Function evaluation 2000 Figure 1. Genetic algorithm optimization approach. Design Variables Population initialization


INTRODUCTION
Using small solid rocket motors (SRM) for space applications became very attractive, because of their advantages compared with liquid propellant driven rocket engines, especially for insertion of small payloads into a circular Low Earth Orbit (LEO).Simplicity, reliability, and easy to fabricate and operate, are some of the main characteristics of the SRMs.In this article, a specific design used for insertion of a payload into a specified orbit commonly known as kick solid rocket motors (KSRM) are analyzed and discussed.
Previous works were performed on analysis of several options of upper stages for different launch vehicles (LV) (McGinnis and Joyner, 2005), and in design and optimization of upper stages for transfer from LEO to geostationary earth orbit (Motlagh and Novinzadeh, 2012).Casalino et al., (2011) carried out an optimization of a hybrid upper stage configuration with a special emphasis on the grain geometry.He Linshu and Murad (2005) developed an improved method for conceptual design of multistage solid rockets based on depleted shutdown condition, which means that the SRM has to burn all-contained propellant, resulting in a reduced gross mass.
The genetic algorithm (GA) global optimization method is increasingly being used in optimization of aerospace and propulsion systems (Kamran et al., 2009;Tedford and Martins, 2010;Riddle et al., 2007;Rafique et al., 2009).Bayley and Hartfield (2007) used GA in multidisciplinary design optimization (MDO).Davis (2001), Goldberg (1989), andFletcher (2000) provided a comprehensive description and step-by-step implementation of the GA in design and optimization of complex systems.
The insertion process of a payload into the orbit of a typical solid propellant LV start after a ballistic flight phase, when the upper stage is in the apogee altitude, which correspond to the required circular orbital altitude.At this point, to reach the required orbital velocity, a kick impulse is necessary, and for this purpose, a KSRM is required.Its propellant should be burnt completely (depleted shutdown condition), avoiding the use of a thrust termination device, resulting in a reduced total mass of the KSRM.The insertion accuracy is maintained by using an attitude control system grouped in a payload deployment module (PDM).
Thus, the objective of this article is to describe the optimization approach of a KSRM for insertion of a small payload into a prescribed circular LEO orbit by using a well-performed GA optimization method as a global optimizer.

OPTIMIZATION STRATEGY
The MDO approach considered aims to find the KSRM optimal design and the corresponding trajectory to successfully accomplish the specified mission, which is an insertion of a small payload into a prescribed circular LEO orbit.To accomplish this task, a widely used GA optimization method is considered.The main advantage of this method relies on its independency of initial point to perform the optimization.
The GA is considered to be a powerful heuristic global optimization tool in solving complex optimization problems that traditionally has been solved using the approximated analysis.An additional advantage of the GA is the ability to solve discrete and continuous variables.
The operation sequence of the GA optimization method is shown in Fig. 1.Here, the procedure starts with the assignment of the design variables, the algorithm performs several operations as population initialization, selection, crossover and mutation until the optimal solution is reached, whereas all the considered constraints are satisfied.
The drawback of the GA is its computational expense, because a large number of function evaluations are required before finding the optimal solution.
The main characteristics of the GA are presented in Table 1.

KiCK solid RoCKet motoR definition Mission analysis
The mission is to deliver a small payload to a circular LEO orbit of h f =600 km altitude.An upper stage composed of a KSRM, a PDM of 100 kg, and, a small payload of 200 kg is considered for this research.
The KSRM should accelerate the upper stage from an initial apogee velocity of V 0 =5000 m/s up to the insertion orbital velocity of V f =7560 m/s, whereas the PDM containing the attitude control system maintains the upper stage witting the allowable insertion accuracy.

KSRM design
There can be many variants of design options; however, for this research effort, a classical configuration was considered.The rocket motor case is considered to be made of high-strength steel, titanium alloy for attachment parts, the ignition system is located in the forward central part of the chamber, and the nozzle has a thrust vector control (TVC) system.

Grain characteristics
For the present analysis, a hydroxyl-terminated polybutadiene (HTPB) grain is selected due to its high performance and suitable for space applications.The composite grain has an internal port and an equivalent length was adopted instead of a complex 3D geometry.This assumption considerably simplifies the model and is acceptable for the conceptual design phase.

PRoPulsion AnAlysis
The propulsion analysis of the KSRM can be calculated using the classical approach.Sutton and Biblarz (2001) and He Linshu (2004a, 2004b) provided a detailed propulsion analysis including the essential parameters, like propellant mass flow rate, burn time, thrust, and nozzle parameters.In this analysis, a constant in time burning surface is considered; a grain geometry shape coefficient k s is used to represent the constant burning surface of the grain S b as a function of the equivalent length of grain L gn and diameter D gn as follows: The mass of the grain is calculated from the design variables and propulsion analysis, and the burn time t b , mass of the grain m gn , and mass flow rate where, u is the burning rate of propellant, ρ gn is the density of grain, L gn =L m +0.3151D m is the equivalent length of the grain, D gn =D m is diameter of the grain, L m is the rocket motor cylindrical length, λ gn is fineness ratio of the grain (grain length/ diameter), and η v is the grain volumetric loading fraction.
The nozzle throat area A t , expansion ratio ε, and nozzle exit area A e are calculated as follows: where, S b is the burning surface of grain, R c =326 J/(kg.K) is gas constant, T c =2790Kº is temperature in the combustion chamber, p e is exit pressure, p c is chamber pressure, p cmax =1.1p c is the maximum value of chamber pressure, and γ=1.21 the specific heat ratio of the gas.
The vacuum specific impulse I vac sp , and the thrust T can be calculated as shown below: where, p a is the atmospheric pressure, I a sp is average specific impulse, g is acceleration due to gravity, m .gn is mass flow, and A e is the nozzle exit area.

mAss AnAlysis
The mass of the upper stage is represented by the following equations: where, m 0 is the upper stage mass, m KSRM is kick solid rocket mass, m PDM is payload deployment module mass, m PAY is payload mass, and m st is the structural mass of the KSRM.
He Linshu (2004b) provided a detailed calculation of the KSRM structural mass, which is highly acceptable for conceptual design phase, and is shown in the following equation: It is composed of the mass of the motor cylinder, m cy , motor dome ends, m c1 and m c2 , forward and aft skirts, m q , forward and aft joints, m j1 and m j2 , forward and aft insulation liners, m in,c1 and m in,c2 , cylindrical section insulation liner, m in,cy , nozzle expansion cone, m noz,ez , nozzle spherical head, m noz,sh , nozzle insulation, m noz,in , igniter, m ig , thrust vector control, m TVC , cables, m cab , and mass of attachment parts, m ap .
The materials selected for the motor case are highstrength steel and titanium alloy, ethylene propylene diene monomer for chamber insulation, and carbon phenolic for the nozzle, whereas the factor of safety was taken as 1.5.
The mass of cylindrical part of the chamber can be calculated by the following relation: where, K cy =1.02 is the ratio of the case cylindrical length to rocket motor equivalent length (L cy /L m ), which was obtained from statistics, f is factor of safety, D ch =D m is diameter of chamber, and σ b is the ultimate strength of the chamber material.
The mass of the forward motor dome end m c1 can be calculated as follows: where, λ e =2 is the chamber ellipsoid ratio, σ is strength ratio σ=σ b /ρ cl , where, ρ cl is the density of closure material and θ 2 taken as 60°.
The relative pressure in the chamber f p is determined as shown below: The aft motor dome end is calculated as follows: where, d n =0.4D ch is the diameter of closure rear end opening for nozzle.
The mass of the forward and aft skirts of the SRM are calculated as follows: where, ρ q is the density of skirt material, δ q is thickness of skirts, and l q1 and l q2 are lengths of forward and aft skirts.
The mass of the forward joint to the igniter m j1 and the aft joint to nozzle assembly m j2 can be calculated as shown below: where, σ j is the strain of joint, ρ j is density of the joint material, and σ bj is ultimate strain of joint.The forward and aft insulation liners m in,c1 , m in,c2 and cylindrical section insulation liner m in,cy can be determined using Eqs.23 to 25.
where, ρ in is the density of insulation material and R a is the rate of ablation of the insulation material.where, L cy is the length of case cylindrical section, ε in is heat transfer coefficient of insulation material, c in is specific heat capacity of insulation, c cy is specific heat capacity of the cylindrical section, α gi is heat transfer coefficient from gas to insulation, ρ cy is density of the case cylindrical section material, T g is temperature of gas, T cy is allowable temperature of the cylindrical section, and T I is the initial temperature of the cylindrical section.
The mass of the nozzle expansion cone m noz,ec , nozzle spherical head m noz,sh , and nozzle insulation m noz,in are expressed using the Eqs.( 27), (28), and (29), respectively.where, ρ ec is the material density of the expansion cone, β n is nozzle expansion half angle, S is submerged coefficient of the nozzle, d e is the nozzle exit diameter, d t is diameter of throat, σ ec is ultimate strength of the expansion cone material, and E ec is the elastic modulus of expansion cone material.where, ρ sh is the material density of spherical head of nozzle.
where, ρ innz is the material density of the nozzle insulation, S nz is surface area of the nozzle, and δ innz is the thickness of the nozzle insulation.
The igniter mass m ig , can be calculated as follows: The nozzle thrust vector control mass m TVC , is considered as shown below: The mass of the cables m cab and attachment parts m ap are expressed in Eqs. ( 32) and (33).
where, ρ c is the linear density of the cable (kg/m).And finally, the mass of the attachment parts is calculated as follows: For trajectory analysis, a 3 degree of freedom (3DOF) model has been developed and implemented using SIMULINK (Zipfel, 2007;Fleeman, 2001).This model uses the upper stage mass, including the KSRM, PDM and payload mass, and thrust as input parameters.The upper stage is considered as a point-mass in two-dimensional coordinate systems (2D) over a spherical and non-rotating earth, where the Coriolis and centrifugal forces are not considered.
The considered state variables are velocity, flight path angle, altitude, and mass, whereas the control variable is the programmed angle of attack.The trajectory analysis computes the state variables by solving the equation of motion presented in Eq. 34, and evaluating the constraint conditions at every phase of flight.
Figure 2 illustrates the forces acting on the upper stage and below a set of governing equations of motion (He Linshu, 2004a;Xiao, 2001).where, V is the velocity, m is vehicle mass, φ is pitch angle, θ is trajectory angle, ϑ is flight path angle, η is range angle, R e is radius of earth, h is height above the ground, and l is the range.
The axial and normal overload coefficients can be calculated in a body centered velocity coordinate systems (0, x, y) and represented as follows: The density variation with altitude can be represented as follows: The variation of gravity with altitude can be represented as follows: where, ρ 0 is the sea level density, β is the density scale height, and μ is the earth gravitational parameter.
The required circular orbital insertion velocity V f for a given final altitude h f can be represented as follows:

oRbit inseRtion PRofile sequenCe
The trajectory optimization was performed considering different constraints that prevent the KSRM, its components, and the payload from failure.The trajectory is modeled considering a depleted shutdown condition at the time of insertion, as follows (He Linshu, 2004b;Qazi and He Linshu, 2005):

Initial Condition
For a typical solid propellant LV orbital payload insertion, the trajectory starts from the time where the upper stage is in the apogee altitude of h 0 =600 km, the initial condition is that the flight path angle should be ϑ 0 =0 degrees, and the initial velocity V 0 =5000 m/s, the angle of attack also should be zero degrees α 0 =0.

Pitch over insertion maneuver
After ignition, the upper stage accelerates by the power of the KSRM and flies following the orbital altitude.This phase  A second order curve is fitted from the initial time to the final insertion time; this curve represents the programmed angle for this phase.

Insertion condition
At this time all parameters should comply with the insertion accuracy requirements.The flight path angle at this time should approach zero degrees and within the boundary conditions, ϑ f ≈0, the programmed angle of attack is constrained to approach zero α f ≈0, the orbital velocity should be V f =7560 m/s corresponding to h f =600 km altitude, and the normal and axial acceleration should be less than its allowable maximum values.

oRbit inseRtion PRofile foRmulAtion
The variation of the flight path angle during insertion flight has substantial influence on the injection accuracy in orbit, acceleration loads, and final orbital velocity.It is influenced by a programmed angle of attack.Figure 3 explains the insertion maneuver, and the angle of attack is programmed using the following relations (He Linshu, 2004a;He Linshu, 2004b;Xiao, 2001): where, α max is the maximum angle of attack, α prog (t) is the programmed angle of attack, k m is the insertion maneuver variable, t is time of flight, t m is time corresponding to maximum angle of attack, t 1 is time of start of insertion maneuver, and t 2 is insertion time, coincident in value with the KSRM burning time t b .

DESIGN OPTIMIZATION PROBLEM objeCtiVe funCtion
For the present research effort, the objective is to minimize the KSRM mass m KSRM .The mathematical description of the objective function is as follows: where, X is the set of variables, X lb is lower bound of variables and X ub is the upper bound of variables.

design VARiAbles
The variables considered in the KSRM design and the insertion maneuver trajectory can be represented in Eq. ( 47) and listed in Table 2  To prevent the failure of the upper stage, KSRM and PDM, several constraints were selected, and are listed in Table 3.

design sequenCe
This section describes a step-by-step sequence for multidisciplinary design of the KSRM as follows:

OPTIMIZATION RESULT
The optimization results show that the optimized KSRM as well as the upper stage insertion trajectory profile successfully reached the objective function.The optimized values of design variables were obtained, and these variables did not violate the considered design constraints.Table 4 shows the lower bound, upper bound, and optimized values of the design variables.
There are several parameters that characterize the KSRM; however, only the most important were calculated from the optimized design variables and are presented in Table 5.     4 it can be seen that the burning time of the KSRM is 45.67 seconds, the end of burning time shows coincidence with the required insertion parameters.Additionally, the figures evidenced that the payload is inserted at the required altitude of 600 km and the obtained circular orbital velocity is 7560 m/s.The insertion flight path angle is within the required accuracy, whereas the axial and normal overloads were maintained within the constraints limits.

sensitiVity AnAlysis
Monte Carlo analysis is widely used in system and early stage of design.It provides a relatively accurate statistical evaluation of the response distribution under input uncertainties.Monte Carlo sensitivity analyses of the main KSRM parameters were conducted to investigate the effect of uncertainties of design variables over the expected result.For this analysis, a ±1% error was added to every optimized value of design variables.The results are presented in Table 6, and the scatter plot is shown in Fig. 6.
The GA optimization method considered successfully reached the optimal solution, a population of 100 with 200 generation was sufficient to perform the present study; however, several trials had been carried out to obtain the desired accuracy of the optimal solution.

Figure 2 .
Figure 2. Forces acting on the upper stage.
the end of burning time (depleted shutdown condition).

Figure 4 .
Figure 4. Insertion trajectory profile of the upper stage.

Figure 5 .
Figure 5. Flight parameters of the upper stage.

Table 4 .
Optimum values of variables.