Multiobjective Optimization Method of Solid Rocket Motor Finocyl Grain Based on Surrogate Model

Abstract

To improve the performance of a solid rocket motor (SRM), a multiobjective optimal design method that can consider the structural integrity, internal ballistic performance, and loading performance of the SRM was proposed based on parametric modeling and surrogate modeling technology. Firstly, the parametric modeling technology was introduced into the field of structural integrity analysis for a high-loading SRM, based on which the influences of load and geometric parameters on the maximum von Mises strain of the SRM grain were analyzed, which effectively improved the sampling speed and prediction accuracy of the surrogate model. Combining the calculation models of the combustion surface area and volume loading fraction of the SRM, the Pareto optimal solution set was obtained based on the NSGA-II algorithm. Under the constraints of the optimization model, the maximum von Mises strain can be reduced by up to 26.72% and the volume loading fraction can be increased by up to 1.83% compared with the original. In addition, the optimal design method proposed in this paper is significantly superior in efficiency, capable of reducing both the single sampling time by more than 95% and the number of numerical simulations from 20,000 to 400, and the average prediction deviation is only 1.87%.

1. Introduction

The new generation of missile weapon systems puts forward the performance requirements of high loading, high-pressure strength, and high reliability for SRM. However, high-loading SRMs are always accompanied by high stress and strain, which brings serious challenges to its grain structural integrity. The multiobjective optimization design of the grain shape of the SRM is an effective method to balance the conflicts between the grain structural integrity, the internal ballistic performance, and the loading performance. However, the process of optimizing the SRM charge shape is repetitive, and obtaining the structural response of the SRM through numerical analysis is an extremely time-consuming job. For this reason, there is an urgent need to explore efficient multiobjective optimization design methods to improve the level of design of complex high-loading SRMs.

At present, the parametrized modeling method has been mainly applied to the internal ballistic calculation models and structural analysis models of some simple SRMs [1,2,3]. Some scholars have carried out a series of works based on CATIA, Pro/E, SolidWorks, and other software for parametrized modeling of SRMs [3,4,5] and realized the calculation and optimization for the internal ballistic performance of SRMs with various grain shapes [6,7,8]. Some scholars have also introduced parametrized modeling methods into the field of structural analysis of SRMs [9,10,11] and applied them to some SRMs with low loading fraction and relatively simple grain shapes.

However, the parametrized modeling method for high-loading SRMs is still lacking in-depth research. The reason for this is that the realization of the analysis process of high-loading SRMs has a drastic conflict between the complex configuration of the grain and the limited preprocessing capability of the finite element software.

As for the surrogate modeling technique, it has been widely used for various structural optimization design problems since it was proposed in the 1960s [12,13], and some scholars have conducted some preliminary studies for SRM. Ye et al. [14] conducted a multidisciplinary optimization design for SRM based on the Kriging model, considering the component cost and the internal ballistic performance; Wu et al. [15] applied the surrogate modeling technique to the grain design field of SRM and obtained the smoothest charge design solution for the thrust.

As for the multiobjective optimal design of SRMs, Miao et al. [16] carried out a multiobjective optimal design of the stress release boot for several performance indicators related to the structural integrity of an SRM. Tola [17] carried out a multiobjective optimal design of a simple two-dimensional star-hole-type grain considering both internal ballistic performance and structural strength criteria.

In summary, for the optimization design of the grain of a high-loading SRM, the complex structural form causes great difficulties in obtaining the structural response data, and there is still a lack of efficient optimization methods that consider the structural integrity, internal ballistic performance, and loading performance of the SRM grain.

In this paper, the implementation of parametrized modeling methods in the field of structural analysis of a high-loading SRM is investigated based on Abaqus software and Python language. On this basis, the influences of geometric parameters and loads on the structural integrity of the grain are analyzed for the SRM with finocyl grain. Then, the applicability of different experimental design methods and surrogate models to the SRM structure optimization problem was analyzed and tested, and the process of establishing the surrogate model for structure optimization was finalized. Finally, the NSGA-II algorithm was used to realize the multiobjective optimization design of the SRM grain by combining the internal ballistic and loading performance calculation models.

2. Model and Methods

In order to implement the multiobjective optimization design of SRM grain, we first need to determine the design variables, objective functions, and constraints, and then establish the mathematical models after transforming it into mathematical problem.

2.1. Optimized Objects

The SRM with finocyl grain is widely used in industry for its advantages of high charge–volume fraction and flexible combustion surface adjustment, etc. [18]. We take a fixed-volume SRM with finocyl grain with an aspect ratio of 5.6 and a diameter of 1.7 m as the research object, whose case dimensions have been fixed. The grain is divided into parts I, II, and III, as shown in Figure 1.

The SRM is in the stage of optimization design of the grain shape after the overall scheme is determined, its internal ballistic performance is basically confirmed, and the case size, grain thickness, and internal hole radius are basically fixed. Therefore, the design variables in this paper are the eighteen configuration parameters associated with finocyl in grain parts I and III, that is,

