A novel aerodynamic modeling method for an axisymmetric missile with tiny units

Abstract

This article investigates the aerodynamic modeling problem of an axisymmetric missile with tiny units. Three aerodynamic force models for an axisymmetric air-to-air standard model missile are established. The mathematical model of trigonometric series model, response surface model, and kriging model for roll moment coefficients are developed. An error reduction approach is presented to determine the configuration of the proposed aerodynamic force models. In this approach, the minimum residual sum of the square criterion is acquired to determine the desired sample data point. Moreover, the modeling accuracy of the three models is further discussed. Finally, a model missile with four control surfaces and a cable cover was applied to test the proposed modeling approach in a transonic wind tunnel named as CARDC FL-24.

Introduction

To achieve great maneuvering capability of larger flight envelope, modern air-to-air flight condition is becoming more and more complex than ever. The large attack angle and the aero-surface deflection lead to serious unsteady flows and nonlinear interference. ¹ Therefore, a number of volumes of aerodynamic data should be provided in the flight control system design, the analysis of flight dynamics, and simulation. Taking missile as an example, its aerodynamic force is affected by the multivariate flight dynamic parameters consisting of the attitude angles, its flight Mach number, altitude, and actuator reflections ² . This lets the establishment of a high-dimensional aerodynamic database, which completely satisfies the demands of control system design and evaluation ³ by wind tunnel (WT) ⁴ and computational fluid dynamics (CFD), ⁵ both being costly and time consuming for missiles.^6,7 Motivated by avoiding this challenge, an appropriate approach is to establishing aerodynamic models using a certain amount of WT tests or using the CFD software. Then, the established models are applied online to predict the aerodynamic force characteristic of the missile at any flight time. Therefore, the problem of establishing aerodynamic model for missiles has received considerable attention.

Although a number of methodologies are available for the mathematical modeling of aerodynamics, most of them were developed based on the linearized aerodynamic force coefficients^8,9 (or their time derivatives). In fact, the dynamics of the modern air-to-air missiles are with strong nonlinearity. These nonlinear characteristic results in that, the modeling methods based on the linearized force coefficients are not suitable to establish a precise aerodynamic model for the modern air-to-air missiles. This problem promotes the development in the modern design of experiments (MDOE).^10,11 This approach can significantly reduce the testing numbers of WT and CFD and shorten the missile research cycle time. That is because a proper mathematic model which can represent the missile’s aerodynamic characteristic well is applied to predict the aero data and to establish the aerodynamic database. ¹²

In MDOE, the missile aerodynamic model is established on the basis of the experimental statistics and the motion mechanism of the missile. The existing missile aerodynamic models can be classified into two categories. ¹³ The first category is the one in which the aerodynamic force and the flight states are mathematically formulated according to the physical phenomenon and the mechanism. The algebraic model, ¹⁴ the series model,^15,16 the state space model, ¹⁷ the step response model, ¹⁸ and so on are the most available results for this category. The other category is knowledge based. The aerodynamic modeling is viewed as a “black box” problem with the complex physical mechanism eliminated. Then, the new cross-disciplinary knowledge is applied to mathematically establish the relation between the aerodynamic force and the flight states. For example, the neural network-based aerodynamic model was developed in Rokhsaz and Steck’s study. ¹⁹ The approximation capability of neural networks to unknown nonlinear functions was applied. In Kong and Wang ²⁰ and Wang et al., ²¹ the aerodynamic model was established using the fuzzy logic theory. In Chen, ²² another model developed using the support vector machine (SVM) was presented for aerodynamics of missiles.

It should be stressed that each aerodynamic model established in the existing literature has its own advantages and disadvantages. Moreover, in general, the model’s flexibility depends on the missile’s aero configuration. Regardless of any category of model employed, however, the following four issues should be addressed during the aerodynamic force modeling:

1. How the configuration of the model is decided?

2. How the unknown coefficients are handled?

3. How the sample data points (or design experiments) are selected?

4. How the accuracy of model prediction is analyzed?

