1.1 Propeller influence on flight dynamics stability

Early studies associated with the modeling and investigation of propeller effects in aircraft stability date from the first decades of the twentieth century. In Lanchester^{(Reference Lanchester1)}, from 1917, a whole appendix is dedicated to notes about the effects of propulsion on a flying machine. Harris^{(Reference Harris2)} develops a mathematical formulation for side forces acting on a propeller as a response to yaw and Glauert^{(Reference Glauert3,Reference Glauert4)} extends that formulation to derive other stability derivatives associated with the propeller.

Following those early investigations, in the work of Katzoff^{(Reference Katzoff5)}, experimental data obtained on the NACA full-scale wind tunnel for eight different aircraft configurations were investigated in terms of the effects of the propeller in longitudinal stability and control. Among the conclusions, it was noticed that the loss of elevator effectiveness at high angles of attack can be largely eliminated by the inclusion of propeller. Also, the rate of increase of effective downwash angle with the angle-of-attack can be considerably increased due to the influence of propellers for certain aircraft configurations, as gull-wing and parasol-wing monoplanes. Ribner^{(Reference Ribner6)} extended some previous formulations to determine the side force of propeller due to yaw by incorporating induction effects and determines an expression for the side forces based on an analogy with fins, with the effective fin area taken as the lateral projected area of the propeller plane. It was observed that besides the side force, a single propeller in yaw also experiences a pitching moment. A dual-rotating propeller develops up to one-third more side force than a single one. Also, the side forces due to angular velocities of pitch or yaw were found to be negligible for typical angular velocities that can be realised in maneuvers, with the exception of spin.

In Butler et al, the slipstream effects in a v/stol aircraft performance and stability were studied. For that, an analytical investigation was proposed and correlated to experimental data. Following a performance investigation, a study on the slipstream influence on stability and control characteristics, as well as the feasibility of use slipstream for stability augmentation. Preliminary results indicated that slipstream can be used as a potential solution for improving the dynamic stability of a tilt-wing aircraft in hovering flight.

Although one century has passed since the early studies, due to the complexity of the problem and new tools available, the investigation and modeling of propeller influence on aircraft stability are in continuous progress in recent years. Some examples of more recent efforts in this topic are the works of Jamison^{(Reference Jamison7)}, Goraj and Cichocka^{(Reference Goraj and Cichocka8)} and Bouquet^{(Reference Bouquet9)}. In Jamison^{(Reference Jamison7)}, flight tests were performed with the E-2C aircraft for two different propeller models, Hamilton-Sundstrand model 54460-1 and model NP2000. The results pointed out that substituting one propeller by the other, with all other test conditions kept the same, influences the static longitudinal aircraft stability significantly.

Goraj and Cichocka^{(Reference Goraj and Cichocka8)} studied the influence of gyroscopic effects on the stability of a light aircraft. Two types of gyroscopic effects were distinguished: weak gyroscopic effects, corresponding to maneuvers with small changes in pitch and yaw rates (classical dynamic stability), and strong gyroscopic effects, corresponding to rapid manoeuvres, with substantial pitching and yaw rates. In this last case, it was found that coupling between lateral and longitudinal degrees of freedom can be significant, potentially causing loss of control.

The influence of propeller slipstream is investigated in the work of Bouquet^{(Reference Bouquet9)}, whose focus was developing, implement, and validate a prediction method for the effects of the propeller on longitudinal stability. The implemented approach is based on Obert’s method, which was found to be computationally inexpensive and with relatively accurate results. It was observed four major effects caused by propeller slipstream: on the longitudinal stability: an additional normal force at propeller disk, the influence of slipstream on the wing, influence in tail-off pitching moment, and change in tail contribution to the pitching moment. In the case of Fokker 50, It was found a decrease in tail effectiveness due to an increase in the downwash angle at the tail with the presence of propellers.

1.2 Propeller influence on aeroelastic stability

1.3 Propeller modeling on nonlinear aeroelastic frameworks

Due to its large deflections under typical loads conditions, VFA is characterised by nonlinear behavior. However, traditional aeroelastic frameworks make use of linear theory for structural and aerodynamics formulations^{(Reference Riso, Di Vicenzo, Ritter, Cesnik and Mastroddi16)}.

Separate research group efforts have developed aeroelastic framework capable of taking into account the nonlinear aspects of the VFA behavior. Examples of state of the art codes in this context cited on literature are UM/NAST^{(Reference Shearer and Cesnik17,Reference Su and Cesnik18,Reference Su and Cesnik19)} from University of Michigan, SHARPy^{(Reference Murua, Palacios and Graham20,Reference Simpson and Palacios21,Reference Hesse, Palacios and Murua22,Reference Del Carre and Palacios23)} from London Imperial College, ASWING^(24,25,26) from Massachusetts Institute of Technology, NATASHA^{(Reference Chang27,Reference Mardanpour and Hodges28)} , from Georgia Institute of Technology, and DLR toolbox^{(Reference Ritter, Cesnik and Kruger29,Reference Ritter, Jones and Cesnik30)} , from Deutsches Zentrum für Luft-und Raumfahrt (DLR).

Concerning the modeling of the propeller effects, usually the aeroelastic frameworks for VFA consider the thrust effect only, modeled as a follower punctual force. Exceptions are, from the best of the authors’ knowledge, for ASWING and the current version of UM/NAST.

ASWING employs an actuator disk model and can model thrust, torque and P-factor loads (loads due to asymmetric loading of blades)^(25,26) .

UM/NAST originally also incorporated just thrust to model the propeller effects. However, developments made by previous works of the authors^{(Reference Teixeira and Cesnik31,Reference Teixeira and Cesnik32,Reference Teixeira and Cesnik33)} made possible the complete inclusion of the propeller effects in a nonlinear aeroelastic framework (composed of all the resultant aerodynamic loads on the hub, the aerodynamic effects of propeller slipstream, and its inertial effects). The unsteady aerodynamics is based on Lifting Line plus Viscous Vortex Particle method for the propellers combined with an Unsteady Vortex Lattice for the lifting surfaces, allowing it to capture propwash effects as well as the interaction of the flow among multiple lifting surfaces. Also presented there is the derivation of the inertial loads associated with the rotating blades and the procedure for the integration of the different involved approaches. Such previous works focused on the development of the framework itself and in the analysis of how the propeller effects impacted static and dynamic simulations of VFA. Investigations on aeroelastic stability were not performed.

1.4 Paper objectives and organization

From the previous sections, one can conclude that although many studies with rigid configurations have demonstrated that propellers can significantly influence an aircraft performance and flight dynamics stability, few investigations have been conducted on the effects of propellers on the aeroelastic stability of very flexible aircraft. Moreover, state-of-the-art nonlinear coupled aeroelastic-flight dynamics frameworks for very flexible configurations usually lack complete modeling of propeller effects, and just the thrust is typically included.

