Computational Determination of Dynamic Stability Derivatives

Sznajder, Janusz

doi:10.2478/tar-2025-0021

Full Article

1.

INTRODUCTION

Aerodynamic derivatives quantify the sensitivity of aerodynamic forces and moments to aircraft motion disturbances (stability derivatives) or control-surface deflections (control derivatives). In linearized aircraft equations of motion, stability derivatives arise from a Taylor-series expansion of aerodynamic forces and moments with respect to the motion variables – namely the components of velocity and acceleration expressed in a chosen reference frame. Retaining only the first-order terms yields an efficient model that represents the aircraft’s aerodynamic response to instantaneous velocity perturbations (the first six terms) as well as its dependence on the history of motion through acceleration-based terms. In the literature, stability derivatives associated with translational velocity disturbances (u, v, w) are commonly referred to as “static” derivatives. Those associated with angular velocities (p, q, r) and with accelerations $(\dot{v}, \dot{w})$ (\dot v,\dot w) are termed “dynamic” derivatives because of their critical role in aircraft dynamic stability, especially in the damping of oscillatory motions triggered by flow disturbances.

The widespread and sustained use of aerodynamic derivatives in flight simulation, stability analysis, and flight-qualities evaluation – extending across many decades – reflects both the efficiency of this modelling approach and its suitability for practical applications. This includes the assessment of flight qualities required for civil and military aircraft certification. Although numerical solutions of the Navier–Stokes equations can provide fully time-accurate aerodynamic loads and thus potentially improve the fidelity of flight-mechanics simulations, such computations remain prohibitively expensive for arbitrary aircraft manoeuvres in everyday engineering practice. Aerodynamic derivatives, by contrast, can be obtained through several complementary methods, including engineering approximations, wind-tunnel experiments, and computational fluid dynamics (CFD). Among these, the determination of dynamic stability derivatives via wind-tunnel testing has become increasingly costly and technically demanding in comparison with modern CFD methods, making CFD-based numerical methods an attractive and often preferred supplement or alternative.

2.

DETERMINATION OF DYNAMIC DERIVATIVES - INERTIAL VS NON-INERTIAL REFERENCE FRAMES

Dynamic derivatives are commonly divided into rotation derivatives and acceleration derivatives, following the classification proposed in [1]. Rotation (or damping) derivatives describe the changes in aerodynamic forces and moments arising from an aircraft’s curvilinear flight at constant angle of attack or sideslip. In the aircraft body-axis system, these effects manifest as a curvature of the incoming streamlines. Acceleration derivatives, in contrast, quantify the sensitivity of aerodynamic loads to the rate of change of angle of attack or sideslip resulting from purely translational motions. Reference [1] provides an even more detailed categorization – including cross derivatives and cross-coupling derivatives – but these lie beyond the scope of the present work.

The focus of this paper is on computational techniques for determining rotation and acceleration dynamic derivatives using CFD simulations, several of which have analogues in experimental practice. These include the Forced Oscillation Method (FOM), in which the aircraft undergoes prescribed harmonic motion and the resulting unsteady aerodynamic forces and moments are extracted, and steady curvilinear-flight simulation using a non-inertial (moving) reference frame (MRF) flow computation. The latter methodology has a corresponding wind-tunnel experimental analogy based on obtaining a curvilinear flow in the wind-tunnel test section [2]. Among these methods, FOM is the most widely used in both computational studies and experimental implementations. A third technique examined in this paper is the Indicial Response Method (IRM), which determines derivatives from the time response to a step perturbation, although this approach currently has no direct experimental counterpart.

2.1

Derivation of the computational methodology exploiting the Forced Oscillations Method

For a model oscillating about its pitch axis, with simulations conducted in an inertial reference frame, the incremental longitudinal force and moment coefficients relative to steady conditions can be expressed, following [3,4], as: 1 $Δ c_{j} = c_{j_{α}} Δ α + \frac{\bar{c}}{2 U_{\infty}} + c_{j_{\dot{α}}} \dot{α} + \frac{\bar{c}}{2 U_{\infty}} c_{j_{q}} q + {(\frac{\bar{c}}{2 U_{\infty}})}^{2} c_{j_{\dot{q}}} \dot{q},$ \Delta {c_j} = {c_{{j_\alpha }}}\Delta \alpha + {{\bar c} \over {2{U_\infty }}} + {c_{{j_{\dot \alpha }}}}\dot \alpha + {{\bar c} \over {2{U_\infty }}}{c_{{j_q}}}q + {\left( {{{\bar c} \over {2{U_\infty }}}} \right)^2}{c_{{j_{\dot q}}}}\dot q, where the index j represents coefficients of longitudinal loads: lift (L), drag (D) (or normal force, axial force) and pitching moment (m), and $\bar{c}$ {\bar c} is the characteristic dimension.

