Numerical Evaluation of Aircraft Aerodynamic Static and Dynamic Stability Derivatives by a Mid-Fidelity Approach

Abstract

The present study aimed to investigate the capability of mid-fidelity aerodynamic solvers in performing a preliminary evaluation of the static and dynamic stability derivatives of aircraft configurations in their design phase. In this work, the mid-fidelity aerodynamic solver DUST, which is based on the novel vortex particle method (VPM), was used to perform simulations of the static and dynamic motion conditions of the Stability And Control CONfiguration (SACCON): an unmanned combat aerial vehicle geometry developed by NATO’s Research and Technology Organisation (RTO), which is used as a benchmark test case in the literature for the evaluation of aircraft stability derivatives. Two different methods were exploited to extract the dynamic stability derivative values from the aerodynamic coefficient time histories that were calculated with DUST. The results for the mid-fidelity approach were in good agreement with the obtained experimental data, as well as with the results obtained using more demanding high-fidelity CFD simulations. This demonstrates its suitability when implemented in DUST for predicting the static and dynamic behavior of airloads in different conditions, as well as in reliably predicting the values of stability derivatives, with the advantage of requiring limited computational effort with respect to classical high-fidelity numerical approaches and the use of wind tunnel tests.

1. Introduction

Aircraft design is characterized by the need to predict the static and dynamic behavior of a vehicle, and this begins with the preliminary stage of the process. Indeed, the ability to obtain reliable and robust quantitative information is useful, as this helps to avoid the need for future redesign operations. The knowledge of stability derivatives allows us, in particular, to create linear aerodynamic models that can be implemented during the preliminary development of aircraft control systems. It also allows for the computation of the modes associated with the free response (i.e., the longitudinal short period, phugoid, roll subsidence, spiral, and dutch roll), as well as to predict the time history of the loads acting on the aircraft in more complex maneuvers. Consequently, different methods have been developed for predicting aircraft control and stability behavior, and this has been achieved using both experiments and numerical simulations.

The experimental methods used for the prediction of aerodynamic coefficients and stability derivatives essentially include wind tunnel tests and flight tests, with the latter representing the most reliable and safe method by which to obtain the required information with a higher accuracy. However, they can only be performed at an advanced stage of design since they require the availability of a working prototype of the aircraft. Wind tunnel tests also provide a high degree of accuracy, but they are also characterized by certain drawbacks such as the difficulty of calculating appropriate corrections for dynamic tests due to the presence of walls and the model supporting a sting, which could alter the airload’s evaluation. Moreover, experimental methods are characterized by the high cost associated with the manufacturing of a model or a prototype, as well as with the realization of a complex test rig (particularly for dynamic tests [1]).

A first review of various experimental methods for estimating stability derivatives was presented by K.J. Orkil-Ruckemann [2]. Starting from the 1970s, there was an increasing demand for the development of accurate methodologies to estimate these quantities. This demand arose in this historical period due to the interest in aeronautics in developing aircraft that are capable of achieving higher angles of attack during extreme maneuvers at faster speeds compared to their predecessors. With the advancement of computational techniques, the information obtained from wind tunnels began to be utilized for a different purpose. It was recognized that the most effective approach for determining stability derivatives could only be achieved through the combined use of computational fluid dynamics (CFD) and experimental techniques. Hence, wind tunnel tests began to be exploited to create large databases that are useful for validating and complementing computational code outputs. The Advisory Group for Aerospace Research and Development (AGARD), for example, started certain activities in the 1990s that aimed at obtaining large amounts of experimental data, which were then to be coupled with numerical results [3]. Schütte et al. [4] presented a collection of dynamic experiments that were useful for the validation of CFD codes. More recently, a huge experimental campaign aimed at determining the ability of computational methods to accurately predict the static and dynamic stability derivatives that were performed in the frame of NATO RTO/STO AVT-161 and AVT 201 [5]. These task groups [6] were focused on providing “best practice” procedures for predicting the static and dynamic behavior of an aircraft, where there was also a focus on non-linear effects. In particular, in the frame of the AVT-161 task group, the unmanned combat aerial vehicle (UCAV) aircraft configuration of the Stability And Control CONfiguration (SACCON) was tested at the DNW-NWB and in NASA’s 14x22-foot Langley tunnel [1]. The experimental results obtained for this aircraft configuration were considered in different studies within international projects (e.g., MEGAFLUG, AeroSUM, SikMA, and SimSAC [7]) for validating the numerical evaluations of dynamic derivatives. Further studies associated with the estimation of the stability and control characteristics of low-aspect-ratio aircraft are provided in [8], where the work of AVT Task Group 189 described how unsteady, reduced aerodynamic models can be obtained from the estimation of stability derivatives. Moreover, Pfnür et al. [9] concluded that the damping derivatives of the SAGITTA flying wing, which were analyzed through the aerodynamic responses to oscillations, revealed non-linearities and varying stability across different freestream conditions, with the vertical tail notably enhancing the yaw damping. Bergmann et al. [10] presented a comprehensive historical overview focusing on the experimental methods employed by the DNW-NWB (the low-speed wind tunnel in Braunschweig, Germany). The article also described the implementation of an oscillatory motion test rig during dynamic tests. In particular, the rolling and spinning derivative support (RDT) system—which allows for a continuous rotation motion around the wind longitudinal axis—and the innovative model-positioning mechanism (MPM)—a six-degrees-of-freedom mechanism with parallel actuation—were described in detail.

Semi-empirical lower-order models were found to be the simplest computational models used when considering computational methods for the calculation of static and dynamic stability derivatives (even if they are characterized by the lowest degree of fidelity). These methods are based on the use of empirical formulas coupled with simple aerodynamic theories (such as strip theory and 2D aerodynamic theories), and they allow for a direct estimation of the stability derivatives of the aircraft. For example, the USAF Stability and Control DATCOM [11] is a database collection of methods for estimating stability and control parameters. Modern computational techniques provide different solutions for the appropriate estimation of stability derivatives, and they are also characterized by different levels of fidelity. These techniques include, for example, the use of high-fidelity computational fluid dynamics (CFD) as URANS solvers or linear methods such as vortex lattice methods (VLMs). Mialon et al. [7] presented an overview of the numerical activity performed within the SimSAC project. In particular, in that work, numerical predictions from the CFD simulations and linear methods (vortex lattice method and doublet lattice method) were compared with the experimental results on the DLR-F12 and Transonic CRuiser (TCR) configurations. The comparison with both static and dynamic test results showed that the linear methods exhibited good correlation with the CFD and experimental data, particularly at relatively low angles of attack prior to the onset of separations and the formation of complex vortical structures. J.F. Roy and S. Morgan [12] performed high-fidelity CFD simulations to validate the use of the elsA software for the calculation of aircraft static and dynamic stability derivatives when using the SACCON UCAV configuration as a benchmark. The results obtained in this latter work will be used as a comparison for the validation of the mid-fidelity approach proposed in the present paper. Moreover, Frink [13] presented a systematic study through which to guide the selection of a numerical solution strategy for the URANS computation of aircraft that are undergoing periodic motion about their body axes using the SACCON aircraft geometry. Da Ronch et al. [14] presented two test cases for calculating stability derivatives using the CFD software PMB of the University of Liverpool. The procedure for computing dynamic stability derivatives, starting from the loads’ time history, is presented in detail in this work; furthermore, two methodologies, in particular, are exploited. The first consists of computing stability derivatives as Fourier coefficients of the temporal history of force and momentum coefficients, while the second consists of performing a regression of the data to obtain the stability derivatives as the solution of a least squares problem. Nevertheless, high-fidelity CFD simulations still require a huge computational effort to accurately predict the flow phenomenology, particularly when considering the dynamic motion conditions of complex aircraft and rotorcraft configurations. These characteristics make this approach unsuitable for the preliminary design stage of a novel aircraft configuration since it requires different geometries, as well as different flight conditions, to be tested.