$$X=(L,R,l_1,l_2,r_1,r_2,\theta ,{\theta }_1,{\theta }_2,L^{\prime },R^{\prime },{l^{\prime }}_1,{l^{\prime }}_2,r_1{}^{\prime },{r^{\prime }}_2,{\theta }^{\prime },{{\theta }^{\prime }}_1,{{\theta }^{\prime }}_2)$$

(1)

where $L,R,\theta ,l_1,l_2,r_1,r_2,{\theta }_1,{\theta }_2$ is the parameter of the finocyl in grain part III; the specific meaning of the parameter is shown in Figure 2; $L^{\prime },R^{\prime },{l^{\prime }}_1,{l^{\prime }}_2,{r^{\prime }}_1,{r^{\prime }}_2,{\theta }^{\prime },{{\theta }^{\prime }}_1,{{\theta }^{\prime }}_2$ is the parameter of the finocyl in grain part I; while the initial value $X_0$ and the range $[X_L,X_U]$ of design variables are shown in

Table 1. Initial values and range of values for the design variables of the SRM grain.

Variable	L/mm	R/mm	l1/mm	l2/mm	r1/mm	r2/mm	θ/°	θ1/°	θ2/°
xL	500.0	3.0	30.0	5.0	40.0	60.0	10.0	80.0	40.0
x0	739.2	18.0	54.3	25.0	60.0	80.0	19.0	90.0	49.0
xU	900.0	27.0	60.0	60.0	80.0	100.0	30.0	90.0	90.0
Variable	L′/mm	R′/mm	${l^{\prime }}_1/\mathrm{mm}$	${l^{\prime }}_2/\mathrm{mm}$	${r^{\prime }}_1/\mathrm{mm}$	${r^{\prime }}_2/\mathrm{mm}$	${\theta }^{\prime }/°$	${{\theta }^{\prime }}_1/°$	${{\theta }^{\prime }}_2/°$
xL	350.0	4.0	5.0	5.0	40.0	40.0	10.0	80.0	40.0
x0	688.0	10.0	42.5	25.0	60.0	60.0	16.1	90.0	40.7
xU	800.0	16.0	60.0	60.0	80.0	80.0	30.0	90.0	90.0

.

2.2. Optimization Objectives and Models

For the SRM, structural integrity, internal ballistic performance, and loading performance are the three key performances for determining its launch capability.

1. Structural integrity

The main failure mode of SRM is the destruction of its structural integrity, especially the structural integrity of its grain. From the beginning of production to the completion of the launch mission, its life history and the corresponding loads are shown in Figure 3. Long-term experience indicates that temperature and pressure loads have the most significant impact on the structural integrity of the SRM, so the load case of the SRM studied in this paper is a combined load consisting of temperature and pressure.

The SRM grain is in a multidirectional stress state under temperature or pressure load, so it is appropriate to use the octahedral shear strain criterion, which is $r_8\le \frac{r_{8m}}{n}$, where $r_{8m}$ is the critical value obtained by the experiment and $n$ is the safety factor.

The expression for the octahedral shear strain is given by

$$r_8=\frac{1}{3}\sqrt{{({\epsilon }_x-{\epsilon }_y)}^2+{({\epsilon }_y-{\epsilon }_z)}^2+{({\epsilon }_z-{\epsilon }_x)}^2+6({\epsilon }_{xy}^2+{\epsilon }_{yz}^2+{\epsilon }_{zx}^2)}$$

(2)

and the expression of von Mises strain is

$${\epsilon }_v=\frac{\sqrt{2}}{3}\sqrt{{({\epsilon }_x-{\epsilon }_y)}^2+{({\epsilon }_y-{\epsilon }_z)}^2+{({\epsilon }_z-{\epsilon }_x)}^2+6({\epsilon }_{xy}^2+{\epsilon }_{yz}^2+{\epsilon }_{zx}^2)}$$

(3)

Therefore, ${\epsilon }_v=\frac{\sqrt{2}}{2}{r}_8$. It can be obtained that the von Mises strain is essentially equivalent to the octahedral shear strain. In addition, for solid propellants with Poisson ratio tending to 0.5, the results of the monotonic tensile test can be applied more simply when using the von Mises strain criterion.

For the SRM, the lower the von Mises strain, the safer the structure of SRM. For this reason, the objective function related to the structural integrity of the SRM was set as

$$\min f_{{\epsilon }_{v,\max }}(X)=\frac{[{\epsilon }_{v,\max }(X)-{\epsilon }_{v,\max }^0]}{{\epsilon }_{v,\max }^0}\times 100\%$$

(4)

where $f_{{\epsilon }_{v,\max }}(X)$ is a function of the maximum von Mises strain ${\epsilon }_{v,\max }$ about $X$, ${\epsilon }_{v,\max }(X)$ is the ${\epsilon }_{v,\max }$ for different $X$, and ${\epsilon }_{v,\max }^0$ is the maximum von Mises strain before optimization.

2. Internal ballistic performance

For the internal ballistic performance of a SRM, the combustion chamber pressure is one of the most critical parameters, because its magnitude and variation law not only directly determine the SRM thrust scheme, but also affect the whole process of combustion and the design of the structural strength of the SRM.