Equation (1) assumes that the aircraft equations of motion are non-dimensionalized, and that definitions of dynamic derivatives are in accordance with the convention followed by Roskam [3]: $c_{j_{α}} = \frac{\partial c_{j}}{\partial α}, c_{j_{\dot{α}}} = \frac{\partial c_{j}}{\partial (\dot{α} \bar{c} / 2 U_{\infty})}, c_{j_{q}} = \frac{\partial c_{j}}{\partial (q \bar{c} / 2 U_{\infty})} .$ {c_{{j_\alpha }}} = {{\partial {c_j}} \over {\partial \alpha }},\quad {c_{{j_{\dot \alpha }}}} = {{\partial {c_j}} \over {\partial \left( {\dot \alpha \bar c/2{U_\infty }} \right)}},\quad {c_{{j_q}}} = {{\partial {c_j}} \over {\partial \left( {q\bar c/2{U_\infty }} \right)}}.

After expressing the change of angle of attack, pitching velocity and their derivatives as harmonic functions of time. $\begin{array}{l} Δ α = α_{A} \sin ω t, \\ q = \dot{α} = ω α_{A} \cos ω t, \\ \dot{q} = - ω^{2} \cos ω t, \end{array}$ \matrix{ {\Delta \alpha = {\alpha _A}\sin \omega t,} \hfill \cr {q = \dot \alpha = \omega {\alpha _A}\cos \omega t,} \hfill \cr {\dot q = - {\omega ^2}\cos \omega t,} \hfill \cr } and substituting these terms into Eq. (1) and regrouping terms multiplied by sine and cosine functions, the equation can be re-written in the form: 2 $Δ c_{j} = α_{A} {\tilde{c}}_{j_{α}} \sin ω t + α A k {\tilde{c}}_{j_{q}} \cos ω t,$ \Delta {c_j} = {\alpha _A}{{\tilde c}_{{j_\alpha }}}\sin \omega t + \alpha Ak{{\tilde c}_{{j_q}}}\cos \omega t, where 3 ${\tilde{c}}_{j_{α}} = (c_{j_{α}} - k^{2} c_{j_{q}}) is the in-phase component of Δ c_{j},$ {{\tilde c}_{{j_\alpha }}} = \left( {{c_{{j_\alpha }}} - {k^2}{c_{{j_q}}}} \right){\rm{ is the in - phase component of }}\Delta {c_j}{\rm{,}} 4 ${\tilde{c}}_{j_{q}} = (c_{j_{\dot{α}}} - c_{j_{q}}) is the out-of-phase component of Δ c_{j}, and$ {{\tilde c}_{{j_q}}} = \left( {{c_{{j_{\dot \alpha }}}} - {c_{{j_q}}}} \right){\rm{ is the out - of - phase component of }}\Delta {c_j}{\rm{, and}} 5 $k = \frac{ω \cdot c}{2 U_{\infty}} is the reduced frequency of the oscillations .$ k = {{\omega \;\cdot\;c} \over {2{U_\infty }}}{\rm{ is the reduced frequency of the oscillations}}{\rm{.}}

A similar formulation can be obtained for lateral forces and moments, with the model oscillating in the XY-plane. In order to obtain the c_{l_p} derivative the model should oscillate around the longitudinal axis.

When plotted as a function of the angle of attack, the coefficient c_j from Eq. (1) forms a hysteresis loop, after a short transition phase, as presented schematically in Fig. 1.

Dynamic derivatives can be evaluated from the out-of-phase components of the unsteady aerodynamic response, which can be extracted from time-dependent experimental or computational data in three different ways, once the solution has converged to a limit cycle:

Computing the difference Δ₁c_j between values at phase angles ωt = 2π and ωt = π (Fig. 1),
applying linear regression of the time-dependent data,
treating Eq. (2) as a Fourier series, truncated at the first term.

The first approach is the simplest one. It can be shown that: 6 $Δ_{1} c_{j} = {c_{j} |}_{ω t = 2 π} - {c_{j} |}_{ω t = π} = 2 (c_{j_{\dot{α}}} + c_{j_{q}}) α_{A} ω \frac{c}{2 U_{\infty}},$ {\Delta _1}{c_j} = {\left. {{c_j}} \right|_{\omega t = 2\pi }} - {\left. {{c_j}} \right|_{\omega t = \pi }} = 2\left( {{c_{{j_{\dot \alpha }}}} + {c_{{j_q}}}} \right){\alpha _A}\omega {c \over {2{U_\infty }}}, from which the out-of-phase term can be obtained.

In the second approach, the linear regression method, a vector function y can be written [6] as: 7 $y = β_{0} + β_{1} x_{1} + β_{2} x_{2} + \dots + β_{k} x_{k} + ε,$ {\bf{y}} = {\beta _0} + {\beta _1}{{\bf{x}}_1} + {\beta _2}{{\bf{x}}_2} + \ldots + {\beta _k}{{\bf{x}}_{\bf{k}}} + \varepsilon , where y is the vector of coefficients c_j·, x_i are independent inputs (vectors of the independent variables in Eq. (1)), β_i are unknown coefficients, and ε is the approximation error.

