Deriving a Control-Oriented Model for an Axisymmetric Vehicle With the Power-Law Revolution Nose

ABSTRACT:

The purpose of this paper is to construct a new general vehicle model as an open fundamental material for the guidance and control research. In this study, parameterized configuration, aerodynamics calculation, control-oriented modeling, stability analysis, and nominal trajectory design are performed for the general vehicle model. First, the aerodynamic configuration is parameterized as an axisymmetric body with a power-law revolution nose. Then, an engineering method considering inviscid flow, base drag and skin friction is used for the aerodynamics calculation, and a control-oriented fitting model of longitudinal aerodynamics is established based on the analysis of the correlation between aerodynamic force and the parameters of Mach number, attack angle, elevator deflection and height. Next, the aerothermodynamic environment prediction of power-law revolution axisymmetric hypersonic vehicle (PRAHV) is discussed, and the nose heating rate formula of PRAHV is established. The stability analysis and nominal trajectory design of PRAHV is performed based on the fitting model and the heating rate formula. The stability analysis shows that both the static stability and dynamic stability of the vehicle are unstable. The nominal trajectory of unpowered longitudinal maneuvering is achieved by the hp-adaptive pseudospectral method, which demonstrated that the availability of the control-oriented model established in this paper. In conclusion, this work provides a fundamental object for further study of vehicle guidance, control, and evaluation.

KEYWORDS:
Kinzhal-like system; Aerodynamic engineering calculation; Control-oriented modeling; Model stability analysis; Nominal trajectory

INTRODUCTION

In 2018, Russia publicly disclosed the Kinzhal (Liao 2018a) air-launched hypersonic missile system. In the same year, the United States media disclosed relevant research projects (Liao 2018b), including advanced hypersonic weapon (AHW), hypersonic conventional strike weapon (HCSW), air-launched rapid response weapon (ARRW), hypersonic aspirated weapon concept (HAWC), long-range hypersonic weapon (LRHW), etc. According to the published pictures and video data, the axisymmetric body with a power-law revolution nose (leading edge) has become one of the typical aerodynamic configurations of hypersonic vehicles. The aerodynamics of the power-law revolution elaborated in the book by Grozovsky et al. (1978) has presented that the aerodynamic configuration with power-law revolution is the optimal shape with the minimum drag and lower heat transfer coefficient. For the research convenience of the hypersonic vehicle of the axisymmetric body with a nose of power-law revolution, the vehicle is named as power-law revolution axisymmetric hypersonic vehicle (PRAHV) in this paper.

According to the reports of the Russian Kinzhal system (Liao 2018a; Military-Today 2019), the vehicle was air-launched and had an initial supersonic speed before the engine started. By means of the rocket propulsion mode, the solid rocket can help the vehicle to obtain a hypersonic speed greater than Mach 5, as well as a large amount of kinetic energy and potential energy. When the fuel is exhausted, the kinetic and potential energy can be further utilized to achieve a larger glide range. However, the unpredictability of the trajectory makes it difficult to intercept due to the continuous maneuver of the glide trajectory. Furthermore, a PRAHV as an incoming target presents a challenge to the air defense system on trajectory estimation, tracking, and interception point of view. Since the establishment of the aerodynamic model of vehicle is the basis of trajectory estimation, control and guidance, the main purpose of this paper is to develop a reasonable PRAHV model as the premise of trajectory estimation and control algorithm research.

Two typical models are widely used for hypersonic control research. One is the generic hypersonic vehicle (GHV) model (Shaughnessy et al. 1990; Keshmiri et al. 2006), which is the first air-breathing hypersonic concept vehicle with winged-cone configuration. The other is the air-breathing hypersonic vehicle (AHV) model (Bolender and Doman 2007; Li et al. 2011a), which is a lift-body with aerodynamics/propulsion integration configuration. The above two models are situated in the research stage, whereas the PRAHV instanced as Kinzhal with a combat capability should have more practical research value for the defense system. The configuration of GHV, AHV, and PRAHV presented in this paper are shown in Fig. 1 for intuitive comparison.

Figure 1
Configuration of two typical models and PRAHV.

The establishment of aerodynamic database (Pamadi et al. 2001; Rufolo et al. 2006) through wind tunnel tests and numerical calculation is a good process for vehicle development, but it is a huge workload and time-consuming work, and it is not suitable for rapid quantitative analysis of models. Therefore, a quick and systematic work is performed around deriving a control-oriented model in this paper. Firstly, a parameterized configuration approach is proposed to develop a three-dimensional (3D) rigid-body profile for PRAHV. On the basis of constructing of 3D rigid-body profile of the Kinzhal-like vehicle, the longitudinal (three degrees of freedom, 3-DOF) aerodynamic data is obtained by using the engineering prediction method. Then, the aerodynamic fitting model is established based on the dimensional correlation of the result data, and the control characteristics of the fitting model are analyzed. From configuration construction to data modeling of PRAHV, this method can quickly and cheaply provide a research object for the research of control related algorithms, especially suitable for the quick judgment of new vehicle or the preliminary research work when lack scientific research funds.

The research significance of this study can be summarized in two aspects. On the one hand, this work can provide theoretical and data support for the research and judgment of the existing Kinzhal or other similar aircraft by establishing a set of complete data sets with lower cost of data acquisition. On the other hand, the open data basis is established for the guidance and control method research of this kind of model.

PARAMETERIZED CONFIGURATION

According to the public images and videos of the hypersonic vehicle, a PRAHV model is constructed by the geometric parameterization method, which is an aerodynamic configuration designed by the integration of rocket and airframe. The parameterization method is shown in Fig. 2. The shape is divided into four parts. Part I is the power-law curve rotating body, and the details of power-law revolution technology is referred to the book of Russian Grozovsky et al. (1978); Part II and III are the truncated vertebral body and the uniform section, respectively; Part IV is the circular section and the rudder.

