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

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. https://doi.org/10.5028/jatm.v12.1164 ORIGINAL PAPER


INTRODUCTION
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. 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. 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

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) 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) 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) 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) 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: x log 10 r*-1.5 , x=log10(Re)-1.5, r*= Re den·z 2.5 , den= T + Sh Z·T+Sh , Sh=110.4, z=1.0+0.115· M 2 ∞ , c fi =0.088/x 2 . 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): 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.
Similar engineering calculation methods have been verified by Che and Tang (2007). Furthermore, to validate the abovementioned 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.

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. 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.
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): 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. 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.

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.
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).
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.

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 =f(C 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. 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). Then u 0 = [7.5972°] is obtained by trimming for longitudinal balance. Linearize the system (Eq. 10) at the equilibrium state x 0 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.
where, A = , and, B = 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 halfplane 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. 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): 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.
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 0 ~60km, v = v f ~10 Mach, γ = -π/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).   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 2 is used to further optimize the heating rate of the trajectory, and the obtained heating rate constraint trajectory is shown in Fig. 13.
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