In zero-dimensional internal ballistics, assuming a uniform pressure in the combustion chamber and ignoring the effect of gas flow [19], the basic differential equation for the time variation of the combustion chamber pressure can be obtained according to the conservation of mass and the gas equation of state as

$$\frac{V_c}{{\mathsf{\Gamma }}^2{c}^{*}{}^2}\frac{d{p}_c}{dt}=p_p{A}_b\alpha p_c^n-\frac{p_c{A}_t}{c^{*}}$$

(5)

where $V_c$ is the free volume of the combustion chamber, $\mathsf{\Gamma }$ is the specific heat ratio function, $c^{*}$ is the characteristic velocity of the propellant, $p_c$ is the gas density, $t$ is the time, $p_p$ is the propellant density, $A_b$ is the area of the combustion surface, $\alpha ,n$ is the calculated correlation coefficient of the combustion velocity, and $A_t$ is the cross-sectional area of the nozzle throat.

For the SRM that has completed the preliminary design of internal ballistics, the combustion chamber shape, propellant type, and nozzle design have been determined, so the change in combustion chamber pressure and combustion surface area with time is basically the same. For this reason, the objective function related to the internal ballistic performance is set as follows: the deviation of the initial combustion surface area from the target combustion surface area is minimized, that is,

$$\min f_{A_b}(X)=\frac{\left|A_{{}_b}(X)-A_b^0\right|}{A_b^0}\times 100\%$$

(6)

where $f_{A_b}(X)$ is a function of the initial combustion surface area $A_{{}_b}$ about $X$, $A_b(X)$ is the initial combustion surface area for different $X$, and $A_b^0$ is the initial combustion surface area before optimization.

3. Loading performance

The volume loading fraction reflects the loading performance of the SRM when the combustion chamber size is constant, and the increase in propellant loading can directly improve the total impulse, so the volume loading fraction should be increased as much as possible during the design process. Therefore, the objective function related to the loading performance is set as follows: maximize the increase in the volume loading fraction relative to that before optimization, that is,

$$\max f_{{\eta }_v}(X)=\frac{[{\eta }_v(X)-{\eta }_v^0]}{{\eta }_v^0}\times 100\%$$

(7)

where $f_{{\eta }_v}(X)$ is a function of the volume loading fraction ${\eta }_v$ about $X$, ${\eta }_v(X)$ is the volume loading fraction for different $X$, and ${\eta }_v^0$ is the volume loading fraction before optimization.

In addition, the optimization objectives should also satisfy certain conditions, namely that the maximum von Mises strain of the SRM grain should not be higher than before the optimization, the volume filling fraction of the SRM should not be lower than before the optimization, and the deviation between the initial combustion surface area and the target combustion surface area should be within 5%.

In summary, the mathematical model of the multiobjective optimization problem of the SRM grain is