The independent inputs can be placed in the matrix: 8 $X = [\begin{matrix} 1 & x_{11} & \dots & x_{k 1} \\ ⋮ & ⋱ & ⋮ \\ 1 & x_{1 n} & \dots & x_{k n} \end{matrix}] .$ X = \left[ {\matrix{ 1 & {{x_{11}}} & \cdots & {{x_{k1}}} \cr \vdots & {} & \ddots & \vdots \cr 1 & {{x_{1n}}} & \cdots & {{x_{kn}}} \cr } } \right].

The vector of unknown coefficients β – the aerodynamic derivatives in Eq. (1) – can be obtained as: 9 $β = {(X X^{T})}^{- 1} (X y) .$ \beta = {\left( {{\bf{X}}{{\bf{X}}^T}} \right)^{ - 1}}({\bf{Xy}}).

Since linear-regression tools are built into most popular spreadsheet software, minimal computational effort is required to obtain the aerodynamic derivatives from computed or measured aerodynamic coefficients c_j.

The third approach is based on the observation that Eq. (2) has the form of a Fourier series truncated at the first harmonic [4]: 10 $f (t) ~ A_{0} + A_{1} \cos (\frac{2 π t}{T}) + B_{1} \sin (\frac{2 π t}{T}),$ f(t)\~{A_0} + {A_1}\cos \left( {{{2\pi t} \over T}} \right) + {B_1}\sin \left( {{{2\pi t} \over T}} \right), and the coefficient A₁ can be expressed in terms of the out-of-phase term ${\tilde{c}}_{j_{q}}$ {{\tilde c}_{{j_q}}}: 11 $A_{1} = α_{A} k {\tilde{c}}_{j_{q}} .$ {A_1} = {\alpha _A}k{{\tilde c}_{{j_q}}}.

Then the ${\tilde{c}}_{j_{q}}$ {{\tilde c}_{{j_q}}} coefficient can be expressed using the formula for Fourier coefficient A₁: 12 ${\tilde{c}}_{j_{q}} = \frac{2}{α_{A} k n_{c} T} \int_{0}^{n_{c} T} Δ c_{j} (t) \cos (ω t) d t .$ {{\tilde c}_{{j_q}}} = {2 \over {{\alpha _A}k{n_c}T}}\int_0^{{n_c}T} \Delta {c_j}(t)\cos (\omega t)dt.

In the author’s experience, all three techniques yield nearly identical results at low and moderate angles of attack, within the attached-flow regime. The relative difference between the values of ${\tilde{c}}_{j_{q}}$ {{\tilde c}_{{j_q}}} derivatives obtained from methods (a)–(c) is typically below 0.1%. Larger discrepancies arise at higher angles of attack, where flow separation develops, the unsteady aerodynamic response becomes stronger, and the hysteresis loops become significantly distorted. This behavior is illustrated later in the paper.

Equation (1) forms the basis for both experimental and computational approaches to determining dynamic stability derivatives. Note, however, that different prescribed motions are required to isolate individual derivatives. For a model oscillating about a fixed pitch axis, the coefficient ${\tilde{c}}_{j_{q}}$ {{\tilde c}_{{j_q}}} is a sum of two derivatives, as shown in Eq. (4). To determine each term separately, distinct motions must be implemented either in the experimental setup or in the numerical simulation. Specifically, evaluation of the $c_{j_{\overset{\cdot}{α}}}$ {c_{{j_{\mathop \alpha \limits^ \cdot }}}} derivatives requires constant-pitch (q=0) oscillations along the vertical axis. Conversely, determination of the C_jq derivative requires oscillations in pitch with constant α. The general rule in designing such tests is to keep all state variables fixed except the one being isolated, which dictates the appropriate translational or rotational motion of the model. More detailed guidance on this methodology is provided in [7]. A simplified numerical implementation is described in [8], where instead of applying harmonic oscillations, steady pitching and rolling motions are used to determine the pitch-damping $(c_{m_{\dot{α}}} + c_{m_{q}})$ \left( {{c_{{m_{\dot \alpha }}}} + {c_{{m_q}}}} \right) and roll damping c_{l_p} derivatives.

2.2

Derivation of the computational methodology exploiting the Indicial Response Method

