Abstract
This study deals with the application of optimization in Finocyl grain design with ballistic objective functions using a genetic algorithm. The classical sampling method is used for space filling; a level-set method is used for simulating the evaluation of a burning surface of the propellant grain. An algorithm is developed beside the level-set code that prepares the initial grain configuration using a computer-aided design (CAD) to export generated models to the level-set code. The lumped method is used to perform internal ballistic analysis. A meta-model is used to surrogate the level-set method in an optimization design loop. Finally, a case study is done to verify the proposed algorithm. Observed results show that the grain design method reduced the design time significantly, and this algorithm can be used in designing any grain type.
1. Introduction
Using solid propellant motors is increasing due to their compactness, simplicity, and low cost [1–5]. Since the solid propellant motor firing static tests are costly during the design or qualification phase, simulation is cheapest and sometimes is the only way to predict and to understand the phenomena which may be essential for design [4]. Grain design is the most obligatory section in perfecting a solid rocket motor (SRM) design in order to satisfy the mission requirements [3]. An effective design of solid propellant grains mostly depends on the specialists’ expertise, and the performance of their design is analyzed by numerical methods or experimental static tests. If this designed grain cannot satisfy the mission requirements, it is essential to modify the geometry of grain configuration. It is clear that this try-and-error method is too expensive and is a very time-consuming method. The first attempt of grain design was performed by Sforzini [6]; other attempts are [2, 3, 7–11] and [5].
Despite development in capacity and processing speed of computers, high computing cost of complicated engineering simulations with high accuracy, is very time consuming process. The strategy is to use approximate models instead of direct use of simulation, which are often referred to as meta-models [12].
In recent years, evolutionary methods are used to solve optimization design problems. Population-based optimization procedures are an attractive choice for this problem because using them is easy and they are effective for intensively nonlinear problems [12–14].
Meta-modeling is widely used for design optimization in many engineering applications; a summary of meta-model applications in aerospace systems can be found in ref. [15]. The last works which have been done for grain design using the level-set code are [3, 5], which first implemented a genetic algorithm and secondly implemented the RBF meta-model and a sequential field approximate optimization algorithm. In this paper, a strategy based on the meta-model technique beside the evolutionary optimization technique is used instead of the direct use of the level-set method in order to decrease the run time and eliminate the divergent effects of numerical methods; besides, the advantages of using the level-set method are its capability to adapt various grain configuration burnbacks and no limitation on the grain configuration. The genetic algorithm (GA) is used as an optimizer because of global search efficiency. Discrete and continuous variables are both supported in GA, which makes it suitable for any major design. GA is able to use historical data from the previous design; GA has neither sensitivity to derivatives nor a starting solution. There are many recent articles that use this method for optimization [1–3, 11, 12, 16], and it has been proved to be a powerful optimization tool.
This paper includes five parts. Section 1 is the introduction. In Section 2, the optimization method is discussed. In Section 3, the grain burnback analysis based on the level-set method is proposed. Section 4 is devoted to using a surrogate model based on the adaptive basis function construction method. In Section 5, the possibility and precision of the proposed algorithm are demonstrated by presenting a case study.
2. Design Optimization
2.1. Objective
Finocyl grain is the most widely used type for 3D grain configuration. Finocyl grain with 13 geometry parameters is presented in Figures 1 and 2.
The objective of this paper is to achieve a thrust history of missile system designer requirement with fixed outside diameter and grain length, in which thrust history, propellant density, grain outer diameter, grain length, and nozzle geometry are design parameters. The mathematical description of design function is as follows:
is the number of design points in the objective thrust history, and the design variable is
2.2. Design Constraints
Design constraints related to Finocyl grain is presented as
Here, is
2.3. Design Boundaries
The upper and lower boundaries of design variables of Finocyl grain [10] are presented as
Here, is
2.4. Optimization Method
The aim of design optimization is to find the optimal solution to a design problem. The genetic algorithm (GA) is the most popular optimization method of the evolutionary algorithms. Exploration properties of GA make it a practical method for aerospace applications. The geometry of optimization of 3D grains requires solving nonlinear equations with numerous design variables. Since gradient-based optimization methods most likely might get caught in the local optima, GA uses a random operator in which its possibility of falling into a local optima is much smaller than the gradient-based methods. Moreover, heuristic methods are capable of applying both continuous and discrete variables suitable for large design problems. GA is able to survey data from a previous design in order to find a pattern in input parameters to produce the desired output.
A surrogate model or meta-model is a simple mathematical approximation of a simulation model that illustrates a correlation between the input and output of the system. Surrogate models are most suitable to design grain geometry, especially when using population-based heuristic methods, since these methods require several function evaluations and each evaluation is associated with a model grid generation and a numerical analysis of the level-set method. So, grain geometry optimization by GA without using surrogate models is very difficult. Using surrogate models decreases computation cost significantly and makes it possible to utilize optimization methods.
To initialize the grain design geometry, grain parameters must be recognized. Then, a design space is defined by design of experiment (DOE) methods. In this paper, a central composite design (CCD) with three levels is used which is one of the classic methods of DOE.
The next step is to develop a code in MATLAB based on the upper and lower bounds of geometry grain parameters and design space defined by the DOE method through an algorithm illustrated in Figure 3. The initial grain geometry is produced based on DOE and computer-aided design (CAD) software. The grid is established for generated geometries to prepare them to be used in the level-set code.
The pseudocode of the training is shown as Pseudocode 1:
| ||||||||||||||||||||||||||||||||||||||||||||||||||||
Pseudocode 1. | ||||||||||||||||||||||||||||||||||||||||||||||||||||
After running the level-set code for generated geometries, the burning surface vs. the web burned is produced then pressure time is predicted by the equilibrium pressure method. The internal pressure of the motor and therefore the thrust history are produced. Now, the database is defined based on geometric grain variables and design space defined by the DOE method. The meta-model is used to train data instead of expensive running of the level-set code and ballistic calculations. To overcome this problem, a surrogate model will be used. The meta-model method used here is Adaptive Basis Function Construction (ABFC). Then, a surrogate model defined by the ABFC meta-model is transformed to GA. The surrogate model in the ABFC method consists of several polynomials. They are defined for any discrete points of thrust-time diagram (pressure vs. time or burning surface vs. web-burned). These polynomials send a function and call anytime in the design loop. Then, the optimization is performed with respect to the defined constraints and objective function. Finally, an optimized geometry of grain is extracted. The resulting geometry through GA would be drawn again. After grid generation, grain burnback analysis is done by the level-set code and internal ballistic analysis is done by the equilibrium pressure method. The results are compared to results of the meta-model method. If results are within the defined tolerance, optimization is finished; otherwise, the design process should be repeated using the last iteration results. Figure 4 shows the 3D grain geometry optimization process which uses the GA based on surrogate models (meta-models).
The pseudocode of the optimization is shown as Pseudocode 2.
| ||||||||||||||||||||||||||||||||||||||||||||||||||||
Pseudocode 2. | ||||||||||||||||||||||||||||||||||||||||||||||||||||
3. Level Set-Based Grain Burnback Analysis
Due to the wide use of solid propellant motors in the aerospace industry, the accurate prediction of a solid propellant motor decreases the number of expensive tests of motors. Therefore, it decreases the development cost of motors [17]. To reach this goal, there is an intense requirement to determine the pressure time and thrust time of the SRMs which greatly depend on the changes of the burning surface area and volume of the combustion chamber through time. There are several methods to describe the solid propellant grain burnback analysis. With respect to the low performance of analytical and drawing methods, it is useful to apply numerical methods. Numerical methods generally survey the moving interface between solid propulsion and combustion gas and is classified by two procedures: frontier following method (Lagrangian method) and frontier capturing method (Eulerian method) [18]. In the frontier following method, the moving front is explicitly followed by tracking the nodes along the path of each fluid particle, and in the frontier capturing method, the grid is kept stationary and the fluid particles on the front are captured [19]. In frontier capturing methods, the front is represented by a grid function . This function is a quick computing technique that uses minimum distance function [20]. Two main equations of the Eulerian method are VOF and level set.
The level-set method is a frontier capturing method defined by a grid function with a different representation of the front; the grid function in this method is an implicit function that represents the zero-level interface to be evolved at any time [18]. In the level-set method, the velocity of the frontier can be either positive or negative. The initial position of the frontier is considered as level zero set from a higher-dimensional function which transforms the governing equations into initial value formulation. It is required that the level-set value of the frontier particle with path must always be zero [18, 19, 21]. Function is defined as a least distance between a point on the motor volume grid and the grain boundary. The level set shows a moving boundary bounding a region of in by Lipschitz function
The value of this function on the boundary (where ) is zero. Therefore,
The motion of can be considered as an evolution function for . After differentiation,
After inserting and into eq. (8), the following equation can be evolved for .
Grain burnback analysis is done with velocity into the burning surface normal direction by using the level-set method based on eq. (9). In the level-set method, the frontier is the boundary between the solid propellant grain phase and gas in the combustion chamber. is the solid propellant grain phase, is negative, and in the gas phase in the combustion chamber, is positive. Therefore, the burning surface is captured by the level-set method which can be used to calculate the burning area and free chamber volume.
The algorithm for grain burnback analysis can be divided into three main categories: grid generation, signed distance function determination, and burnback parameter calculation.
In the grid generation category, two grid types are being generated: one for the initial grain port and one for the volume of the motor. The initial grain port grid is generated using elements of the tetrahedral unstructured grid. The motor volume grid is a Cartesian grid that is generated by cubic elements. This grid includes information of coordinates of the vertices in the grid and how to connect the nodes to construct nearby elements as a numerical file called in the main program. The motor volume grid is created by a code in such a way to include the entire motor volume. In determining the signed distance function section, the level-set function is defined as a minimum distance between a point from the motor volume grid and the grain interface. Calculating the signed distance function has three steps: (1)Calculating the distance from each node of the motor volume grid and all vertices of the initial port grid and finding the minimum of this value which is considered as signed distance function for each individual node in a Cartesian grid. This process is performed for all points of the motor volume grid(2)Determining the location of each point in the motor volume grid with respect to the port interface using tetrahedral elements of the initial port grid based on the normal vector method. Then demonstrating their locations by identifying the internal and external nodes(3)Multiplying the value of the signed distance function multiplied by -1 for each point inside of the interface.
Since determining the location of each point in the Cartesian grid is too time-consuming, to reduce run time, the Gambit output file is used in which the output file surface element is determined. So, this is enough to determine the surface element location in the Cartesian grid; the run time decreases intensely.
After computing the minimum distance function in each time step, it is time to calculate the grain burnback parameters such as burning surface, port area, and port volume. Based on previous research, there are four methods to calculate the grain burnback parameters: (1)Capturing cell method [22, 23](2)Section method [24](3)Dirac Delta and Heaviside function methods [18, 25](4)Cut cell method [26]
In this paper, the cut cell method is used. In this method, elements that are cut by the interface are accounted for. This method is more precise with a shorter run-time in comparison with other methods [27].
Grain burnback analysis is performed by a level-set method and compared with a drawing method, using CAD software (SolidWorks), for a Finocyl grain represented in Figure 5. This comparison is done to validate the level-set code developed by this method.
4. Using the Meta-Model Based on Adaptive Basis Function Construction
Meta-modeling or approximation techniques were used as surrogate of expensive simulation to increase the efficiency of calculations. Now, they are used as a valuable tool in order to support a wide domain of activities in engineering, especially in design optimization [28, 29]. A surrogate model or meta-model is a simple mathematical approximation from a simulation model that is used to illustrate the correlation between the input and output of the system. In other words, finding a correlation between the system input and output reduces the number of simulations that are needed to be evaluated [28]. The meta-models can be used in engineering applications as surrogates for the detailed model when a large number of evaluations are needed, due to the high computational cost [30], as in grain geometry optimization, or when it is too time-consuming to run the detailed model for each evaluation. The mathematical model, i.e., meta-model type, suitable for the approximation could vary depending on the intended use or the underlying physics that the model should be captured.
Haftka et al. [31] discussed the relation between experiments and optimization. The advantages of meta-model-based design optimization (MBDO) are as follows: synchronizing different simulation codes, enabling parallel simulations in multiple design points, filtering noises more efficiently than gradient-based methods, and representing a general view of the design space. And since the entire domain of the design is analyzed, it is easier to detect errors in simulation.
Meta-modeling requires four steps. The first step is to choose a DOE method to generate data. The second step is to choose a model to demonstrate data. The third step is to fit the model. The fourth step is to validate the observed model data based on DOE [32].
Different sampling techniques might be used for building different meta-models. The action of determining the place of the design points in the design space is called DOE [16]. In this paper, the location of data points in the design space where true responses are evaluated is determined by a classical DOE called central composite design (CCD) [33]. The CCD method for experimental design in three variables is illustrated in Figure 6.
When running a detailed simulation model, a vector of input (design variable values), , results in a vector of output (response values), , where a vector of input contains grain geometry parameters and a vector of output contains grain burnback analysis vs. web increment.
There are different approaches to build the meta-models. Parametric techniques are based on a priori correlation between the design variables and the response. The meta-model is fitted to the dataset of design variables and corresponding responses from the detailed model by determining the coefficients of the chosen function. Nonparametric techniques are used to build different types of neural network models [16, 30].
One of the parametric techniques is regression. Most of the regression modeling methods use basis function. In nonadaptive modeling methods, the degree of basis function is fixed, and the model always contains a fixed number of basis functions. However, in adaptive modeling methods, basis functions are adapted to data by applying search methods, and this method does not have any limitations on degree of polynomials [34–36]. However, increasing the model’s polynomial degree increases the number of basis functions exponentially which increases the number of data samples making the meta-model impossible.
In this paper, the “Adaptive Basis Function Construction” (ABFC) method is used. This method requires basis functions to be constructed based on heuristic search methods. Moreover, this method makes adaptive models to be used without any restriction on the degree of the model and execution time. The number of required input variables and the degree and complexity of the target model are similar to polynomial methods, but the run time does not increase with time exponentially. Finally, it can be performed without the need to repeat the model-making process [34].
The thrust time (pressure time or burning surface vs. web-burned) diagram should be discrete to a finite number of points. Each point has two coordinates: (1) thrust and (2) time. Therefore, for each point, there are two training polynomial models, one for thrust and one for time.
The general form of a polynomial model is presented in eq. (12). In this equation, is a basis function, and coefficient can be calculated using ordinary least squares (OLS) [35].
In Equation (12), are coefficients of basis functions, and are basis functions. is a vector of grain geometry parameters (input variables); for Finocyl grain, there are 13 input variables.
The basis function for the ABFC method can be defined by input variables each with a single exponent as shown in eq. (13). In this equation, is an matrix with positive integer degrees, is the number of input variables, is the number of rows in the matrix, and is the degree of the th variable in the th basis function [36].
For example, for a discrete point (one point of the thrust-time diagram), the model is as follows:
is a vector of grain geometry parameters. Here, only four geometry parameters () are used to define one point of the thrust time diagram. The corresponding matrix with (number of basis functions) and (number of design variables used to train polynomial) is The basis function is as follows:
Finding the optimum set of basis functions () results in finding the optimum matrix with the optimum combination of positive integer values:
Using the ABFC method, polynomial models of arbitrary complexity can be generated, for each input variable with an arbitrary number of basis functions and an arbitrary power. Therefore, the process of constructing the model is extremely flexible.
By applying this regression model (ABFC) [35], to use meta-modeling/surrogate modeling is constructed to describe a relation between input variables which are geometry parameters of grain and output variables which are specific points that define thrust time (burning surface vs. web plot).
The meta-model accuracy is measured by the following criteria that would be shown in the analysis of variance (ANOVA) table: (i)Average absolute error: averages of differences between real and approximated values(ii)Maximum error: maximum of differences between real and approximated values(iii)Variance the expectation of the squared deviation of a random variable from its mean(iv)Standard deviation : the square root of the variance(v)Sum of squared errors of prediction (SSE): the sum of the squares of residuals (deviations predicted from actual empirical values of data). It is a measure of the discrepancy between the data and an estimation model. A small RSS indicates a tight fit of the model to the data(vi)Relative root mean squared error (RRMSE): RMSE divided by STD(vii)Mean square error (MSE): measures the average of the squares of the errors(viii)Root mean square error (RMSE): squared differences between real and approximated values (ix)Coefficient of determination or -squared (): the proportion of the variance in the dependent variable that is predictable from the independent variable, , is varied between 0 and 1; means that there is no error between real and approximated values. where , , and are the observed values, mean of the observed value, and approximated values, respectively.
5. Internal Ballistic Analysis
The steady-state lumped-parameter method [37] was used for predicting the pressure-time curve based on equal variations of the mass inside the combustion chamber and that exiting through the throat. After computing the burning area vs. web-burned based on the level-set method, and using propellant parameters (such as burning rate, density), nozzle throat area and characteristic velocity in combustion chamber, it is time to calculate pressure from
Thrust time is calculated using eq. (21), and the thrust coefficient is calculated using eq. (22):
This method is used when the velocity inside of the chamber is low in comparison with local sound speed.
6. Results
In this paper, to validate the grain design process, a Finocyl grain configuration with 13 independent geometric parameters based on the lowest number of input datasets using design of experiment (DOE) method has been investigated. In this method, the process of work is that the initial grain configuration should be generated in CAD software; then, input settings are performed for data. To do this, the response surface method (RSM) from the DOE method is used. Different sampling techniques can be implemented for generating different meta-models. Two sampling methods for a star grain with 7 independent geometry parameters were tested: a classical DOE method, namely, the central composite design (CCD) method, and a space-filling method called the Latin hypercube sampling (LHS) method. The number of samplings was 40 for the CCD method and the LHS method. Due to the existence of analytic methods for the star-shaped grain burnback analysis, the burning area vs. web-burned diagram was obtained from analytic relations, then trained with different meta models based on two DOE methods.
The results obtained from the CCD and the LHS methods are compared for a sampling size of 40. It is observed that, although the LHS method randomly selects the parameters from the start, the middle, and the end of the parameter range, the CCD sampling method provides better results as compared with the LHS method. As can be seen in Table 1 for different meta-models, the values of RMSE for the central composite design are lower than those of the LHS method, and values of are more than those of the LHS method.
In this step, since the CCD sampling method is selected, and Finocyl grain has 13 geometry parameters (design variables), the minimum number of samples for this method is 124. If after training with 124 samples using the Adaptive Basis Function Construction (ABFC) meta-model, the mean difference between the thrust time (pressure time or burning surface-web-burned) diagram obtained from the level set and values obtained from the meta-model was less than 0.1%, a convergence occurs. Here, with 124 samples, the average absolute error is 0.03%, so there is no need to add any sample to the database and start the training from the first.
According to Figure 3, the geometries of grains are generated automatically according to the DOE file. By entering these files into grid generator software, the user can generate a nonstructure grid with tetrahedral elements. Now the file is ready to be exported to the level-set code, to do grain burnback analysis. The grid size for this Finocyl grain is . The level-set code run for each variant on an Intel server with Pentium Core i7 2.0 GHz and 6 GB RAM takes 14 hours. ANOVA results of the grain burnback analysis by the ABFC method are displayed in Table 2. A comparison of burning surface vs. web-burned from the level-set method and training data based on ABFC is illustrated in Figures 7 and 8.
Since each run is independent of others, one can run more than 1 variant on a system depending on the value of the system’s RAM. Now the design objective is to design a grain geometry in such a way that the thrust vs. time is according with the desired diagram. The genetic algorithm method is used to find design variables of the grain geometry in order to fit the thrust history with the desired thrust history. The purpose of optimization is to minimize the difference between the grain burnback diagram with a specific basis diagram. In other words, the aim of optimization is to find values of geometric grain parameters, in a way that burning surface vs. time of grain which is trained by the ABFC method has the least difference with the assumed base diagram. The range of input geometry parameters (design variables) and constraints used for optimization was discussed before and are illustrated in eqs. 3 and 5.
The population size of individuals is considered equal to 100, and the algorithm was set to stop if upon 10 consecutive iterations, the mean difference between the objective function values obtained during the 10th iteration and the previous iteration was less than 0.001. Here, individuals are considered as grain designs, with each design containing design parameters and a thrust time (pressure time or burning surface-web-burned) diagram. Convergence occurred after 51 iterations, and the corresponding value obtained for the objective function was 54. Optimal values of the design variables for the Finocyl grain based on the desired burning surface vs. web-burned are presented in Table 3. Figure 9 illustrates the burning area vs. web-burned; in order to validate this method, grain burnback analysis is done by SolidWorks. Figures 10 and 11 show the pressure and thrust histories for the optimal grain design and the desired pressure and thrust histories based on propellant and nozzle parameters from Table 4. Grain burnback analysis of optimum grain variant is illustrated in Figure 12.
Here, is number of Fins, which is an integer design variable. So, is a discrete variable; in order to use it in GA, it is mapped to a continuous variable then mapped to a discrete variable using floor function.
7. Conclusions
In this paper, a procedure is proposed to perform a solid propellant motor grain design. The sequence of the task consists of several sections: the first step of this method was to create a module that creates initial grain geometry. To do this, first a DOE program is used to appoint design points in the design space. A central composite sampling method is employed for this purpose. Second, a MATLAB code is used to create geometric parameters of each grain configuration, and third, an application programming interface (API) toolkit code is used to create initial grain configuration after the grain geometry model is generated. Now, it is time for grid generation and then meshed geometry is ready to perform grain burnback analysis by the level-set code. A dataset will be created by these level-set code runs. Using the level-set code is one of the important parts of this research because of its accuracy and capability to do the burnback analysis for any grain configuration. The training grain burnback analysis is done, using a meta-model based on Adaptive Basis Function Construction. This method in comparison with the level-set method greatly decreases the run time, and it is possible to perform burnback analysis for any grain configuration in this design space. The level-set code takes 14 hours to execute a single variant compared with the meta-model which only takes 0.3 seconds for each variant. A genetic algorithm is used for optimization, and a surrogate model is used to evaluate objective function instead of the level-set code. The most important part of this research was to train grain burnback analysis, using the surrogate model based on adaptive basis function construction, especially when using population-based optimization methods. Run time decreases greatly, and one can rapidly use this surrogate model for any optimization objective. The advantages of the meta-model in this research are its ability to connect expensive simulation codes (level set, CAD, grid generators, internal ballistic codes, etc.), ability in parallel simulations in multiple design points (level set codes can be run independently and parallel and the meta-model can be used in optimization), ability to make filter noises better than gradient-based methods, ability to represent a general view of the design space, and, since the entire domain of design is analyzed, ease in detecting errors in simulation. Finally, a Finocyl grain geometry is used to validate this method. Modeling results confirm that the proposed algorithm can be used instead of level-set grain burnback analysis.
Glossary
| N: | Number of fins |
| : | Motor front bore |
| : | Grain radius |
| : | Motor rear bore |
| : | Middle bore |
| : | Fin height |
| : | Fin thickness |
| : | Grain length |
| : | Front web |
| : | Front cone |
| : | Rear cone |
| : | Rear cylinder |
| : | Fin length |
| [deg.]: | Fin angle |
| : | Fin fillet |
| : | Maximum chamber pressure |
| : | Chamber pressure |
| : | Nozzle exit pressure |
| : | Ambient pressure |
| : | Propellant mass |
| : | Time step |
| : | Velocity orthogonal to boundary |
| : | Level set function |
| : | Cartesian grid steps |
| : | Basis function |
| : | Throat area |
| : | Burning area |
| : | Propellant density |
| : | Burning rate coefficient |
| : | Burning rate index |
| : | Expansion ratio |
| : | Characteristic velocity |
| : | Specific heat ratio |
| : | Sum of squared error |
| : | Mean squared error |
| : | Root mean squared error |
| : | Standard deviation |
| : | Relative root mean squared error |
| : | Variance |
Data Availability
The initial data used to support the findings of this study are included within the article. Other data based on request for data will be considered by the corresponding author.
Conflicts of Interest
The authors of this manuscript declare that there is no conflict of interest regarding the publication of this paper.