Motivated by solving the above four challenges, a novel and practical aerodynamic modeling approach is presented for an axisymmetric air-to-air missile. In this approach, three aerodynamic force models including the response surface model, the trigonometric series model, and the kriging model are established for the missile. The main contributions of the article are as follows:

1. Comparing with the existing results, the selection of modeling sample data points and the optimal required sample points are both closely related with the aerodynamic mathematic model.

2. This article also shows that the trigonometric series model is more precise with limited sample points size than the response surface and the kriging models, because it reveals the characteristic of the aero configuration of axisymmetric missile with tiny units. Hence, the related results can be used to design wind tunnel test and CFD simulation, which would save lots of time and costs.

3. An error reduction method is used to determine the configuration of the aerodynamic force models, and the minimum residual sum of square criterion is acquired to determine the desired sample data points. The rest of the article is organized as follows: a brief description of an air-to-air missile and aerodynamic force data samples is introduced in Section 2; the mathematical equations of three aerodynamic force models are presented in Section 3; the model configuration determination and the method to handle unknown coefficients are presented in Section 4; the method for modeling data samples selection is shown in Section 5; some conclusions are given in Section 6.

Missile configuration and data samples

The body configuration of most air-to-air missiles is axisymmetric. However, the overall shape is asymmetric due to some protrusive units such as cable cover, antenna, hoop, and ramjet inlet. Figure 1 shows an air-to-air missile standard model. It is tested in FL-24 wind tunnel at China Aerodynamic Research and Develop Center (CARDC). This missile model is also chosen as a plant whose aerodynamic force will be mathematically modeled in this work. For missiles, the longitudinal size of those protrusions along the missile body may be close to the length of the missile. The protrusion height is far less than the missile’s diameter. Although those tiny units may have no substantial effect on the actuator efficiency, they have a significant influence on the aerodynamic force and even destroy the axisymmetric aerodynamic force functional form. These two issues let the modeling problem be very difficult.

Figure 1.

The configuration of an air-to-air missile standard model.

The missile standard model consists of a cone nose and a cylinder body with four cruciform canards and four tail actuators. A thin cable cover is also mounted on the body. For convenience, the roll moment coefficients are chosen for mathematical modeling. The missile’s roll moment coefficients acquired in WT tests are partly shown in Figure 2. The rolling angle starts at −90° and ends at 90°, with the intervals being 5.625°. Hence, each curve that corresponds to a total angle of attack has 33 data points.

Figure 2.

The roll moment coefficients of the missile standard model.

The aerodynamic modeling approach

In this section, a novel and practical aerodynamic modeling approach is presented for an axisymmetric air-to-air missile. Three aerodynamic force models, that is, the response surface model, the trigonometric series model, and the kriging model are established for the missile.

The development of the response surface model

The mathematical equation of the response surface model is the function of an output response variable F and a set of input variables $x_1,x_2,\dots ,x_k$ . This function is specified as

$$F=f(x_1,x_2,\dots ,x_k)$$

(1)

Generally, this function F, in equation (1), can be approximated by polynomial. Especially, when the function F denotes a curve or a surface, its mathematical model can be written as

$$F={\beta }_0+\sum _{i=1}^k{\beta }_i{x}_i+\sum _{i=1}^k{\beta }_{1i}{x_i}^2+\sum \sum {\beta }_{1j}{x}_i{x}_j$$

(2)

where ${\beta }_0$ , ${\beta }_i$ , ${\beta }_{1i}$ , and ${\beta }_{\mathrm{ij}}$ are coefficients.

Because the roll moment coefficient mainly depends on the roll angle, the missile’s response surface model of roll moment coefficient can be established using the following equation

$$m_x=\sum _{i=0}^n{m}_{\mathrm{xi}}{\varphi }_n^i$$

(3)

where $m_{\mathrm{xi}}$ is regression coefficient and ${\varphi }_n$ is the air flow roll angle.

In practical engineering, the response surface model (equation (3)) cannot use a much higher-order polynomial due to its ill-conditioned fitting characteristics.