Another approach that exploits unsteady-flow simulations in an inertial reference frame is the Indicial Response Method (IRM), originally presented in [9] and further detailed in [6] and [7]. In this method, the unsteady motion of the model consists of a step change in either angle of attack or pitch rate. The approach is based on representing the unsteady normal-force and pitching-moment responses to arbitrary aircraft motion through a superposition (convolution) integral: 13 $c_{j} (t) = c_{j} (0) + \int_{0}^{t} c_{j_{α}} (t - τ) \frac{d}{d τ} α (τ) d τ + \int_{0}^{t} c_{j_{q}} (t - τ) \frac{d}{d τ} q (τ) d τ .$ {c_j}(t) = {c_j}(0) + \int_0^t {{c_{{j_\alpha }}}} (t - \tau ){d \over {d\tau }}\alpha (\tau )d\tau + \int_0^t {{c_{{j_q}}}} (t - \tau ){d \over {d\tau }}q(\tau )d\tau . where the c_jα (t – τ) and c_jq (t – τ) are the unsteady indicial responses to unit step change in α and q at t = τ. As the time approaches infinity, these functions approach steady-state values of derivatives c_jα and c_jq. It is also convenient to formulate them in terms of the deficiency functions: 14 $c_{j_{α}} (t - τ) = c_{j_{α}} (\infty) - F_{j_{α}} (t - τ),$ {c_{{j_\alpha }}}(t - \tau ) = {c_{{j_\alpha }}}(\infty ) - {F_{{j_\alpha }}}(t - \tau ), 15 $c_{j_{α}} (t - τ) = c_{j_{α}} (\infty) - F_{j_{α}} (t - τ),$ {c_{{j_\alpha }}}(t - \tau ) = {c_{{j_\alpha }}}(\infty ) - {F_{{j_\alpha }}}(t - \tau ), where F_jα (t – τ) and F_jq (t – τ) are the differences between final steady-state responses and their unsteady values. Under the additional assumption of quasi-steady character of the flow, the Eq. (13) can be written in the form [6]: 16 $c_{j} (t) = c_{j} (0) + c_{j_{α}} [α - α (0)] + c_{j_{q}} (\infty) q (t) - [\int_{0}^{t} F_{j α} (t - τ) d τ] \dot{α} (t)$ {c_j}(t) = {c_j}(0) + {c_{{j_\alpha }}}[\alpha - \alpha (0)] + {c_{{j_q}}}(\infty )q(t) - \left[ {\int_0^t {{F_{j\alpha }}} (t - \tau )d\tau } \right]\dot \alpha (t)

Equation (16) is the basis for determining the dynamic derivatives c_jq and $c_{j \dot{α}}$ {c_{j\dot \alpha }} using the Indicial Response Method. The rotational derivatives c_jq can be obtained as the final, asymptotic values of the unsteady response to a unit step increase of the pitching velocity from its initial, zero value, and the $c_{j \dot{a}}$ {c_{j\dot a}} derivatives can be obtained from the deficiency function in simulation of the unsteady response to a unit step increase of the angle of attack: 17 $c_{j_{α}} = - \int_{0}^{t} F_{j_{α}} (t - τ) d τ,$ {c_{{j_\alpha }}} = - \int_0^t {{F_{{j_\alpha }}}} (t - \tau )d\tau , where index j can be used for lift, drag, normal, or axial force and pitching moment.

For the forms of motion of the model aircraft mentioned above – oscillating, translational, or with steady rotation – the flow solver must take into account the motion of the computational domain. In the present work, the ANSYS Fluent solver with its Dynamic Mesh model [10] was used. The Dynamic Mesh model allows simulation of flows in which the domain boundaries move in time. The motion may be prescribed or coupled to the flow solution. The generic transport equation for a flow variable, represented by a generic scalar ϕ and arbitrary control volume V with moving boundaries, is: 18 $\frac{d}{d t} \int_{V} ρ ϕ d V + \int_{\partial V} ρ ϕ (\vec{u} - {\vec{u}}_{g}) \cdot d \vec{A} = \int_{\partial V} Γ \nabla ϕ \cdot d \vec{A} + \int_{V} S_{ϕ} d V,$ {d \over {dt}}\int_V \rho \phi dV + \int_{\partial V} \rho \phi \left( {\vec u - {{\vec u}_g}} \right)\;\cdot\;d\vec A = \int_{\partial V} \Gamma \nabla \phi \;\cdot\;d\vec A + \int_V {{S_\phi }} dV, where:

ρ is the fluid density,

$\vec{u}$ {\vec u} is the flow velocity vector,

${\vec{u}}_{g}$ {{\vec u}_g} is the mesh velocity of the moving mesh,

Γ is the diffusion coefficient,

S_ϕ is the source term of ϕ,

∂V is the boundary of control volume V.

In this study, the Dynamic Mesh model was applied to the entire computational domain surrounding the aircraft configuration. The motion was prescribed using a user-defined function (UDF) written in C, employing the macro definitions provided by the Fluent solver.

2.3

Brief description of approaches involving flow simulations in a non-inertial reference frame