In view of such issues, this paper has the following goals:

1. (i) Investigate the impact of propeller effects on the aeroelastic stability of VFA;

2. (ii) Get insight on the isolated contribution of the different kinds of loads transmitted by the propeller.

To accomplish these objectives, an enhanced aeroelastic framework developed in previous works by the authors^{(Reference Teixeira and Cesnik31,Reference Teixeira and Cesnik32,Reference Teixeira and Cesnik33)} , and capable of incorporating the different propeller effects, is employed to generate snapshots data for a HALE aircraft model with propellers. Also, a method to extract frequencies, damping, and modes from time-series data, based on a combination of Proper Orthogonal Decomposition (POD) and system identification (sys ID), is proposed and verified for a purely structural case, for which reference data is available.

The remaining of this paper is organised as follows: Section 2 describes key aspects of the enhanced aeroelastic framework with propellers previously developed by the authors^{(Reference Teixeira and Cesnik31,Reference Teixeira and Cesnik32,Reference Teixeira and Cesnik33)} . Section 3 describes the motivation and strategy of an alternative approach to extract dynamic information from simulations of VFA with propellers. The method is then verified and applied for numerical investigations in Section 4. Finally, concluding remarks and key contributions of this work are summarised in Section 5.

2.0 AEROELASTIC FRAMEWORK WITH PROPELLER

The point of departure for the aeroelastic framework used in this work is the University of Michigan’s Nonlinear Aeroelastic Simulation Toolbox (UM/NAST)^{(Reference Shearer and Cesnik17,Reference Su and Cesnik18,Reference Su and Cesnik19)} . UM/NAST is a multidisciplinary tool developed for the simulation of very flexible aircraft. It features solutions for the modal characterisation about different steady state conditions, static, trim, and nonlinear transient simulations and aeroelastic stability analyses. For that, the equations for structural dynamics, nonlinear 6-DoF vehicle dynamics and aerodynamics are solved simultaneously.

In order to take into account the different propeller effects, the original UM/NAST framework was enhanced, as described in previous works by the authors^{(Reference Teixeira and Cesnik31,Reference Teixeira and Cesnik32,Reference Teixeira and Cesnik33)} , with an Unsteady Vortex Lattice for the lifting surfaces and a Lifting Line and Viscous Vortex Particle (LL/VVP) methods to model the propeller aerodynamics. Furthermore, inertia effects associated with the rigid rotating blades were also incorporated.

Figure 1 presents an overview of the enhanced coupled aeroelastic-flight dynamics framework, for the case of dynamic simulation. For each time step, based on current values of loads and control inputs, the new geometric configuration of the structure is determined. Beam coordinate deformations are then converted to panel grid point deformations, and aerodynamic loads are calculated with uVLM coupled to LL/VVP. Those loads are converted to concentrated loads at beam nodes, and the dynamic process continues until the predetermined simulation time is reached. It is important to note that the contribution of the elastic deformation of the body in the effective free-stream speed is considered in both uVLM and LL/VVPs.

Figure 1. Enhanced coupled aeroelastic-flight dynamics framework with propeller effects.

In the next sections, key aspects of the various components of this framework are described, focusing on the parts added by the authors in previous works^{(Reference Teixeira and Cesnik31,Reference Teixeira and Cesnik32,Reference Teixeira and Cesnik33)} . For more details about the formulation, numerical implementation, as well as verification cases of the different parts included, the reader can report to Refs (Reference Teixeira and Cesnik33,Reference Teixeira34).

2.1 Structural model

The structural model is based on a geometrically nonlinear formulation using a strain-based nonlinear finite element model of Su and Cesnik^{(Reference Su and Cesnik35)}. In this approach, constant strain in extension, twist, and in- and out-of-plane bending is assumed inside each beam element. Nonlinear equations of motion are also solved in terms of those strain values, and the related displacements are post-processed.

2.2 Lifting surfaces aerodynamics

The original aerodynamic modeling in UM/NAST framework employed a corrected strip theory based on Peters’ inflow theory^{(Reference Peters and Johnson36,Reference Peters and Cao37)} . However, as concluded in the investigation performed by Ritter et al^{(Reference Ritter, Teixeira and Cesnik38)}, although strip theory has been a popular method for the aeroelastic simulation of very flexible aircraft^{(Reference Hodges, Patil and Chae13,Reference Su and Cesnik18,Reference Shearer and Cesnik17,Reference Kima, Leeb, Leeb and Hanb39,Reference Patil, Hodges and Cesnik40,Reference Georgiou, Vio and Cooper41,Reference Weisshaar42)} , its inability to take into account mutual lifting surfaces influence and other 3D effects limit its accuracy and the possible range of applications. In order to capture the interaction of propellers and aircraft lifting surfaces, a method capable of taking into account that mutual aerodynamic influence was necessary. For that purpose, the Unsteady Vortex Lattice method (UVLM) was selected, and a UVLM code developed by Ritter^{(Reference Ritter43)} was adapted and integrated to UM/NAST framework, as an additional aerodynamic option.

Vortex Lattice has the advantage to be capable of modeling the aircraft undergoing large translations and rotations and has become a popular method in aeroelastic frameworks for the simulation of VFA^{(Reference Ritter, Teixeira and Cesnik38,Reference Murua, Palacios and Graham20,Reference Huang and Su44,Reference Teixeira, Palacios and Cesnik45,Reference Chao and Yi46,Reference Xie, Wang, Yang and Liu47,Reference Guimaraes Neto48,Reference Hesse and Palacios49,Reference Nguyen, Reynolds, Trinh and Frost50)} . Also, for conditions within the limitations and assumptions considered in its formulation, Vortex Lattice results demonstrate remarkable agreement with CFD and experimental data^{(Reference James51,Reference Ritter, Dillinger and Meddaikar52,Reference Falkner53)} . Its range of validity is typically the case of very flexible aircraft, whose usual velocities are within the incompressible regime, but with high enough Reynolds number to assume that the viscous effects are confined in a thin boundary layer.

Assuming potential flow and applying vortex ring elements distributed along panels and wake panels of lifting surfaces, the circulation at each panel is solved satisfying the no-penetration boundary condition:

(1) \begin{equation}[AIC]\,{\Gamma} = - {u_n},\end{equation}

where [AIC] is the aerodynamic influence coefficient matrix (which gives influence of one panel on another one), ${\Gamma}$ is the vector of circulation at each panel, and ${u_n}$ represents the normal component of total velocity at panel’s control points due to free-stream, body elastic deformation, mutual influence between panels, influence of the wake and any other influence of velocity field, as the presence of the propeller slipstream.

The wake is updated using a time-stepping approach. In this work, although any influence on the velocity flow field can be included, just the onflow velocity (velocity due to rigid body motion and atmospheric perturbation) was considered for the wake convection. This reduces the computational cost (and possible numerical instabilities related to wake roll-up) while still retaining enough information to capture wake-surfaces interaction. Also, it was observed in simulations of representative VFA that the effects of wake roll-up are usually not significant^{(Reference Murua, Palacios and Michael54,Reference Murua, Palacios and Michael55)} .