The establishment of the trigonometric series model

The trigonometric series model is different from the response surface model (equation (3)). It is closely related to the missile aero configuration. In the normal condition, the aerodynamic force coefficient function is decided by the velocity vector, the control surface deflection, and the angular velocity relative to the center of gravity. Hence, neglecting the aero damp effect, the following equation can be established as the aerodynamic force model

$$F_{\mathrm{dc}}=F({\alpha }_n,{\varphi }_n,\delta ,{\varphi }_{\delta })$$

(4)

where ${\alpha }_n$ is the total angle of attack, $\delta$ is the actuator deflection, and ${\varphi }_{\delta }$ is the actuator’s installation angle related to body.

According to the space geometry relation of the four control surfaces illustrated in Figure 1, the orthographic equivalent deflection angles ${\delta }_x$ , ${\delta }_y$ , and ${\delta }_z$ of each control surface can be computed by

$${\delta }_x=\delta$$

(5)

$${\delta }_z=\delta \cos ({\varphi }_{\delta })$$

(6)

$${\delta }_y=\delta \sin ({\varphi }_{\delta })$$

(7)

Once ${\alpha }_n$ and $\delta$ are determined, the aerodynamic force model in equation (7) can be converted to a particular kind of mathematical function of ${\varphi }_n$ and ${\varphi }_{\delta }$ . Due to the axisymmetric configuration of the missile, the aerodynamic force is periodic with a period angle $2\pi /n$ . It is seen in Figure 2 that n is the number of control surfaces displaced azimuthally at the equal angle intervals. In this case, the model of aerodynamic force coefficient can be expressed as the following double trigonometric series

$$F_{\mathrm{dc}}=\sum _{i,j=-\infty }^{\infty }(a_{i,j}\sin (p)+b_{i,j}\cos (p))$$

(8)

where $p=i{\varphi }_n+j{\varphi }_{\delta }$ .

In practice, however, the missile configuration is not perfectly symmetric because of various protrusions such as cable cover, hook, and antenna. Although these units are much smaller than the missile body, they indeed have an influence on the aerodynamic force. More specifically, they may destroy the periodicity. Therefore, equation (8) must be improved to eliminate the effect of the asymmetrical units. In this article, the practical correcting method reported in Tang et al. ¹² is applied to achieve that improvement. Moreover, the following equation shows the revised roll moment coefficient trigonometric series model of the missile

$$m_x=\sum _{i=0}^{\infty }{m}_{\mathrm{xi}}\sin (i(n{\varphi }_n))+\sum _{j=0}^{\infty }{m}_{\mathrm{xj}}\sin (j({\varphi }_n-{\phi }_r))$$

(9)

where ${\phi }_r$ is the installation azimuthal angle of cable cover in Figure 1.

The design of the kriging model

The kriging model contains a regression model F and a random function (stochastic process) ²³

$${\hat{y}}_l(x)=F({\beta }_{:,l},x)+z_l(x),l=1,2,\dots ,q$$

(10)

where ${\beta }_{:,l}$ are the coefficients ${\beta }_0$ , ${\beta }_i$ , ${\beta }_{1i}$ , and ${\beta }_{\mathrm{ij}}$ .

In this work, a regression model with a linear combination of p chosen functions $f_j$ can be used for equation (10). Then, ${\hat{y}}_l(x)$ can be rewritten as

$${\hat{y}}_l(x)={\beta }_{1,l}{f}_1(x)+\cdots +{\beta }_{p,l}{f}_p(x)$$

(11)

where the coefficients ${\beta }_{k,l}$ are regression parameters. The random process z is assumed to have zero mean, and its covariance between $z_l(\omega )$ and $z_l(x)$ is given by

$$\mathrm{Cov}[z_l(\omega ),z_l(x)]={\sigma }_l^{2}R(\theta ,\omega ,x),l=1,2,\dots ,q$$

(12)

where ${\sigma }_l^2$ is the process variance for the lth component of the response. The term $R(\theta ,\omega ,x)$ is the correlation model with parameter $\theta$ .

