Research Article  Open Access
Unsteady Vibration Aerodynamic Modeling and Evaluation of Dynamic Derivatives Using Computational Fluid Dynamics
Abstract
Unsteady aerodynamic system modeling is widely used to solve the dynamic stability problems encountering aircraft design. In this paper, single degreeoffreedom (SDF) vibration model and forced simple harmonic motion (SHM) model for dynamic derivative prediction are developed on the basis of modified Etkin model. In the light of the characteristics of SDF time domain solution, the free vibration identification methods for dynamic stability parameters are extended and applied to the time domain numerical simulation of blunted cone calibration model examples. The dynamic stability parameters by numerical identification are no more than 0.15% deviated from those by experimental simulation, confirming the correctness of SDF vibration model. The acceleration derivatives, rotary derivatives, and combination derivatives of ArmyNavy Spinner Rocket are numerically identified by using unsteady NS equation and solving different SHV patterns. Comparison with the experimental result of Army Ballistic Research Laboratories confirmed the correctness of the SHV model and dynamic derivative identification. The calculation result of forced SHM is better than that by the slender body theory of engineering approximation. SDF vibration model and SHM model for dynamic stability parameters provide a solution to the dynamic stability problem encountering aircraft design.
1. Introduction
In aircraft design and development, it is difficult and costly to achieve dynamic stability within the flight envelop. Shock waveinduced separation, swirling motion and rupture, and the interaction between them lead to strong unsteady and nonlinear behaviors of fluid motion, which generally results in aerodynamic effect beyond expectation and even some overwhelming consequences [1]. In the aircraft design stage, it is hard to foresee the boundaries of dynamic operation stability and the severity level of dynamic problems, most of which are not discovered until in the flight test stage. This directly leads to multiplication of the aircraft design cycle, significant increase of design cost, and unavoidable aircraft performance loss in connection with local reshaping. Since the beginning of this century, dynamic stability problems have been a common phenomenon [2]. To better understand the causes of flight stability deterioration from the fluidity mechanism horizon, unsteady aerodynamic mathematic model has been widely applied to the studies on dynamic flight quality issues [3]. A mature, reliable package of unsteady aerodynamic mathematic models enables us to evaluate and select different aircraft design plans so as to minimize the aircraft design cost and mitigate risk exposure [4].
In aircraft flight, atmospheric turbulence, asymmetric transition, base flow interference, maneuver implementation, and other disturbances can result in highly complex aerodynamic consequences such as asymmetric vortex separation, shedding or breakup, shock motion, and vortex or boundary layer disturbances, which eventually give rise to dynamic stability problems of the aircraft [5]. Longitudinal or lateral vibration amplitude and frequency vary linearly with the motion attitude angle and its variation rate: the vibration process can be either convergent or divergent. The vibration amplitude can span from a few to dozens of degrees or even be divergent [6]. These greatly affect the stability control and trajectory design of a hypersonic aircraft and even the imaging of the airborne camera. Severe dynamic behavior problems could cause the aircraft to disintegrate [2].
In the hypersonic or supersonic stage, when a spacecraft reenters the atmosphere, it is often challenged by the rollinghelical modal coupling problem. In its X33 research, the United States discovered that the lifting body it designed is noticeably laterally unstable in construction. All the validation projects implemented by the United States Defense Advanced Research Projects Agency (DARPA) and the United States Air Forces on Falcon hypersonic technology vehicle (HTV) ended up with failure [7], and all these failures are closely associated with the stability control of the aircraft attitude. The stability analysis result of University of Maryland suggests that the waverider aircraft tends to be longitudinally unstable under three hypersonic states [8]. The analysis report on the first flight failure of X43A points out that the specially structured aerodynamic forces of the aircraft have far greater than expected nonlinearity, and the inaccurate aerodynamic modeling is considered to be one of the critical failure aspects [9].
Dynamic unsteady aerodynamic system modeling is widely used to solve these dynamic stability problems encountering aircraft design [4]. An unsteady aerodynamic model is used to determine the motion state variables that decide aerodynamics and the relation among these variables. In this model, the aerodynamic load to which the aircraft is subject is expressed as the function of the motion state parameter and the transient value of its derivative [10]. The coefficients of the orders in the function are called dynamic derivatives, which directly decide the divergence feature of the aircraft openloop system when oscillating under disturbance [11]. For a free flight aircraft, dynamic derivative prediction is indispensable and directly decides the flight stability. Figure 1 gives the five motion states of a free flight aircraft after disturbance. The stability of the aircraft closedloop system is not only decided by aerodynamic stability; it is also associated with the control system. The control system can change the system stability by feeding back the attitude plane and angular speed. Dynamic derivatives decide the gain factor of the control system. Excessive adjustment of the feedback gain can degrade the aircraft maneuverability and deteriorated controllability and flight qualities though good flight stability is guaranteed.
(a) Steady linear flight
(b) Static instability
(c) Static neutral stability
(d) Static stability, dynamic stability
(e) Static stability, dynamic neutral stability
(f) Static stability, dynamic instability
The idea of stability and dynamic derivative was first proposed by Bryan [12] in 1911 and is still used today in some engineering problems as a conventional derivative model for aerodynamic load [10]. Over the hundred years’ history of unsteady aerodynamic mathematic model, one after another, modified models targeting at the shortfalls of Bryan model have been introduced. These include state space model [13] (differential equation model), indicating function model [14] (integral equation model), Etkin fullorder model, and reducedorder model having gained rapid development over the past few years [15–17]. Table 1 lists the dynamic stability derivatives of the moment coefficient, where , , and represent rolling, pitching, and yawing moment coefficients, respectively; and are the angle of attack and angle of sideslip; , , and are the angular speed components of the rolling axis, pitching axis, and yawing axis.

