# Unsteady Vibration Aerodynamic Modeling and Evaluation of Dynamic Derivatives Using Computational Fluid Dynamics

Authors: Xu Liu; Wei Liu; Yunfei Zhao
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

Unsteady aerodynamic system modeling is widely used to solve the dynamic stability problems encountering aircraft design. In this paper, single degree-of-freedom (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 Army-Navy Spinner Rocket are numerically identified by using unsteady N-S 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.


---

## Body

## 1. Introduction

In aircraft design and development, it is difficult and costly to achieve dynamic stability within the flight envelop. Shock wave-induced 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 rolling-helical modal coupling problem. In its X-33 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 X-43A 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 open-loop 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 closed-loop 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.Figure 1
Five motion states of aircraft after disturbance.
(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 instabilityThe 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 full-order model, and reduced-order model having gained rapid development over the past few years [15–17]. Table 1 lists the dynamic stability derivatives of the moment coefficient, where C
l, C
m, and C
n represent rolling, pitching, and yawing moment coefficients, respectively; α and β are the angle of attack and angle of sideslip; p, q, and r are the angular speed components of the rolling axis, pitching axis, and yawing axis.Table 1
Dynamic stability derivatives of the moment coefficient.
Rotation derivatives
Acceleration derivatives
Vibration around a fixed axis
Damping derivatives
Cross derivatives
Cross-coupling derivatives
Roll moment derivatives
C
l
p
,
C
l
q
,
C
l
r
C
l
α
˙
,
C
l
β
˙
C
l
p
+
C
l
β
˙
sin
⁡
α
C
l
r
-
C
l
β
˙
cos
⁡
α
C
l
q
+
C
l
α
˙
Pitch moment derivatives
C
m
p
,
C
m
q
,
C
m
r
C
m
α
˙
,
C
m
β
˙
C
m
q
+
C
m
α
˙
C
m
r
-
C
m
β
˙
cos
⁡
α
C
m
p
+
C
m
β
˙
sin
⁡
α
Yaw moment derivatives
C
n
p
,
C
n
q
,
C
n
r
C
n
α
˙
,
C
n
β
˙
C
n
r
-
C
n
β
˙
cos
⁡
α
C
n
p
+
C
n
β
˙
sin
⁡
α
C
n
q
+
C
n
α
˙In our study, single degree-of-freedom (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 Army-Navy Spinner Rocket, dynamic derivatives at different positions of center of mass are identified by solving three-dimensional unsteady N-S 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 three-dimensional dimensionless Navier-Stokes equation for perfect gases with ignored mass forces is as follows:(1)
∂
Q
∂
t
+
∂
E
∂
x
+
∂
F
∂
y
+
∂
G
∂
z
=
∂
E
v
∂
x
+
∂
F
v
∂
y
+
∂
G
v
∂
z
,where Q is the vector of conserved variables, Q
=
ρ
,
ρ
u
,
ρ
v
,
ρ
w
,
ρ
e
T, ρ is the density, and u, v, and w are the components of velocity given by the Cartesian velocity vector U
=
u
,
v
,
w
T. The total energy per unit mass is e. The inviscid flux vectors, E, F, and G, are given by(2)
E
=
ρ
u
ρ
u
2
+
p
ρ
u
v
ρ
u
w
ρ
h
u
,
F
=
ρ
v
ρ
u
v
ρ
v
2
+
p
ρ
v
w
ρ
h
v
,
G
=
ρ
w
ρ
u
w
ρ
v
w
ρ
w
2
+
p
ρ
h
w
,while the viscous flux vectors, F
v, G
v, and H
v, contain terms of the heat flux and viscous forces exerted on the body and can be represented by(3)
E
v
=
1
R
e
0
,
τ
x
x
,
τ
x
y
,
τ
x
z
,
u
τ
x
x
+
v
τ
x
y
+
w
τ
x
z
+
k
′
∂
T
∂
x
T
F
v
=
1
R
e
0
,
τ
y
x
,
τ
y
y
,
τ
y
z
,
u
τ
y
x
+
v
τ
y
y
+
w
τ
y
z
+
k
′
∂
T
∂
y
T
G
v
=
1
R
e
0
,
τ
z
x
,
τ
z
y
,
τ
z
z
,
u
τ
z
x
+
v
τ
z
y
+
w
τ
z
z
+
k
′
∂
T
∂
z
T
.Here, τ
i
j is the stress tensor components. k
′ is the heat transfer coefficient. Consider(4)
k
′
=
1
γ
-
1
M
∞
2
μ
L
P
r
L
+
μ
T
P
r
T
.Here, γ is the specific heat ratio. P
r
L and P
r
T are the laminar and turbulent Prandtl number. M
∞ and R
e are the freestream Mach number and Reynolds number, respectively. The Spalart-Allmaras (SA) model is employed to evaluate the turbulent viscosity μ
T. The laminar viscosity, μ
L, is evaluated using Sutherland’s law:(5)
μ
L
=
T
3
/
2
1
+
c
T
+
c
,
c
=
110.4
K
T
~
∞
.Second-order upwind NND scheme finite volume [18, 19] was used to discretize the spatial derivative term of the flow-controlled equation group. The LU-SGS 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 right-hand rule. Record the centroid position at the initial time as (
x
o
1
,
y
o
1
,
z
o
1
), the aircraft centroid position at time t as (
x
o
2
,
y
o
2
,
z
o
2
), and the Euler angle as (
θ
,
ψ
,
γ
). Then, for any grid point at the initial time (
x
0
,
y
0
,
z
0
), the coordinate (
x
1
,
y
1
,
z
1
) at time t is(6)
x
1
-
x
o
2
y
1
-
y
o
2
z
1
-
z
o
2
=
M
1
ϕ
M
2
ψ
M
3
θ
x
0
-
x
o
1
y
0
-
y
o
1
z
0
-
z
o
1
,where M
i 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, far-field 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 Degree-of-Freedom 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:(7)
I
z
z
θ
¨
=
C
m
+
μ
z
θ
˙
.Here, I
z
z is the dimensionless moment of inertia; μ
z 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 C
m can be expressed as(8)
C
m
t
=
C
m
α
t
,
α
˙
t
,
α
¨
t
,
…
,
q
t
,
q
˙
t
,
q
¨
t
,
…
.Here, α
,
α
˙
,
α
¨
,
… are the derivatives of orders of the angle of attack and its time; q
,
q
˙
,
q
¨
,
… 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 high-order terms):(9)
C
m
=
C
m
α
0
θ
+
C
m
q
+
C
m
α
˙
0
θ
˙
.Here, (
C
m
α
)
0 and (
C
m
q
+
C
m
α
˙
)
0 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 second-order differential equation for aircraft SDF vibration:(10)
I
z
z
θ
¨
-
C
m
q
+
C
m
α
˙
0
+
μ
z
θ
˙
-
C
m
α
0
θ
=
0
.When giving the analytic solution to this differential equation, we regard(
C
m
q
+
C
m
α
˙
)
0 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(11)
θ
=
A
e
λ
t
cos
⁡
ω
t
+
φ
.From (11), the motion of an object is a vibration over the time cycle T
=
2
π
/
ω. However, different from simple harmonic vibration (SHV), its amplitude A
e
λ
t varies exponentially with time t as shown in Figure 2. The envelop line equation connected by angular amplitudes is(12)
θ
=
θ
0
e
λ
t
=
θ
0
e
C
m
q
+
C
m
α
˙
0
+
μ
z
t
/
2
I
z
z
.As in numerical calculation, μ
z
=
0, when(i)
(
C
m
q
+
C
m
α
˙
)
0
<
0, the amplitude attenuates representing dynamic stability as shown in Figure 2(a);
(ii)
(
C
m
q
+
C
m
α
˙
)
0
=
0, the amplitude stays unchanged representing dynamic neutral stability as shown in Figure 2(b);
(iii)
(
C
m
q
+
C
m
α
˙
)
0
>
0, the amplitude increases representing dynamic instability as shown in Figure 2(c).Figure 2
Vibration at small damping modal.
(a)
Amplitude attenuated
(b)
Amplitude unchanged
(c)
Amplitude divergedBy numerically simulating the free pitch process of aircraft, we get the time-related function curve of angle of pitch as shown in Figure3. The function expression of this curve has been obtained from the analysis above (11).Figure 3
Dynamic derivative identification curve of pitch angular vibration.Select angles of pitch from the curve one time cycle apart from each other,θ
1 and θ
2 (corresponding to t
1 and t
2, and t
2
=
t
1
+
T):(13)
θ
1
=
A
e
λ
t
1
cos
⁡
ω
t
1
+
φ
,
θ
2
=
A
e
λ
t
2
cos
⁡
ω
t
2
+
φ
=
A
e
λ
t
2
cos
⁡
ω
t
1
+
φ
+
2
π
.By dividing them, we get(14)
C
m
q
+
C
m
α
˙
0
=
2
I
z
z
T
ln
⁡
θ
2
θ
1
-
μ
z
.


### 3.2. Large Damping Modal

For the large damping modal, its characteristic equation has two unequal real roots. The general solution to the equation is(15)
θ
=
c
1
e
r
1
t
+
c
2
e
r
2
t
.Here, arbitrary constants c
1 and c
2 can be determined by the initial conditions. If static stability is included, we have r
1,2
<
0. From (15), there is a maximum of one t that makes θ
=
0. 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.Figure 4
Vibration at large damping modal.
(a)
Amplitude attenuated
(b)
Amplitude divergedDynamic stability parameters are solved using a linear regression model. Define independent variablex and dependent variable y:(16)
x
=
C
¯
m
θ
C
¯
m
θ
˙
,
y
=
C
m
t
1
-
C
m
0
C
m
t
2
-
C
m
0
⋮
C
m
t
N
-
C
m
0
.MatrixA that associates the unknown parameters and dependent variable is defined as(17)
A
=
Δ
α
1
Δ
α
˙
1
Δ
α
2
Δ
α
˙
2
⋮
⋮
Δ
α
N
Δ
α
˙
N
.Using linear regression process for solution, we get(18)
A
T
A
x
=
A
T
y
.Normally, matrixA 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 A
T
A, Gaussian elimination method, Moore-Penrose generalized inversion method, and QR factorization though Moore-Penrose generalized inversion method and QR factorization are more accurate than direct transposition and Gaussian elimination. In terms of calculation quantity, it is O
(
n
2
) for QR factorization and O
(
n
3
) for Moore-Penrose generalized inversion method.


### 3.3. Critical Damping Modal

For the critical damping modal, its characteristic equation has two unequal real roots:(19)
r
1,2
=
M
α
˙
2
I
z
z
.The general solution to the equation is(20)
θ
=
e
λ
t
c
1
+
c
2
t
,
λ
=
C
m
q
+
C
m
α
˙
0
+
μ
z
2
I
z
z
.Arbitrary constants c
1 and c
2 can be determined by the initial conditions. From (20), there is a maximum of one t that makes θ
=
0. Hence the object does not vibrate any more.


## 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(11)
θ
=
A
e
λ
t
cos
⁡
ω
t
+
φ
.From (11), the motion of an object is a vibration over the time cycle T
=
2
π
/
ω. However, different from simple harmonic vibration (SHV), its amplitude A
e
λ
t varies exponentially with time t as shown in Figure 2. The envelop line equation connected by angular amplitudes is(12)
θ
=
θ
0
e
λ
t
=
θ
0
e
C
m
q
+
C
m
α
˙
0
+
μ
z
t
/
2
I
z
z
.As in numerical calculation, μ
z
=
0, when(i)
(
C
m
q
+
C
m
α
˙
)
0
<
0, the amplitude attenuates representing dynamic stability as shown in Figure 2(a);
(ii)
(
C
m
q
+
C
m
α
˙
)
0
=
0, the amplitude stays unchanged representing dynamic neutral stability as shown in Figure 2(b);
(iii)
(
C
m
q
+
C
m
α
˙
)
0
>
0, the amplitude increases representing dynamic instability as shown in Figure 2(c).Figure 2
Vibration at small damping modal.
(a)
Amplitude attenuated
(b)
Amplitude unchanged
(c)
Amplitude divergedBy numerically simulating the free pitch process of aircraft, we get the time-related function curve of angle of pitch as shown in Figure3. The function expression of this curve has been obtained from the analysis above (11).Figure 3
Dynamic derivative identification curve of pitch angular vibration.Select angles of pitch from the curve one time cycle apart from each other,θ
1 and θ
2 (corresponding to t
1 and t
2, and t
2
=
t
1
+
T):(13)
θ
1
=
A
e
λ
t
1
cos
⁡
ω
t
1
+
φ
,
θ
2
=
A
e
λ
t
2
cos
⁡
ω
t
2
+
φ
=
A
e
λ
t
2
cos
⁡
ω
t
1
+
φ
+
2
π
.By dividing them, we get(14)
C
m
q
+
C
m
α
˙
0
=
2
I
z
z
T
ln
⁡
θ
2
θ
1
-
μ
z
.


## 3.2. Large Damping Modal

For the large damping modal, its characteristic equation has two unequal real roots. The general solution to the equation is(15)
θ
=
c
1
e
r
1
t
+
c
2
e
r
2
t
.Here, arbitrary constants c
1 and c
2 can be determined by the initial conditions. If static stability is included, we have r
1,2
<
0. From (15), there is a maximum of one t that makes θ
=
0. 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.Figure 4
Vibration at large damping modal.
(a)
Amplitude attenuated
(b)
Amplitude divergedDynamic stability parameters are solved using a linear regression model. Define independent variablex and dependent variable y:(16)
x
=
C
¯
m
θ
C
¯
m
θ
˙
,
y
=
C
m
t
1
-
C
m
0
C
m
t
2
-
C
m
0
⋮
C
m
t
N
-
C
m
0
.MatrixA that associates the unknown parameters and dependent variable is defined as(17)
A
=
Δ
α
1
Δ
α
˙
1
Δ
α
2
Δ
α
˙
2
⋮
⋮
Δ
α
N
Δ
α
˙
N
.Using linear regression process for solution, we get(18)
A
T
A
x
=
A
T
y
.Normally, matrixA 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 A
T
A, Gaussian elimination method, Moore-Penrose generalized inversion method, and QR factorization though Moore-Penrose generalized inversion method and QR factorization are more accurate than direct transposition and Gaussian elimination. In terms of calculation quantity, it is O
(
n
2
) for QR factorization and O
(
n
3
) for Moore-Penrose generalized inversion method.


## 3.3. Critical Damping Modal

For the critical damping modal, its characteristic equation has two unequal real roots:(19)
r
1,2
=
M
α
˙
2
I
z
z
.The general solution to the equation is(20)
θ
=
e
λ
t
c
1
+
c
2
t
,
λ
=
C
m
q
+
C
m
α
˙
0
+
μ
z
2
I
z
z
.Arbitrary constants c
1 and c
2 can be determined by the initial conditions. From (20), there is a maximum of one t that makes θ
=
0. Hence the object does not vibrate any more.


## 4. Forced Simple Harmonic Motion Modeling

Figures5–7 show the aircraft flight attitudes and motion profiles under three types forced SHV (pitching/plunging/flapping).Figure 5
Forced pitching vibration of aircraft.Figure 6
Forced plunging vibration of aircraft.Figure 7
Forced flapping vibration of aircraft.Forced pitching SHV is used to calculate the pitch direct damping derivativeC
m
q
+
C
m
α
˙, the vibration equation of which is(21)
α
t
=
α
0
+
θ
t
,
θ
t
=
θ
m
sin
⁡
k
t
.Here, α
0 is the initial angle of attack; θ
m is the pitch angle amplitude; k is the dimensionless reduced frequency (k
=
ω
~
L
r
e
f
/
2
V
∞, where ω
~ is the circular frequency; L
r
e
f is the reference length; V
∞ is the inflow velocity).Forced plunging SHV is used to calculate the acceleration derivative of pitch momentC
m
α
˙, the vibration equation of which is(22)
α
t
=
α
0
+
α
m
sin
⁡
k
t
,
θ
t
=
0
.Forced flapping SHV is used to calculate the rotary derivative of pitch momentC
m
q, the vibration equation of which is(23)
α
t
=
α
0
,
θ
t
=
θ
m
sin
⁡
k
t
.The superposition of plunging vibration and flapping vibration equals pitching SHV. In the next paragraphs, the pitch direct damping coefficientC
m
q
+
C
m
α
˙ 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 t as the explicit state-dependent variable. In such case, the state variables of the rolling round pitch axis motion are α and q. As the base flight state is symmetric straight linear flight and subject to very modest disturbance amplitude, only first-order dynamic derivatives are included in the calculation. By performing Taylor expansion on (8) and ignoring small quantities of higher orders, we have(24)
C
m
t
=
C
m
0
+
C
¯
m
θ
Δ
θ
+
C
¯
m
θ
˙
Δ
θ
˙
.Here(25)
C
¯
m
θ
=
C
m
α
-
k
2
C
m
α
¨
+
C
m
q
˙
+
k
4
C
m
α
⃛
+
C
m
q
⃛
-
⋯
,
C
¯
m
θ
˙
=
C
m
α
˙
+
C
m
q
-
k
2
C
m
α
⃛
+
C
m
q
¨
+
⋯
.
C
¯
m
θ and C
¯
m
θ
˙ are the in-phase and out-phase components of Δ
C
m [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 derivesC
¯
m
θ and C
¯
m
θ
˙ by calculating the first Fourier coefficient of Δ
C
m over the time history of the first cycle T:(26)
C
¯
m
θ
=
k
π
α
m
∫
t
s
t
s
+
T
C
m
t
sin
⁡
k
t
d
t
C
¯
m
θ
˙
=
1
π
α
m
∫
t
s
t
s
+
T
C
m
t
cos
⁡
k
t
d
t
.From (25) we can see that the frequency effect of the in-phase component C
¯
m
θ is significantly greater than the out-phase component C
¯
m
θ
˙. As the order of magnitude of the high-order coefficient C
m
α
¨
+
C
m
q
˙
,
C
m
α
⃛
+
C
m
q
¨
,
C
m
α
⃛
+
C
m
q
⃛
,
… quickly reduces with the increase of the order of derivative. By ignoring the influence of the high-order coefficient and frequency effect, we can approximate to(27)
C
m
θ
≈
C
m
α
,
C
m
θ
˙
≈
C
m
α
˙
+
C
m
q
.The in-phase/out-phase component approximation (27) applies to the extent that the aircraft flies at a small angle of attack under linear aerodynamic conditions, when the in-phase and out-phase 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 high-order term coefficient is not ignorable [24, 25]. However, it is noted that when high-order coefficient is needed to express aerodynamics, as the orthogonal nature of model sequence itself keeps the fundamental harmonic unchanged, in-phase and out-phase 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 in-phase/out-phase components [29]. The finite Fourier integer of a continuous time scalar function x
t within a definite time interval t
∈
0
,
t is defined as(28)
F
x
t
≡
x
-
i
ω
=
∫
0
T
x
t
e
-
i
2
π
f
t
d
t
.Here i is an imaginary unit and f is the spatial frequency.Finite Fourier integer can be calculated at the discrete frequencyf
k, which is equal spaced from zero frequency through to Nyquist frequency f
N. 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 Z 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](29)
G
i
ω
=
F
C
m
t
F
α
t
=
R
ω
e
i
ϕ
ω
.


### 4.3. Phase Method

When forcing the aircraft to make small-amplitude, low-frequency SHV, perform Fourier expansion on the unsteady pitch momentC
m:(30)
C
m
t
=
C
m
0
+
C
-
m
sin
⁡
k
t
+
τ
+
u
t
.Here C
-
m 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; u
t is the high-order harmonic component.To the extent that the unsteady calculation time is long enough and the calculation cyclen is large enough, the influence of the initial effect is ignorable and the aerodynamic force/moment reaches a periodical steady-state value.Orderingk
t
=
2
n
π
+
0.5
π
n
=
1,2
,
3
,
…, we have(31)
C
-
m
θ
=
C
-
m
sin
⁡
τ
+
1
/
2
π
θ
m
=
C
m
4
n
+
1
π
/
2
k
-
C
m
0
θ
m
.Orderingk
t
=
2
n
π
n
=
1,2
,
3
,
…, we have(32)
C
-
m
θ
˙
=
C
-
m
sin
⁡
τ
k
θ
m
=
C
m
2
n
π
/
k
-
C
m
0
k
θ
m
.


## 4.1. Integral Method

Integral method derivesC
¯
m
θ and C
¯
m
θ
˙ by calculating the first Fourier coefficient of Δ
C
m over the time history of the first cycle T:(26)
C
¯
m
θ
=
k
π
α
m
∫
t
s
t
s
+
T
C
m
t
sin
⁡
k
t
d
t
C
¯
m
θ
˙
=
1
π
α
m
∫
t
s
t
s
+
T
C
m
t
cos
⁡
k
t
d
t
.From (25) we can see that the frequency effect of the in-phase component C
¯
m
θ is significantly greater than the out-phase component C
¯
m
θ
˙. As the order of magnitude of the high-order coefficient C
m
α
¨
+
C
m
q
˙
,
C
m
α
⃛
+
C
m
q
¨
,
C
m
α
⃛
+
C
m
q
⃛
,
… quickly reduces with the increase of the order of derivative. By ignoring the influence of the high-order coefficient and frequency effect, we can approximate to(27)
C
m
θ
≈
C
m
α
,
C
m
θ
˙
≈
C
m
α
˙
+
C
m
q
.The in-phase/out-phase component approximation (27) applies to the extent that the aircraft flies at a small angle of attack under linear aerodynamic conditions, when the in-phase and out-phase 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 high-order term coefficient is not ignorable [24, 25]. However, it is noted that when high-order coefficient is needed to express aerodynamics, as the orthogonal nature of model sequence itself keeps the fundamental harmonic unchanged, in-phase and out-phase 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 in-phase/out-phase components [29]. The finite Fourier integer of a continuous time scalar function x
t within a definite time interval t
∈
0
,
t is defined as(28)
F
x
t
≡
x
-
i
ω
=
∫
0
T
x
t
e
-
i
2
π
f
t
d
t
.Here i is an imaginary unit and f is the spatial frequency.Finite Fourier integer can be calculated at the discrete frequencyf
k, which is equal spaced from zero frequency through to Nyquist frequency f
N. 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 Z 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](29)
G
i
ω
=
F
C
m
t
F
α
t
=
R
ω
e
i
ϕ
ω
.


## 4.3. Phase Method

When forcing the aircraft to make small-amplitude, low-frequency SHV, perform Fourier expansion on the unsteady pitch momentC
m:(30)
C
m
t
=
C
m
0
+
C
-
m
sin
⁡
k
t
+
τ
+
u
t
.Here C
-
m 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; u
t is the high-order harmonic component.To the extent that the unsteady calculation time is long enough and the calculation cyclen is large enough, the influence of the initial effect is ignorable and the aerodynamic force/moment reaches a periodical steady-state value.Orderingk
t
=
2
n
π
+
0.5
π
n
=
1,2
,
3
,
…, we have(31)
C
-
m
θ
=
C
-
m
sin
⁡
τ
+
1
/
2
π
θ
m
=
C
m
4
n
+
1
π
/
2
k
-
C
m
0
θ
m
.Orderingk
t
=
2
n
π
n
=
1,2
,
3
,
…, we have(32)
C
-
m
θ
˙
=
C
-
m
sin
⁡
τ
k
θ
m
=
C
m
2
n
π
/
k
-
C
m
0
k
θ
m
.


## 5. Pitching Motion Stability of Blunted Cone

### 5.1. Computational Model and Conditions

For the 10-degree blunted cone calibration model, the bottom diameter isd
=
133.35 mm and bluntness ratio is r
n
/
r
b
=
0.4. Inflow conditions are as follows: M
∞
=
4.0, Reynolds number based on bottom diameter R
e
d
=
1.0
×
1
0
6, α
0
=
4
°
,
8
°
,
1
2
°, and T
~
∞
=
56.7 K.Three groups of examples are included. The first group uses different angular amplitudes (α
0
=
4
°, α
0
=
8
°, and α
0
=
1
2
°) to examine the effect of initial angular amplitude. The second group uses different centroid positions (X
c
g
=
0.6040, X
c
g
=
0.6336, and X
c
g
=
0.6618) to examine the effect of vibration axis position. The third uses different moments of inertia (I
z
z
=
0.01, I
z
z
=
0.1, and I
z
z
=
1) 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 Figure8. The number of cells is 99
×
73
×
57 (flow direction × circumferential × normal). The minimum cell distance is 0.92
×
1
0
-
4
R
n (R
n 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 197
×
145
×
113. The total number of grids is eight times that of the previous grids.Figure 8
Mesh of a 10-degree blunted cone calibration model.
(a)
3D mesh
(b)
Rear view of meshThe 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 Figures9 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.Figure 9
Free vibration time history of a blunted cone calibration model.Figure 10
Free vibration phase curve of a blunted cone calibration model.


### 5.3. Flow Field Characteristics of Blunt Cone Configuration

Figures11 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.Figure 11
Pressure distribution cloud of symmetric plane.Figure 12
Pressure spatial distribution and surface limiting streamline.


### 5.4. Dynamic Stability of Blunt Cone Configuration

Figure13 shows the effect of initial angular amplitude α
0 on pitching damping derivative. Here, “+,” which is the computational point of our thesis, is obtained assuming I
z
z
=
0.1 and X
c
g
=
0.6040. Experiment [32] indicates that, for large bluntness ratio, 10-degree cone configuration (e.g., r
n
/
r
b
=
0.4), 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 (
C
m
α
˙
+
C
m
q
)
0 with initial angular amplitude is insignificant.Figure 13
Effect of initial angular amplitude on obtaining pitching damping derivative by free vibration approach.In Figure14, the lateral axis is the axis position of free vibration, taken as X
c
g
=
0.6040
l, X
c
g
=
0.6336
l, and X
c
g
=
0.6618
l, respectively. The longitudinal axis is the pitching damping derivative. “+” is the result of our calculation, assuming I
z
z
=
0.1, α
0
=
1
2
°. We can see that as the axis moves backward, the dynamic stability reduces. Our numerical calculation agrees well with the experiment.Figure 14
Effect of vibration axis position on obtaining pitching damping derivative by free vibration approach.Assuming the maximum angular amplitudeα
0
=
1
2
° and X
c
g
=
0.6040
l, the free pitching process of the 10-degree blunt cone is numerically simulated at I
z
z
=
0.1
,
0.01
,
0.1 and 1. The corresponding angles of pitch variations are shown in Figures 15(a)~15(d). From the profile of the curves, as I
z
z increases, the amplitude attenuation becomes slower and slower, and the vibration period extends gradually. At this time, λ and ω in (11) are(33)
λ
=
C
m
q
+
C
m
α
˙
0
+
μ
z
2
I
z
z
ω
=
C
m
q
+
C
m
α
˙
0
+
μ
z
2
I
z
z
2
+
C
m
α
0
I
z
z
.Figure 15
Free vibration configuration under increased moment of inertia.
(a)
I
z
z
=
0.001
(b)
I
z
z
=
0.01
(c)
I
z
z
=
0.1
(d)
I
z
z
=
1We can see that whenI
z
z increases, the absolute of λ reduces and causes the amplitude to attenuate more slowly. The reduction of ω also leads to a longer vibration period. When I
z
z 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 Figures15(b)–15(d); that is, I
z
z
=
0.01, I
z
z
=
0.1, and I
z
z
=
1 correspond with frequencies k
=
6.91, k
=
2.33, and k
=
0.75. 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.Table 2
Comparison of calculation results between free vibration and forced vibration.
M
∞
=
4, R
e
d
=
1.0
×
1
0
6, α
0
=
1
2
°, X
c
g
=
0.6040
Free/forced vibration
-
(
C
m
q
+
C
m
α
˙
)
0
-
(
C
m
α
)
0
Frequency 1
I
z
z
=
0.01
0.651
0.665
k
=
6.91
0.581
0.617
Frequency 2
I
z
z
=
0.1
0.696
0.563
k
=
2.33
0.621
0.569
Frequency 3
I
z
z
=
1
0.707
0.556
k
=
0.75
0.641
0.556


## 5.1. Computational Model and Conditions

For the 10-degree blunted cone calibration model, the bottom diameter isd
=
133.35 mm and bluntness ratio is r
n
/
r
b
=
0.4. Inflow conditions are as follows: M
∞
=
4.0, Reynolds number based on bottom diameter R
e
d
=
1.0
×
1
0
6, α
0
=
4
°
,
8
°
,
1
2
°, and T
~
∞
=
56.7 K.Three groups of examples are included. The first group uses different angular amplitudes (α
0
=
4
°, α
0
=
8
°, and α
0
=
1
2
°) to examine the effect of initial angular amplitude. The second group uses different centroid positions (X
c
g
=
0.6040, X
c
g
=
0.6336, and X
c
g
=
0.6618) to examine the effect of vibration axis position. The third uses different moments of inertia (I
z
z
=
0.01, I
z
z
=
0.1, and I
z
z
=
1) 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 Figure8. The number of cells is 99
×
73
×
57 (flow direction × circumferential × normal). The minimum cell distance is 0.92
×
1
0
-
4
R
n (R
n 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 197
×
145
×
113. The total number of grids is eight times that of the previous grids.Figure 8
Mesh of a 10-degree blunted cone calibration model.
(a)
3D mesh
(b)
Rear view of meshThe 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 Figures9 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.Figure 9
Free vibration time history of a blunted cone calibration model.Figure 10
Free vibration phase curve of a blunted cone calibration model.


## 5.3. Flow Field Characteristics of Blunt Cone Configuration

Figures11 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.Figure 11
Pressure distribution cloud of symmetric plane.Figure 12
Pressure spatial distribution and surface limiting streamline.


## 5.4. Dynamic Stability of Blunt Cone Configuration

Figure13 shows the effect of initial angular amplitude α
0 on pitching damping derivative. Here, “+,” which is the computational point of our thesis, is obtained assuming I
z
z
=
0.1 and X
c
g
=
0.6040. Experiment [32] indicates that, for large bluntness ratio, 10-degree cone configuration (e.g., r
n
/
r
b
=
0.4), 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 (
C
m
α
˙
+
C
m
q
)
0 with initial angular amplitude is insignificant.Figure 13
Effect of initial angular amplitude on obtaining pitching damping derivative by free vibration approach.In Figure14, the lateral axis is the axis position of free vibration, taken as X
c
g
=
0.6040
l, X
c
g
=
0.6336
l, and X
c
g
=
0.6618
l, respectively. The longitudinal axis is the pitching damping derivative. “+” is the result of our calculation, assuming I
z
z
=
0.1, α
0
=
1
2
°. We can see that as the axis moves backward, the dynamic stability reduces. Our numerical calculation agrees well with the experiment.Figure 14
Effect of vibration axis position on obtaining pitching damping derivative by free vibration approach.Assuming the maximum angular amplitudeα
0
=
1
2
° and X
c
g
=
0.6040
l, the free pitching process of the 10-degree blunt cone is numerically simulated at I
z
z
=
0.1
,
0.01
,
0.1 and 1. The corresponding angles of pitch variations are shown in Figures 15(a)~15(d). From the profile of the curves, as I
z
z increases, the amplitude attenuation becomes slower and slower, and the vibration period extends gradually. At this time, λ and ω in (11) are(33)
λ
=
C
m
q
+
C
m
α
˙
0
+
μ
z
2
I
z
z
ω
=
C
m
q
+
C
m
α
˙
0
+
μ
z
2
I
z
z
2
+
C
m
α
0
I
z
z
.Figure 15
Free vibration configuration under increased moment of inertia.
(a)
I
z
z
=
0.001
(b)
I
z
z
=
0.01
(c)
I
z
z
=
0.1
(d)
I
z
z
=
1We can see that whenI
z
z increases, the absolute of λ reduces and causes the amplitude to attenuate more slowly. The reduction of ω also leads to a longer vibration period. When I
z
z 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 Figures15(b)–15(d); that is, I
z
z
=
0.01, I
z
z
=
0.1, and I
z
z
=
1 correspond with frequencies k
=
6.91, k
=
2.33, and k
=
0.75. 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.Table 2
Comparison of calculation results between free vibration and forced vibration.
M
∞
=
4, R
e
d
=
1.0
×
1
0
6, α
0
=
1
2
°, X
c
g
=
0.6040
Free/forced vibration
-
(
C
m
q
+
C
m
α
˙
)
0
-
(
C
m
α
)
0
Frequency 1
I
z
z
=
0.01
0.651
0.665
k
=
6.91
0.581
0.617
Frequency 2
I
z
z
=
0.1
0.696
0.563
k
=
2.33
0.621
0.569
Frequency 3
I
z
z
=
1
0.707
0.556
k
=
0.75
0.641
0.556


## 6. CFD Predictions of Dynamic Derivatives for Army-Navy Spinner Rocket

### 6.1. Computational Model and Conditions

The validation example is the Army-Navy 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 5-caliber and 9-caliber ANSR calibration models, where the length of one caliber is 20 mm.Figure 16
ANSR dimensions.In our study, the acceleration derivative, rotary derivative, and combination derivative of 5-caliber and 9-caliber 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 Table3) 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).Table 3
Position of ANSR center of mass.
Length of ANSR
Position of center of mass
5
2.5
3.0
3.5
9
4.0
5.038
5.885


### 6.2. Computational Grids

The three-dimensional calculation mesh is generated by rotating the symmetric plane mesh. Figure17 shows the symmetric plane mesh near the rocket and topological structure. The numbers of cells of the 5-caliber and 9-caliber ANSR calibration meshes are 275,000 and 296,000. The scale of the first layer of normal mesh is 10
-
5
D, where D 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.Figure 17
ANSR symmetric plane mesh topology.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. Figure18 shows the time history and hysteresis loop curves of pitching moment coefficient of 5-caliber 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 4-process 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.Figure 18
Time history curve and hysteresis loop of the pitching moment coefficient.
(a)
Time history curve
(b)
Hysteresis loop


### 6.3. Dynamic Derivative Simulation Result of ANSR

Figure19 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 C
m
a
˙
,
C
N
a
˙ and combination derivative C
m
a
˙
+
C
m
q
,
C
N
a
˙
+
C
N
q. The forced vibration amplitude is one degree and the dimensionless reduced frequency k is taken as 0.1. Quasi-steady method is used to obtain the rotary derivative C
m
q
,
C
N
q. By directly adding C
m
q and acceleration derivative C
m
a
˙, we get C
m
a
˙
+
C
m
q
¯. 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.Table 4
Comparison of ANSR combination derivatives.
c.g.
C
m
a
˙
+
C
m
q
¯
C
m
a
˙
+
C
m
q
diff.
2.5
−25.71
−25.85
0.54%
3.0
−23.78
−23.79
0.04%
3.5
−24.62
−25.02
1.62%
4.0
−138.54
−138.22
0.23%
5.038
−118.94
−119.39
0.38%
5.885
−114.80
−114.92
0.10%Figure 19
ANSR model longitudinal dynamic derivative variation as a function of center of mass.
(a)
5-caliber ANSR model normal force derivative
(b)
5-caliber ANSR model pitch moment derivative
(c)
9-caliber ANSR model pitch moment derivative
(d)
9-caliber ANSR model normal force derivativeIn Figure19, 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 N-S 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 angle-of-attack 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.5-caliber 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.


## 6.1. Computational Model and Conditions

The validation example is the Army-Navy 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 5-caliber and 9-caliber ANSR calibration models, where the length of one caliber is 20 mm.Figure 16
ANSR dimensions.In our study, the acceleration derivative, rotary derivative, and combination derivative of 5-caliber and 9-caliber 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 Table3) 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).Table 3
Position of ANSR center of mass.
Length of ANSR
Position of center of mass
5
2.5
3.0
3.5
9
4.0
5.038
5.885