For a set m of sample data sites, $S=[\begin{array}{cccc}{s}_1\hfill & s_2\hfill & \cdots \hfill & s_m\hfill \end{array}{]}^{\mathrm{T}}$ and the responses $y=[\begin{array}{cccc}{y}_1\hfill & y_2\hfill & \cdots \hfill & y_m\hfill \end{array}{]}^{\mathrm{T}}$ , at an untried point x. Let $r(x)$ be defined as

$$r(x)=[\begin{array}{cccc}R(\theta ,s_1,x)\hfill & R(\theta ,s_2,x)\hfill & \cdots \hfill & R(\theta ,s_m,x)\hfill \end{array}{]}^{\mathrm{T}}$$

(13)

The predictor equation at x is given by

$$\hat{y}(x)=f^{\mathrm{T}}{\beta }^{*}+r^{\mathrm{T}}(x){\gamma }^{*}$$

(14)

Then, for the regression problem, it has

$$F\beta =[\begin{array}{cccc}f(s_1)\hfill & f(s_2)\hfill & \cdots \hfill & f(s_m)\hfill \end{array}{]}^{\mathrm{T}}\beta =Y$$

(15)

To this end, the generalized least squares (with respect to R) can be obtained from equation (15) as

$${\beta }^{*}=(F^{\mathrm{T}}{R}^{-1}F{)}^{-1}{F}^{\mathrm{T}}{R}^{-1}Y$$

(16)

In accordance, using equation (16), ${\gamma }^{*}$ can be computed via the following residual

$${\gamma }^{*}=Y-F{\beta }^{*}$$

(17)

More specifically, the form of the correlation model $R(\theta ,\omega ,x)$ in equation (12) usually has seven choices. ²⁴ These choices are the EXP, the EXPG, the GAUSS, the LIN, the SPHERICAL, the CUBIC, and the SPLINE. Conversely, the aerodynamic force for a missile tends to show a parabolic behavior near the basic point. This means that the Gaussian, the cubic, and the spline should be chosen. Therefore, let a Gaussian process be assumed to be such that

$$R(\theta ,\omega ,x)=\sum _{j=1}^n\exp (-{\theta }_j{({\omega }_j-x_j)}^2)$$

(18)

Then, solving the correlation function yields the optimal coefficient ${\theta }^{*}$ having the form

$$R(\theta ,\omega ,x)=\underset{\theta }{max}\{\varphi (\theta )=|R{|}^{1/m}{\sigma }^2\}$$

(19)

where $|R|$ is the determinant of $R(\theta ,\omega ,x)$ .

The selection of the aerodynamic force model configuration

For the response surface model and the trigonometric series model, the roll moment coefficients’ mathematical equations about the air flow roll angle contain infinite terms. This infinity property allows the equations not to be applied directly. However, once the roll moment coefficients model is selected, then the number of terms should be determined. Therefore, this work will be done in this section.

The selection of the aerodynamic force model configuration can be done according to two steps. The first step is to build candidate set. Then, the second step is to select the specific form of model by the stepwise regression analysis ²⁵ or the orthogonal least square algorithm ²⁶ with sample data.

Taking the response model as an example, the error reduction method is applied to determine the roll moment coefficients model of the missile. Suppose that equation (3) has n terms for the $k\mathrm{th}$ sample data, then the square of the real roll moment coefficient is

$$\overline{m_x^2(k)}=\sum _{i=1}^n{m}_{\mathrm{xi}}^2\overline{{\theta }_n^{2i}(k)}+\overline{{\phi }^2(k)}$$

(20)

where $\phi$ is the prediction error.

In the specific case of n = 0, the maximum prediction error is equal to $\overline{m_x^2(k)}$ . Whenever the term, for example, the $i\mathrm{th}$ term $m_{\mathrm{xi}}\overline{{\varphi }_n^i}$ , is added into the model, the square of prediction error will be reduced as $m_{\mathrm{xi}}^2\overline{{\theta }_n^{2i}(k)}$ . Hence, the ratio illustrated in equation (21) affords a useful standard. Following this standard, whether $m_{\mathrm{xi}}{\varphi }_n^i$ should be selected into the response model or not can be decided