Figure 2
Parameterized vehicle of slender body with power-law rotating.

Referring to the public data of the Kinzhal system and Iskander missile data, a PRAHV model instance can be obtained, as shown in Fig. 3. According to the public reports (Military-Today 2019; Wikipedia 2019a; Liao 2018a), the Kinzhal resembles an air-launched version of the Iskander short-range ballistic missile. This paper assumes that the Kinzhal has the same appearance as Iskander, and generate a PRAHV model of Kinzhal-like vehicle with the appearance of Iskander. The total length of the vehicle is 7.2 meters, and other parameters can be obtained by the pixel ratio of the photo. Fig.3 (a) is a high-definition photo downloaded from Wikipedia (Wikipedia 2019b), and an obtained 3D PRAHV model instance is shown in Fig. 3 (b). The detail parameters of the 3D PRAHV model instance are shown in Table 1.

Figure 3
A PRAHV model of Kinzhal-like vehicle. (a) Iskander measured in photoshop (pixels); (b) 3D model of a PRAHV same size with Iskander

Table 1
The detail parameters of the 3D PRAHV model instance.

Objects	L 1	L 2	L 3	Rn	y1	y2
Pixels in photo of Fig. 3(a)	91.4	368.6	743	~	50.15	101.2
Meters in 3D model of Fig. 3(b)	0.8856	3.5712	7.2	0.0136	0.486	0.9806

AERODYNAMICS CALCULATION

According to the characteristics of engineering methods analyzed by Huang (1995) and Li (2012), a comprehensive method combining inviscid flow, base drag, and skin friction was adopted for aerodynamics calculation of PRAHV. A similar calculation for quasi-waverider vehicle has been proposed by Che and Tang (2007). The impact angles between the incoming flow of each windward component and surface are quite different for hypersonic vehicles under different flight states. Therefore, the aerodynamic characteristics of the PRAHV as a whole cannot be simply described by the Newton method or other single calculation methods. Let the pressure coefficient be C_P=0 for the flow field of leeward surface. Dahlem-Buck method is used for the windside surface flow field. The Dahlem-Buck method (Li and Gao 2014) is a simple combination of Newtonian and tangent wedge/cone methods. The pressure coefficient is calculated by the Newtonian method when the impact angle is greater than 22.5°; otherwise, it shall be computed by the tangent wedge/cone method, as in Eq. 1.

(1) ${\mathit{C}}_{\mathit{P}}=\left\{\begin{array}{cc}\left[\frac{1.0}{\sin {\left(4\mathit{\delta }\right)}^{\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$4$}\right.}}+1\right]{\sin }^2\mathit{\delta }& \mathit{\delta }\le 22.5°\\ \mathit{K}{\sin }^2\mathit{\delta }& \mathit{\delta }>22.5°\end{array}\right.$

where δ is the impact angle, and the surface is windside when δ > 0. K is determined by a given Mach number and is derived from the positive shock wave relation and isentropic condition (Moore and Williams 1989), which is expressed in Eq. 2.

(2) $\mathit{K}=\frac{2}{\mathit{\gamma }{\mathit{M}}_{\infty }^2}\left[{\left(\frac{\mathit{\gamma }+1}{2\mathit{\gamma }{\mathit{M}}_{\infty }^2-\mathit{\gamma }+1}\right)}^{\frac{1}{\mathit{\gamma }-1}}{\left(\frac{\mathit{\gamma }+1}{2}{\mathit{M}}_{\infty }^2\right)}^{\frac{\mathit{\gamma }}{\mathit{\gamma }-1}}-1\right]$

The aerodynamic resultant force of the upwind surface flow field can be obtained by integrating the pressure coefficient of the surface element. The normal force coefficient, axial force coefficient and pitch moment coefficient with the reference point of (0,0,0) can be obtained as shown in Eq. 3.

(3) $\begin{array}{c}{\mathit{C}}_{\mathit{N}}=\left(1/{\mathit{S}}_{\mathit{ref}}\right)\sum {\mathit{C}}_{\mathit{p}}·{\mathit{n}}_{\mathit{z}}·d\mathit{s}\hfill \\ {\mathit{C}}_{\mathit{A},\mathit{w}}=\left(1/{\mathit{S}}_{\mathit{ref}}\right)\sum {\mathit{C}}_{\mathit{p}}·{\mathit{n}}_{\mathit{x}}·d\mathit{s}\hfill \\ {\mathit{C}}_{\mathit{m},0}=\left(1/{\mathit{L}}_{\mathit{ref}}\right)\sum \left({\mathit{C}}_{\mathit{p}}·{\mathit{n}}_{\mathit{z}}·d\mathit{s}·\mathit{z}-{\mathit{C}}_{\mathit{p}}·{\mathit{n}}_{\mathit{x}}·d\mathit{s}·\mathit{x}\right)\hfill \end{array}$

where S_ref is the reference area, and L_ref is the reference length. The ds is the area of a triangle surface. The n_X and n_Z are x and z components of the normal vector N, respectively.

The bottom of a vehicle flying at hypersonic speed can be generally considered as a vacuum. However, due to the influence of real gas viscosity and other factors, the bottom will not be a complete vacuum and there exists partial pressure. Considering the influence of Mach number on the bottom vacuum, a modified formula (Eq. 4) is used to calculate base drag.