$$\left\{\begin{array}{l}\hspace{1em}\hspace{1em}\hspace{1em} \min f_{{\epsilon }_{v,\max }}(X)\\ \hspace{1em}\hspace{1em}\hspace{1em} \min f_{A_b}(X)\\ \hspace{1em}\hspace{1em}\hspace{1em} \max f_{{\eta }_v}(X)\\ s.t. f_{{\epsilon }_{v,\max }}(X)\le 0.00\%\\ \hspace{1em} f_{A_b}(X)\le 5.00\%\\ \hspace{1em} f_{{\eta }_v}(X)\ge 0.00\%\\ \hspace{1em} {\rm X}\in \left[{{\rm X}}_L,{{\rm X}}_U\right],X=(L,R,\theta ,l_1,l_2,r_1,r_2,{\theta }_1,{\theta }_2,L^{\prime },R^{\prime },{\theta }^{\prime },{l^{\prime }}_1,{l^{\prime }}_2,r_1{}^{\prime },{r^{\prime }}_2,{{\theta }^{\prime }}_1,{{\theta }^{\prime }}_2)\end{array}\right.$$

(8)

3. Grain Structural Integrity Analysis Based on Parametric Modeling Technology

The key and difficulty of the multiobjective optimization of the grain shape is how to obtain the structural response under different design variables. The optimization process requires repeated invocations of the numerical simulation model to obtain the structural response data, and the high-precision simulation analysis model is always accompanied by high computational costs, resulting in long optimization cycles and low solution efficiency.

In order to improve the optimization efficiency, the parametric modeling method was first introduced into the field of structural analysis of a high-loading SRM, and it was further combined with the surrogate modeling technique to establish the mapping relationship between input and output data through limited sample data, which can effectively reduce the total computational effort of the optimization process.

3.1. Parametric Modeling Technology for SRM

The structural integrity analysis of the grain of the SRM is a really complex and time-consuming process, and the manual modeling time is about tens of hours. For this reason, a Python script is used to replace the user’s pre- and postprocessing operations, which reduces the cost of the structural analysis to tens of minutes; the basic principle is shown in Figure 4.

The Python script for the SRM mainly contains six functions: geometric modeling, meshing, material definition, boundary setting, load setting, and result extraction. The content and order of the script are shown in Figure 5.

The specific features of the six functions are:

1. Function of geometric modeling: since the SRM grain has nine holes with the same circumferential position on the front and rear sides, combined with its cyclic symmetry, the 1/18 geometric model, as shown in Figure 6a, will be created by the function of geometric modeling.

2. Function of meshing: dividing the geometric model automatically into the finite element model shown in Figure 6b, which has the number of hexahedral elements at around 150,000.

3. Function of material definition: assigning the material properties of each part to the corresponding geometry. The material parameters of the case, insulation, cladding, and propellant are shown in

Table 2. Material parameters of the SRM.

Material Parameters	Elastic Modulus/MPa	Poisson’s Ratio	Density/(kg/m3)	Expansion Coefficient/°C−1
Case	1.86 × 105	0.300	7.9 × 103	1.10 × 10−5
Insulator	60.00	0.498	2.1 × 103	2.20 × 10−5
Cladding	Prony series	0.498	1.86 × 103	8.60 × 10−5
Propellant	Prony series	0.498	1.86 × 103	8.60 × 10−5

. The propellant and cladding are the same viscoelastic materials, and their relaxation modulus is characterized using the Prony series as shown in Equation (9), and the first five order parameters are shown in

Table 3. The first 5 order parameters of Prony series.

i	0	1	2	3	4	5
τi		0.0037	0.1099	3.2977	98.9295	2967.8852
Ei	18.2082	8.3728	4.2322	2.3112	0.9354	0.8212

.

$$E\left(t\right)=E_0-{\displaystyle \sum _{i=1}^n{E}_i\left(1-e^{-t/{\tau }_i}\right)}$$

(9)

The time–temperature equivalence equation of propellant and cladding is shown in Equation (10), where ${\alpha }_T$ is the translational factors at different temperatures; $C_1$, $C_2$ is the material constant, determined by the material’s own characteristics, $C_1=5.0474$, $C_2=144.9207$; $T_0$ is the reference temperature, $T_0=293.15$.

$$\lg {\alpha }_T=\frac{-C_1(T-T_0)}{C_2+(T-T_0)}$$

(10)

4. Function of boundary setting: applying symmetric displacement constraints on both sides of the SRM according to the cyclic symmetry, and then applying axial displacement constraints on the intersection of the tail head and the straight section of the case.

5. Function of load setting: setting the analysis step and loading the load to the finite element model of the SRM. When the temperature load is applied, a linear decrease in temperature from 58 °C to 20 °C in 86,400 s is placed on the SRM; when the internal pressure load is applied, a linear increase in pressure to 10 MPa in 0.3 s is placed on the internal surface of the SRM [20].

3.2. Analysis of Influence of Load

Elasticity theory indicates that for the plane strain model, the stress concentration coefficients caused by uniform temperature load and uniform internal pressure load are the same. For the linear viscoelastic three-dimensional grain model used in this paper, the grain will shrink under the temperature load, which was similar to the surrounding pressure environment in which the pressure load was applied, and the von Mises stress cloud of the grain when the temperature and pressure load are applied separately is shown in Figure 7.

From Figure 7, it can be obtained that the stress cloud diagrams of the grain are extremely similar when the temperature load and the pressure load are applied separately, and the locations of the danger points of each part of the grain are exactly the same, in which the danger points of grain part I and III appear near the junction point between the tail of the finocyl and the inner hole, while the danger points of grain part I appear on the central surface of the inner hole, and the maximum von Mises strains of three parts are described as ${\epsilon }_{v,\max }^{\mathrm{I}},{\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}},{\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}\mathrm{I}}$, respectively.

The ${\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}\mathrm{I}}$ was taken as the object of study, keeping the rest parameter fixed and making $L$ and $\theta$ change in their respective ranges, while the ${\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}\mathrm{I},t}$ under temperature load, the ${\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}\mathrm{I},p}$ under pressure load and the ${\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}\mathrm{I},t+p}$ under combined load were extracted, respectively, as shown in Figure 8.

From Figure 8, the effects of temperature and pressure load are almost directly superimposable for ${\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}\mathrm{I},t}+{\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}\mathrm{I},p}\approx {\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}\mathrm{I},t+p}$. From von Mises strain Equation (3), it is obtained that

$${\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}\mathrm{I},t+p}=\frac{\sqrt{2}}{3}\sqrt{{({\epsilon }_1^{\mathrm{I}\mathrm{I}\mathrm{I},t+p}-{\epsilon }_2^{\mathrm{I}\mathrm{I}\mathrm{I},t+p})}^2+{({\epsilon }_1^{\mathrm{I}\mathrm{I}\mathrm{I},t+p}-{\epsilon }_3^{\mathrm{I}\mathrm{I}\mathrm{I},t+p})}^2+{({\epsilon }_2^{\mathrm{I}\mathrm{I}\mathrm{I},t+p}-{\epsilon }_3^{\mathrm{I}\mathrm{I}\mathrm{I},t+p})}^2}$$

(11)