The unsteady-flow approaches described above – required for the Indicial Response Method and the Forced Oscillation Method – are computationally expensive because the solution must evolve in time until the aerodynamic coefficients reach a limit cycle or approach a steady state after an initial step change. CFD methodology also enables the determination of dynamic derivatives in a non-inertial reference frame, where the undisturbed flow is steady, offering significantly reduced computation time. Three such approaches have been described in the literature:

Steady-flow solution in a rotating (non-inertial) reference frame, in which aerodynamic coefficients are obtained for prescribed rotational velocities. This approach applies to the evaluation of rotational derivatives.
The Reduced Frequency Approach (RFA), which assumes harmonic flow excitation and response. In this method, the conservative flow variables and the terms in the flux-conservation equations are represented by a finite Fourier series in time. A forced-oscillation motion of the model is prescribed, and a limited set of frequencies is used to represent the flowfield. The transformed time-dependent conservation equations are then solved as a steady problem using conventional CFD techniques [11].
Adjoint-based sensitivity analysis, in which the adjoint method is used to compute sensitivities of the aerodynamic coefficients with respect to rotational velocity as part of a design-optimization framework [12].

In approach (a), the flow around an aircraft undergoing steady rotation (pitching, rolling, or yawing) can be computed using steady-flow techniques formulated in a non-inertial reference frame associated with the rotating model. Such capability is available, for example, in the ANSYS Fluent solver through the Moving Reference Frame (MRF) module. The basic flow equations for a steadily moving reference frame take the following form:

- continuity: 19 $\frac{\partial ρ}{\partial t} + \nabla \cdot ρ {\vec{v}}_{r} = 0,$ {{\partial \rho } \over {\partial t}} + \nabla \cdot\rho {{\vec v}_r} = 0

- conservation of momentum: 20 $\frac{\partial}{\partial t} ρ \vec{v} + \nabla \cdot (ρ \vec{v_{r}} \vec{v}) + ρ [\vec{ω} \times (\vec{v} - {\vec{v}}_{t})] = - \nabla p + \nabla \cdot \bar{\bar{τ}} + \vec{F}$ {\partial \over {\partial t}}\rho \vec v + \nabla \;\cdot\;\left( {\rho \overrightarrow {{v_r}} \vec v} \right) + \rho \left[ {\vec \omega \times \left( {\vec v - {{\vec v}_t}} \right)} \right] = - \nabla p + \nabla \;\cdot\;\bar \bar \tau + \vec F

- conservation of energy: 21 $\frac{\partial}{\partial t} ρ E + \nabla \cdot (ρ {\vec{v}}_{r} H + p {\vec{u}}_{r}) = \nabla \cdot (k \nabla T + \overset{=}{τ} \cdot \vec{v}) + S_{h},$ {\partial \over {\partial t}}\rho E + \nabla \;\cdot\;\left( {\rho {{\vec v}_r}H + p{{\vec u}_r}} \right) = \nabla \;\cdot\;(k\nabla T + \mathop \tau \limits^ = \;\cdot\;\vec v) + {S_h}, where

${\vec{v}}_{r}$ {{\vec v}_r} is relative velocity, observed in the moving reference system: ${\vec{v}}_{r} = \vec{v} - {\vec{u}}_{r},$ {{\vec v}_r} = \vec v - {{\vec u}_r},

$\vec{v}$ {\vec v} is the absolute velocity, viewed from the stationary reference system, ${\vec{u}}_{r}$ {{\vec u}_r} is the velocity of the moving reference system relative to the inertial system, ${\vec{v}}_{t}$ {{\vec v}_t} is translation velocity, and $\vec{ω}$ {\vec \omega } is angular velocity: 22 ${\vec{u}}_{r} = {\vec{v}}_{t} + ω x \vec{r} .$ {{\vec u}_r} = {{\vec v}_t} + \omega x\vec r.

In this formulation, the Coriolis and centripetal accelerations are simplified into a single term $[\vec{ω} \times (\vec{v} - {\vec{v}}_{t})] . \vec{F}$ \left[ {\vec \omega \times \left( {\vec v - {{\vec v}_t}} \right)} \right].\vec F is the external body force, e.g. from interaction with another phase, while S_h is the enthalpy source term. In the cases investigated herein, both these terms were equal to zero.

The aerodynamic derivatives with respect to rotational velocity can then be obtained using finite difference quotients of aerodynamic coefficients with respect to the rotational velocities: 23 $c_{i_{x}} = \frac{\partial c_{i}}{\partial \frac{x \cdot d}{2 U_{\infty}}} \approx \frac{Δ c_{i}}{\frac{Δ x \cdot d}{2 U_{\infty}}},$ {c_{{i_x}}} = {{\partial {c_i}} \over {\partial {{x\cdotd} \over {2{U_\infty }}}}} \approx {{\Delta {c_i}} \over {{{\Delta x\cdotd} \over {2{U_\infty }}}}}, where i is the index of an aerodynamic coefficient, x is rotational velocity (p, q, r), and d is the characteristic dimension of the configuration.