(4) ${\mathit{C}}_{\mathit{P}0}=\left\{\begin{array}{cc}\frac{0.771{\left(1.0-0.011{\mathit{M}}_{\infty }^2\right)}^{3.5}-1.43}{{\mathit{M}}_{\infty }^2}& {\mathit{M}}_{\infty }\le 9.5\\ -\frac{1}{\frac{\mathit{\gamma }}{2}{\mathit{M}}_{\infty }^2}& {\mathit{M}}_{\infty }>9.5\end{array}\right.$

The formula for calculating turbulent friction is derived from the formula for calculating turbulent friction on a flat plate (Li 2012). Turbulent friction c_f is shown in Eq. 5:

(5) ${\mathit{c}}_{\mathit{f}}=\left\{\frac{{\left({\mathit{x}}^{\prime }\right)}^2·{\mathit{c}}_{\mathit{fi}}/\mathit{z},\mathit{Re}>1.0\times {10}^6}{1.328/\sqrt{\mathit{Re}},\mathit{Re}\le 1.0\times {10}^6}\right.$

where, $x^{\prime }=\frac{\mathit{x}}{{\log }_{10}\mathit{r}*-1.5}$, $\mathit{r}*=\frac{\mathit{Re}}{\mathit{den}·{\mathit{z}}^{2.5}}$, $\mathit{den}=\frac{\mathit{T}+S\mathit{h}}{\mathit{Z}·\mathit{T}+S\mathit{h}}$, Sh=110.4, z=1.0+0.115·M²_∞, c_fi=0.088/_x². T and Re represent the temperature and Reynolds of the incoming flow, respectively, both of which can be calculated by the formula of atmospheric parameters.

Then, the friction coefficient C_Df can be obtained by the following formula (Eq. 6):

(6) ${\mathit{C}}_{\mathit{Df}}={\mathit{c}}_{\mathit{f}}·{\mathit{S}}_{\mathit{sin}\mathit{k}}/{\mathit{S}}_{\mathit{ref}}$

where S_sink is the infiltrating area and S_ref is the reference area.

Finally, the drag coefficient C_D, lift coefficient C_L and pitching moment coefficient C_m under the specific center of mass (x_cg, 0, z_cg) are calculated by Eq. 7. The center of mass is dimensionless parameters in this paper, and the unit and dimension for the center of mass electric field are expressed by the full length (L_em = 7.2m) of the vehicle.

(7) $\begin{array}{c}{\mathit{C}}_{\mathit{A}}={\mathit{C}}_{\mathit{A},\mathit{w}}+{\mathit{C}}_{\mathit{p}0}\hfill \\ {\mathit{C}}_{\mathit{D}}={\mathit{C}}_{\mathit{A}}·\mathbf{cos}\left(\mathit{\alpha }\right)+{\mathit{C}}_{\mathit{N}}·\mathbf{sin}\left(\mathit{\alpha }\right)+{\mathit{C}}_{\mathit{Df}}\hfill \\ {\mathit{C}}_{\mathit{L}}=-{\mathit{C}}_{\mathit{A}}·\mathbf{sin}\left(\mathit{\alpha }\right)+{\mathit{C}}_{\mathit{N}}·\mathbf{cos}\left(\mathit{\alpha }\right)\hfill \\ {\mathit{C}}_{\mathit{m}}={\mathit{C}}_{\mathit{m},0}-\left({\mathit{C}}_{\mathit{A}}·{\mathit{z}}_{\mathit{cg}}-{\mathit{C}}_{\mathit{N}}·{\mathit{x}}_{\mathit{cg}}\right)·{\mathit{L}}_{\mathit{en}}\hfill \end{array}$

Similar engineering calculation methods have been verified by Che and Tang (2007). Furthermore, to validate the above-mentioned method for the PRAHV, several symbolic computational fluid dynamics (CFD) calculations (solving the Reynolds Average N-S equations by k-ε model, which is reliable accuracy but relatively high-cost method) are carried out. Three grids of PRHAV for CFD computation were made corresponding to the rudder deflection degree of 0, 5, and 10 respectively. The aerodynamics is computed at flight altitude of 30 km, flight speed of M_∞ = 10,12, and attack angle of α = 0° ~ 8°. As shown in Fig. 4, the results show that the engineering method is basically consistent with CFD calculation results. The limitation of this method is that the errors become slightly bigger in lower Mach number, but fortunately all force trends are consistent and studies of control characteristics and methods should be effective. The low Mach errors may be caused by insufficient leeward consideration. It is feasible for the preliminary aerodynamics analysis of PRAHV because the PRAHV is considered as a hypersonic vehicle which could fly at more than Mach 10. The axis coordinates are same in the six charts, the horizontal axis is attack angle, and the vertical axis is values of lift coefficient and drag coefficient.

Figure 4
Lift and drag curve of PRAHV for contrast verification.

DATA SET ANALYSIS

Sufficient aerodynamic data sets are obtained for the control analysis through the above algorithms. The correlation analysis is carried out on the calculated result data. For the convenience of analysis, the definitions of variables in this section are declared as follows. First, the symbols and units are declared in the parentheses following the variable name, as follows: attack angle (α, degree), Mach number (M_∞, the local velocity of sound), elevator angle (δ_p, degree) and flight altitude (H, km). Second, some expressions like M_∞ = 3 and H = 40 km in the figure are abbreviated as ‘M3’ and ‘H40km’, respectively.

Figure 5 shows the relationship between the calculated lift coefficient, drag coefficient, pitching moment coefficient, attack angle, Mach number, elevator angle, and flight altitude. Since the center of mass of the model in this paper is within the range of 0.75L ~ 0.65L, the fixed moment reference point is set as (0.7L, 0, 0) to quantitatively analyze aerodynamic characteristics.

Figure 5
Relationship curves of aerodynamic coefficients.

The results show that the lift coefficient raises with the increase of the angle of attack and the elevator, which presents the characteristic of centrosymmetric cubic curve. The Mach number has a small effect on the lift coefficient, while the height has little effect on the lift coefficient. The drag coefficient decreases first and then raises with the increase of the angle of attack, and it has the same change with the elevator. However, the drag coefficient decreases with the increase of the Mach number, the pitch moment coefficient raises with the increase of the height, the pitch angle raises with the increase of the angle of attack, and the lift decreases with the increase of the angle of the elevator, which presents the characteristic of centrosymmetric cubic curve. The influence of Mach number and height on pitch moment coefficient is very limited.

As the above analysis shows that the Mach number and height have little influence on the aerodynamic coefficient, the influence is further analyzed in Fig. 6. In order to focus on the influence of the height and Mach number, only the analysis results of the angle of attack at 4 degrees are shown in the figure, and the results are similar within other angles of attack. It shows that the lift coefficient, the drag coefficient and the pitching moment present a nonlinear relationship on Mach, which can be approximated as a quadratic curve. The effect of height on the lift coefficient is negligible. The relationship between height and drag coefficient can be approximately linear. The effect of height on the pitching moment is shown as a small nonlinear relationship, which can be approximately linear.

Figure 6
Relationship between aerodynamic coefficients and Mach number and height.

In summary, in order to facilitate the analysis of the control characteristics and the design of the controller, the relationship between force coefficients and parameters should be expressed as the continuous function approximately. The lift coefficient can be approximated as the quadratic function of the Mach number, the linear function of the attack angle and the cubic function of the elevator angle under the condition that the fitting precision is satisfied. The drag coefficient is a quadratic function of Mach number, attack angle and elevator deflection, as well as a linear function of height. The pitching moment is approximated as the linear function of the attack angle and height, the quadratic function of the Mach number and the cubic function of the elevator angle.

AERODYNAMIC MODEL

Through the analysis in the previous section, the basic fitting function of each coefficient can be expressed as Eq. 8.

(8) $\begin{array}{c}{\mathit{C}}_{\mathit{L}}={\mathit{a}}_0{\mathit{M}}_{\infty }^2+{\mathit{a}}_1{\mathit{M}}_{\infty }+{\mathit{a}}_2\mathit{\alpha }+{\mathit{a}}_3{\mathit{\delta }}_{\mathit{p}}^3+{\mathit{a}}_4{\mathit{\delta }}_{\mathit{p}}^2+{\mathit{a}}_5{\mathit{\delta }}_{\mathit{p}}+{\mathit{a}}_6{\mathit{\alpha }\mathit{\delta }}_{\mathit{p}}+{\mathit{a}}_7\hfill \\ {\mathit{C}}_{\mathit{D}}={\mathit{b}}_0{\mathit{M}}_{\infty }^2+{\mathit{b}}_1{\mathit{M}}_{\infty }+{\mathit{b}}_2{\mathit{\alpha }}^2+{\mathit{b}}_3\mathit{\alpha }+{\mathit{b}}_4{\mathit{\delta }}_{\mathit{p}}^2+{\mathit{b}}_5{\mathit{\delta }}_{\mathit{p}}+{\mathit{b}}_6{\mathit{\alpha }\mathit{\delta }}_{\mathit{p}}+{\mathit{b}}_7\mathit{H}+{\mathit{b}}_8\hfill \\ {\mathit{C}}_{\mathit{m}}={\mathit{c}}_0{\mathit{M}}_{\infty }^2+{\mathit{c}}_1{\mathit{M}}_{\infty }+{\mathit{c}}_2\mathit{\alpha }+{\mathit{c}}_3{\mathit{\delta }}_{\mathit{p}}^3+{\mathit{c}}_4{\mathit{\delta }}_{\mathit{p}}^2+{\mathit{c}}_5{\mathit{\delta }}_{\mathit{p}}+{\mathit{c}}_6{\mathit{\alpha }\mathit{\delta }}_{\mathit{p}}+{\mathit{c}}_7\mathit{H}+{\mathit{c}}_8\hfill \end{array}$

The vector [M_∞, α, δ_p,H] is the parameters of the fitting model (Eq. 8), where the dimension units of each component are the local velocity of sound, degree, degree, and kilometer (km), respectively. The undetermined coefficient in the formula can be solved by the least square method (Eq. 9):

(9) $\mathbf{a}=\left[\begin{array}{c}-3.78046238064082e-06\hfill \\ 0.000127277324253248\hfill \\ 0.0374144130261035\hfill \\ 9.80441547014421e-06\hfill \\ -2.15046985750635e-07\hfill \\ 0.00460072818180406\hfill \\ -5.00461946588615e-06\hfill \\ 0.00726766542956915\hfill \end{array}\right],\mathbf{b}=\left[\begin{array}{c}0.000628445665793072\hfill \\ -0.0173353487394285\hfill \\ 0.000966644232057383\hfill \\ 0.000632302167887976\hfill \\ 0.000235495442102079\hfill \\ -8.50034940855780e-06\hfill \\ 0.000366887347544208\hfill \\ 3.660866135459709e-05\hfill \\ 0.149767920614365\hfill \end{array}\right],\mathbf{c}=\left[\begin{array}{c}-1.15155410551467e-05\hfill \\ 0.000424997398866636\hfill \\ 0.0532370917577559\hfill \\ -1.83350642091320e-05\hfill \\ -5.12209384489568e-07\hfill \\ -0.00829204196445136\hfill \\ 7.889540588874551e-06\hfill \\ 5.27361418925611e-09\hfill \\ 0.0145120939976004\hfill \end{array}\right]$

As shown in Table 2, the correlation coefficients of all aerodynamic coefficients are greater than 0.98 with a small error and a large F-statistical value, which means the fitting model can meet the regression requirements.

Table 2
Statistics of the fitting model.

Coefficients	R2-statistic	F-statistic	p-value	Error variance
C L	0.9948	1.305e6	0	2.015e-4
CD	0.9858	4.107e5	0	3.065e-5
CM	0.9929	8.259e5	0	6.317e-4

According to Fig. 7, when the Mach number is 10, it can be intuitively analyzed that each coefficient surface is basically consistent with the trend of the original data points, and there are errors in some data points (results are similar in other Mach numbers). Through the results of quantitative and qualitative analysis, it can be seen that the fitting accuracy of each aerodynamic coefficient meets the design requirements of the controller.

Figure 7
Direct qualitative error analysis of fitting model on M_∞ =10.

AEROTHERMODYNAMIC ENVIRONMENTS PREDICTION

A hypersonic aerothermodynamic environment platform has been developed to predict the complex configuration vehicles in an early thesis. The thermal prediction method of the platform is the engineering method based on flow field topology and has been validated by several typical vehicles in the thesis of Meng (2019). Here, the thermal environment of PRAHV is calculated by the platform directly. It can be known by trial calculation or research experience that the thermal environment at the nose and rudder is the most serious. When setting the conditions of M_∞ = 10, H = 20 km, α = 0, and δ_p = 0, the aerothermodynamic environment of PRAHV has been given in Fig. 8. It is shown that the heating rate is distinctly serious at the nose and the wing rudder, and the maximum heating rate occurs at the nose.

Figure 8
Aerothermodynamic environment (M_∞ = 10, H = 20 km, α = 0, δ_p = 0). (a) thermal environment over PRAHV surface; (b) heating rate on a slice along center line

The nose is the most serious location of heating rate because it is a stagnation point in the flow field. Furthermore, the heating rate of the nose usually sets as the constraint of designing in lots of research. So, a heat flow data set are calculated along the flight corridor, and then a heating rate formula of nose is established by parameter estimation (Eq. 10).

(10) ${\mathit{q}}_{\mathit{w}}=1.52\mathit{e}-6{\mathit{\rho }}^{0.5}{\mathit{V}}^3$

where ρ is atmospheric density related to altitude H, which can be obtained by a standard atmosphere program.

In Fig. 9(a), the results calculated from Eq. 10 and the results of the platform are contrastively presented in the same flight corridor. It can be seen that Eq. 10 is feasible for PRAHV’s nose heating rate prediction and can be used as the data basis for trajectory calculation. The heating rate distribution mapped on the altitude and Mach number is shown in Fig. 9(b). The logarithm of heating rate at base 10 is depicted by contour lines to clearly define the available flight conditions for each heating rate. Flight trajectory design should refer to this heating rate distribution as well.

Figure 9
Heating rate of nose along the flight corridor of PRAHV.

STABILITY ANALYSIS OF MODEL

In engineering, the partial derivative C_m^CL = dC_m/dC_L is generally used to measure the static stability of a model, which is called static stability. The aerodynamic model established above is sufficient to support the static stability analysis. First, the function c_m=fC_L) is established for the pitching moment coefficient and lift coefficient, and then the derivative C_m^CL is calculated.

The instantaneous moment balance is adopted to assume that the flight path is stationary or composed of several tiny stationary flights. In order to make the vehicle always keep the longitudinal balance state at a certain flight attack angle, the elevator must be deflected at the corresponding angle, namely the pitch trim angle. Figure 10(a) shows the cure of trim angle function δ_p = f(α) in the range of -10 ~ 10 degrees of attack at each Mach number corresponded height. Static stability is the balance maintenance ability of the vehicle itself. Some factors may make the attitude angle of attack change and generate additional lift, and then the additional lift will cause the change of pitching moment. The relationship between lift and pitching moment of the model established in this paper is shown in Fig. 10(b). Whether the moment has the ability to recover the attitude is the static stability. In order to measure the static stability quantitatively, the stability function C_m^CL = f(α) is further calculated. As shown in Fig. 10(c), C_m^CL>0 for all the attitude angle, that is, the rise of moment is generated to further increase the angle of attack away from the balanced state when the angle of attack increases, indicating that the longitudinal stability of the vehicle is unstable.

Figure 10
Trim angle and static stability of model.

The dynamic stability was further analyzed, and the flight motion dynamic system was established and derived from dynamics and kinematics details in reference book (Zong et. al. 2016), as shown in Eq. 11. The dynamic system has been applied in a large number of literatures for guidance and control research, such as model analysis of air-breathing hypersonic vehicle by Bolender of U.S. Air Force Research Laboratory (Bolender and Doman 2007), and algorithm research of output feedback control by Li et al. (2011b).

(11) $\begin{array}{c}\frac{\mathit{dh}}{\mathit{dt}}=\mathit{\nu }\sin \mathit{\gamma }\hfill \\ \frac{\mathit{d}\mathit{\nu }}{\mathit{dt}}=\frac{1}{\mathit{m}}\left(-\mathit{D}-\mathit{mg}\sin \mathit{\gamma }\right)\hfill \\ \frac{\mathit{d}\mathit{\delta }}{\mathit{dt}}=\frac{1}{\mathit{m}\mathit{\nu }}\left(\mathit{L}-\mathit{mg}\cos \mathit{\gamma }\right)\hfill \\ \frac{\mathit{d}\mathit{\alpha }}{\mathit{dt}}=\mathit{q}-\frac{1}{\mathit{m}\mathit{\nu }}\left(\mathit{L}-\mathit{mg}\cos \mathit{\gamma }\right)\hfill \\ \frac{\mathit{dq}}{\mathit{dt}}=\frac{\mathit{M}}{\mathit{I}}\hfill \end{array}$

where D = 0.5ρV²S_refC_D, _L = 0.5ρV²S_refC_L, M = 0.5ρV²S_refC_m. Let l = 1e7kg·m² and m = 2900 kg. Let x = [h, v, γ, α, q]^T, u = δ_p and X0 = [30km, 10Mach, 0rad, 1°, 0rad/s]^T.

Then u₀ = [7.5972°] is obtained by trimming for longitudinal balance. Linearize the system (Eq. 10) at the equilibrium state x₀ by the Jacobi linearization method to obtain the linear control system, as shown in Eq. 12. The dimension units of variables in the system (Eq. 12) are exactly the same as those in the system (Eq. 11), where each component of x corresponds to m, m/s, rad, degree, and rad/s, respectively.

(12) $\begin{array}{c}\dot{\mathbf{x}}=\mathit{A}\mathbf{x}+\mathit{B}\mathbf{u}.\hfill \\ \mathrm{where},\mathit{A}=\left[\begin{array}{ccccc}0& 0& 3017.8& 0& 0\\ -1.033e-3& -2.117e-2& -9.800& -0.1511& 0\\ 0& 1.340e-06& 0& 3.496e-4& 0\\ 0& -1.340e-06& 0& -3.496e-4& 1\\ 4.317e-11& 5.281e-09& 0& 4.363-4& 0\end{array}\right],\mathrm{and},\mathit{B}=\left[\begin{array}{c}\hfill 0\\ \hfill 0.1111\\ \hfill 0.1775\\ \hfill -0.1775\\ \hfill -9.386e-5\end{array}\right]\hfill \end{array}$

The poles of the linearized system are obtained as follows: 0.0027 + 0.0122i, 0.0027 - 0.0122i, -0.0268 + 0.0000i, -0.0208 + 0.0000i, 0.0207 + 0.0000i. According to the linear control theory, the model is unstable at the equilibrium point because it has the right half-plane pole. Therefore, the appropriate controller design must be carried out to achieve stable flight or command tracking for the vehicle.

APPLICATION OF CONTROL-ORIENTED MODEL

All the above works should be the basis of guidance and control research. In general, Fig. 11 shows that three parts woks, i.e., nominal trajectory design, guidance law design, and attitude controller design, should be further handled as an open problem based on the control-oriented model (Tian and Zong 2011). For the sake of descriptive integrality, this section takes PRAHV as an application example and presents a nominal trajectory design in the case of the unpowered longitudinal maneuvering based hp-adaptive pseudospectral method. The nominal trajectory is regarded as the particle trajectory without considering the attitude control of the vehicle.

Figure 11
A full research framework of guidance and control.

Since the static longitudinal stability of the vehicle is unstable as analyzed in the previous section, the pitch angle should keep varying to maintain longitudinal balance of the vehicle at a certain flight attack angle. Therefore, the unpowered maneuver trajectory design problem is considered as a maximum range problem as follow (Eq. 13):

(13) $\begin{array}{c}\min \mathit{J}\left[\mathbf{x},\mathbf{u}\right]=-{\int }_0^{{\mathit{t}}_{\mathit{f}}}\mathit{\nu }\mathit{dt}\hfill \\ \begin{array}{cc}\frac{\mathit{s}.\mathit{t}.}{}\hfill & \frac{\mathit{dh}}{\mathit{dt}}=\mathit{\nu }\sin \mathit{\gamma }\hfill \\ \hfill & \frac{\mathit{d}\mathit{\nu }}{\mathit{dt}}=\frac{1}{\mathit{m}}\left(-\mathit{D}-\mathit{mg}\sin \mathit{\gamma }\right)\hfill \\ \hfill & \frac{\mathit{d}\mathit{\gamma }}{\mathit{dt}}=\frac{1}{\mathit{m}\mathit{\nu }}\left(\mathit{L}-\mathit{mg}\cos \mathit{\gamma }\right)\hfill \\ \hfill & \mathit{M}=0\hfill \\ \hfill & {\mathit{q}}_{\mathit{w}}<{\mathit{Q}}_{\max }\hfill \end{array}\hfill \end{array}$

where the drag D, the lift L and the moment M are the same as in Eq. 11.

The mass and mass center of the vehicle are constant (m= 900 kg, x_cg = 0.7 L_en) for the unpowered maneuvering which corresponded to the thrust shut-off by fuel exhaustion. The state and control are x = [h, v, γ] and u = α respectively in the problem. The heating rate q_w is calculated by Eq. 10, and the Q_max is a user-specified parameter to limit the heating rate of the trajectory.

The maneuver trajectory design is a typical nonlinear optimal control problem with two-point boundary values. And the hp-adaptive pseudospectral method proposed by Darby et al. (2011) is suitable for solving the problem. Several mesh refinement algorithms (include hp-adaptive) for the pseudospectral method have been discussed in previous research (Huang et al. 2019). Therefore, the hp-adaptive pseudospectral method is directly used in this paper to solve the problem.

Figure 9 shows that the heating rate restricts the speed and altitude of the vehicle. The values of two-point boundary in the trajectory design problem should be given, but it directly affects the heating rate in the trajectory. To find a nominal trajectory with a low heating rate, several pairs of boundaries are tried to calculate the problem (Eq. 12) without the constraint of heating rate. Therefore, the values of the two-point boundary are given in Table 3.

Table 3
The values of two-point boundary for the unpowered trajectory problem.

State boundary	Initial state			Terminal state
	h0(km)	v0(Mach)	γ0(rad)	hf(km)	vf(Mach)	γf(rad)
B(30,10-5,5)	30	10	-π/2~π/2	5	5	-π/2~π/2
B(30,10-10,5)	30	10	-π/2~π/2	10	5	-π/2~π/2
B(30,10-5,3)	30	10	-π/2~π/2	5	3	-π/2~π/2
B(50,10-5,3)	50	10	-π/2~π/2	5	3	-π/2~π/2

And the boundaries of the interior variables in trajectory are given as h = h₀ ~60km, v = v_f ~10Mach, γ= -π/2 ~π/2rad, α= -10 ~ 10°. The hp-adaptive pseudospectral method performed on each boundary of Table 3, and the results of trajectories are shown in Fig. 12. Both numerical solutions achieved by hp-adaptive pseudospectral and the simulation results obtained by ode45 are presented in the Fig. 12. The discrete numerical solutions solved by hp-pseudospectral method are presented by discrete circles, and every variable includes several intervals generated by hp mesh refinement. The one discrete point (named as the node of pseudospectral) of variable is presented as a circle, and one interval is presented by the same one color. The simulation results are plotted as the solid line with the same color as the corresponding intervals. All the chart styles of the figures in this section are totally same as previous research (Huang et al. 2019).

Figure 12
The trajectories under different two-point boundaries without heating constraint.

After iterative steps of mesh refinement, the numerical solutions are obtained. For the four settings of two-point boundary, every final iterative solution constructed by several intervals and each interval contains a number of nodes. As shown in Fig. 12, the charts (a)-(c) show the state variables changed with time, and the chart (d) shows the change of the control variable with time. The continue attack angles of time are interpolated by piecewise cubic Hermite interpolating polynomial (PCHIP) function on the discrete numerical solutions which are drawn with the solid line in the figure. Then let the continue control as the input of system (Eq. 13) and ode45 function be used for simulation, and the continuous state variables obtained by ode45 simulation are shown in (a)-(c). The plot (e) shows the pitch angles corresponding to attack angles to keep the longitudinal balance of vehicle. It is shown that the feasible trajectories are optimized by pseudospectral method because the optimal solution of the discrete numerical solution satisfies the simulation results completely. The heating rates along trajectories are presented in the chart (f), and the values of two-point boundary directly affect the heating rate of trajectory.

The high heating rates occur on the fixed states of initial point and terminal point cannot be optimized in the trajectory designing problem. But the heating rates of the interior trajectory should be optimized for the situation of B(50,10-5,3). Hence, the Q_max = 2.8e6 w/m² is used to further optimize the heating rate of the trajectory, and the obtained heating rate constraint trajectory is shown in Fig. 13.

Figure 13
The trajectory constrained of heating rate in situation of b (50,10-5,3).

After adding the heating rate constraint, the maximum value of the trajectory is less than Q_max, but the range and time of trajectory become smaller. Therefore, different computational settings are needed to balance the specific task requirements with the constraints of the thermal protection system. Anyway, the availability of the aerodynamic model and the heating rate formula in this paper is demonstrated by the above trajectory optimization calculation. Lots of control-related method researches can be performed based on the aerodynamic model and the heating rate formula, such as the researches of guidance law and attitude controller for the unstable dynamics.

CONCLUSION

In this paper, a 3D model of Kinzhal-like described as PRAHV (a hypersonic vehicle of axisymmetric body with a power-law revolution nose) is constructed based on data of the public literature. The aerodynamic data set of the 3D model was obtained by engineering prediction methods, with the range of Mach number 2 ~ 20, angle of attack -10° ~ 10°, elevator angle -20° ~ 20° and altitude 0 ~ 80 km. The aerodynamic fitting model based on dimensional correlation is established to reduce the complexity of the model and maintain the characteristics of the original data. The aerothermodynamic environment of PRAHV is discussed, and the heating rate formula of the nose is established as the basis of trajectory designing.

According to the theories of flight mechanics and motion control, the longitudinal static stability and the longitudinal dynamic stability of the model are quantitatively analyzed. The results show that the model is unstable on both static stability and dynamic stability.

Finally, to validate the availability of the aerodynamic model and the heating rate formula, the nominal trajectory design problem is solved by the hp-adaptive pseudospectral method for the case of unpowered longitudinal maneuvering. The work of this paper provides theoretical and data support for the research and judgment of Kinzhal-similar vehicle and model basis for the research of guidance and control.

• FUNDING
  Sichuan Provincial Open Foundation of Civil-Military Integration Research Institute
  Grants #2017SCII0219; #2017SCII0220

ACKNOWLEDGMENTS

Editors and authors are thankful to Fundação Conrado Wessel for providing the financial support for publishing this article.

REFERENCES

• Bolender MA, Doman DB (2007) Nonlinear longitudinal dynamical model of an air-breathing hypersonic vehicle. J Spacecr Rockets 44(2):374-387. https://doi.org/10.2514/1.23370
  » https://doi.org/10.2514/1.23370
• Che J, Tang S (2007) The engineering calculation of aerodynamics for quasi-waverider vehicle. Acta Aerodyn Sin 25(3):381-385.
• Darby CL, Hager WW, Rao AV (2011) An hp‐adaptive pseudospectral method for solving optimal control problems. Optim Control Appl Meth 32(4):476-502. https://doi.org/10.1002/oca.957
  » https://doi.org/10.1002/oca.957
• Li HF (2012) Hypersonic vehicle guidance and control technology (upper and lower) (Chinese Edition). Beijing: China Astronautic Publishing House.
• Grozovsky GL, An XF, Yu DB (1978) The supersonic aerodynamics of power law rotation (in Chinese). China: National Defense Industry.
• Huang J, Liu ZG, Liu ZQ, Wang QF, Fu J (2019) A pk-adaptive mesh refinement for pseudospectral method to solve optimal control problem. IEEE Access 7:161666-161679. https://doi.org/10.1109/ACCESS.2019.2952139
  » https://doi.org/10.1109/ACCESS.2019.2952139
• Keshmiri S, Colgren R, Mirmirani M (2006) Six-DOF modeling and simulation of a generic hypersonic vehicle, for control and navigation purposes. Paper presented AIAA Guidance, Navigation, and Control Conference and Exhibit. AIAA; Keystone, United States.
• Li P, Gao Z (2014) An engineering method of aerothermodynamic environments prediction for complex reentry configurations. Paper presented AIAA SPACE 2014 Conference and Exposition. AIAA; San Diego, California. https://doi.org/10.2514/6.2014-4414
  » https://doi.org/10.2514/6.2014-4414
• Li HF, Lin P, Xu DJ (2011a) Control-oriented modeling for air-breathing hypersonic vehicle using parameterized configuration approach. Chinese J Aeronaut 24(1):81-89. https://doi.org/10.1016/S1000-9361(11)60010-1
  » https://doi.org/10.1016/S1000-9361(11)60010-1
• Li XD, Xian B, Diao C, Yu YP, Yang KY, Zhang Y (2011b) Output feedback control of hypersonic vehicles based on neural network and high gain observer. Sci China Inf Sci 54(3):429-447. https://doi.org/10.1007/s11432-011-4194-y
  » https://doi.org/10.1007/s11432-011-4194-y
• Liao MH (2018a) Review and judgement for “Kinzhal” air-launched hyper-missile of Russian (in Chinese). China Aviation News, 2018-04-10(005). Bejing: Aviation Industry Corporation of China. [accessed May 10 2019]. https://kns.cnki.net/KCMS/detail/detail.aspx?dbcode=CCND&dbname=CCNDLAST2018&filename=CHQB201804100050
  » https://kns.cnki.net/KCMS/detail/detail.aspx?dbcode=CCND&dbname=CCNDLAST2018&filename=CHQB201804100050
• Liao MH (2018b) Hypersonic attack and defense: More important details disclosure about US air-sea-land hypersonic boost glide missiles (in Chinese). China Aviation News, 2018-11-06(007). Bejing: Aviation Industry Corporation of China. [accessed May 10 2019].https://kns.cnki.net/KCMS/detail/detail.aspx?dbcode=CCND&dbname=CCNDLAST2018&filename=CHQB201811060071
  » https://kns.cnki.net/KCMS/detail/detail.aspx?dbcode=CCND&dbname=CCNDLAST2018&filename=CHQB201811060071
• Military-Today (2019) Kh-47M2 Kinzhal: Air-launched ballistic missile [Internet]. [accessed May 10 2019]. http://www.military-today.com/missiles/kh_47m2_kinzhal.htm
  » http://www.military-today.com/missiles/kh_47m2_kinzhal.htm
• Moore M, Williams J (1989) Aerodynamic prediction rationale for analyses of hypersonic configurations. Paper presented 27th Aerospace Sciences Meeting. AIAA; Reno, United States. https://doi.org/10.2514/6.1989-525
  » https://doi.org/10.2514/6.1989-525
• Pamadi BN, Brauckmann GJ, Ruth MJ, Fuhrmann HD (2001) Aerodynamic characteristics, database development, and flight simulation of the X-34 vehicle. J Spacecr Rockets 38(3):334-344. https://doi.org/10.2514/2.3706
  » https://doi.org/10.2514/2.3706
• Zong Q, Zeng FL, Zhang XB, You M, (2016) Hypersonic vehicle modeling and model validation (Chinese Edition). Beijing: Science.
• Rufolo G, Roncioni P, Marini M, Votta R, Palazzo S (2006) Experimental and numerical aerodynamic data integration and aerodatabase development for the PRORA-USV-FTB_1 reusable vehicle. Paper presented 14th AIAA/AHI Space Planes and Hypersonic Systems and Technologies Conference. AIAA; Canberra, Australia. https://doi.org/10.2514/6.2006-8031
  » https://doi.org/10.2514/6.2006-8031
• Shaughnessy JD, Pinckney SZ, McMinn JD, Cruz CI, Kelley M-L (1990) Hypersonic vehicle simulation model: winged-cone configuration. NASA Technical Memorandum 102610. Hampton: NASA Langley Research Center. [accessed May 10 2019]. https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19910003392.pdf
  » https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19910003392.pdf
• Tian B, Zong Q (2011) Optimal guidance for reentry vehicles based on indirect Legendre pseudospectral method. Acta Astronaut 68(7-8):1176-1184. https://doi.org/10.1016/j.actaastro.2010.10.010
  » https://doi.org/10.1016/j.actaastro.2010.10.010
• Wikipedia: the free encyclopedia (2019a) 9K720 Iskander. [Internet]. [accessed Jun 20 2019]. https://en.wikipedia.org/w/index.php?title=9K720_Iskander&oldid=898253956
  » https://en.wikipedia.org/w/index.php?title=9K720_Iskander&oldid=898253956
• Wikipedia: the free encyclopedia (2019b) Kh-47M2 Kinzhal [Internet]. [accessed Jun 20 2019]. https://en.wikipedia.org/w/index.php?title=Kh-47M2_Kinzhal&oldid=900203756
  » https://en.wikipedia.org/w/index.php?title=Kh-47M2_Kinzhal&oldid=900203756
• Meng XL (2019) Research and Implementation of thermal environment prediction platform for complex shape aircraft based on surface flow field topology (Thesis) MianYang: Southwesh University of Science and Technology. [accessed May 01 2019]. https://kns.cnki.net/KCMS/detail/detail.aspx?dbcode=CMFD&dbname=CMFD201902&filename=1019171321.nh
  » https://kns.cnki.net/KCMS/detail/detail.aspx?dbcode=CMFD&dbname=CMFD201902&filename=1019171321.nh
• Huang ZC (1995) Aerodynamics of hypersonic vehicle (Chinese Edition). Beijing: National Defense and Industry.

Edited by

Section Editor: T John Tharakan

Publication Dates

• Publication in this collection
  16 Sept 2020
• Date of issue
  2020

History

• Received
  11 May 2019
• Accepted
  08 Mar 2020