Combined with the principal strain situation of the grain, for linear viscoelastic materials, the position of ${\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}\mathrm{I}}$ and the direction of the principal strain under the effect of temperature or pressure load are completely same. For this reason, the principal strain under combined load is equal to the sum of the principal strains when temperature and pressure loads are applied separately, that is,

$${\epsilon }_i^{\mathrm{I}\mathrm{I}\mathrm{I},t+p}={\epsilon }_i^{\mathrm{I}\mathrm{I}\mathrm{I},t}+{\epsilon }_i^{\mathrm{I}\mathrm{I}\mathrm{I},p}(i=1,2,3)$$

(12)

where ${\epsilon }_i^{\mathrm{I}\mathrm{I}\mathrm{I},t},{\epsilon }_i^{\mathrm{I}\mathrm{I}\mathrm{I},p},{\epsilon }_i^{\mathrm{I}\mathrm{I}\mathrm{I},t+p}$ are the principal strains at the danger point under the effect of temperature or pressure load alone or in combination, respectively. Obviously, when the principal strains satisfy

$${\epsilon }_i^{\mathrm{I}\mathrm{I}\mathrm{I},t}/{\epsilon }_i^{\mathrm{I}\mathrm{I}\mathrm{I},p}=K$$

(13)

then

$${\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}\mathrm{I},t}=\frac{\sqrt{2}}{3}K\sqrt{{({\epsilon }_1^{\mathrm{I}\mathrm{I}\mathrm{I},p}-{\epsilon }_2^{\mathrm{I}\mathrm{I}\mathrm{I},p})}^2+{({\epsilon }_1^{\mathrm{I}\mathrm{I}\mathrm{I},p}-{\epsilon }_3^{\mathrm{I}\mathrm{I}\mathrm{I},p})}^2+{({\epsilon }_2^{\mathrm{I}\mathrm{I}\mathrm{I},p}-{\epsilon }_3^{\mathrm{I}\mathrm{I}\mathrm{I},p})}^2}=K{\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}\mathrm{I},p}$$

(14)

Therefore,

$${\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}\mathrm{I},t}+{\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}\mathrm{I},p}=\frac{\sqrt{2}}{3}(1+K)\sqrt{{({\epsilon }_1^{\mathrm{I}\mathrm{I}\mathrm{I},p}-{\epsilon }_2^{\mathrm{I}\mathrm{I}\mathrm{I},p})}^2+{({\epsilon }_1^{\mathrm{I}\mathrm{I}\mathrm{I},p}-{\epsilon }_3^{\mathrm{I}\mathrm{I}\mathrm{I},p})}^2+{({\epsilon }_2^{\mathrm{I}\mathrm{I}\mathrm{I},p}-{\epsilon }_3^{\mathrm{I}\mathrm{I}\mathrm{I},p})}^2}=(1+K){\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}\mathrm{I},p}$$

(15)

and from Equations (11) and (12),

$${\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}\mathrm{I},t+p}=\frac{\sqrt{2}}{3}(1+K)\sqrt{{({\epsilon }_1^{\mathrm{I}\mathrm{I}\mathrm{I},p}-{\epsilon }_2^{\mathrm{I}\mathrm{I}\mathrm{I},p})}^2+{({\epsilon }_1^{\mathrm{I}\mathrm{I}\mathrm{I},p}-{\epsilon }_3^{\mathrm{I}\mathrm{I}\mathrm{I},p})}^2+{({\epsilon }_2^{\mathrm{I}\mathrm{I}\mathrm{I},p}-{\epsilon }_3^{\mathrm{I}\mathrm{I}\mathrm{I},p})}^2}=(1+K){\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}\mathrm{I},p}$$

(16)

Combining Equations (15) and (16), we can obtain

$${\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}\mathrm{I},t}+{\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}\mathrm{I},p}={\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}\mathrm{I},t+p}$$

(17)

The numerical analysis result shows that under the material parameters and load conditions of this paper, the ratio of principal strain under temperature and pressure load is approximately equal to 0.39, that is,

$${\epsilon }_i^{\mathrm{I}\mathrm{I}\mathrm{I},t}/{\epsilon }_i^{\mathrm{I}\mathrm{I}\mathrm{I},p}\approx K=0.39(i=1,2,3)$$

(18)

Therefore,

$${\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}\mathrm{I},t}+{\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}\mathrm{I},p}\approx {\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}\mathrm{I},t+p}\approx (1+K){\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}\mathrm{I},p}$$

(19)

For this reason, when it is required to obtain ${\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}\mathrm{I},t+p}$ only the numerical simulation of temperature load or pressure load can be performed, which can save half of the time for the obtaining of the structural response data.

3.3. Analysis of the Influence of Parameters

The Saint-Venant principle indicates that when the magnitude of the force is constant, the distribution of the load has little effect on the area far from the effect of the load. For this reason, grain part III shown in Figure 2 was studied, and the parameters were grouped into those close to the danger point—$L,d,R,\theta$—and those far from the danger point—$l_1,l_2,r_1,r_2,{\theta }_1,{\theta }_2$.