In our study, single degreeoffreedom (SDF) vibration model and forced simple harmonic motion (SHM) model for dynamic derivative prediction is developed on the basis of modified Etkin model [5]. To validate the predictions, two test cases are used. For blunted cone calibration model, the coupling computation between aircraft motion and flow field motion equations is studied using the unsteady SDF vibration model. For ArmyNavy Spinner Rocket, dynamic derivatives at different positions of center of mass are identified by solving threedimensional unsteady NS equation and compared with the computational result of the United States Army Research Laboratory and experimental result of the United States Army Ballistic Research Laboratories.
2. Computational Fluid Dynamics Formulation
Under the Cartesian inertial coordinate system, the threedimensional dimensionless NavierStokes equation for perfect gases with ignored mass forces is as follows:where is the vector of conserved variables, , is the density, and , , and are the components of velocity given by the Cartesian velocity vector . The total energy per unit mass is . The inviscid flux vectors, , , and , are given bywhile the viscous flux vectors, , , and , contain terms of the heat flux and viscous forces exerted on the body and can be represented byHere, is the stress tensor components. is the heat transfer coefficient. ConsiderHere, is the specific heat ratio. and are the laminar and turbulent Prandtl number. and are the freestream Mach number and Reynolds number, respectively. The SpalartAllmaras (SA) model is employed to evaluate the turbulent viscosity . The laminar viscosity, , is evaluated using Sutherland’s law:
Secondorder upwind NND scheme finite volume [18, 19] was used to discretize the spatial derivative term of the flowcontrolled equation group. The LUSGS containing dual time steps [20, 21] was used to obtain higher efficiency and precision of time solution.
Rigid grid generation method where grids are fixedly connected to the object is used for unsteady dynamic calculation. When the object is in motion, the coordinate system remains unchanged while the grids, fixedly connected to the object, make rigid motion. The aircraft motion is decomposed into translation with the centroid and rotation with the centroid. Assume that the Euler angle sequence during aircraft rotation is pitching , then yawing , and then rolling . The Euler angle direction is positive along the positive direction of the coordinate axis according to the righthand rule. Record the centroid position at the initial time as , the aircraft centroid position at time as , and the Euler angle as . Then, for any grid point at the initial time , the coordinate at time iswhere is the elementary transformation matrix of the coordinate system. The greatest merit of dynamic grids generated by this method is that they are simple and applicable for unsteady motions of any configuration and any amplitude. Dynamic grids are as good as static grids and conform to geometric conservation law.
For dynamic unsteady calculation, farfield boundaries were treated with nonreflecting boundary conditions (NRBC) based on Riemann variable and applicable for dynamic movement. For wall boundary, the velocity was accounted for nonslip conditions, the temperature was accounted for diathermic wall, and the pressure conditions accounted for centrifugal forces.
3. Single DegreeofFreedom Vibration Modeling
The free pitch motion of aircraft around the fixed axis is expressed by the equation for rigid body rotation with a fixed axis:Here, is the dimensionless moment of inertia; is the dimensionless mechanical damping coefficient; is the angle of pitch starting from the balance point.
When flight aerodynamics discusses the dynamic stability of aircraft under disturbance, the aerodynamic forces on the aircraft are expressed as the linear function of the motion dynamic parameters [22]. For SDF pitch motion, under the hypothesis of Etkin’s aerodynamic model [5], the pitch moment coefficient can be expressed asHere, are the derivatives of orders of the angle of attack and its time; are the derivatives of orders of the pitch angular speed and its time. For (8), Taylor series expansion is performed at the balance angle of attack (omitting the highorder terms):Here, and are the pitch rigidity and pitch damping derivatives at the balance angle of attack. From (9), the free vibration of aircraft is not only subject to aerodynamic recovery moment due to deviation from the balance angle of attack, but also subject to the air damping moment caused by changes in the angle of attack or the angle of pitch [10]. Substitute (9) into (7) and we get the secondorder differential equation for aircraft SDF vibration:
When giving the analytic solution to this differential equation, we regard as a constant and (10) as a constant coefficient linear differential equation.
3.1. Small Damping Modal
For the small damping modal, the roots of its characteristic equation are a conjugate complex root pair. Hence the general solution to (10) is
From (11), the motion of an object is a vibration over the time cycle . However, different from simple harmonic vibration (SHV), its amplitude varies exponentially with time as shown in Figure 2. The envelop line equation connected by angular amplitudes isAs in numerical calculation, , when(i), the amplitude attenuates representing dynamic stability as shown in Figure 2(a);(ii), the amplitude stays unchanged representing dynamic neutral stability as shown in Figure 2(b);(iii), the amplitude increases representing dynamic instability as shown in Figure 2(c).
(a) Amplitude attenuated
(b) Amplitude unchanged
(c) Amplitude diverged
By numerically simulating the free pitch process of aircraft, we get the timerelated function curve of angle of pitch as shown in Figure 3. The function expression of this curve has been obtained from the analysis above (11).
Select angles of pitch from the curve one time cycle apart from each other, and (corresponding to and , and ):By dividing them, we get
3.2. Large Damping Modal
For the large damping modal, its characteristic equation has two unequal real roots. The general solution to the equation isHere, arbitrary constants and can be determined by the initial conditions. If static stability is included, we have . From (15), there is a maximum of one that makes . That is, the object can cross over the balance point one time at most. Hence the object does not vibrate anymore. The function graph is as shown in Figure 4.
(a) Amplitude attenuated
(b) Amplitude diverged
Dynamic stability parameters are solved using a linear regression model. Define independent variable and dependent variable :
Matrix that associates the unknown parameters and dependent variable is defined as
Using linear regression process for solution, we get
Normally, matrix is not a square array and has more rows than columns. Many numerical methods can solve this least square (LS) problem [23]. These include direct transposition , Gaussian elimination method, MoorePenrose generalized inversion method, and QR factorization though MoorePenrose generalized inversion method and QR factorization are more accurate than direct transposition and Gaussian elimination. In terms of calculation quantity, it is for QR factorization and for MoorePenrose generalized inversion method.
3.3. Critical Damping Modal
For the critical damping modal, its characteristic equation has two unequal real roots:The general solution to the equation isArbitrary constants and can be determined by the initial conditions. From (20), there is a maximum of one that makes . Hence the object does not vibrate any more.
4. Forced Simple Harmonic Motion Modeling
Figures 5–7 show the aircraft flight attitudes and motion profiles under three types forced SHV (pitching/plunging/flapping).
Forced pitching SHV is used to calculate the pitch direct damping derivative , the vibration equation of which isHere, is the initial angle of attack; is the pitch angle amplitude; is the dimensionless reduced frequency (, where is the circular frequency; is the reference length; is the inflow velocity).
Forced plunging SHV is used to calculate the acceleration derivative of pitch moment , the vibration equation of which is
Forced flapping SHV is used to calculate the rotary derivative of pitch moment , the vibration equation of which is
The superposition of plunging vibration and flapping vibration equals pitching SHV. In the next paragraphs, the pitch direct damping coefficient is calculated using pitching SHV.
For aircraft under SHM, when the dynamic flow field reaches the periodical solution, it means the disturbance started infinitely long time ago, and modified Etkin model [5] does not include time as the explicit statedependent variable. In such case, the state variables of the rolling round pitch axis motion are and . As the base flight state is symmetric straight linear flight and subject to very modest disturbance amplitude, only firstorder dynamic derivatives are included in the calculation. By performing Taylor expansion on (8) and ignoring small quantities of higher orders, we haveHere and are the inphase and outphase components of [24].
The resulting time domain data are postprocessed using formula (25) to identify the dynamic derivatives. Methods to realize this include integral method, domain frequency transformation, and phase method. In our example, the postprocessing results by integral method, frequency domain transformation, and phase method do not show significant divergence. Under a small reduced frequency, however, domain frequency transformation provides more accurate result than other methods. Phase method, which uses only two points of the time domain data, provides lower calculation accuracy and is thereby less frequently used.
4.1. Integral Method
Integral method derives and by calculating the first Fourier coefficient of over the time history of the first cycle :
From (25) we can see that the frequency effect of the inphase component is significantly greater than the outphase component . As the order of magnitude of the highorder coefficient quickly reduces with the increase of the order of derivative. By ignoring the influence of the highorder coefficient and frequency effect, we can approximate to
The inphase/outphase component approximation (27) applies to the extent that the aircraft flies at a small angle of attack under linear aerodynamic conditions, when the inphase and outphase components are able to approximate the static derivative and combination derivative.
When the flight is at a large angle of attack or subject to significant frequency effect, the influence of the highorder term coefficient is not ignorable [24, 25]. However, it is noted that when highorder coefficient is needed to express aerodynamics, as the orthogonal nature of model sequence itself keeps the fundamental harmonic unchanged, inphase and outphase components are still very important terms in aerodynamic expressions.
When the flight is in an area subject to high nonlinear, unsteady aerodynamics, there can be a problem with the applicability of linearized model. In such case, a more advanced mathematic model [26–28] is needed to express the nonlinear, unsteady aerodynamics.
4.2. Frequency Domain Transformation
Frequency domain analysis has many merits and is applied to different areas. By transforming the frequency domain that acts on one cycle of stable harmonic, it is possible to obtain the frequency spectral characteristics of the aerodynamics and its inphase/outphase components [29]. The finite Fourier integer of a continuous time scalar function within a definite time interval is defined asHere is an imaginary unit and is the spatial frequency.
Finite Fourier integer can be calculated at the discrete frequency , which is equal spaced from zero frequency through to Nyquist frequency . Its frequency resolution is the reciprocal of the time step length [30]. The solution of any frequency on a given wave band can be obtained by linear frequency modulation transformation [31]. The amplitude ratio and phase angle closely related to input signals in the system response can be determined by the transfer function of input and output. For pitch moment coefficient, the transfer function of pitch moment coefficient is [1]
4.3. Phase Method
When forcing the aircraft to make smallamplitude, lowfrequency SHV, perform Fourier expansion on the unsteady pitch moment :Here is the aerodynamic force/moment amplitude of the fundamental harmonic component; is the phase angle between the displacement and aerodynamic moment of the aircraft when under forced vibration; is the highorder harmonic component.
To the extent that the unsteady calculation time is long enough and the calculation cycle is large enough, the influence of the initial effect is ignorable and the aerodynamic force/moment reaches a periodical steadystate value.
Ordering , we have
Ordering , we have
5. Pitching Motion Stability of Blunted Cone
5.1. Computational Model and Conditions
For the 10degree blunted cone calibration model, the bottom diameter is mm and bluntness ratio is . Inflow conditions are as follows: , Reynolds number based on bottom diameter , , and K.
Three groups of examples are included. The first group uses different angular amplitudes (, , and ) to examine the effect of initial angular amplitude. The second group uses different centroid positions (, , and ) to examine the effect of vibration axis position. The third uses different moments of inertia (, , and ) to examine the effect of moment of inertia and also compares the results between free vibration and forced vibration.
5.2. Computational Grids
The calculation mesh is shown in Figure 8. The number of cells is (flow direction × circumferential × normal). The minimum cell distance is ( is the head radius). To examine the effect of grid density on the calculation result, the grids are refined onefold in all the three directions into fine grids. The number of flow direction × circumferential × normal grids is . The total number of grids is eight times that of the previous grids.
(a) 3D mesh
(b) Rear view of mesh
The free pitch process of the cone is numerically simulated. By solving the SDF vibration model, the time history and phase curve of the free vibration of the blunted cone calibration model are derived as shown in Figures 9 and 10. From these charts we can see that the blunted cone calibration model belongs to a small damping modal. It takes 22 hours to calculate one state of coarse grids using single process method. The calculation result from fine grids in the chart is completely the same as that from coarse grids, suggesting that the density of coarse grids is enough for calculating unsteady aerodynamic and aerodynamic damping derivatives.
5.3. Flow Field Characteristics of Blunt Cone Configuration
Figures 11 and 12 give the pressure distribution and surface limiting streamline at the angle of attack of 12 degrees. The shock waves on the windward side moved toward the wall while those on the leeside elevated as a result of the angle of attack and the exposure of their strength to expansion waves gradually decreased. After the shock, the direction and speed of the flow field changed. From the pressure cloud at the angle of attack of 12 degrees, we can clearly see the distribution of bow shocks and Mach waves.
5.4. Dynamic Stability of Blunt Cone Configuration
Figure 13 shows the effect of initial angular amplitude on pitching damping derivative. Here, “+,” which is the computational point of our thesis, is obtained assuming and . Experiment [32] indicates that, for large bluntness ratio, 10degree cone configuration (e.g., ), the difference of initial angular amplitude does not make much difference to dynamic stability. As can be observed from the chart, the calculation result agrees well with the experiment. The variation of with initial angular amplitude is insignificant.
In Figure 14, the lateral axis is the axis position of free vibration, taken as , , and , respectively. The longitudinal axis is the pitching damping derivative. “+” is the result of our calculation, assuming , . We can see that as the axis moves backward, the dynamic stability reduces. Our numerical calculation agrees well with the experiment.
Assuming the maximum angular amplitude and , the free pitching process of the 10degree blunt cone is numerically simulated at and 1. The corresponding angles of pitch variations are shown in Figures 15(a)~15(d). From the profile of the curves, as increases, the amplitude attenuation becomes slower and slower, and the vibration period extends gradually. At this time, and in (11) are
(a)
(b)
(c)
(d)
We can see that when increases, the absolute of reduces and causes the amplitude to attenuate more slowly. The reduction of also leads to a longer vibration period. When is too small, overdamping occurs in Figure 15(a), preventing the curve from vibrating.
After investigating the effect of moment of inertia, we compared free vibration with forced vibration. The amplitude of forced vibration is 1°. The vibration frequency is provided by free vibration. Vibration frequencies are read out from Figures 15(b)–15(d); that is, , , and correspond with frequencies , , and . The calculation results are listed in Table 2. The two methods agree well with each other in terms of pitching rigidity, but the pitching damping derivatives vary a little bit significantly. However, as the vibration frequency reduces, the two calculations are closer to each other.