$$[\mathrm{ERR}{]}_i=\frac{{\hat{m}}_{\mathrm{xi}}^2\overline{{\varphi }_n^{2i}(k)}}{m_x^2(k)},i=1,2,\dots ,n$$

(21)

For the sample data points shown in Figure 2, the candidate response model of the roll moment coefficients can be assumed as

$$\begin{array}{cc}\multicolumn{1}{c}{m_x=}& m_{x0}+m_{x1}{\varphi }_n+m_{x2}{\varphi }_n^2+m_{x3}{\varphi }_n^3\hfill \\ \multicolumn{1}{c}{}& +m_{x4}{\varphi }_n^5+m_{x5}{\varphi }_n^7+m_{x6}{\varphi }_n^9\hfill \end{array}$$

(22)

where the coefficients $m_{x0}~m_{x6}$ can be calculated using the least square algorithm. Moreover, the relative error for each term of equation (22) is listed in Table 1.

Table 1.

The values of the relative error in the response surface model.

Model terms	The value of the relative errors
$m_{x0}$	0.000011284218419
$m_{x1}{\varphi }_n$	0.010988767419702
$m_{x2}{\varphi }_n^2$	0.000001007544944
$m_{x3}{\varphi }_n^3$	0.237446922247964
$m_{x4}{\varphi }_n^5$	0.540477570729898
$m_{x5}{\varphi }_n^7$	0.202095810013662
$m_{x6}{\varphi }_n^9$	0.008978637825410

According to the above listed relative error values, the third term $m_{x2}{\varphi }_n^2$ has the minimum impact on the response model, and then the first term $m_{x0}$ . If the allowable modeling accuracy is given, the configuration of the response model can be decided by deleting the minimum impact term alternately. The configuration of the trigonometric model can also be decided by the same procedure.

The requirements of the sample data points

The sample data points are usually acquired by the WT tests or the CFD software scheduled via the design of experiments. According to the MDOE, the sample data points are required to satisfy two factors: (1) a sufficient number of data to decide model’s configuration and (2) an extra number of data to evaluate modeling accuracy. Hence, the six-term trigonometric series model given in equation (23), the nine-order response model given in equation (24), and the kriging model with the Gaussian correlation model are used to analyze the effect of the sample data points selection on the aerodynamic force model

$$\begin{array}{c}{m}_x=m_{x0}+m_{x1}\sin (4{\varphi }_n)+m_{x2}\sin (8{\varphi }_n)\hfill \\ +m_{x3}\sin ({\varphi }_n-135^{°})+m_{x4}\sin (3({\varphi }_n-135^{°}))\hfill \\ +m_{x5}\sin (5({\varphi }_n-{135}^{°}))\hfill \end{array}$$

(23)

$$\begin{array}{cc}\multicolumn{1}{c}{m_x=}& m_{x0}+m_{x1}{\varphi }_n+m_{x3}{\varphi }_n^3\hfill \\ \multicolumn{1}{c}{}& +m_{x4}{\varphi }_n^5+m_{x5}{\varphi }_n^7+m_{x6}{\varphi }_n^9\hfill \end{array}$$

(24)

In Tang et al., ¹² a saturated D-best criterion was used to determine the sample data points. In this article, the minimum residual sum of square criterion in Tang et al. ¹² is adopted to eliminate the low modeling information value points in sequence. The residual sum of squares is defined as

$$S{S}_E=\sum _{i=1}^n{(y_i-{\hat{y}}_i)}^2$$

(25)

where $y_i$ and ${\hat{y}}_i$ are the actual value and model predicting value of the $i\mathrm{th}$ verification data point. The results of the eliminated sample data points shown in Figure 2 for three roll moment coefficients model are listed in Tables 2–4, respectively, where the point having the lowest modeling information is eliminated first.

Table 2.

A series of eliminated data points for the trigonometric series model.