Also, induced and profile drags are included. In order to estimate the profile drag, the local lift coefficient, $c_l$ , of each strip is first calculated. Then, the corresponding drag coefficient, $c_d$ , is determined based on a polar table (e.g. XFOIL), matching the local value of $c_l$ with the corresponding value of $c_d$ . Rotation matrices based on UM/NAST formulation are used to rotate forces from wind axis to body system for each deformed configuration.

2.3 Propeller blade aerodynamics

In this work, the Lifting Line (LL) approach was chosen for the blades aerodynamics due to its simplicity, small computational cost, and good capability to represent unsteady behaviour of blade circulation, as observed in Abedi et al^{(Reference Abedi, Davidson and Voutsinas56)}. It should be noted that, in this context, LL is applied to model the bound vorticity only and will be coupled to another approach for the wake model, as will be discussed soon.

For each section, the local blade twist ( $\beta$ ), chord, aerofoil type, and associated polar tables (for a range of Reynolds number and angle-of-attack) are provided. In this work, data from XFOIL was used to generate the polar tables, but any other data source, such as experiment or CFD, could be used as well. XFOIL employs the $e^N$ method for transition prediction, which depends on the choice of the parameter $N_{\rm crit}$ , corresponding to the amplification factor of the most amplified frequency that causes the transition. The default value is 9, and it was applied in this work unless otherwise stated.

At each time step, based on the local flow velocity, the sectional Reynolds number and effective angle-of-attack, $\alpha$ , are calculated and used to interpolate the local aerodynamic coefficients, $c_l$ , $c_d$ and $c_m$ from the polar tables provided. The loads at each section are then calculated (using the component of total velocity aligned with the section), and the resultant loads acting on the propeller hub can be determined. Those loads are then transferred to the structural nodes where the propellers are attached.

The bound vortex segments circulation is calculated from the combination of 2D Kutta–Jukowski theory and the definition of lift coefficient:

(2) \begin{equation}\Gamma = 1/2 c V_{\rm section} c_l\end{equation}

where $\Gamma$ is the circulation at the considered control point, c is the local chord, $V_{\rm section}$ is the component of total velocity at the control point that is in the plane of aerofoil section considered, and $c_l$ is the local lift coefficient.

2.4 Propeller wake aerodynamics

There are many possibilities of propeller wake modeling available on literature, ranging from high-cost, high-fidelity solutions, as CFD approaches, to low-cost, low-fidelity solutions, as Momentum Theory, vortex filaments or panels (prescribed and free-wake), and semi-empirical approaches.

While CFD can be used to capture complex effects, as viscous effects, dynamic stall, and flow separation, and has demonstrated good capability to reproduce experimental results, it is usually too expensive, in particular in the case of multiple surfaces and propellers.

On the other hand, although the lower-fidelity approaches provide a much quicker evaluation of propeller slipstream velocities, their accuracy is limited: many of them assume a pre-determined wake shape, isolated propeller configuration, and usually capture just averaged effects. Although those methods can be enough for some applications, this is not the case of a dynamically deforming VFA, which may contain multiple surfaces and propellers.

2.4.2 Basic equation and assumptions

In the VVP, both viscosity and rotational effects are taken into account. The assumption of incompressible flow, however, is retained. From the classical Navier–Stokes equations for incompressible flow, the momentum conservation for a differential fluid element is:

(3) \begin{equation}\frac{\partial \vec{u}}{\partial t} + \left( \vec{u} \cdot \nabla \right) \vec{u} = - \frac{1}{\rho} \nabla p + \nu \nabla ^2 \vec{u},\end{equation}

where $\vec{u}(\vec{x},t)$ is the velocity field, $\omega(\vec{x},t)$ is the vorticity field, $p(\vec{x},t)$ is the pressure field and $\rho$ and $\nu$ are the fluid density and kinematic viscosity, respectively.

Taking the curl of Equation (3), using the fact that the fluid is incompressible (then, from continuity equation, $\nabla \cdot \vec{u} = 0$ ) and remembering the definition of vorticity, ( $\omega = \nabla \times \vec{u})$ , Equation (3) can be re-written as:

(4) \begin{equation}\frac{D\vec{\omega} }{dt} = [\vec{\omega}\cdot \nabla] \vec{u} + \nu \nabla^2 \vec{\omega},\end{equation}

where $D()/dt = \partial()/\partial t + \vec{u} \cdot \nabla()$ is the material derivative.

Equation (4) is the vorticity-velocity in the Lagrangian form, and it is the fundamental equation for solving the transport of vorticity in the VVP.

2.4.3 Vorticity field discretization

The key idea behind the VVP is that the vorticity field can be discretised into vortex particles, corresponding to influencing elements characterised by a volume, a position, and strength. The global vorticity field is then approximated by the sum of each particle’s individual vorticity field.

In the classical proposition of the Vortex Particle method, the vorticity field is assumed to be concentrated at discrete points, called singular particle or point vortex^{(Reference Winckelmans and Leonard65,Reference Speck66)} and can be written as:

(5) \begin{equation}\vec{\omega}(\vec{x},t) = \sum_{i=1}^{N} \delta(\vec{x} - \vec{x}_i)\vec{\alpha}_i(t),\end{equation}

where $\delta(x)$ is the 3D $\delta$ -function, and $\vec{x}_i$ is the particle position and $\vec{\alpha}_i$ corresponds to its strength, which is given by the particle volume times the vorticity, $\vec{\alpha}_i = \vec{{\omega}}_i vol_p$ . Then, in this approach, the region outside of the vortex particles is irrotational.

The velocity field associated with such particles representation can be demonstrated^{(Reference Winckelmans and Leonard65)} to be:

(6) \begin{equation}\vec{u}_{\Psi}(\vec{x},t) = \nabla \times \vec{\Psi}(\vec{x},t) = \sum_{i=1}^{N} \vec{K}(\vec{x} - \vec{x}_i) \times \vec{\alpha}_i(t),\end{equation}

where $\vec{K}(\vec{x} - \vec{x}_i)$ is the Biot–Savart kernel, given by:

(7) \begin{equation}\vec{K}(\vec{x} - \vec{x}_i) = - \frac{1}{4\pi} \frac{\vec{x} - \vec{x}_i}{|\vec{x} - \vec{x}_i|^3}\end{equation}

From a quick examination of the Biot–Savart kernel, one can conclude that it is a singular function, leading to numerical instabilities when particles move toward each other. To overcome this problem, the concept of “regularised particles,” also known as “vortex blob” was introduced^{(Reference Beale and Majda67)} and is widely applied by studies using VVP.