6. CFD Predictions of Dynamic Derivatives for ArmyNavy Spinner Rocket
6.1. Computational Model and Conditions
The validation example is the ArmyNavy Spinner Rocket (ANSR). A number of dynamic stability derivative experiments and calculations have been carried out targeting at the shape of ANSR calibration models [33–36]. Figure 16 shows the dimensions of the 5caliber and 9caliber ANSR calibration models, where the length of one caliber is 20 mm.
In our study, the acceleration derivative, rotary derivative, and combination derivative of 5caliber and 9caliber ANSR calibration models along the pitch direction at the Mach number of 2.5 are identified. The correctness of the algorithm is validated by comparing calculations and experiments. Besides, the influence of the position of center of mass on the damping derivative is investigated by selecting three different positions of center of mass (see Table 3) for the calibration models. Here, one caliber equals 20 mm. For all the examples, the initial angles of attack and angles of sideslip are zero degrees. The free flow conditions are taken as the sea level conditions under standard atmosphere (101.325 kPa, 288 K).

6.2. Computational Grids
The threedimensional calculation mesh is generated by rotating the symmetric plane mesh. Figure 17 shows the symmetric plane mesh near the rocket and topological structure. The numbers of cells of the 5caliber and 9caliber ANSR calibration meshes are 275,000 and 296,000. The scale of the first layer of normal mesh is , where is the rocket diameter. A small cell Reynolds number of around 10 is maintained. The mesh extension ratio near the head is limited to 1.1 maximum.
To examine the effect of grid density on the calculation result, the grids are refined threefold in all the three directions into fine grids. The total number of grids is eight times that of the previous grids. Figure 18 shows the time history and hysteresis loop curves of pitching moment coefficient of 5caliber ANSR model when its centroid position is 3.5. When the initial transient disturbance of forced pitching SHM deteriorates, the flow parameters and aerodynamic load come to harmonic, as shown in Figure 18(a). In Figure 18(b), the hysteresis loop under dynamic conditions appears to be smoothly inclined ellipse. The hysteresis loop area indicates the level of contribution of environment during the vibration of the projectile. The rotation direction of the hysteresis loop of plunging motion is also marked. The rotation direction of a hysteresis loop represents the direction of such work. Anticlockwise rotation means the work of the projectile on the environment, which consumes the dynamic energy of the projectile and contributes to dynamic stability. It takes 3 hours to calculate one state of coarse grids using 4process parallel method. In the chart, the calculation result from fine grids fully agrees with that from coarse grids, suggesting that the density of coarse grids is enough for calculating unsteady aerodynamic and aerodynamic damping derivatives.
(a) Time history curve
(b) Hysteresis loop
6.3. Dynamic Derivative Simulation Result of ANSR
Figure 19 compares the calculation result of our study (CFD), the referenced [36] result (PNS, SBT), and the referenced [34] experimental result (EXP) for the longitudinal derivatives of the ANSR calibration models. In our study, forced vibration is used to obtain the acceleration derivative and combination derivative . The forced vibration amplitude is one degree and the dimensionless reduced frequency is taken as 0.1. Quasisteady method is used to obtain the rotary derivative . By directly adding and acceleration derivative , we get . Table 4 compares it with the combination derivative by forced vibration method, where the difference is not more than 2%, suggesting that the aerodynamics under these conditions satisfy the assumption of linear superposition. The main aerodynamic parametric variations are linear to the disturbance quantity.