No.	${\varphi }_n(°)$	No.	${\varphi }_n(°)$	No.	${\varphi }_n(°)$
1	11.25	10	28.125	19	−11.25
2	67.5	11	78.75	20	−56.25
3	−67.5	12	−50.625	21	−73.125
4	−61.875	13	−45	22	5.625
5	−22.5	14	−28.125	23	−16.875
6	16.875	15	84.375	24	50.625
7	0	16	90	25	45
8	−29.375	17	39.375	26	−84.375
9	33.75	18	61.875	27	−5.625

Table 3.

A series of eliminated data points for the response surface model.

No.	${\varphi }_n(°)$	No.	${\varphi }_n(°)$	No.	${\varphi }_n(°)$
1	56.25	10	61.875	19	−84.375
2	84.375	11	73.125	20	33.75
3	−67.5	12	−22.5	21	−28.125
4	−11.255	13	−11.25	22	78.75
5	16.875	14	50.625	23	67.5
6	−45	15	0	24	−39.375
7	−16.875	16	−50.625	25	−56.25
8	28.125	17	5.625	26	−90
9	−61.875	18	39.375	27	−5.625

Table 4.

A series of eliminated data points for the kriging model.

No.	${\varphi }_n(°)$	No.	${\varphi }_n(°)$	No.	${\varphi }_n(°)$
1	−28.125	10	39.375	19	−50.625
2	61.875	11	−90	20	−5.625
3	67.5	12	−61.875	21	5.625
4	22.5	13	−73.125	22	−22.5
5	33.75	14	−67.5	23	−78.75
6	−33.75	15	84.375	24	56.25
7	11.25	16	78.75	25	−84.375
8	45	17	−16.875	26	0
9	50.625	18	16.875	27	90

As listed in Tables 2–4, the selection of the sample data points is closely related to the modeling methods. Moreover, the modeling information contained in each data point is different. The information heavily depends on the different mathematic models. Figures 3–6 show the prediction results of the roll moment coefficients obtained from the proposed three aerodynamic models and the experimental testing with 32, 20, 10, and 6 sample data points.

Figure 3.

The prediction results with 32 sample data points.

Figure 4.

The prediction results with 20 sample data points.

Figure 5.

The prediction results with 10 sample data points.

Figure 6.

The prediction results with 6 sample data points.

Because the kriging model has no modeling error at the sample data points, its prediction results are more accurate than the results provided by the other two models, when all 32 points are used for modeling. However, less sample data points will decrease the prediction accuracy of all the presented three models. If the least necessary modeling points, that is, 6 sample data points, are chosen, then it is observed that the kriging model and the response surface model will have a definite offset from the WT test data, as shown in Figure 7. More specifically, it can be seen in Figure 7 that when the reduced number of sample data points is more than 20, the prediction relative errors of the kriging and the response surface model are increased significantly. In Figure 8, the optimal 6 sample data points are given for the three roll moment coefficients models.

Figure 7.

The prediction relative errors of different aerodynamic models.

Figure 8.

The optimal 6 sample points for different aerodynamic models.

Conclusion

The investigation of this article concludes that the selection of the modeling sample data points is closely related with the aerodynamic mathematic model, and the optimal required sample points are different for different models. This means that the aerodynamic model is the main influence factor in MDOE. Suitable aerodynamic model can greatly improve the effectiveness of the WT tests and CFD. Comparing with the three models above, the kriging model has the most accurate modeling results when the sample is large enough. The trigonometric series model is more precise when the sample is small. The trigonometric series model has significant advantages in the modeling precision and in the modeling sample size because it reveals the characteristic of the aero configuration of axisymmetric missile with tiny units. The aero configuration of the missile is another important factor for aerodynamic modeling, especially for approximately plane symmetric missile. Solving problems is one of the future works. Moreover, based on the mathematical model established and validated in this article, the development of nonlinear controller for the missile considered with high performance ensured should also be done in the future. More specifically, this could be achieved by applying the advanced control approaches in literature.^27–31