Firstly, $l_1,l_2,r_1,r_2,{\theta }_1,{\theta }_2$ were analyzed, keeping $L,d,R,\theta$ constant, and adjusting $l_1,l_2,r_1,r_2$ and ${\theta }_1,{\theta }_2$, respectively, the variation of ${\epsilon }_{v,\max }^{\mathrm{I},p},{\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I},p},{\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}\mathrm{I},p}$ with $l_1,l_2,r_1,r_2$ and ${\theta }_1,{\theta }_2$ are shown in Figure 9 and Figure 10, in which the left and right cloud maps are the configuration and von Mises stress cloud maps of the grain with the first and last group of parameters, respectively.

From Figure 9 and Figure 10, it can be obtained that ${\epsilon }_{v,\max }^{\mathrm{I},p},{\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I},p},{\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}\mathrm{I},p}$ fluctuate within a very small range when $l_1,l_2,r_1,r_2,{\theta }_1,{\theta }_2$ are varied with constant $L,d,R,\theta$, and the positions of the danger points remain basically unchanged. On the one hand, it indicates that the change in $l_1,l_2,r_1,r_2,{\theta }_1,{\theta }_2$ does not affect the structural response of grain part I and II, and on the other hand, it also indicates that the change in $l_1,l_2,r_1,r_2,{\theta }_1,{\theta }_2$ does not affect the von Mises strain response at the danger point of grain part III.

To analyze the influences of $L,d,R,\theta$, ${\epsilon }_{v,\max }^{\mathrm{I}},{\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}},{\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}\mathrm{I}}$ at different $L,d,R,\theta$ are shown in Figure 11.

From Figure 11, it can be obtained that:

1. ${\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}\mathrm{I},p}$ shows a monotonic decreasing trend with the increase in $d,R$, and the relative increase in $R$ can make ${\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}\mathrm{I},p}$ decrease when the sum of $d$ and $R$ is constant;

2. ${\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}\mathrm{I},p}$ gradually decreases with the increase in $L$, and the decreasing trend gradually becomes slower;

3. ${\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}\mathrm{I},p}$ shows the trend of decreasing and then increasing with the increase in $\theta$, and reaches the minimum value near $\theta =30°$, and ${\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}\mathrm{I},p}$ does not change with $\theta$.

4. ${\epsilon }_{v,\max }^{\mathrm{I},p},{\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I},p}$ also does not change with $L,d,R,\theta$.

In summary, based on the analysis of the effects of loads and parameters, for the acquisition of ${\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}\mathrm{I},t+p}$ with different parameters, based on the superposition analysis of load effects in Section 3.2, we can firstly obtain that

$${\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}\mathrm{I},t+p}=(1+K){\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}\mathrm{I},p}$$

(20)

After partitioning the SRM as shown in Figure 1, Equation (20) is transformed into

$${\epsilon }_{v,\max }^{t+p}=(1+K)\max ({\epsilon }_{v,\max }^{\mathrm{I},p},{\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I},p},{\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}\mathrm{I},p})$$

(21)

Combined with the conclusions obtained from the parameter influence analysis in Section 3.3, we can obtain

$$\left\{\begin{array}{c}{\epsilon }_{v,\max }^{\mathrm{I},p}=f_{\mathrm{I}}^p(d^{\prime },R^{\prime },L^{\prime },{\theta }^{\prime })\\ {\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}}={\epsilon }_0^p\\ {\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I},p}=f_{\mathrm{I}\mathrm{I}\mathrm{I}}^p(d,R,L,\theta )\end{array}\right.$$

(22)

where ${\epsilon }_0^p$ is the maximum von Mises strain of part II of grain and ${\epsilon }_0^p$ is a constant.

Then, according to conclusion (2) of Section 3.3, controlling $d^{\prime }=2.0 \mathrm{mm}$ and $d=3.0 \mathrm{mm}$ in the subsequent analysis, Equation (22) is simplified to

$$\left\{\begin{array}{c}{\epsilon }_{v,\max }^{\mathrm{I},t+p}=(1+K)f_{\mathrm{I}}^p(R^{\prime },L^{\prime },{\theta }^{\prime })\\ {\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}}=(1+K){\epsilon }_0^p\\ {\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I},t+p}=(1+K)f_{\mathrm{I}\mathrm{I}\mathrm{I}}^p(R,L,\theta )\end{array}\right.$$

(23)

4. Establishment and Validation of Strain Surrogate Model

The multiobjective optimization design of SRM grain requires repeated invocation of numerical simulation models to obtain structural response data. High-precision simulation analysis models are always accompanied by high computing costs, resulting in long optimization cycles and low solution efficiency. For this reason, the surrogate model technology is introduced to establish the mapping relationship between ${\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}\mathrm{I},t+p}$ and $X$ to improve the optimization efficiency and reduce the optimization cost.

Based on the conclusion of Section 3.2 and Equations (21) and (23), the establishment and validation process of the surrogate model for ${\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}\mathrm{I},t+p}$ is shown in Figure 12.

4.1. Acquisition of Sample Data