The basic idea of regularised particles is to substitute the 3D $\delta$ -function, $\delta(\vec{x} - \vec{x}_i)$ , in the vorticity field representation (Equation (5)) by a distribution function (known as regularisation function, regularised smoothing kernel or also cutoff function), $\xi_\sigma(\vec{x} - \vec{x}_i)$ , such that the point vortex particles are transformed into vortex “blobs’’ with a finite core $\sigma$ .

Using this regularisation function, the vorticity field is now written as:

(8) \begin{equation}\vec{\omega}_\sigma(\vec{x},t) = \xi_\sigma(\vec{x} - \vec{x}_i) \circledast \vec{\omega}(\vec{x},t) = \sum_{i=1}^{N} \xi_\sigma(\vec{x} - \vec{x}_i) \vec{\alpha}_i (t),\end{equation}

where $\vec{\omega}_\sigma(\vec{x},t)$ is the regularised representation of the vorticity field and $\vec{\omega}(\vec{x},t)$ is the vorticity field represented by the singular particles (Equation (5)).

There are many possibilities for the choice of the regularised function (regularised smoothing kernel). Examples of two and three-dimensional kernels can be found in Refs (Reference Beale and Majda68,Reference Hald and Del Prete69,Reference Nordmark70,Reference Winckelmans and Leonard65). In this work, Gaussian smoothing kernel, with order $r = 2$ , was used. Gaussian smoothing regularisation function is a common choice for regularised vortex particles, as in Refs (Reference Singh and Friedmann58,Reference He and Zhao59,Reference Rajmohan, Zhao and He71,Reference Zhao and He72,Reference Zhao and He73).

Detailed mathematical demonstration for the convergence and stability of the vorticity field represented by regularised vortex particles has been developed in many studies^{(Reference Hald and Del Prete69,Reference Beale and Majda67,Reference Cottet74,Reference Anderson and Greengard75,Reference Beale76)} . It has been shown that, for stability and accuracy, the particles need to satisfy the overlap criteria, which couples inter-particle distances and core sizes:

(9) \begin{equation}\frac{\sigma}{h} > 1,\end{equation}

where h is the typical distance between neighbour particles.

2.4.4 Evolution equations

Substituting the vorticity representation, given by Equation (8), into the incompressible vorticity-velocity equation in Lagrangian form, Equation (4), results in the evolution equation for the strength of each particle i, as given by Equation (10). Also, vortex particles are considered fluid elements and, as so, travel with the local flow velocity (taking into account free-stream, velocity induced by particles, velocity induced by blades and any other influences in the velocity field), leading to the particle convection equation (Equation (11)).

(10) \begin{equation}\frac{D\vec{\alpha}_i}{dt} = [\vec{\alpha}_i\cdot \nabla] \vec{u} + \nu \nabla^2 \vec{\alpha}_i\end{equation}

(11) \begin{equation}\frac{d\vec{x}_i}{dt} = \vec{u}(\vec{x}_i,t).\end{equation}

Equations (10) and (11) are the viscous vortex particle evolution equations, governing the update of particles strength and positions at each time step. The first term in the right-hand side of Equation (10), $[\vec{\alpha}_i\cdot \nabla] \vec{u}$ , corresponds to the stretching effect, corresponding to vortex stretching and rotation due to the velocity field gradient. The second term in the right-hand side of that equation, $\nu \nabla^2 \vec{\alpha}_i$ , corresponds to the viscous diffusion, corresponding to the vorticity diffusion due to the viscous effects.

2.4.5 Generating new particles

After each time step, new particles are generated, satisfying the conservation of vorticity.

(12) \begin{equation}\vec{\Gamma}_{\rm wake} = - \frac{d\vec{\Gamma}}{dt} + \vec{u}_b(\nabla \cdot \vec{\Gamma})\end{equation}

where $\vec{\Gamma}_{\rm wake}$ is the circulation shed into the wake, $\vec{\Gamma}$ is the bound circulation and $\vec{u}_b$ is the resultant relative velocity between the airflow and the blade. The first term in the right-hand side of Equation (12) corresponds to the vorticity shed into the wake due to the time-varying circulation, called shed vorticity (associated with shedding particles). The second term corresponds to the vorticity released into the wake due to the gradient of circulation (in this case, the spanwise distribution of circulation) and is called trailing vorticity (associated with trailing particles).

2.5 Propeller inertial effects

The derivation of propeller inertia formulation is described in detail in previous work by the authors^{(Reference Teixeira and Cesnik33)}, and for completeness, it is summarised here.

The propeller is modeled as a rigid, massless rod with a distribution of concentrated masses. Also, it is assumed the propeller configuration is such that: (i) all blades have the same geometry and mass distribution, and (ii) the same angle between adjacent blades. Note that from (i) and (ii) one has:

(13) \begin{equation}\sum_{i=1}^{N_{\rm blades}} \vec{r}_{{\rm blade\,CG},{\it i}} = \vec{0},\end{equation}

where $N_{\rm blades}$ is the number of blades and $\vec{r}_{{\rm blade\,CG},{\it i}}$ is the CG position of the i-th blade relative to the propeller hub, written in any desired frame. Both assumptions above are reasonable for most practical propellers.

Additionally, the following set of coordinate systems is considered (as presented in Fig. 2): (i) Inertial frame (I): as the name implies, a fixed inertial global frame; (ii) Body frame (B): a frame attached to the body, that may translate and rotate with relation to the inertial frame; (iii) Local beam frame (W): a frame attached to a structural beam node and it may translate and rotate in relation to the body frame due to elastic deformations; (iv) Propeller frame (P): a frame whose origin is attached at the propeller hub position and has no rotation or translation with relation to local frame W; (v) Blade frame (b): a frame that is attached to one of the propeller blades, has the same origin as the propeller frame, and rotates with relation to it with the propeller angular velocity.

Figure 2. Representation of various reference systems used on modeling blade inertial effects. For simplicity, just one concentrated mass per blade is shown.

Based on the coordinate systems defined, the position of one of the concentrated masses i defined with relation to the inertial frame (I) can be written in terms of the body frame (B) as:

(14) \begin{equation}\vec{p}_{m_i} = \vec{p}_B + \vec{p}_W + \vec{p}_P + \vec{p}_{i},\end{equation}

where, as illustrated in Fig. 2, $\vec{p}_B$ is the position of the origin of the body frame with relation to the inertial frame, $\vec{p}_W$ is defined as the position of the origin of the local frame with relation to the body frame, $\vec{p}_P$ is the position of the origin of the propeller frame with relation to the local frame, and $\vec{p}_{i}$ is the position of the point mass i with relation to the propeller frame.

Recalling rigid-body dynamics theory, the time derivatives of a vector $\vec{r}$ expressed with respect to two different frames, (I) and (B), rotating with respect to each other is given by:

(15) \begin{equation}^I\dot{\vec{r}} = ^B\dot{\vec{r}} + ^I\vec{\omega}^B \times \vec{r},\end{equation}