On the other hand, mid-fidelity methods, which have been widely investigated in the literature, have been made suitable for the calculation of stability derivatives with low computational effort. For instance, the vortex lattice method was used in the work by Tomac et al. [15] to calculate the aerodynamic characteristics of the SACCON configuration, and they compared their obtained results with higher fidelity methods. Another example is the work by Boschetti et al. [16], where an unsteady vortex lattice method was used for evaluating the ground effect in aerodynamic coefficients. Moreover, the work by Green et al. [17] showed the comparisons of the numerical and experimental evaluations of static and dynamic aerodynamic coefficients for the F-16XL aircraft (which is characterized by extreme geometries and a delta wing). In particular, in this work, three different methods for obtaining separate stability derivatives using a low-order, time-dependent panel methodand CFD were compared.

Recently, Politecnico di Milano developed a novel, flexible mid-fidelity aerodynamic computational tool called DUST, which aims at representing a fast and reliable asset for the simulation of the aerodynamics of complex aircraft configurations such as electrical vertical take off and landing (eVTOL) aircraft and tiltrotors. DUST is an open-source code that was released under an MIT license and that integrates different aerodynamic models for solid bodies as thick surface panels, thin vortex lattice elements, and lifting line elements. Moreover, a vortex particle method (VPM) was implemented for wake modeling, whereby a stable Lagrangian description of a free-vorticity flow field was provided, which is suitable for the numerical simulations of configurations that are characterized by strong aerodynamic interactions. Details about the mathematical formulation of the solver can be found in [18]. The DUST tool enables one to consider the simulations of complex geometries, ones that are characterized by the presence of multiple rotors and that create arbitrarily imposed complex motions. These features, when combined with a robust VPM for wake modeling, allowed for easier and faster aerodynamic simulations of aircraft or rotorcraft dynamic conditions when compared to CFD (this was even the case when also accurately considering the interactional aerodynamic phenomena that are typical of eVTOLs and tiltrotors [19,20]).

The present work, therefore, aims at investigating the capability of the DUST mid-fidelity aerodynamic solver in performing a preliminary evaluation of the static and dynamic stability derivatives of aircraft configurations in their design phase, as well as to understand the limits of the application of this new tool. In particular, numerical simulations were performed in this framework for both the static attitudes and dynamic motions over the SACCON configuration, which is considered a benchmark aircraft, as was conducted in [12]. The results obtained with the DUST simulations were used to compute the static and dynamic derivatives of the SACCON geometry using two different methods that have been described in the literature [14]. The accuracy of the static and dynamic derivatives obtained using the mid-fidelity numerical tool were then validated by comparison with the experimental data from [1] and with the results from Roy and Morgan [12], which were obtained using high-fidelity CFD simulations. The novelty of the present work is in showing the robustness and capability of the modern mid-fidelity numerical solver DUST with respect to evaluating a key feature in the design process of an aircraft as the calculation of aerodynamic dynamic derivatives. In particular, the advantages and limits of the present approach will be discussed via a thorough comparison with the experimental and CFD results.

This paper is organized as follows. Section 2 provides a summary of the mathematical model used to compute the dynamic stability derivatives, including both the methods described in the literature and a review of the oscillatory motions to be implemented for their calculation. Section 3 provides a brief description of DUST’s mathematical formulation. Section 4 provides a description of the numerical model that was built in DUST for the simulations of the SACCON aircraft configuration. Section 5 presents a comparison and discussion of the results that were obtained from the DUST simulations with respect to the experiments and high-fidelity CFD simulations’ evaluation available in the literature for the same test case. The conclusions are drawn in Section 6.

2. Dynamic Stability Derivative Calculation Methodology

The computation of dynamic stability derivatives was performed by imposing a sinusoidal oscillatory motion on the aircraft in its partial or complete configuration around a material point that belongs to its volume (typically, the center of gravity). The definition of the main stability derivatives along with the aircraft motions to be considered for their calculation are briefly summarized in

Table 1. Definition of the stability derivatives in the aircraft body axes (Figure 7b) and oscillatory motions to be simulated for their calculation.

	X Axis	Y Axis	Z Axis
Rolling	$C_{l_p}+C_{l_{\dot{\beta }}}sin\alpha$	$C_{Y_p}+C_{Y_{\dot{\beta }}}sin\alpha$	$C_{n_p}+C_{n_{\dot{\beta }}}sin\alpha$
Pitching	$C_{X_q}+C_{X_{\dot{\alpha }}}$	$C_{m_q}+C_{m_{\dot{\alpha }}}$	$C_{Z_q}+C_{Z_{\dot{\alpha }}}$
Yawing	$C_{l_r}-C_{l_{\dot{\beta }}}cos\alpha$	$C_{Y_r}-C_{Y_{\dot{\beta }}}cos\alpha$	$C_{n_r}-C_{n_{\dot{\beta }}}cos\alpha$
Phugoid	$C_{X_q}$	$C_{m_q}$	$C_{Z_q}$
Plunging osc.	$C_{X_{\dot{\alpha }}}$	$C_{m_{\dot{\alpha }}}$	$C_{Z_{\dot{\alpha }}}$

.

The hypothesis underlying the definition of stability derivatives is that the aerodynamic forces and moments acting on the aircraft follow a linear behavior with respect to variations of both the angle of attack and its derivative, i.e., $\alpha ,\dot{\alpha }$, the angle of sideslip and its derivative, that is, $\beta ,\dot{\beta }$, and the angular rates and accelerations in the body axes (i.e., $p,q,r,\dot{p},\dot{q}$, and $\dot{r}$). Taking this hypothesis into account, it is then possible to use a first-order Taylor series development for the description of the aerodynamic coefficients. Non-linear, higher order, frequency-dependent, and time-dependent terms were neglected. In the following mathematical description, for the sake of consistency, only longitudinal aerodynamic force coefficients were considered, while lateral and directional directions will be described in detail in Section 2.3.

As demonstrated and described in [14], the following equation can be written for longitudinal aerodynamic coefficients $C_j$ (j can be X, Z or m):

$$\begin{array}{c}\hfill \Delta C_j=C_j-C_{j_0}=C_{j_{\alpha }}\Delta \alpha +\frac{l}{v}{C}_{j_{\dot{\alpha }}}\dot{\alpha }+\frac{l}{v}{C}_{j_q}q+{\left(\frac{l}{v}\right)}^2{C}_{j_{\dot{q}}}\dot{q},\end{array}$$

(1)

where l is the characteristic dimension of the aircraft (as well as the chord for the longitudinal motions and the semi-span for the lateral-directional motions in the present work) and v is the freestream velocity. Under the validity of Equation (1), the behavior of the longitudinal aerodynamic coefficient $C_j$ with respect to the angle of attack, i.e., the $C_j-\alpha$ curve, follows an elliptical path, as shown in Figure 1. The following kinematic relations were obtained by assuming an imposed sinusoidal motion with respect to the pitch axis:

$$\left\{\begin{array}{c}\Delta \alpha ={\alpha }_{A}sin\left(\omega t\right),\hfill \\ \dot{\alpha }=q=\omega {\alpha }_{A}cos\left(\omega t\right),\hfill \\ \ddot{\alpha }=\dot{q}=-{\omega }^2{\alpha }_{A}sin\left(\omega t\right).\hfill \end{array}\right.$$

(2)

Then, by substituting the terms of (2) in Equation (1), we obtained

$$\begin{array}{c}\hfill \Delta C_j={\alpha }_A(C_{j_{\alpha }}-k^2{C}_{j_{\dot{q}}})sin\left(\omega t\right)+{\alpha }_{A}k(C_{j_{\dot{\alpha }}}+C_{j_q})cos\left(\omega t\right),\end{array}$$

(3)

where $(C_{j_{\alpha }}-k^2{C}_{j_{\dot{q}}})$ is the “in-phase” stability derivative of the aerodynamic coefficient, while $(C_{j_{\dot{\alpha }}}+C_{j_q})$ is the “out-of-phase” derivative, which can also be defined as

$$\begin{array}{c}\hfill {\overline{C}}_{j_{\alpha }}=(C_{j_{\alpha }}-k^2{C}_{j_{\dot{q}}}),\\ \hfill {\overline{C}}_{j_q}=(C_{j_{\dot{\alpha }}}+C_{j_q}).\end{array}$$

(4)

This led to obtaining the following expression for the aerodynamic coefficient:

$$\begin{array}{c}\hfill \Delta C_j={\alpha }_A{\overline{C}}_{j_{\alpha }}sin\left(\omega t\right)+{\alpha }_{A}k{\overline{C}}_{j_q}cos\left(\omega t\right).\end{array}$$

(5)

Equation (5) provides the starting point for the two methods that are used to calculate stability derivatives, i.e., the single-point method and the Fourier coefficient method, which are briefly described in the following for the sake of completeness. For a more detailed description of the mathematical methods, the readers are referred to [14].

2.1. Single-Point Method

By inverting Equation (5), when $\omega t=0$, the out-of-phase stability derivative can be derived. This operation coincides with measuring the thickness of the aerodynamic coefficient hysteresis ellipse around ${\alpha }_0$ (see Figure 1). Substituting $\omega t=0$ in Equation (5) yields

$$\begin{array}{c}\hfill {\left.\Delta C_j\left(t\right)\right|}_{t=t_{{\alpha }_0}}={\alpha }_{A}k{\overline{C}}_{j_q}.\end{array}$$

(6)

Consequently, the dynamic stability derivative can be obtained as

$$\begin{array}{c}\hfill {\overline{C}}_{j_q}=\frac{{\left.\Delta C_j\left(t\right)\right|}_{t=t_{{\alpha }_0}}}{{\alpha }_{A}k}.\end{array}$$

(7)

Figure 1 shows that the extent of the ellipse area around ${\alpha }_0$ was proportional to the combined out-of-phase stability derivative.

2.2. Fourier Coefficient Method

A second method, which is based on considering Equation (5) as a Fourier series truncated at the first term, was then investigated.

The sine–cosine form of the Fourier series, from Ref. [21], reads as follows:

$$f\left(t\right)\sim A_0+\sum _{n=1}^{\infty }\left(A_{n}cos\left(\frac{2\pi nt}{T}\right)+B_{n}sin\left(\frac{2\pi nt}{T}\right)\right)$$

(8)

A Fourier series allows us to represent a general periodic function as a sum of sinusoidal functions (sine and cosine). The Fourier series coefficients are defined as

$$\begin{array}{cc}\hfill A_0& =\frac{1}{T}{\int }_{-T/2}^{T/2}f\left(t\right)dt\hfill \\ \hfill A_n& =\frac{2}{T}{\int }_{-T/2}^{T/2}f\left(t\right)cos\left(\frac{2\pi nt}{T}\right)dt\hfill \\ \hfill B_n& =\frac{2}{T}{\int }_{-T/2}^{T/2}f\left(t\right)sin\left(\frac{2\pi nt}{T}\right)dt.\hfill \end{array}$$

(9)

The Fourier series when truncated to the first harmonic can be written as

$$f\left(t\right)\sim A_0+A_{1}cos\left(\frac{2\pi t}{T}\right)+B_{1}sin\left(\frac{2\pi t}{T}\right).$$

(10)

However, the Taylor series expansion of the aerodynamic force coefficient $C_j=f\left(t\right)$, from Reference [14], can be written as follows:

$$f\left(t\right)=C_{j_0}+{\alpha }_A{\overline{C}}_{j_{\alpha }}sin\left(\omega t\right)+{\alpha }_{A}k{\overline{C}}_{j_q}cos\left(\omega t\right).$$

(11)

This equation represents a good approximation for the aerodynamic force or moment coefficient as a function of time. It can be rewritten as

$$f\left(t\right)=C_{j_0}+{\alpha }_A{\overline{C}}_{j_{\alpha }}sin\left(\frac{2\pi t}{T}\right)+{\alpha }_{A}k{\overline{C}}_{j_q}cos\left(\frac{2\pi t}{T}\right).$$

(12)

By comparing (10) and (12), it is easy to see that the out-of-phase dynamic stability derivatives ${\overline{C}}_{j_q}$ can be calculated starting from the following:

$${\alpha }_{A}k{\overline{C}}_{j_q}=A_1.$$

(13)

When explicitly writing the expression of the Fourier coefficient $A_1$, Equation (13) becomes

$${\alpha }_{A}k{\overline{C}}_{j_q}=\frac{2}{T}{\int }_{-T/2}^{T/2}f\left(t\right)cos\left(\omega t\right)dt.$$

(14)

This can be rewritten as

$${\overline{C}}_{j_q}=\frac{2}{{\alpha }_{A}kT}{\int }_{-T/2}^{T/2}f\left(t\right)cos\left(\omega t\right)dt,$$

(15)

which can be demonstrated to be equivalent to the following expression (when considering $n_c$ the number of periods of oscillation):

$${\overline{C}}_{j_q}=\frac{2}{{\alpha }_{A}k{n}_{c}T}{\int }_0^{n_{c}T}\Delta C_j\left(t\right)cos\left(\omega t\right)dt.$$

(16)

Equation (16) is the expression used within the context of the present dissertation to compute the out-of-phase dynamic stability derivatives.

2.3. Simulated Motions for Dynamic Derivatives

As already mentioned, to calculate dynamic stability derivatives, the harmonic motions of the aircraft have to be simulated to evaluate the time histories of the aerodynamic coefficients acting on the aircraft. A detailed description of the dynamic motions simulated in the present work using the DUST software for SACCON is provided in the following. On the other hand, static stability derivatives ($C_{j_{\alpha }},C_{j_{\beta }}$) can be computed via performing aerodynamic simulations of the aircraft at various angles of attack or sideslip angles, which is then followed by calculating the slope of the $C_j-\alpha$ or $C_j-\beta$ curves. The linearity of the relationship must be verified within the context of linear aerodynamics.

2.3.1. Rolling Oscillation

It was assumed that the harmonic roll maneuver (Figure 2) was performed at a zero angle of attack and at an angle of sideslip ($\alpha ,\beta =0$). The Taylor series expansion of the lateral-directional aerodynamic coefficients, when coupled with the kinematic law of motion, reads as follows:

$$\left\{\begin{array}{c}\Delta C_j=C_j-C_{j_0}=\frac{l}{v}{C}_{j_p}p+{\left(\frac{l}{v}\right)}^2{C}_{j_{\dot{p}}}\dot{p}\hfill \\ \varphi \left(t\right)={\varphi }_{A}sin\left(\omega t\right),\hfill \\ \dot{\varphi }=p=\omega {\varphi }_{A}cos\left(\omega t\right),\hfill \\ \ddot{\varphi }=\dot{p}=-{\omega }^2{\varphi }_{A}sin\left(\omega t\right),\hfill \end{array}\right.$$

(17)

which can be rewritten as

$$\begin{array}{c}\hfill \Delta C_j=-{\varphi }_A{k}^2\left(C_{j_{\dot{p}}}\right)sin\left(\omega t\right)+{\varphi }_{A}k\left(C_{j_p}\right)cos\left(\omega t\right).\end{array}$$

(18)

Therefore, the methods described in Section 2.1 and Section 2.2 can be applied to Equation (18). Figure 2 shows a schematic of the rolling motion.

2.3.2. Pitching Oscillation

It was assumed that a harmonic pitch maneuver (Figure 3) is performed at a zero angle of sideslip ($\beta =0$). Both the angle of attack and the pitch angle were found to vary sinusoidally. The Taylor series expansion of the longitudinal aerodynamic coefficients, which was coupled with the kinematic law of motion, reads as follows:

$$\left\{\begin{array}{c}\Delta C_j=C_j-C_{j_0}=C_{j_{\alpha }}\Delta \alpha +\frac{l}{v}{C}_{j_{\dot{\alpha }}}\dot{\alpha }+\frac{l}{v}{C}_{j_q}q+{\left(\frac{l}{v}\right)}^2{C}_{j_{\dot{q}}}\dot{q}\hfill \\ \Delta \alpha =\theta \left(t\right)={\alpha }_{A}sin\left(\omega t\right),\hfill \\ \dot{\alpha }=\dot{\theta }=q=\omega {\alpha }_{A}cos\left(\omega t\right),\hfill \\ \ddot{\alpha }=\ddot{\theta }=\dot{q}=-{\omega }^2{\alpha }_{A}sin\left(\omega t\right),\hfill \end{array}\right.$$

(19)

which becomes

$$\begin{array}{c}\hfill \Delta C_j={\alpha }_A(C_{j_{\alpha }}-k^2{C}_{j_{\dot{q}}})sin\left(\omega t\right)+{\alpha }_{A}k(C_{j_{\dot{\alpha }}}+C_{j_q})cos\left(\omega t\right).\end{array}$$

(20)

to which the methods described in Section 2.1 and Section 2.2 may be applied.

It is important to underline that, from this simulation, it was possible to only obtain the combined out-of-phase dynamic derivative ($C_{j_{\dot{\alpha }}}+C_{j_q}$). To separate these two derivatives, either a plunging oscillation (Section 2.3.5) or a phugoid motion (Section 2.3.4) is required, even if (since the calculations are performed within the context of a linearized theory) the value obtained from pitching is not perfectly coincident with the sum of the values obtained from the phugoid and plunging oscillations. Figure 3 shows a schematic of the pitching motion.

2.3.3. Yaw Oscillation

A harmonic yaw motion, at a zero angle of attack ($\alpha =0$), was considered. The Taylor series expansion, which was coupled with the kinematic law of motion, in this case, reads as

$$\left\{\begin{array}{c}\Delta C_j=C_j-C_{j_0}=C_{j_{\beta }}\Delta \beta +\frac{l}{v}{C}_{j_{\dot{\beta }}}\dot{\beta }+\frac{l}{v}{C}_{j_r}r+{\left(\frac{l}{v}\right)}^2{C}_{j_{\dot{r}}}\dot{r}\hfill \\ \psi =-\Delta \beta ={\psi }_{A}sin\left(\omega t\right),\hfill \\ \dot{\psi }=-\dot{\beta }=r=\omega {\psi }_{A}cos\left(\omega t\right),\hfill \\ \ddot{\psi }=-\ddot{\beta }=\dot{r}=-{\omega }^2{\psi }_{A}sin\left(\omega t\right),\hfill \end{array}\right.$$

(21)

which can be rewritten as

$$\begin{array}{c}\hfill \Delta C_j={\psi }_A(k^2{C}_{j_{\dot{r}}}-C_{j_{\beta }})sin\left(\omega t\right)+{\psi }_{A}k(C_{j_r}-C_{j_{\dot{\beta }}})cos\left(\omega t\right).\end{array}$$

(22)

These were applied once again similarly to the methods described in Section 2.1 and Section 2.2. The simulation enables one to obtain only the combined out-of-phase dynamic derivative ($C_{j_r}-C_{j_{\dot{\beta }}}$). To separately obtain these two derivatives’ values, either a lateral oscillation simulation or a lateral phugoid motion simulation needs to be performed. Figure 4 shows a schematic of the yawing motion.

2.3.4. Phugoid Oscillation

It was assumed that a phugoid motion is performed at a zero angle of attack ($\alpha =0$). The phugoid motion is a longitudinal oscillatory motion of an aircraft. It is characterized by alternating cycles of pitch-up and pitch-down movements. During the phugoid motion, the aircraft experiences changes in altitude and pitch angle while maintaining a constant angle of attack; hence, the only derivatives that appear in the Taylor series expansion are the ones associated with the angular rate q as follows:

$$\left\{\begin{array}{c}\Delta C_j=C_j-C_{j_0}=\frac{l}{v}{C}_{j_q}q+{\left(\frac{l}{v}\right)}^2{C}_{j_{\dot{q}}}\dot{q}\hfill \\ \theta \left(t\right)={\theta }_{A}sin\left(\omega t\right),\hfill \\ \dot{\theta }=q=\omega {\theta }_{A}cos\left(\omega t\right),\hfill \\ \ddot{\theta }=\dot{q}=-{\omega }^2{\theta }_{A}sin\left(\omega t\right),\hfill \end{array}\right.$$

(23)

which becomes

$$\begin{array}{c}\hfill \Delta C_j={\theta }_A(-k^2{C}_{j_{\dot{q}}})sin\left(\omega t\right)+{\theta }_{A}k\left(C_{j_q}\right)cos\left(\omega t\right).\end{array}$$

(24)

It is to this that the methods described in Section 2.1 and Section 2.2 may be applied. This method permits one to directly evaluate the separated dynamic stability derivative $C_{j_q}$. Figure 5 shows a schematic of the phugoid motion.

It is important to note that, in the calculations derived from Equation (23), it is assumed that the pitch angle time history $\theta \left(t\right)$ is known. Usually, it is easier to express the description of the phugoid motion in terms of the amplitude of the oscillation in z, which is strictly related to the pitch angle $\theta$ by the fact that the angle of attack ${\alpha }_0$ must remain constant. If it is assumed that

$$\begin{array}{c}\hfill z\left(t\right)=z_{A}sin(\omega t+\pi /2)\\ \hfill \dot{z}\left(t\right)=\omega z_{A}cos(\omega t+\pi /2),\end{array}$$

(25)

then, if it is imposed that $\alpha ={\alpha }_0=const.$, the relation between ${\theta }_A$ and $z_A$ is

$$\begin{array}{c}\hfill {\theta }_A=-\frac{z_A\omega }{V_{inf}}.\end{array}$$

(26)

2.3.5. Plunge Oscillation

A harmonic plunging maneuver (one that features a rigid vertical translation and is, therefore, $q=0$) is defined as being performed with a zero angle of sideslip ($\beta =0$). The Taylor series expansion, which was coupled with the kinematic law of motion, was, thus, simplified as

$$\left\{\begin{array}{c}\Delta C_j=C_j-C_{j_0}=C_{j_{\alpha }}\Delta \alpha +\frac{l}{v}{C}_{j_{\dot{\alpha }}}\dot{\alpha }\hfill \\ \Delta \alpha ={\alpha }_{A}sin\left(\omega t\right),\hfill \\ \dot{\alpha }=\omega {\alpha }_{A}cos\left(\omega t\right),\hfill \\ \ddot{\alpha }=-{\omega }^2{\alpha }_{A}sin\left(\omega t\right),\hfill \end{array}\right.$$

(27)

which becomes

$$\begin{array}{c}\hfill \Delta C_j={\alpha }_A\left(C_{j_{\alpha }}\right)sin\left(\omega t\right)+{\alpha }_{A}k\left(C_{j_{\dot{\alpha }}}\right)cos\left(\omega t\right).\end{array}$$

(28)

The methods described in Section 2.1 and Section 2.2 could then again be applied. From this simulation, it was possible to obtain the separate out-phase dynamic derivative $C_{j_{\dot{\alpha }}}$. This information can be used together with the combined results of Section 2.3.2 $C_{j_{\dot{\alpha }}}+C_{j_q}$ to derive the stability derivative $C_{j_q}$. Alternatively, it was possible to calculate $C_{j_q}$ with a phugoid simulation. Figure 6 shows a schematic of the phugoid motion.

It is common to impose the value of the vertical translation $\Delta z$ when approaching a plunging simulation. Once $\Delta z$ is fixed, it is easy to verify that the angle of attack amplitude of the oscillation caused by the plunging motion is

$$\begin{array}{c}\hfill {\alpha }_A=\frac{\Delta z\omega }{V_{inf}}.\end{array}$$

(29)

3. DUST Mid-Fidelity Aerodynamic Software

DUST is an open-source software that was developed using the object-oriented paradigms of the FORTRAN standards. It was designed to solve complex aerodynamics problems with a flexible and reliable approach. The vortex particle method (VPM) is a computational technique that involves discretizing the vorticity field using Lagrangian particles. In the VPM, the velocity field is decomposed into a potential component and a rotational component. The vorticity equation is resolved using a Lagrangian approach (which will be briefly described later). In this formulation, a lifting body sheds a wake behind itself, and it starts from an identified trailing edge. Initially, the wake is modeled with panels, where each panel represents a vortex ring. After each time step, the wake is advected, and a new set of wake panels is shed from the lifting body. These wake panels can then be converted into vortex particles, thereby providing a more robust formulation of the wake.

The governing equation for the mathematical problem associated with the vortex particle method implemented in DUST was based on the assumption of an incompressible fluid. It was from the Stokes hypothesis (when applied to the constitutive equation of a Newtonian fluid) that the following form of the Navier–Stokes equations was obtained:

$$\nabla ·\mathbf{u}=0,$$

(30)

$$\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}·\nabla \mathbf{u}=-\frac{1}{\rho }\nabla p+\nu {\nabla }^2\mathbf{u}.$$

(31)

The Helmholtz decomposition theorem was then applied, with which the velocity field can be written as a sum of two different terms, i.e., an irrotational vector field ${\mathbf{u}}_{\varphi }$ and a solenoidal vector field ${\mathbf{u}}_{\psi }$, as follows:

$${\mathbf{u}}_{\varphi }=\nabla \varphi ,$$

(32)

$${\mathbf{u}}_{\psi }=\nabla \psi ,$$

(33)

where ${\mathbf{u}}_{\varphi }$ is the potential velocity, ${\mathbf{u}}_{\psi }$ is the rotational field, $\varphi$ is the scalar potential, and $\psi$ is the vector potential. Thanks to their definitions, as well as by taking the divergence and the curl of the decomposed velocity field, it was possible to obtain the Laplace equation and the Poisson equation as follows:

$$\Delta \varphi =0,$$

(34)

$$-\Delta \psi =\omega ,$$

(35)

where $\omega$ is the vorticity vector field, and under the incompressible fluid hypothesis, its form reads as

$$\frac{\partial \omega }{\partial t}+\mathbf{u}·\nabla \omega =\omega ·\nabla \mathbf{u}+\nu {\nabla }^2\omega .$$

(36)

In the context of the vortex particle method (VPM), the vorticity field was approximated by a Lagrangian discrete method, and it was discretized through vortex particles of intensity $\alpha \left(t\right)$ in position $\mathbf{r}\left(t\right)$. The Lagrangian elements were called vortex particles. The approximation of the vorticity field can be written as

$$\omega (\mathbf{r},t)=\sum _{i=1}^{N_p}{\alpha }_{\mathit{i}}\left(t\right)\xi (\mathbf{r}-{\mathbf{r}}_i\left(t\right)),$$

(37)

where $\xi \left(\mathbf{r}\right)$ is a cut-off function of the vorticity induced by the vortex particles. The evolution equations for the position and intensity of a generic particle i can be derived by substituting Equation (37) into Equation (36), thus obtaining

$$\frac{d\mathbf{r}}{dt}=\mathbf{u}({\mathbf{r}}_i\left(t\right),t),$$

(38)

$$\frac{d{\alpha }_i}{dt}=\nabla \mathbf{u}({\mathbf{r}}_i\left(t\right),t){\alpha }_{\mathit{i}}+{\nu }^{″}\Delta {\alpha }_{\mathit{i}}^{″}.$$

(39)

The first equation describes the convection of the particle via the means of the local flow velocity $\mathbf{u}\left({\mathbf{r}}_i\right)$, while the second describes the variation in time of the intensity ${\alpha }_{\mathit{i}}$, which was determined with a vortex stretching–tilting term and a vortex diffusion term. Equation (38) is the governing equation of the problem.

The rotational part of the velocity ${\mathbf{u}}_{\psi }(\mathbf{r},t)$ was produced from the vortex particles and from the line vortices of the vortex panels. Starting from Equations (35) and (38), as well as using Green’s function method, the rotational velocity can be calculated from the vorticity field as follows:

$${\mathbf{u}}_{\psi }({\mathbf{r}}_i,t)=\sum _{i=1}^{N_p}{\mathbf{K}}_{\xi }\times {\alpha }_{{\mathit{i}}_{\mathit{p}}},$$

(40)

where $\mathbf{K}(\mathbf{r},{\mathbf{r}}_0)$ is the consistent velocity kernel—that is, the Biot–Savart kernel regularized by the cut-off function $\xi \left(r\right)$—which eliminates the singularity of $\mathbf{K}(\mathbf{r},{\mathbf{r}}_0)$ when r is equal to $r_0$. This equation represents the contribution of the velocity induced by the particles. More details about the formulation can be found in [18].

The potential velocity field ${\mathbf{u}}_{\varphi }(\mathbf{r},t)$ is instead generated by the superposition of surface elementary singularities, ones that are associated with the aerodynamic element types, surface panels, vortex rings, and lifting lines, as well as the part of the wake that is modeled by panels.

4. DUST Numerical Model of the SACCON Aircraft

This section describes the numerical model built to perform the dynamic simulations with DUST for the SACCON aircraft.

SACCON is an unmanned combat aerial vehicle (UCAV), which is characterized by the configuration of a flying wing with a sweep angle of 53 degrees, and it is not equipped with horizontal and vertical empennages. The geometrical scheme, with the dimensions of the aircraft geometry, can be observed in Figure 7a. The chord at the aircraft centerline that characterizes the fuselage was $1.608$ $\mathrm{m}$, while the chord at the wing root of the wing ($c_{ref}$, which is used for the calculation of the aerodynamic coefficients for the longitudinal motions) was $0.479$ $\mathrm{m}$. The wing semi-span $b/2$ was used for the computation of the aerodynamic coefficients for lateral and directional motions. The aerodynamic model of the aircraft for the DUST simulations was built using surface panel elements.

Surface panels were used in DUST to model solid bodies with a non-negligible thickness. Triangular or quadrilateral surface panel elements were used to discretize these surfaces. A piecewise-uniform distribution of doublets and sources was associated with these elements according to a Morino-like formulation for the velocity potential. In using a uniform-intensity panel discretization and defining the intensity of the doublets and the sources on the surface panels, the discrete counterpart of the problem was written as a linear system, and it could be solved for the singularities’ intensities (implicit scheme). The resulting formulation is shown in detail in [18]. Surface panel elements do not consider viscous effects, and they are more costly with respect to the lifting line and vortex lattice elements in terms of computational effort; however, having said that, they do not consider thickness effects.

The reference system employed for the DUST simulations is depicted in Figure 7b. All the oscillating motions simulated in the present work were performed around the point of rotation, which is represented in Figure 7a, which also shows the moment reduction reference point (MRP).

A spatial sensitivity analysis of the mesh was performed by considering static simulations at a constant angle of attack. The selected flight condition for this study was characterized by a freestream velocity of 50 $\mathrm{m} {\mathrm{s}}^{-1}$ and an angle of attack of $\alpha =5^{\circ }$. Figure 8a illustrates the variation in the aerodynamic coefficients $C_z$ and $C_m$, which were computed by changing the number of the surface panel elements in the chord while the spanwise discretization was adjusted proportionally. The simulation results showed a particularly small variation in the aerodynamic coefficients between the two finer configurations; thus, the mesh with 20 elements along the wing chord was chosen as it provided the best compromise between computational time and accuracy. Figure 8b shows the final mesh considered in the present activity.

A sensitivity analysis was also performed to assess the influence of the length of the computational domain. This was conducted as DUST allows one to select the perimeter within which the vortex particles are considered to be influencingthe aerodynamic problem. The results shown in Figure 9 demonstrate that increasing the length of the domain by more than 5 m provides a negligible variation in the computed aerodynamic coefficients. Therefore, a domain length of 5 m was chosen for all the simulations performed in the present activity.

A time step sensitivity analysis was also performed by considering the simulations of the harmonic rolling motion of the aircraft. As shown by the DUST results presented in Figure 10 (which were obtained by varying the time discretization of the simulations), the solution was time step-independent, even for the especially small time steps. Consequently, the time step $\Delta t=0.00625$ s (160 time steps in one period of dynamic oscillation) was chosen for all the simulations performed in the present activity.

The computational time required to perform DUST simulations for each harmonic motion condition using the numerical model parameters selected from these studies was approximately a few minutes when using a single commercial laptop.

5. Results and Discussion

The test conditions selected for the DUST simulations included both static tests. These tests were performed by varying the angle of attack and the angle of sideslip, finalizing them to obtain static stability derivatives, as well as for the rolling, pitching, yawing, phugoid, and plunging motions (which allowed us to obtain a wide set of dynamic stability derivatives for the aircraft object of this study). The flight parameters and freestream conditions considered for the DUST simulations on the SACCON’s geometry are summarized in

Table 2. The flight parameters considered for the DUST simulations on the SACCON’s geometry.

Motion	${\mathit{\alpha }}_0$	k	$\mathit{Amp}.$	$\mathit{\varphi }$	$\mathit{\theta }$	$\mathit{\psi }$
Static $\mathit{\alpha }$	0, 3, 6, 9, 12, 15°	/	/	0°	0, 3, 6, 9, 12, 15°	0°
Static  $\mathit{\beta }$	10°	/	/	0°	10°	±0, 3, 6, 10°
Rolling	0, 5, 10°	0.070	5°	−5 to 5°	0, 5, 10°	0°
Pitching	0, 5, 10°	0.060	5°	0°	−5° to 15°	0°
Yawing	0, 5, 10°	0.060	5°	0°	0, 5, 10°	−5 to 5°
Phugoid	5°	0.060	0.6 m	0°	±4.32°	0°
Plunging	10°	0.150	0.05 m	0°	0°	0°

and

Table 3. Freestream conditions considered for the DUST simulations on the SACCON’s geometry.

Motion	${\mathit{V}}_{inf}$ [m/s]	${\mathit{\rho }}_{inf}$ [kg/m3]	$\mathit{Ma}$	$\mathit{Re}$	${\mathit{P}}_{inf}$ [Pa]
Static  $\mathit{\alpha }$	50	1.225	0.147	1.63 $\times {10}^6$	101,325
Static  $\mathit{\beta }$	50	1.225	0.147	1.63 $\times {10}^6$	101,325
Rolling	43	1.225	0.126	1.40 $\times {10}^6$	101,325
Pitching	50	1.225	0.147	1.63 $\times {10}^6$	101,325
Yawing	50	1.225	0.147	1.63 $\times {10}^6$	101,325
Phugoid	50	1.225	0.147	1.63 $\times {10}^6$	101,325
Plunging	50	1.225	0.147	1.63 $\times {10}^6$	101,325

.

The DUST simulation results were validated in the present work using experimental data, which are available from the tests at the DNW-NWB and the tests in NASA’s 14 × 22 ft Langley subsonic tunnel [1]. Finally, some of the flight conditions considered by the DUST simulations were compared with the high-fidelity simulations that were performed with the elsA CFD software and that were used by ONERA in their contribution to the AVT161 RTO group [12].

5.1. Static Stability Derivative Analysis

In this section, the static simulation results are analyzed to validate the numerical model built with DUST for the calculation of static stability derivatives. Each test point was obtained by performing a time marching simulation of 2 $\mathrm{s}$, and the integral load values were extracted at the end of the transient after the solution reached a stationary condition. The freestream velocity for all simulations was set to 50 $\mathrm{m} {\mathrm{s}}^{-1}$. For the simulations performed by varying the sideslip angle, the angle of attack was kept constant at 10°.

Figure 11 shows the aerodynamic coefficients calculated by DUST for static test conditions as a function of the angle of attack when compared to the experiments and CFD data. A good agreement with the data available in the literature was found for both the computed vertical and horizontal force coefficients in the body axes’ reference frames (i.e., $C_z$ and $C_x$), as well as for the pitching moment coefficient (i.e., $C_m$). The offset of the experimental curves from the numerical data, which was particularly evident for the pitching moment curve, can be justified by the presence of the sting, which was not corrected in the wind tunnel data. This was visible in [22], where the wind tunnel rear and ventral mountings were compared, as well as in [6,23], where the CFD models did not include the sting and the gap in the pitching moment was visible. Figure 12 shows the trend of the static stability derivatives with respect to the angle of attack computed by extracting the slope of the previous curves at each angle of attack that was tested. The static stability derivatives evaluated by the experiments and CFD simulations were, again, particularly well reproduced by the data provided by the DUST simulations, thus validating the suitability of the numerical model for this study within the linear aerodynamic region. Further validation of the suitability of the mid-fidelity numerical tool is provided by Figure 13, wherein good agreement with the pressure coefficient distributions calculated by DUST on the surface of the SACCON aircraft with respect to the experimental data from the DNW-NWB T2373 and high-fidelity numerical data was obtained from the CFD elsA. The pressure distributions were evaluated on the measurement lines, which are visible in Figure 14b. Moreover, the flow topology over the aircraft surface was also captured well by DUST with respect to CFD [12,13], as shown by the visualization that was created of the surface streamlines and pressure coefficient contours that are presented in Figure 14a and Figure 14b, respectively. Figure 14c shows the flow field and pressure contours for $\alpha ={16.83}^{\circ }$, thereby showing that DUST is not capable of capturing the complex flow topology that is generated by the configuration, which is a consequence of the leading edge separation beginning and is typical of delta wing aircraft.

Figure 15 shows the aerodynamic coefficients that were calculated via DUST for the static test conditions as a function of the sideslip angle in comparison with the experiments.

The roll moment coefficient $C_l$, lateral force coefficient $C_Y$, and yawing moment coefficient $C_n$ curves were found to be in good agreement with the experimental data that were obtained by the DNW-NWB (test T2373 [1]) experiment, even if a slight discrepancy of the curves’ slope was also found from the DUST simulation results. A partial justification for these discrepancies could be provided by the fact that SACCON is a flying wing aircraft without any vertical empennages, thus providing especially low values of lateral-directional loads; nevertheless, the errors produced by the DUST simulations for all the considered force and moment coefficients can be considered acceptable.

Figure 16 shows the static stability derivatives with respect to the sideslip angle, which were computed by extracting the slope of the previous curves at each sideslip angle that was tested. Considering the low values obtained due to the lack of vertical empennages, the agreement of the calculations that are based on the DUST simulations may provide a good preliminary indication of the stability derivative values.

5.2. Dynamic Stability Derivatives Analysis

This section presents the results that were obtained using the DUST simulation outcomes to compute the dynamic stability derivatives of the SACCON aircraft along the rotation axes. All of the DUST simulations presented in this section were conducted with a total marching time of 2 $\mathrm{s}$.

5.2.1. Roll Dynamic Oscillation

In order to calculate the dynamic stability derivatives around the roll axis, the harmonic rolling motions of the aircraft with an oscillation frequency were set to 1 $\mathrm{Hz}$, and a freestream velocity of 43 $\mathrm{m} {\mathrm{s}}^{-1}$ was considered for the DUST simulations. In particular, different sinusoidal oscillation cycles with a roll angle amplitude fixed to $\varphi$ = 5° around an angle of attack of 0°, 5°, and 10° were simulated. A visualization of the solution computed by DUST for a roll oscillation case is shown in Figure 17; there, the dissymmetry that occurred due to the rolling-motion-induced velocity on the pressure coefficient distribution is visible.

Figure 18, in particular, shows the hysteresis ellipses that were obtained for the aerodynamic coefficients of the oscillation cycles, which were calculated by DUST, with 0° and 5° mean angles of attack. A satisfactory agreement with the experimental data obtained by NASA at Langley (test T134, Run 15 [1]) was clearly visible, thus highlighting the capabilities of the mid-fidelity solver in reproducing the unsteady airloads that characterized the dynamic test conditions.

When starting from the computed hysteresis ellipses that are shown in Figure 18, the aircraft stability derivatives were extracted by using both of the methods described in Section 2.1 and Section 2.2. The results presented in Figure 19 show that, since the hypothesis of the linear behavior turned out to be valid, the aerodynamic coefficients followed an elliptical path; therefore, both methods used to calculate the dynamic derivatives provided analogous outcomes. A particularly good agreement with the data was obtained by DUST with respect to the experimental evaluation, as well as for the dynamic derivative values (particularly for those oscillations around 0° and 5° mean angles of attack). A limited discrepancy with respect to the experimental data can be observed for the roll oscillation at a mean angle of attack equal to 10° (see Figure 19 and Figure 20). This was most likely due to the onset of a leading edge separation of the delta wing, which—in the dynamic case—occurred at a lower angle of attack compared to the static case [24]. Despite this, the good correspondence between the DUST and NASA Langley experimental results confirmed the effectiveness of the mid-fidelity approach that was implemented in DUST, as well as the preliminary evaluation of the aircraft dynamic derivatives (which were obtained with low effort with respect to the experimental test activities).

5.2.2. Pitch Dynamic Oscillation

In order to calculate the dynamic stability derivatives around the pitch axis, the sinusoidal pitching motions of the aircraft that had an oscillation frequency were set to 1 $\mathrm{Hz}$, and a freestream velocity of 50 $\mathrm{m} {\mathrm{s}}^{-1}$ was considered for the DUST simulations. Different sinusoidal pitching cycles were analogously reproduced in the rolling motion simulations with an angle of attack amplitude that was fixed to $\alpha$ = 5° around a mean angle of attack of 0°, 5°, and 10°.

Figure 21 shows the hysteresis ellipses of the vertical force and pitching moment coefficients that were obtained via DUST for the pitching cycles with mean angles of attack of 5° and 10° compared with the experimental curves measured under the same motion conditions at the DNW-NWB [1].

The airload experimental hysteresis loop amplitudes were well captured via the DUST simulations for both presented pitching motion test cases. The offset of the numerical curves with respect to the experimental ones was again due to the presence of the supporting sting during the wind tunnel tests, albeit that these were not corrected in the result representation (see the wind tunnel mounting dependency study in [22] and the CFD computations without the sting in [6]).

Figure 22 shows the comparison of the dynamic stability derivatives that were computed from the DUST pitching cycle simulations using both the single-point and Fourier coefficient methods with the experimental data obtained from the DNW-NWB tests [1]. In particular, the stability derivatives presented in Figure 22 were combined derivatives, as discussed in Section 2.3.2. Indeed, due to the nature of the pitching motion, it was not possible to directly derive the stability derivatives of the force coefficient $C_j$ separately from a motion involving $\dot{\alpha }$ and q.

A good resemblance with the experimental data of the values obtained from the DUST simulations can be observed from this comparison by using both of the methods. In particular, the dynamic behavior of the vertical force coefficient was especially limited, while the pitching moment presented more accentuated dynamic characteristics. As observed in the rolling motion test case, a limited discrepancy with respect to the experimental data could be observed for the pitch oscillation with a mean angle of attack $\alpha =$ 10°. This underlined the limits of the mid-fidelity approach for the evaluation of aircraft aerodynamic performances at high angles of attack that could be characterized by local or extended flow separation regions.

5.2.3. Yaw Dynamic Oscillation

In order to calculate the dynamic stability derivative around the yaw axis, sinusoidal motions of the aircraft with an oscillation frequency set to 1 $\mathrm{Hz}$ and a freestream velocity of 50 $\mathrm{m} {\mathrm{s}}^{-1}$ were considered for the DUST simulations. The simulated motions were characterized by oscillations with fixed amplitudes of $\Delta \beta$ = 5° while keeping the angle of attack constant at ${\alpha }_0$ = 0°, 5°, and 10°.

Figure 23 shows the airload hysteresis ellipses that were obtained via the DUST simulations for the yawing motions around angles of attack of 0°, 5°, and 10°. The comparison of the experimental curves for the aerodynamic coefficients is not reported due to the lack of the tests’ data in [1]. The numerical airload curves exhibited relatively small magnitudes and hysteresis amplitudes, as was expected due to the absence of vertical empennages. However, the slopes of the curves demonstrated consistent trends when compared to the static curves shown in Section 5.1.

Figure 24 shows the yawing stability derivatives that were obtained from the DUST simulations when compared to the data obtained from the experiments at the DNW-NWBT2373 [1] under the same flight conditions. The experimental data are herein presented for the yawing motions. The lack of aircraft vertical empennages led to an increased dispersion in the results. The yawing stability derivatives obtained from the DUST simulations remained consistent with the experimental findings as their trend for all the motions considered were inside the experimental data dispersion.

5.2.4. Phugoid Dynamic Oscillation

In order to calculate the separated pitching dynamic stability derivatives, i.e., those influenced solely by the pitch rate q, a simulation of the phugoid motion (which was described in detail in Section 2.3.4) was performed using DUST. In particular, the simulation was conducted with a freestream velocity of 50 $\mathrm{m} {\mathrm{s}}^{-1}$, an oscillation frequency of 1 $\mathrm{Hz}$, a plunge oscillation amplitude of 0.6 m, and a pitch oscillation amplitude of 4.32°. The pitch oscillation amplitude was calculated to maintain a constant 5° angle of attack for the aircraft.

The aerodynamic coefficient time histories for the vertical force and pitching moment coefficients are presented in Figure 25, while Figure 26 presents the separated stability derivatives extracted from the airload coefficient hysteresis curves. The values obtained using the two previously described methods, starting from the DUST simulation outcomes, were compared with the results obtained from the high-fidelity CFD simulations (elsA [12]). However, the high-fidelity CFD simulations were performed over a quasi-steady pitching motion, which allows one to calculate (similar to the phugoid motion) the stability derivative $C_{j_q}$. Indeed, a quasi-steady pitching motion refers to an aircraft’s motion that is characterized by a constant pitch attitude over a specific period. The aerodynamic forces acting on the aircraft were assumed to be primarily determined by the instantaneous pitch angle, and the effects of acceleration and transient phenomena were ignored. This approximation allows for a simplified analysis and prediction of the aircraft’s behavior in steady-state or quasi-steady conditions during pitching maneuvers. In particular, the quasi-steady pitching flight condition tested with elsA was characterized by $M=0.17$, $Re=1.92\times {10}^6$ and a constant angle of attack $\alpha =5^{\circ }$.

Figure 26 shows that, although the simulated motions were different, the methodologies led to especially similar results for the separate pitch stability derivative. This result highlighted the potential of the mid-fidelity approach in obtaining a faster calculation of the dynamic stability derivatives with the same degree of accuracy with respect to the high-fidelity CFD.

5.2.5. Plunge Dynamic Oscillation

Finally, in order to calculate the separated dynamic stability derivatives, i.e., those influenced solely by $\dot{\alpha }$, a plunge motion of the aircraft was simulated using DUST. The considered flight condition was characterized by a freestream velocity of 50 $\mathrm{m} {\mathrm{s}}^{-1}$ and an oscillation frequency of 2.5 Hz. By imposing a vertical translation motion with an amplitude of 0.05 m, it was possible to calculate the corresponding change in the angle of attack during the period—which amounted to 0.90°. A direct comparison of the computed airloads with the experiments was possible as the same flight condition was tested at the DNW-NWB T2373 [1]. Figure 27 shows the comparison of the vertical force and pitching moment coefficients as a function of the angle of attack, which was obtained from both static and dynamic (plunging) tests. Once again, a shift of the numerical curves with respect to the experimental ones can be observed due to the presence of the supporting sting in the wind tunnel tests. Nevertheless, DUST captured the slope and the hysteresis amplitude of the experimental vertical force coefficient particularly well, while the pitching moment hysteresis amplitude was slightly underestimated.

The comparison of the separated $\dot{\alpha }$ dynamic stability derivative that was obtained from the numerical and experimental aerodynamic coefficients hysteresis curves is presented in Figure 28. The vertical force coefficient stability derivative that was computed from the DUST simulations was well inside the experimental data dispersion related to the different leading edge geometry that was used for the wind tunnel model. Meanwhile, the pitch moment coefficient stability derivative was slightly underestimated by DUST with respect to the experiments. This was probably due to the fact that, as already observed, one is close to a condition of separation at $\alpha$ = 10° in the dynamic case. Nevertheless, the correlation with the experimental values obtained from the DNW-NWB test campaign was again good. This outcome represented a particularly interesting result considering the fact that the DUST simulations only required low effort to obtain the dynamic stability derivatives (which also applied for the complex dynamic motions of the aircraft).

6. Conclusions

The present activity aimed at investigating the potential in the use of a mid-fidelity aerodynamic software for the calculation of the static and dynamic stability derivatives of an aircraft. With this aim, a well-defined geometry as the SACCON UCAV was used as the benchmark test case, which was possible thanks to the availability in the literature of a comprehensive experimental database for the evaluation of the stability derivatives of this aircraft. Simulations with the mid-fidelity solver DUST were performed to evaluate the aircraft aerodynamic coefficients for both the static attitudes and harmonic motions of the aircraft, as well as their outcomes. These were used to, respectively, compute the static and dynamic stability derivatives with two different mathematical methods.

The results that were obtained using the mid-fidelity numerical approach showed good agreement with the experimental data, and they also provided a similar level of accuracy with respect to the stability derivative values that were calculated with a high-fidelity CFD approach (which is available in the literature). A limitation of the mid-fidelity approach was found when considering the calculation of the dynamic stability derivatives for harmonic motions around a high angle of attack: a higher discrepancy within the experimental evaluation most likely occurs due to the occurrence of a consistent flow separation, which could not be accurately reproduced by a VPMpanel.

Nevertheless, the low computational effort required for this approach and the good accuracy obtained with respect to both the experiments and high-fidelity CFD open up new scenarios in the evaluation of aircraft stability derivatives. Indeed, this numerical activity was completely performed using the computing performance of a single commercial laptop for all the static and harmonic motion simulations. Thus, the lesson learned from this activity was that mid-fidelity aerodynamic solvers are in a mature state for improving the level of accuracy regarding the preliminary design of novel aircraft configurations. This is the case thanks to their capabilities of performing fast and accurate simulations, thereby allowing one to characterize the aerodynamics and flight mechanics performance parameters of complex aircraft architectures.