In order to enable the surrogate model to effectively predict the response in the design space, a certain number of initial sample points need to be obtained first. The acquisition of sample points is based on the experimental design method, and to make the surrogate model effectively replace the true objective function, a suitable experimental design method should be selected so that the initial sample points can be distributed as uniformly as possible in the design space.

The most widely used experimental design method is optimal Latin hypercube sampling (OPLHS). OPLHS is based on Latin hypercube sampling (LHS), in which the sampling process will evaluate and optimize different combinations based on the total spacing between sample points using the optimization algorithm, and finally derive the most uniformly distributed combination of sample points.

For this purpose, 200 data samples were collected in the sampling space shown in Figure 1 for grain parts I and III, respectively, and the spatial distribution of the sampling points is shown in Figure 13.

4.2. Comparison and Validation of Surrogate Model

The commonly used surrogate modeling techniques include the response surface methodology (RSM) [21], Kriging [22] and radical basis function (RBF) [23] models, etc. The basic mathematical expression of the RBF model is

$$\tilde{f}\left(x\right)={\displaystyle \sum _{i=1}^N{w}_i}\phi (‖x-x_i‖)$$

(24)

where $w_i$ is the weight coefficient of the $i$th basis function.

To compare and verify the prediction accuracy of the three models, a cross-validation method was used, where cross-validation means that for the initial sampling data obtained, most of them are used to build a surrogate model, and a small part is used to compare the deviation between the actual and predicted values, and then evaluate the prediction accuracy of the model. There are four commonly used indicators for model accuracy evaluation as follows.

• Root-mean-square error, RMSE

$$RMSE=\sqrt{\frac{1}{m^{\prime }}{\displaystyle \sum _{i=1}^{m^{\prime }}{(y_i-{\widehat{y}}_i)}^2}}$$

(25)

2

$R^2$

$$R^2=1-\frac{{\displaystyle \sum _{i=1}^{m^{\prime }}{(y_i-{\widehat{y}}_i)}^2}}{{\displaystyle \sum _{i=1}^{m^{\prime }}{(y_i-\overline{y})}^2}}$$

(26)

where $\overline{y}$ is the mean response value of the $m^{\prime }$ sample points for validation.

3

Average absolute error, $AAE$

$$AAE=\frac{\frac{1}{m^{\prime }}{\displaystyle \sum _{i=1}^{m^{\prime }}\left|y_i-{\widehat{y}}_i\right|}}{\max (y_i)-\min (y_i)}$$

(27)

4

Maximum absolute error, $MAE$

$$MAE=\frac{\max (\left|y_i-{\widehat{y}}_i\right|)}{\max (y_i)-\min (y_i)}$$

(28)

In which, $RMSE$, $R^2$, and $AAE$ are the global prediction accuracy evaluation indexes, and $MAE$ is the local accuracy evaluation index. When $RMSE$ and $AAE$ are closer to 0 and $R^2$ is closer to 1, the global prediction accuracy of the model is higher, and when $MAE$ is closer to 0, the local prediction accuracy of the model is higher.

Based on the sampling point data obtained by OPLHS, the RSM, Kriging and RBF models for ${\epsilon }_{v,\max }^{\mathrm{I}\mathrm{I}\mathrm{I},p}$ were cross-validated five times, and the results are shown in

Table 4. Validation results and partial Pareto optimal solutions based on NSGA-II algorithm.

	RSM				Kriging				RBF
	RMSE	R2	AAE	MAE	RMSE	R2	AAE	MAE	RMSE	R2	AAE	MAE
Test 1	0.012	0.999	0.010	0.021	0.027	0.991	0.020	0.068	0.005	1.000	0.004	0.015
Test 2	0.013	0.998	0.010	0.034	0.016	0.995	0.011	0.054	0.006	1.000	0.005	0.011
Test 3	0.006	1.000	0.005	0.013	0.037	0.987	0.024	0.083	0.004	1.000	0.003	0.008
Test 4	0.014	0.998	0.012	0.027	0.022	0.992	0.016	0.051	0.007	0.999	0.005	0.021
Test 5	0.009	0.999	0.007	0.024	0.022	0.994	0.017	0.042	0.008	0.999	0.005	0.026
Average value	0.011	0.998	0.009	0.024	0.025	0.992	0.018	0.059	0.006	1.000	0.004	0.016

.

It was observed that the applicability of different surrogate models was ranked as follows: RBF model > RSM model > Kriging model, and the $R^2$ of the RBF model was always higher than 0.998. It was concluded that the prediction accuracy of the RBF model established based on OPLHS met the requirements, and the subsequent optimization of the SRM grain could be carried out after the surrogate model for ${\epsilon }_{v,\max }^{\mathrm{I},p}$ was established by the same method.

5. Result of Multiobjective Optimization for SRM Grain

After obtaining the required surrogate model, combined with the computational models of $A_b$ and ${\eta }_v$, the optimization model shown in Equation (1) is solved based on the NSGA-II algorithm according to the flowchart shown in Figure 14.

NSGA-II is one of the most widely used multiobjective optimization algorithms. The prototype was NSGA, proposed by Srimivas and Deb in 1994; and then K. Deb and S. Agrawal improved NSGA based on NSGA in 2001 and proposed the second-generation of NSGA with elite strategy, namely NSGA-II.