where $^I\dot{\vec{r}}$ is the time derivative of $\vec{r}$ in frame I, $^B\dot{\vec{r}}$ is the time derivative of $\vec{r}$ with respect to B, and $\vec{^I\omega^B}$ is the angular velocity vector of B with relation to I.

The inertial force that each mass $m_i$ transfers to the structure is:

(16) \begin{equation}\vec{F}_{m_i} = - m_i \; ^I\ddot{\vec{p}}_{m_i}.\end{equation}

Applying Equation (15) into Equation (14) to calculate $^I\ddot{\vec{p}}_{m_i}$ , substituting in Equation (16) and adding the contributions of all concentrated masses, one can find the following expression for the total force acting on the propeller hub due to inertial effects ( $\vec{F}_{\rm hub}$ ):

(17) \begin{eqnarray}\vec{F}_{\rm hub} = & m_p\,(\ddot{\vec{p}}_B + \ddot{\vec{p}}_W + \ddot{\vec{p}}_P + 2\,^I\vec{\omega}^B\times\,(\dot{\vec{p}}_B + \dot{\vec{p}}_W + \dot{\vec{p}}_P) \nonumber\\&+ ^I\dot{\vec{\omega}}^B\times\,(\vec{p}_B + \vec{p}_W + \vec{p}_P) + ^I\vec{\omega}^B\times\,^I\vec{\omega}^B\times\,(\vec{p}_B + \vec{p}_W + \vec{p}_P)\,).\end{eqnarray}

Note that if one concentrates the mass $m_p$ at the propeller hub location and follows a similar derivation, the same Equation (17) would be obtained. Then, as expected, the inertial forces acting on the propeller hub are equivalent to the force acting in a punctual mass $m_p$ located at that location.

Now, consider the propeller inertial moments transferred to the body. The objective is to determine the moment due to the inertial forces acting on each rotating punctual mass transferred to the local node position, where the origin of the W frame is situated. The moment of those forces in relation to the origin of the W frame can be determined by:

(18) \begin{equation}\vec{M}_W = \vec{p}_P\times\vec{F}_{\rm hub} + \sum_{k=1}^{N_{\rm seg}} \sum_{i=1}^{N_{\rm blades}} [\,\vec{p}_i\times \vec{F}_{mi}\,],\end{equation}

where $N_{\rm blades}$ and $N_{\rm seg}$ are the number of blades and blades segments, respectively and the index k is related to the radial position along the blade i.

Then, the moment transferred to the origin of the W frame due to the inertial forces from the rotating masses is equivalent to the moment with relation to the origin of the W frame due to a concentrated mass $m_P$ at the propeller hub ( $\vec{p}_P\times\vec{F}_{\rm hub}$ ) plus a contribution due to the fact that masses are rotating around the hub ( $\vec{M}_{\rm rot} = \sum_{k=1}^{N_{\rm seg}} \sum_{i=1}^{N_{\rm blades}} [\,\vec{p}_i\times \vec{F}_{mi}\,]$ ).

Substituting the individual contributions $\vec{F}_{mi}$ of each concentrated mass (calculated as described previously), one can find:

(19) \begin{eqnarray} \vec{M}_{\rm rot} = & \sum_{k=1}^{N_{\rm seg}} \sum_{i=1}^{N_{\rm blades}} - m_k [\,\vec{p}_{i,k} \times (\ddot{\vec{p}}_{i,k} + 2 ^I\vec{\omega}^B\times\,\dot{\vec{p}}_{i,k} \nonumber\\& +\, ^I\dot{\vec{\omega}}^B\times\,\vec{p}_{i,k} + ^I\vec{\omega}^B\times\,^I\vec{\omega}^B\times\,\vec{p}_{i,k})],\end{eqnarray}

with:

(20) \begin{equation}\dot{\vec{p}}_{i,k} = (^B\vec{\omega}^P + ^P\vec{\omega}^b) \times \vec{p}_{i,k}\end{equation}

(21) \begin{eqnarray}\ddot{\vec{p}}_{i,k} = & ^B\dot{\vec{\omega}}^P\,^b\vec{p}_{i,k} + (^B\vec{\omega}^P + ^P\vec{\omega}^b) \nonumber\\& \times(^B\vec{\omega}^P +\, ^P\vec{\omega}^b)\times \vec{p}_{i,k},\end{eqnarray}

where $m_k$ is the mass of a concentrated mass located at blade segment k, $\vec{p}_{i,k}$ is the position of the concentrated mass with relation to a frame attached to the propeller hub, $^I\vec{\omega}^B$ is the angular velocity of the body with relation to an inertial frame, $^B\vec{\omega}^P$ is the angular velocity of propeller frame with relation to body frame, and $^P\vec{\omega}^b$ is the angular velocity of blade frame with relation to propeller frame (given by propeller rotation per unit time). Also, all derivatives are defined in the body frame.

Then, the propeller inertial loads are incorporated by defining a concentrated mass $m_p$ at each propeller hub and calculating an additional moment term given by Equations (19), (20) and (21).

2.6 Interface structure and aerodynamics

Since the structural model employs a 1D finite element representation, and the uVLM uses 2D non-planar grid, an interface between the two solvers is required to transfer loads and displacements.

Assuming rigid behavior in the chordwise direction (no change in local camber), the local beam coordinates and orientation given by the UM/NAST structural solver are converted to a non-planar panel grid by re-writing the camber line coordinates of each section along the span from aerofoil coordinates to body coordinates.

The transfer of loads from panels to beam structural nodes is currently accomplished in two steps: (i) energy-consistent transfer of loads from the panels to the neighbouring points along the beam given by panel discretization (based on equivalence of virtual work between the forces at panels and forces transferred to beams) and (ii) linear transfer of loads from those points to the structural nodes, in case of different structural and aerodynamic spanwise discretisations. This approach is in accordance with the assumption of constant strains along each beam element.

2.7 Interface propeller and lifting surfaces aerodynamics

To integrate the uVLM for lifting surfaces with propeller modeling, the mutual influence between them needs to be taken into account. Figure 3 illustrates how the process of integration occurs.

Figure 3. uVLM and propeller model integration approach.

At a given time step, the wing panel and wake panel circulations are updated by uVLM considering the induced velocity due to current particle distribution (calculated in the previous time step) and blade bound circulation. Then, based on the new wing and wake panel circulations, the new circulation of the blades are calculated, the existing particles have positions and strength updated, and new particles are generated, defining new particle distribution for the next time step. It is important to notice that the time step used in the propeller solver should be smaller than the propeller period divided by the number of blades. To not restrict the time step of the dynamic solver by the propeller solver requirements, one should consider the propeller time step smaller than the global dynamic solver time step. In this work, and unless otherwise stated, the propeller time step was half of the time step on the dynamic solver.

3.0 ALTERNATIVE APPROACH TO EXTRACT DYNAMIC INFORMATION