(a) 5caliber ANSR model normal force derivative
(b) 5caliber ANSR model pitch moment derivative
(c) 9caliber ANSR model pitch moment derivative
(d) 9caliber ANSR model normal force derivative
In Figure 19, PNS is the longitudinal dynamic stability derivative obtained by a British name Qin and Paul Weinacht from Army Research Laboratory when solving coning motion and helical motion in a noninertia system using parabolic NS equation [36]. The merit of this method lies in its high efficiency that directly obtains the component forms of the combination derivative (acceleration derivative and rotary derivative). Its shortfall, however, lies in that it only applies to the small angleofattack linear state of an axial symmetric or planar symmetric shape and can only calculate the pitching and yawing channels. SBT is the result by the slender body theory of engineering approximation [36]. Except where the position of center of mass is 2.5caliber at which time the pitch moment derivatives are the same as the CFD calculations, there are order of magnitude and even symbol differences in all the other states. EXP is the experimental result of the pitch moment combination derivative by the United States Army Ballistic Research Laboratories [34]. Engineering method SBT conforms to the experiment in terms of the dynamic derivative prediction trend, but the numerical difference is significant. PNS method has the same results as our CFD calculation and experiment, validating the correctness of our method.
For the pitch moment derivative of ANSR models, the proportion of the acceleration derivative representing flowing time lag effect is highly variable and can be as high as 40%–50%. However, as the center of mass shifts backwards, this proportion gradually reduces and is even opposite to the combination derivative symbol, contributing to dynamic instability. For the normal force derivative of ANSR models, the rotary derivative is linear to the position of center of mass while the acceleration derivative is not subject to the position of center of mass. The proportion of the acceleration derivative in the combination derivative is larger than that of the rotary derivative.
7. Conclusions
A forced SHM model and SDF vibration model are developed on the basis of modified Etkin model. The coupling computation between aircraft motion and flow field motion equations is studied using the computational fluid dynamics formulations. The forced SHM model and SDF vibration model we developed fully applies to the dynamic derivative calculation of hypersonic missiles.
The acceleration derivatives, rotary derivatives, and combination derivatives of 5caliber and 9caliber ANSR models for ArmyNavy Spinner Rocket along the pitch direction at different positions of center of mass at the Mach number of 2.5 are numerically identified by using threedimensional unsteady NS equation and solving different motion patterns. Comparison with the experimental result of United States Army Ballistic Research Laboratories confirmed the correctness of the motion model and dynamic derivative identification. Our calculation result is significantly better than that by the slender body theory of engineering approximation.
Free vibration identification methods for dynamic stability parameters are developed and applied to the time domain numerical simulation of blunt cone calibration model examples. According to the shape of a 10degree blunt cone calibration model, the time domain numerical simulation of free vibration at the Mach number of 4 and angle of attack of 12 is carried out. The dynamic stability parameters by numerical identification are no more than 0.15% deviated from those by experimental simulation, confirming the correctness of our model. The pressure distribution and surface limiting streamline of a blunt cone calibration model clearly indicate the bow shock and Mach wave distribution. Our computational methods for dynamic stability parameters provide a solution to the dynamic stability problem encountering aircraft design.
Conflict of Interests
The authors declare that they have no conflict of interests.
Acknowledgment
The research presented in this document is supported by the National Natural Science Foundation of China under Grant nos. 90716015 and 11172325.
References
 A. D. Ronch, D. Vallespin, M. Ghoreyshi, and K. J. Badcock, “Evaluation of dynamic derivatives using computational fluid dynamics,” AIAA Journal, vol. 50, no. 2, pp. 470–484, 2012. View at: Publisher Site  Google Scholar
 J. R. Chambers and R. M. Half, “Historical review of uncommanded lateraldirectional motions at transonic conditions,” Journal of Aircraft, vol. 41, no. 3, pp. 436–447, 2004. View at: Publisher Site  Google Scholar
 S. H. Woodson, B. E. Green, J. J. Chung, D. V. Grove, P. C. Parikh, and J. R. Forsythe, “Understanding Abrupt Wing Stall with computational fluid dynamics,” Journal of Aircraft, vol. 42, no. 3, pp. 578–585, 2005. View at: Publisher Site  Google Scholar
 J. R. Forsythe, C. M. Fremaux, and R. M. Hall, “Calculation of static and dynamic stability derivatives of the F/A18E in abrupt wing stall using RANS and DES,” in Computational Fluid Dynamics, pp. 537–542, Springer, Berlin, Germany, 2006. View at: Google Scholar
 B. Etkin, Dynamics of Atmospheric Flight, Courier Dover Publications, 2012.
 K. OrlikRückemann, “Techniques for dynamic stability testing in wind tunnels,” AGARD Dyn. Stability Parameters (SEE N 7915061 0608), 1978. View at: Google Scholar
 S. Walker, M. Tang, S. Morris, and C. Mamplata, “Falcon HTV3X—a reusable hypersonic test bed,” in Proceedings of the 15th AIAA International Space Planes and Hypersonic Systems and Technologies Conference, 2008, paper no. AIAA20082544. View at: Publisher Site  Google Scholar
 L. A. Marshall, G. P. Corpening, and R. Sherrill, “A chief engineer's view of the NASA X43A scramjet flight test,” AIAA Paper 3332, AIAA, 2005. View at: Google Scholar
 C. Bahm, E. Baumann, J. Martin, D. Bose, R. E. Beck, and B. Strovers, “The X43A hyperX mach 7 flight 2 guidance, navigation, and control overview and flight test results,” in Proceedings of the AIAA/CIRA 13th International Space Planes and Hypersonics Systems and Technologies Conference, AIAA 20053275, Capua, Italy, May 2005. View at: Google Scholar
 T. Boyd, “One hundred years of GH Bryan's stability in aviation,” Journal of Aeronautical History, paper 2011/4, 2011. View at: Google Scholar
 J. Nielsen, Missile Aerodynamics, 1988.
 G. Bryan, “Stability in aviation: an introduction to dynamical stability as applied to the motions of aëroplanes,” Nature, vol. 88, no. 2204, pp. 406–407, 1912. View at: Publisher Site  Google Scholar
 M. Goman and A. Khrabrov, “Statespace representation of aerodynamic characteristics of an aircraft at high angles of attack,” Journal of Aircraft, vol. 31, no. 5, pp. 1109–1115, 1994. View at: Publisher Site  Google Scholar
 W. R. Sears, “Some aspects of nonstationary airfoil theory and its practical application,” Journal of the Aeronautical Sciences, vol. 8, no. 3, pp. 104–108, 1941. View at: Publisher Site  Google Scholar
 M. Balajewicz, F. Nitzsche, and D. Feszty, “Application of multiinput volterra theory to nonlinear multidegreeoffreedom aerodynamic systems,” AIAA Journal, vol. 48, no. 1, pp. 56–62, 2010. View at: Publisher Site  Google Scholar
 M. Balajewicz, F. Nitzsche, and D. Feszty, “Reduced order modeling of nonlinear transonic aerodynamics using a pruned volterra series,” in Proceedings of the 50th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, Palm Springs, Calif, USA, May 2009. View at: Google Scholar
 M. Balajewicz, “Reduced order modeling of transonic flutter and limit cycle oscillations using the pruned Volterra series,” in Proceedings of the 51st AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference, April 2010. View at: Google Scholar
 H. Zhang and F. Zhuang, “NND schemes and their applications to numerical simulation of twoand threedimensional flows,” Advances in Applied Mechanics, vol. 29, pp. 193–256, 1991. View at: Publisher Site  Google Scholar
 W. Liu, H. Zhang, and H. Zhao, “Numerical simulation and physical characteristics analysis for slender wing rock,” Journal of Aircraft, vol. 43, no. 3, pp. 858–861, 2006. View at: Publisher Site  Google Scholar
 L. P. Zhang and Z. J. Wang, “A block LUSGS implicit dual timestepping algorithm for hybrid dynamic meshes,” Computers & Fluids, vol. 33, no. 7, pp. 891–916, 2004. View at: Publisher Site  Google Scholar
 A. Jameson, “Time dependent calculations using multigrid, with applications to unsteady flows past airfoils and wings,” AIAA Paper 911596, 1991. View at: Google Scholar
 V. Klein and E. A. Morelli, Aircraft System Identification: Theory and Practice, American Institute of Aeronautics and Astronautics, Reston, Va, USA, 2006.
 A. Quarteroni, R. Sacco, and F. Saleri, Numerical Mathematics, vol. 37, Springer, Berlin, Germany, 2007. View at: MathSciNet
 K. Vladislav, C. M. Patrick, J. C. Timothy, and M. B. Jay, “Analysis of wind tunnel longitudinal static and oscillatory data of the F16XL aircraft,” Tech. Rep., 1997. View at: Google Scholar
 V. Klein and K. D. Noderer, “Modeling of aircraft unsteady aerodynamic characteristics. Part 2 parameters estimated from wind tunnel data,” NASA TM 110161, 1995. View at: Google Scholar
 M. Tobak and G. T. Chapman, Nonlinear Problems in Flight Dynamics Involving Aerodynamic Bifurcations, National Aeronautics and Space Administration (NASA), 1985.
 P. H. Reisenthel, “Development of a nonlinear indicial model for maneuvering fighter aircraft,” in Proceedings of the 34th Aerospace Sciences Meeting and Exhibit, AIAA Paper 896, p. 1996, 1996. View at: Google Scholar
 W. Silva, “Identification of nonlinear aeroelastic systems based on the volterra theory: progress and opportunities,” Nonlinear Dynamics, vol. 39, no. 12, pp. 25–62, 2005. View at: Publisher Site  Google Scholar
 E. A. Morelli, “High accuracy evaluation of the finite Fourier transform using sampled data,” NASA Technical Report TME110340, 1997. View at: Google Scholar
 E. O. Brigham and R. E. Morrow, “The fast Fourier transform,” IEEE Spectrum, vol. 4, no. 12, pp. 63–70, 1967. View at: Publisher Site  Google Scholar
 L. R. Rabiner, R. W. Schafer, and C. M. Rader, “The chirp $z$transform algorithm and its application,” The Bell System Technical Journal, vol. 48, pp. 1249–1292, 1969. View at: Publisher Site  Google Scholar  MathSciNet
 R. H. Prislin, “High amplitude dynamic stability characteristics of blunt 10degree cones,” in Proceedings of the 3rd and 4th Aerospace Sciences Meeting, New York, NY, USA, 1966. View at: Google Scholar
 J. DeSpirito, S. I. Silton, and P. Weinacht, “Navierstokes predictions of dynamic stability derivatives: evaluation of steadystate methods,” Journal of Spacecraft and Rockets, vol. 46, no. 6, pp. 1142–1154, 2009. View at: Publisher Site  Google Scholar
 C. Murphy and L. Schmidt, “The effect of length on the aerodynamic characteristics of bodies of revolution in supersonic flight,” DTIC Document, 1953. View at: Google Scholar
 L. Schmidt and C. Murphy, “The aerodynamic properties of the 7Caliber ArmyNavy Spinner Rocket in transonic flight,” Tech. Rep. BRLMR775, US Army Ballistics Research Laboratory, Aberdeen Proving Ground, Md, USA, 1954. View at: Google Scholar
 P. Weinacht, W. B. Sturek, and L. B. Schiff, “Navierstokes predictions of pitch damping for axisymmetric projectiles,” Journal of Spacecraft and Rockets, vol. 34, no. 6, pp. 753–761, 1997. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2015 Xu Liu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.