The approach described in point (b) takes advantage of representing the aerodynamic response to a forced, harmonic motion of the aircraft model using a small, predetermined number of frequency components. It has been shown that a single frequency component is sufficient to simulate angle-of-attack oscillations of a NACA 0012 airfoil in inviscid flow [11]. The computational cost of the Reduced Frequency Approach (RFA) scales approximately with N times the cost of a steady, static solution, where N is the number of samples of the flow variables per oscillation period. The method also requires storing N copies of each flow variable in the numerical scheme. In the NACA 0012 example, N = 3 for one frequency component and N = 5 for two frequency components. The RFA has been applied to determine the pitch-damping sum $(c_{m_{\dot{α}}} + c_{m_{q}})$ \left( {{c_{{m_{\dot \alpha }}}} + {c_{{m_q}}}} \right) and the roll-damping derivative for the Finner missile model and for the Standard Dynamic Model (a configuration loosely resembling the F-16). In both cases, the agreement with experimental measurements was satisfactory.

The approach described in point (c) was developed to enable efficient determination of rotational derivatives within computational frameworks originally designed for optimization problems. In this formulation, the rotational derivatives can be evaluated directly within the optimization loop of the design process, thereby improving overall efficiency. The CFD solver equations are modified to compute sensitivities with respect to rotational velocities alongside sensitivities to the design variables. To account for aircraft rotation, a Moving Grid methodology is employed to solve the Euler equations formulated in the inertial reference frame. The approach leverages efficient automatic differentiation of the solver using specialized numerical tools, allowing for rapid and accurate sensitivity generation.

3.

RESULTS OF COMPUTATIONS

In the present work, both steady and unsteady approaches for determining dynamic derivatives were tested on two model configurations. The first was the Basic Finner missile model, a widely studied configuration representative of high-velocity vehicles. The second was the SZD-9 Bocian glider, representative of low-velocity aircraft.

The geometry of the Finner model is shown in Fig. 2. It is a simple, easily reproducible cone–cylinder body with rectangular fins, making it suitable for scaling and systematic investigation. Detailed geometric parameters, including the fin radii, are provided in [13]. The center of gravity is located at 6.1d from the model nose. In the experimental studies reported in [6] and [14], the reference diameter was d = 30 mm. In missile aerodynamics, the diameter d serves as the characteristic length used in the definition of aerodynamic coefficients.

Two unsteady-flow approaches – the Forced Oscillation Method and the Indicial Response Method – and one steady-flow approach were applied. The steady-flow method employed the solution of the Navier–Stokes equations in a non-inertial reference frame attached to the rotating body.

Figure 3 presents a comparison of the values of the (pitch-damping) derivatives c_{m_q}, computed in the presented work with the experimental values reported in [14] over a wide range of Mach numbers, from 1.0 to 4.5. The experimental methodology of [14] was based on identification of aerodynamic characteristics from free-flight trajectories on aeroballistics test range, where projectiles were launched using explosive charges.

The computational methods applied in the present work included steady-flow simulations in the non-inertial reference frame (MRF in Ansys Fluent) and, for selected Mach numbers, also unsteady-flow simulations with the oscillatory method using Dynamic Mesh module in Ansys Fluent. All stability derivatives were computed in the body-fixed coordinate system with the origin at the model’s center of gravity, the X-axis directed forward along the fuselage centerline, the Y-axis in the right direction, and the Z-axis completing the system, ensuring its right-handedness. The reduced frequency of the oscillations (Eq. 5) was set at k = 0.025, consistent with values encountered in the literature (e.g., [6]).

The difference between results of the two computational approaches applied is approximately 2.5% of the MRF value at Mach numbers above 2.5, and approximately 4% at lower Mach numbers, close to 1. Such differences are expected because the FOM introduces flow unsteadiness, whereas the MRF solution represents a steady rotating flowfield. The unsteady simulations employed a nondimensional time step $\bar{Δ t} = Δ t \cdot \frac{U_{\infty}}{d} = 0.05$ \overline {\Delta t} = \Delta t\cdot{{{U_\infty }} \over d} = 0.05 (corresponding to a motion of 0.05d per time step), as in reference [6]. The steady MRF simulations involved two rotational velocities: 1 rad/s and -1 rad/s. The significant scatter in the experimental data arises from several factors, including differences in the type of projectile trajectory (with and without roll), small variations due to manufacturing tolerances, and trajectory measurement errors influencing the final results.

Figures 4 and 5 show the convergence of the FOM method with respect to decreasing simulation time step. Two rotation derivatives – c_{m_q} and c_{N_q} – were evaluated with this method at flow conditions of M = 1.0 and α = 0°. In practice, the results obtained with relative time step $\bar{Δ t} = 0.1$ \overline {\Delta t} = 0.1 can be considered converged, as halving the time step from $\bar{Δ t} = 0.1$ \overline {\Delta t} = 0.1 resulted in differences of only 0.15% for c_{m_q} and 0.12% for c_{N_q}.