Traditionally, two main approaches are used to extract dynamic information (frequencies, damping, and modes) and perform a stability analysis of an aeroelastic system: frequency and time-domain methods. The frequency-domain methods are primarily developed for linear aeroelastic systems, as they assume the superposition of signals with different constant frequencies. Therefore, they are more suitable for problems with small deformations. Typical methods for stability analysis based on frequency domain are the k method^{(Reference Bisplinghoff and Ashley77)}, p-k method^{(Reference Hassig78)} and p method^{(Reference Abel79)}.

For very flexible aircraft undergoing large deformations, however, strong nonlinearities can arise, and time-domain methods are usually preferred^{(Reference Hallissy and Cesnik80)}. Examples of time-domain methods for stability analysis are reviewed in McNamara and Friedmann^{(Reference McNamara and Friedmann81)}: moving-block approach (MBA)^{(Reference Hammond and Dogget82)}, least squares curve-fitting method (LSCFM)^{(Reference Bennett and Desmarais83)} and system identification using the autoregressive moving-average (ARMA) model^{(Reference Onoda84)}.

Even for linear systems, the complexity of the problem can make the use of frequency-domain techniques a cumbersome task, as the equations are not easily manipulated. In this case, time-domain simulations may also be preferred. In Silva^{(Reference Silva85)}, for example, the dynamic information of a linear transonic aeroelastic problem based on modal solutions is obtained by the system identification of CFD simulations.

3.1 Why an alternative approach?

In UM/NAST, the flutter solution is obtained based on the linearisation of the governing equations of the system. The model is first linearised about its trimmed condition, and then the dynamic information is extracted based on the eigendecomposition of the so-called state matrix A. This was developed by the original UM/NAST aerodynamic formulation, which used Peters’ inflow theory.

Incorporating the UVLM aerodynamics for the lifting surfaces and LL/VVP formulation for the propellers in the UM/NAST linearised formulation would require an in-depth review of the linearisation formulation. Also, while the analytical expressions of Peters’ inflow theory allow a direct coupling with the structural equations, the same is not the case for the UVLM and LL/VVP approaches, which are solved separately from the structural solver at each time step. This option was, then, out of the scope for the present work. Instead, an alternative approach for the stability analysis of VFA was proposed and applied. The procedure makes use of system identification combined with Proper Orthogonal Decomposition (POD) modes (Sys ID + POD) in order to reduce the number of degrees of freedom (outputs) of the problem. This is important, as the system identification is not well suited for a high number of outputs. Another advantage is that it can be applied to any time signal (e.g. experimental data), and a modal structural approach is not required (the method can be used for deformed configurations far from the undeformed condition).

3.2 Combining POD and Sys ID for stability analysis

In order to extract dynamic information (frequencies, damping, and modes) from the coupled aeroelastic-flight dynamics framework with propellers, it is first necessary to find a linearised mathematical representation of the system containing all the coupled structural, flight dynamics, lifting surfaces aerodynamics and propeller effects (inertial and aerodynamics).

One natural choice for this is to apply system identification and find a state-space representation of the global system. However, one difficulty in directly applying system identification in this context is that the VFA model often has a high number of degrees of freedom (the model considered later, for example, has more than 500 degrees of freedom), resulting in a high number of output states. However, system identification is not well suited for a high number of outputs, and the prediction performance is deteriorated.

One option to tackle this problem would be to write the output data as a function of a linear combination of the aircraft linear modes. This was done, for example, by Silva^{(Reference Silva85)} and studies for aeroelastic problems involving small deformations. However, for a VFA, a representation of the deformed equilibrium state with linear modes is no longer adequate.

Another alternative, used in this work, is to reduce the output matrix dimension by using POD modes. Then, instead of directly providing the output matrix containing the snapshots of displacements, the snapshots containing the coefficients of a much smaller number of POD modes are provided, and the original degrees of freedom are later recovered. This has the additional advantage of providing a better basis representation for an aircraft whose equilibrium state exhibits large deformation compared to the undeformed shape. Figure 4 illustrate how this combination of POD + Sys ID works.

Figure 4. POD + Sys ID workflow.

First, the original data snapshots are collected by performing a time simulation disturbing the aircraft from its equilibrium condition. It is worthwhile to note that, for the extraction of dynamic information, it suffices to apply small disturbance around the deformed equilibrium condition, making possible the use of linear system identification approaches. Then, the input matrix consists of the snapshots of the values of loads (disturbances) applied. The original output matrix consists of the snapshots of the displacements (in x, y and z directions) of each structural nodes with relation to an initial equilibrium condition.

Next, the POD modes of the original output need to be determined and selected. For this, SVD is first applied to the original output matrix M:

(22) \begin{equation}M_{m\times n} = U_{m\times m}\,\Sigma_{m\times n}\,V^T_{n\times n}\end{equation}

where m is the number of d.o.f., n is the number of snapshots, U and V are orthogonal matrices containing left-singular vectors and right-singular vectors of M, with the columns of U corresponding to the POD modes, and $\Sigma$ is a diagonal matrix of non-negative real number ordered in descending order. Also, the product $\Sigma_{m\times n}\,V^T_{n\times n}$ is associated with the coefficients of the POD modes which reconstitute the matrix M.

Usually, the first few POD modes are responsible for more than 99% of the trace of Matrix $\Sigma$ . Reducing the number of POD modes to p, one obtains the new output matrix containing snapshots of the coefficients of the p POD modes which approximate the matrix M where $p << m$ . For this, a reduced singular value decomposition representation of M can be found by using just the p POD modes selected as:

(23) \begin{equation}M_{m\times n} = U_{r,m\times p}\,\Sigma_{r,p\times p}\,V^T_{r,p\times n}\end{equation}

The new reduced output matrix N to be provided to the system identification is then given by:

(24) \begin{equation}N_{p\times n} = \Sigma_{p\times p}\,V^{T}_{p\times n}\end{equation}

Once the reduced output matrix is determined, an input-output based Sys ID is applied and an equivalent mathematical model given by the system of Equations (25) and (26) is determined:

(25) \begin{equation}\dot{x}(t) = A\,x(t) + B\,u(t) + K\,e(t)\end{equation}

(26) \begin{equation}y(t) = C\,x(t) + K\,e(t),\end{equation}

where A, B, C, and D are the estimated matrices of state space, K a matrix which accounts by disturbances, u(t) is the input, y(t) is the output and x(t) is a vector with $n_x$ states associated with this mathematical representation. The identification was performed using the tool n4sid available on Matlab and more details about the method can be found in Van Overschee and De Moor^{(Reference Van Overschee and De Moor86)}.

Finally, the dynamic information of the system can be found by the eigenvalue decomposition of matrix A, in Equation (25). As the state-space representation is in a continuous time domain, the frequencies and damping can be directly obtained from the eigenvalues $\lambda_i$ of A, as:

(27) \begin{equation}freq = \frac{Re(\lambda_i)}{2\pi}\end{equation}

(28) \begin{equation}damp = \frac{Imag(\lambda_i)}{2\pi},\end{equation}

where $Re(\lambda_i)$ and $Imag(\lambda_i)$ are the real and imaginary parts of $\lambda_i$ , respectively.

From the eigenvectors of A, and remembering that they are related to coefficients of the chosen POD modes, the mode shapes can be obtained. Considering E as the matrix containing the eigenvectors of A and by Equation (26), the modes in terms of displacements at structural nodes can be recovered by:

(29) \begin{equation}N_{p\times n} = U_{r,m\times p} C_{p\times k} E_{k\times k}\end{equation}

where k is the order chosen for the system identification method.

Due to nonlinearities, non-proportional damping, and possible noise, the modes obtained by the system identification can be complex. In Rainieri and Fabbrocino^{(Reference Rainiere and Fabbrocino87)}, a discussion of those complex modes from system identification is presented, and an approach, used in the present work, to convert them to real mode shapes was proposed. The mode shapes can be obtained by adding the real modes of displacements with the initial equilibrium condition.

Finally, it was observed that, for better accuracy, the number of frequencies determined should be no more than the number of POD modes selected.

4.0 NUMERICAL STUDIES

4.1 Simulation details

A verification of the POD combined with system identification (sys ID) is first performed by a purely structural case, and then the method is applied for the stability analysis of the aeroelastic cases. For these studies, two models are used: (i) the University of Michigan’s X-HALE UAS^{(Reference Cesnik, Senatore, Su, Atkins and Shearer88)} for the purely structural case and (ii) the University of Michigan’s X-HALE UAS with an added tip mass of 0.5kg located 0.36m behind the wing trailing edge for the aeroelastic cases. This tip mass was added such that the model could present unstable aeroelastic behaviour in a range of velocities within the aircraft flight envelope. As shown in Fig. 5, the X-HALE has a wingspan of 6m with a 0.2-m chord, five pods along the wing, five tails, three fins and five electric motor-propeller combinations located in front of each pod at spanwise locations $y = -2, -1, 0, 1$ , and 2m. The wing is mounted with an incidence angle of 5 degrees. For the simulations, the propellers are located 20cm ahead and 2.8cm below of the wing leading edge, with a pitch angle of 5 degrees between the propeller axis and the X-HALE wing, as the wing has an incidence angle with respect to the pods. In order to reduce asymmetric loads, the propellers on the right wing have a different direction of rotation with relation to propellers at the left wing. A summary of the propeller parameters used here is presented in Table 1 and more details about the choice of blade discretisation and time step (based on a parametric study) and particle core radius ( $\sigma$ ) can be found in Teixeira and Cesnik^{(Reference Teixeira and Cesnik33)}. A cut-off distance of two radii after the end of the X-HALE vertical tails was applied to save computational cost. In the aeroelastic solver, a stiffness-proportional damping coefficient of 0.005s and a time step of 0.002s (for 6000RPM) and 0.0017 (for 7000RPM) were considered (with a sub-time step of half structural step). The dynamic viscosity is $\mu = 1.7855 \; 10^5$ N.s/m $^2$ and the air density is $\rho _\infty = 1.225$ kg/m $^3$ . For the aeroelastic cases, gravity effects are also considered.

Figure 5. Undeformed panel model for the X-HALE UAS vehicle (units: meters).

Table 1 Two-bladed APC 11X5.5E propeller parameters

4.2 Verifying POD + Sys ID for purely structural case

To verify the capability of the POD + Sys ID method to extract frequencies, damping, and modes based on the snapshots generated by the simulation of a VFA, a purely structural case was used, for which a reference from the UM/NAST modal solver was available about its undeformed shape. For the POD + Sys ID, a time simulation was performed for the clamped model from the undeformed configuration disturbed after 1.5s by a 5Nm torsion moment and a 1N step force in the vertical and chordwise direction, all applied at the right wing tip of the clamped model and chosen to excite different modes. The snapshots were then provided to the POD + Sys ID method, 5 POD modes were used (contributing for more than 99% of the snapshot energy), and an order of 18 was chosen for the system identification method based on the fitting quality of the POD coefficients. As the reference was for the undeformed, undamped case, no gravity or damping effects were included. Also, as the identification order of 18 provides nine frequencies, but just five POD modes are used (and up to five frequencies can be more accurately identified), the main frequencies were determined by choosing the five higher norms of the corresponding coefficients identified for the POD modes. Therefore the frequencies associated with the five columns of $C_{p\times k} E_{k\times k}$ with higher Euclidean norms were kept and shown in Fig. 6.

Figure 6. Comparison of mode shapes and natural frequencies for the X-HALE UAS vehicle about its undeformed configuration. (a) First bending (UM/NAST): 0.5943Hz. (b) First bending (POD + SysID): 0.5923Hz. (c) First torsion (UM/NAST): 2.5747Hz. (d) First torsion (POD + SysID): 2.5299Hz. (e) Second bending (UM/NAST): 3.6986Hz. (f) Second bending (POD + SysID): 3.6631Hz. (g) First in-plane bending (UM/NAST): 4.4491Hz. (h) First in-plane bending (POD + SysID): 4.4900Hz. (i) Second torsion (UM/NAST): 6.5696Hz. (j) Second torsion (POD + SysID): 6.5286Hz.

Table 2 presents the comparison between frequencies as well as a comparison of the corresponding modes by applying Model Assurance Criteria (MAC). For this purely structural case, one can see the method had an excellent agreement with the data calculated by the modal solver in UM/NAST with a maximum error in frequency prediction smaller than 2% and a MAC value of about 0.99 for all five modes.

Table 2 Comparison of natural frequencies for purely structural case

4.3 Influence of the inclusion of aerodynamic and gyroscopic propeller effects on the aeroelastic stability

In order to check the additional influence of the propeller slipstream and gyroscopic effects, besides the influence of thrust, transient solutions for a clamped aeroelastic case were conducted. Beginning from an equilibrium state and after 1.5s of simulation, a step perturbation of 5Nm in the torsional moment and a 1N step force in vertical and chordwise directions were applied. In this case, gravity and a stiffness-proportional damping coefficient of 0.005s were considered, to have a more realistic response. Three cases with different speeds: $12.5$ , 13 and $13.5$ m/s were simulated for a model with just thrust at the propeller location (equivalent to the thrust produced by the isolated propeller at each speed for the considered RPM) and a complete propeller model, including thrust (and other loads at hub), slipstream, and gyroscopic effects, with a rotation of 6000RPM.