## 6.2. Computational Grids

The three-dimensional calculation mesh is generated by rotating the symmetric plane mesh. Figure17 shows the symmetric plane mesh near the rocket and topological structure. The numbers of cells of the 5-caliber and 9-caliber ANSR calibration meshes are 275,000 and 296,000. The scale of the first layer of normal mesh is 10
-
5
D, where D 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.Figure 17
ANSR symmetric plane mesh topology.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. Figure18 shows the time history and hysteresis loop curves of pitching moment coefficient of 5-caliber 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 4-process 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.Figure 18
Time history curve and hysteresis loop of the pitching moment coefficient.
(a)
Time history curve
(b)
Hysteresis loop


## 6.3. Dynamic Derivative Simulation Result of ANSR

Figure19 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 C
m
a
˙
,
C
N
a
˙ and combination derivative C
m
a
˙
+
C
m
q
,
C
N
a
˙
+
C
N
q. The forced vibration amplitude is one degree and the dimensionless reduced frequency k is taken as 0.1. Quasi-steady method is used to obtain the rotary derivative C
m
q
,
C
N
q. By directly adding C
m
q and acceleration derivative C
m
a
˙, we get C
m
a
˙
+
C
m
q
¯. 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.Table 4
Comparison of ANSR combination derivatives.
c.g.
C
m
a
˙
+
C
m
q
¯
C
m
a
˙
+
C
m
q
diff.
2.5
−25.71
−25.85
0.54%
3.0
−23.78
−23.79
0.04%
3.5
−24.62
−25.02
1.62%
4.0
−138.54
−138.22
0.23%
5.038
−118.94
−119.39
0.38%
5.885
−114.80
−114.92
0.10%Figure 19
ANSR model longitudinal dynamic derivative variation as a function of center of mass.
(a)
5-caliber ANSR model normal force derivative
(b)
5-caliber ANSR model pitch moment derivative
(c)
9-caliber ANSR model pitch moment derivative
(d)
9-caliber ANSR model normal force derivativeIn Figure19, 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 N-S 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 angle-of-attack 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.5-caliber 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 5-caliber and 9-caliber ANSR models for Army-Navy 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 three-dimensional unsteady N-S 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 10-degree 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.



---

*Source: 813462-2015-08-06.xml*