Compared with NSGA, the major difference of NSGA-II is the introduction of crowding degree ranking, which makes the distribution of individuals in each generation of the population more uniform, combined with the Pareto nondominance ranking to determine the fitness of individuals in the population.

Due to the high number of design variables, the number of populations was set to 200 and the number of iterations was set to 50 in order to ensure the diversity of populations. The distribution of the obtained Pareto optimal solution in the space of objective functions was shown in Figure 15.

From Figure 15, it can be obtained that the optimization method based on the NSGA-II algorithm is able to obtain the complete Pareto front for this problem. In addition, from the distribution of Pareto optimal solutions in the two-dimensional plane of the objective function, it can be seen that there is an obvious constraint relationship between the ${\epsilon }_{v,\max }$ and ${\eta }_v$, and the multiobjective optimization design method proposed in this paper can effectively balance the contradiction.

To verify the accuracy of the multiobjective optimization results of the drug shape based on the surrogate model, the numerical analysis of some of the Pareto optimal solutions was performed, and the comparison results are shown in

Table 5. Partial Pareto optimal solutions and validation results of multiobjective optimization for SRM grain.

	${\mathit{f}}_{{\mathit{\epsilon }}_{\mathit{v},\mathbf{max}}}(\mathit{X})/\mathit{\%}$	${\mathit{f}}_{{\mathit{\eta }}_{\mathit{v}}}(\mathit{X})/\mathit{\%}$	${\mathit{f}}_{{\mathit{A}}_{\mathit{b}}}(\mathit{X})/\mathit{\%}$	${\mathit{\eta }}_{\mathit{v}}/\mathit{\%}$	${\mathit{\epsilon }}_{\mathit{v},\mathbf{max}}^{\mathit{t}+\mathit{p}}/\mathit{\%}$	${\overline{\mathit{\epsilon }}}_{\mathit{v},\mathbf{max}}^{\mathit{t}+\mathit{p}}/\mathit{\%}$	$\frac{({\mathit{\epsilon }}_{\mathit{v},\mathbf{max}}^{}-{\overline{\mathit{\epsilon }}}_{\mathit{v},\mathbf{max}}^{})}{{\mathit{\epsilon }}_{\mathit{v},\mathbf{max}}^{\mathit{t}+\mathit{p}}}/\mathit{\%}$
Original	0.00	0.00	0.00	84.84	35.36	35.24	0.33
Partial Pareto optimal solutions	2.11	0.17	−26.72	84.98	26.27	25.91	1.36
	1.26	0.13	−26.59	84.95	26.29	25.96	1.26
	3.24	0.40	−26.36	85.18	26.64	26.04	2.26
	4.98	1.42	−12.58	86.04	31.32	30.91	1.30
	0.00	1.05	−1.91	85.73	33.65	34.69	−3.10
	0.02	0.22	−26.24	85.03	26.79	26.08	2.62
	0.06	0.88	−17.68	85.59	29.00	29.11	−0.40
	0.16	0.23	−25.94	85.03	26.73	26.19	2.03
	0.25	1.22	−13.11	85.87	31.33	30.73	1.93
	0.36	0.84	−20.17	85.55	28.92	28.23	2.40

, where ${\epsilon }_{v,\max }^{}$ is the actual calculated value and ${\overline{\epsilon }}_{v,\max }^{}$ is the maximum von Mises strain value of the SRM grain predicted by the surrogate model.

The validation results show that the maximum prediction deviation of the surrogate model for ${\epsilon }_{v,\max }^{t+p}$ was 3.10% and the average prediction deviation was 1.87%, which confirm the effectiveness and accuracy of the multiobjective optimization method of SRM grain based on the surrogate model.

6. Conclusions

In this paper, a multiobjective optimization design of grain is achieved for the SRM with finocyl grain based on parametric modeling and surrogate modeling methods, and the following conclusions are obtained:

(1) The multiobjective optimization of the SRM grain considering the structural integrity, internal ballistic performance, and loading performance was realized, and the complete Pareto optimal solution set was obtained based on NSGA-II with significant optimization effect. Under the constraints, the maximum von Mises strain can be reduced by up to 26.72% and the volume loading fraction can be increased by up to 1.83% compared with that before the optimization, and a series of balanced Pareto optimal solutions were also obtained for the decision maker to choose.

(2) The optimal design method based on parametric modeling and surrogate modeling techniques is significantly superior in efficiency, which can reduce the single calculation time by more than 95% and the number of calculations from 20,000 to 400, with an average prediction deviation of only 1.87%.

(3) Load and geometric parameters were analyzed for the SRM grain under combined temperature and pressure loads, and it was found that the effects of temperature and pressure loads on the von Mises strain at the danger point are superimposable and the ratio between them is fixed; at the same time, the only key geometric parameter affecting the von Mises strain at the danger point is $R,L,\theta$ or $R^{\prime },L^{\prime },{\theta }^{\prime }$. The resulting conclusions can significantly reduce the work required to build the surrogate model.