First, the same set of parameters determined in the verification case was tried (five POD modes and Sys ID order of 18), but for the cases with aerodynamic and, for some of them, gyroscopic moment, the adjusting of the fifth POD coefficient using the same set of parameters for all cases was hard to get (especially for cases after the flutter boundary). Then, just four POD modes were considered (representing more than 99% of the snapshot energy), and an order of 18 was again included, improving the fitting of the first four POD coefficients. For the purely structural case, the effect of reducing the POD modes to four, keeping the same order for the Sys ID, was a reduction in the accuracy of higher frequencies. But the first four frequencies were kept with good accuracy (less than 2% error). Then, for the cases in this section, just four frequencies and modes are compared, corresponding to the most important modes for the system response. Also, the sensitivity of the results with the number of snapshots provided before the perturbation starts was found to be higher for cases with higher numerical perturbations, as in the cases after flutter. All cases use the same set of parameters, and it is expected that the results can capture frequency variations due to the different effects being modeled.

Figure 7 presents a comparison of the frequencies, dampings, and dynamic responses in terms of wing tip vertical and angular displacements for each simulated speed. From the dynamic responses, it is possible to see that for this model and propeller RPM, the flutter boundary is about 13m/s for both cases, i.e. thrust only, and including complete propeller modeling. However, as one can also observe from the dynamic responses, for all speeds increasing differences in phase arise between cases with thrust only and with all propeller effects for all speeds. For the same set of POD + Sys ID parameter choice, one can also note differences in frequencies and especially damping for modes farther away from the stability boundary. Considering that the variations of frequency and damping are captured, this indicates propellers can affect the dynamic response.

Figure 7. Frequencies, dampings, and response after perturbation with step loads of 5N in vertical and chordwise direction and 1Nm in torsion (6000RPM). (a) $v = \text{12.5}$ m/s. (b) $v = \text{12.5}$ m/s. (c) $v = \text{12.5}$ m/s. (d) $v = \text{13}$ m/s (e) $v = \text{13}$ m/s. (f) $v = \text{13}$ m/s. (g) $v = \text{13.5}$ m/s. (h) $v = \text{13.5}$ m/s. (i) $v = \text{13.5}$ m/s.

4.4 Influence of increasing propeller RPM

In order to investigate the effect of propeller RPM in the aeroelastic stability, a transient solution similar to the one presented in Fig. 7 was simulated for the clamped X-HALE model with tip masses at $v = 13$ m/s and with a higher propeller rotation, 7000RPM. Increasing the propeller RPM with all other parameters constant means a stronger influence on the velocity flowfield behind the propeller plane of rotation, as well as a higher thrust and gyroscopic loads. It can be interpreted as the degree of propeller effects intensity, which could also be affected by other parameters like propeller mass, geometry, velocity flowfield, etc.

Figure 8 presents a comparison of the wing tip dynamic responses. Now, more noticeable differences can be observed in the stability behaviour; while the case with 6000RPM is yet stable (although close to flutter), the case with 7000RPM shows an unstable behaviour with increasing amplitudes, then a lower flutter speed. This destabilising effect with the increase in RPM is possibly related to the increase of dynamic pressure proportioned by the propeller slipstream, causing higher lift forces for the same local angle-of-attack. This example indicates that depending on the intensity of propeller effects the influence on stability boundary can be more significant. Damping and frequency were not compared, as the same set of parameters used in $v = 13$ m/s for 6000RPM do not provide a good fitting for the case of 7000RPM.

Figure 8. Wing tip response after perturbation with step loads of 5N in vertical and chordwise direction and 1Nm in torsion, including propeller aerodynamics and gyroscopic effects. (a) Tip vertical displacement. (b) Tip twist.

4.5 Contribution of different propeller effect components

Figure 9 compares the dynamic response as well as frequency and damping for the same disturbance, as in Fig. 8 for the clamped X-HALE model with tip masses at $v = 13$ m/s. This case is already in the unstable regime. Due to more difficulties in find a common set of parameters with good fitting and in order to reduce the perturbation due to numerical noise, just three POD modes were incorporated (representing more than 97% of the snapshot energy), keeping the identification order as 18, and the snapshots were provided with 0.5s after beginning of perturbations.

Figure 9. Frequencies, dampings, and response after perturbation with step loads of 5N in vertical and chordwise direction and 1Nm in torsion (7000RPM). (a) Root locus. (b) Wing tip vertical displacement. (c) Wing tip twist.

In order to observe the isolated propeller effect components (mainly thrust, slipstream, and gyroscopic moment), three different modelings of propeller effects are considered:

1. (i) with thrust only (equivalent to thrust produced by isolated propeller at $v = 13$ m/s and 7000RPM);

2. (ii) with thrust (and other loads at the hub) and slipstream (propeller aero), and;

3. (iii) with thrust (and other loads at the hub), slipstream, and gyroscopic effects (propeller aero + gyro).

From the dynamic response, it is clear that the inclusion of propeller aerodynamic effects can cause a non-negligible difference in the aeroelastic stability, with the model including just thrust having a smoother amplitude increase, suggesting a higher flutter boundary. As previously commented, this destabilising effect may be related to the increased dynamic pressure due to the propeller slipstream, causing higher lift (and consequently higher structural deformations) for same freestream conditions. An increasing difference of phase with time is also noticeable between the case with just thrust and other two cases. The inclusion of gyroscopic effects causes negligible differences for vertical displacements but more visible differences in the amplitude of angular displacement, which is smaller for the case including gyroscopic effects. This is in accordance with the root locus, where it is clear that the positive damping for the unstable modes in the case with slipstream and no gyroscopic effects is higher, suggesting the gyroscopic effect has a stabilising effect.

5.0 FINAL REMARKS

This work investigated the influence of propeller aerodynamics and gyroscopic effects on the aeroelastic stability of very flexible aircraft. For that, an enhanced aeroelastic framework with propellers developed in previous efforts by the authors^{(Reference Teixeira and Cesnik31,Reference Teixeira and Cesnik33)} was applied to extract time-data snapshots of the clamped modified X-HALE UAS representing a HALE aircraft. A method based on POD plus system identification to extract frequencies, dampings and modes from the time-series data was successfully verified for a purely structural case for which a reference solution was available. The method was then used to analyse the effect of the propeller on the aeroelastic stability of the very flexible aircraft based on a set of generated snapshots. From the analysis of aeroelastic cases and the clamped model considered, it was observed that the presence of propeller aerodynamics and gyroscopic effects influences the values of damping and frequencies of some modes and can influence stability boundary. Moreover, an increase in phase delay and differences in amplitude in response to a perturbation close to flutter were shown as compared to the case with just thrust. A reduction of flutter boundary was found by an increase in propeller RPM. This destabilising effect may be associated with the impact of slipstream, increasing the dynamic pressure and local lift. Also, a slight stabilising effect due to the gyroscopic moment was noticed, suggesting this effect may be more important for cases with higher gyroscopic moment loads, i.e. higher RPM, and propeller mass. Finally, the method of POD plus system identification showed to be effective to extract frequencies, dampings, and modes for a purely structural case. Further investigations are necessary to make the method more robust for aeroelastic cases.