Table 1 presents the convergence behavior of the FOM method, in terms of percentage change in the computed derivatives as the time step is reduced. The relative errors of derivatives serve as a useful indication of how much accuracy is lost with increasing simulation time step. Simulations with $\bar{Δ t} = 0.05$ \overline {\Delta t} = 0.05 required 2514 time steps per oscillation period, whereas simulations with $\bar{Δ t} = 0.3$ \overline {\Delta t} = 0.3 required only 420 time steps per period, illustrating the trade-off between computational cost and desired accuracy.

Table 1.

Percentage change in dynamic derivatives after decreasing the time step.

Δ t* U_∞/d	Δ% C_{m_q}	Δ% C_{N_q}
0.3	-2.47%	-1.93%
0.2	-1.08%	-0.91%
0.1	-0.15%	-0.12%
0.05	0.00%	0.00%

Figure 6 presents a comparison of the results of the FOM and IRM methods implemented in the present work for determining the $c_{m_{\dot{α}}}$ {c_{{m_{\dot \alpha }}}} derivative. Both simulations were performed using the same nondimensional time step $\bar{Δ t} = 0.05$ \overline {\Delta t} = 0.05. Very good agreement between the results of both methods is visible, except at a single point at M = 1.2, where both methods predict a sharp, discontinuous rise of the derivative. This discontinuity also appears in the experimental data shown in Figure 7, where the sum of derivatives $c_{m_{\dot{α}}}$ {c_{{m_{\dot \alpha }}}} and c_{m_q} computed in the present work is compared with the experimental results.

More insight into the properties of the computational methods can be obtained by examining an object operating at lower speeds but over a wider range of angles of attack. The second model used in this study was the two-seat glider SZD-9 Bocian, a classic glider configuration featuring high-aspect-ratio wings and typical flight speeds in the range of 60–200 km/h. Its geometry is described in detail in [15] and [16]. For this aircraft, the rolling-moment derivative c_{l_p} (the derivative of the rolling-moment coefficient with respect to roll rate) was evaluated using two approaches: the steady-flow solution in a non-inertial reference frame with steady undisturbed flow, and unsteady-flow simulations in an inertial reference frame.

For determining the rolling moment in a non-inertial reference frame under the assumption of quasi-steady flow at a nonzero angle of attack, the undisturbed flowfield can be viewed as the superposition of two components: the flow induced by rolling motion about the aircraft X-axis, and the flow associated with steady flight at the given angle of attack. Such a combined flowfield can be imposed by prescribing a rotational velocity about a point offset along the Y-axis, as illustrated schematically in Fig. 9. In ANSYS Fluent, this was implemented by specifying the appropriate rotational velocity and center-of-rotation coordinates in the Cell Zone Conditions menu after activating the Moving Reference Frame option. The axes shown in Fig. 9 correspond to the default coordinate system of the solver.

In the present work the c_{l_p} derivative was computed as a one-sided difference quotient: $c_{l_{p}} \approx \frac{Δ c_{l}}{(\frac{b}{2 U_{\infty}}) Δ p},$ {c_{{l_p}}} \approx {{\Delta {c_l}} \over {\left( {{b \over {2{U_\infty }}}} \right)\Delta p}}, where p₁ = 0.2 rad/s and p₂ = 0.

For an aircraft oscillating around its longitudinal axis at zero-angle of attack, the rolling moment coefficient can be described as: 24 $c_{l} = c_{l 0} + (\frac{l}{2 V}) c_{l_{p}} p + {(\frac{1}{2 V})}^{2} c_{l_{p}} \dot{p},$ {c_l} = {c_{l0}} + \left( {{l \over {2V}}} \right){c_{{l_p}}}p + {\left( {{1 \over {2V}}} \right)^2}{c_{{l_p}}}\dot p, where p is the rolling angular velocity and c_l0 is the rolling moment coefficient at zero rolling velocity (equal to zero for a symmetric configuration in symmetric flow).

In the present simulations, the velocity p was harmonic in time, as a result of imposing oscillation with the roll angle φ varying according to: $φ = φ_{A} \sin ω t .$ \varphi = {\varphi _A}\sin \omega t. with amplitude φ_A = 5°.

If the angle of attack is different from zero, the oscillations around the aircraft’s longitudinal axis involve the appearance of sideslip velocity, as shown in Figure 10. In this figure V_Z is a component of a free-stream velocity along the aircraft’s Z-axis, due to the non-zero angle of attack, at zero roll angle. The sideslip velocity is a component of V_Z along the aircraft’s Y_ac-axis. For Eq. (24) to describe the rolling moment properly, a compensating side-slipping motion of the model aircraft must be introduced to zero out the sideslip velocity, as postulated in [7]. In the present case, the compensating motion involved motion of the entire mesh in each time step in the direction of the Y_ac-axis: $\begin{matrix} V_{Y mesh} = V_{Z} \sin φ \cdot \cos φ, \\ V_{Z mesh} = V_{Z} \sin^{2} φ . \end{matrix}$ \matrix{ {{V_{Y{\rm{ mesh }}}} = {V_Z}\sin \varphi \;\cdot\;\cos \varphi ,} \cr {{V_{Z{\rm{ mesh }}}} = {V_Z}{{\sin }^2}\varphi .} \cr }

The reference aircraft velocity was adjusted in each time step with the components of the mesh velocity along the domain X- and Z-axes. The compensating motion was accomplished using a macro from the Dynamic Mesh module of ANSYS Fluent.

The simulations of harmonic rolling oscillations yield a hysteresis loop of the rolling-moment coefficient as a function of roll angle, and the derivative c_{l_p} derivative can be obtained using the previously introduced Equations (6), (7) and (12) for the Forced Oscillations Method.

The first step was to investigate the influence of oscillation frequency on the computed value of c_{l_p}. Simulations performed at an angle of attack of 4° showed that this influence is relatively small, as illustrated in Figure 11. In the range of reduced frequencies decreasing from 0.5 to 0.1 (corresponding to 0.6 to 0.12 Hz), the value of the derivative decreases from -2.445 to -2.392, a change of approximately 2.2%.

A comparison of the values of the c_{l_p} derivative obtained based on c_l values determined in the non-inertial reference frame (using the MRF method) with those computed in the inertial reference frame (using the oscillatory method) was performed for a range of angle-of-attack values of 0° – 20°. This range covers both attached-flow conditions and flow with increasing separation, eventually leading to the appearance of autorotation moment, which is an important phenomenon for flight dynamics and control. The results of the comparison are presented in Fig. 12. For low angles of attack, in the attached-flow regime, the results of the two computational methods and the results of the three result-processing techniques of the oscillatory method are identical. Differences between results of the two methods appear at angle-of-attack values where flow separation regions increase, and the results converge once again when the flow separation covers almost entire upper wing surface, as illustrated in Figure 16.

The hysteresis characteristics of the rolling-moment coefficient c_l = f (φ) at three angles of attack, 4°, 15°, 20°, shown in Figures 13, 14, and 15, help explain the differences of the result-processing methods at α = 15°. In conditions of partially separated flow, the hysteresis characteristics become highly distorted, as the local loads diverge from linear dependence on local a and $\dot{α}$ {\dot \alpha }. In the presented calculations, the phase-angle results were averaged over the number of oscillation periods simulated in the test run for given α, and the two other techniques were applied accordingly for the entire time span of the test runs.

In the absence of reference data from wind-tunnel or flight tests for the c_lp derivative at high angle of attack, it is impossible to state which of the two computational methods is more accurate in this regime. The approach based on flow simulations in non-inertial reference frame is more in agreement with the assumption of quasi-steady flow conditions than the approach involving forced simulations of the model. The results obtained with the FOM method in the α range from 10° to 18° could change if more oscillation periods were simulated, as suggested by the hysteresis behavior seen in Fig. 14, but at the cost of increasing the computational time. Despite this uncertainty, the results clearly indicate the flow conditions where caution is required when applying either of these methods, and, where the results may be accepted with greater confidence. Importantly, both methods predict the occurrence of an autorotative moment in the a range where flow separation develops.

4.

CONCLUSIONS

In this study, three computational approaches for determining the dynamic stability derivatives of aircraft and missile configurations were implemented and evaluated. The first approach, based on solving the flow equations in a non-inertial reference frame, is particularly well suited for determining rotational derivatives. The second approach employs forced oscillations of the model using the Dynamic Mesh capability of ANSYS Fluent and can be applied to both rotational and acceleration derivatives. The third approach models the indicial response following a rapid change in angle of attack and was used here to determine the acceleration derivatives of the Finner missile model.

For the damping-sum derivative $(c_{m_{\hat{α}}} + c_{m_{q}})$ \left( {{c_{{m_{\hat \alpha }}}} + {c_{{m_q}}}} \right), the computational results were within the scatter of the available experimental data. The only significant discrepancy between the Indicial Response and Forced Oscillation methods occurred at M = 1.2, where both the simulations and the experimental data exhibit a discontinuity.

Two of the methods – the non-inertial reference frame (MRF) approach and the Forced Oscillation Method – were applied to evaluate the roll-damping derivative c_{l_p}. The results agree closely at low and moderate angles of attack and begin to diverge as flow separation develops. For fully separated flow, the results of both methods converge again. The MRF-based simulations predict c¡p values consistent with a stronger autorotative moment than those obtained using the Forced Oscillation Method.

Computational Determination of Dynamic Stability Derivatives

Full Article

Paradigm

My account