Nonlinear Analysis and Bifurcation Characteristics of Whirl Flutter in Unmanned Aerial Systems

Abstract

Aerial drones have improved significantly over the recent decades with stronger and smaller motors, more powerful propellers, and overall optimization of systems. These improvements have consequently increased top speeds and improved a variety of performance aspects, along with introducing new structural challenges, such as whirl flutter. Whirl flutter is an aeroelastic instability that can be affected by structural or aerodynamic nonlinearities. This instability may affect the prediction of potentially dangerous behaviors. In this work, a nonlinear reduced-order model for a nacelle-rotor system, considering quasi-steady aerodynamics, is implemented. First, a parametric study for the linear system is performed to determine the main aerodynamic and structural characteristics that affect the onset of instability. Multiple polynomial nonlinearities in the two degrees of freedom nacelle-rotor model are tested to simulate possible structural nonlinear effects including symmetric cubic hardening nonlinearities for the pitch and yaw degrees of freedom; purely yaw nonlinearity; purely pitch nonlinearity; and a combination of quadratic, cubic, and fifth-order nonlinearities for both degrees of freedom. Results show that the presence of hardening structural nonlinearities introduces limit cycle oscillations to the system in the post-flutter regime. Moreover, it is demonstrated that the inclusion of quadratic nonlinearity introduces asymmetric oscillations and subcritical behavior, where large and potentially dangerous deformations can be reached before the predicted linear flutter speed.

1. Introduction

Drones, or unmanned vehicles, are becoming more common for a variety of applications due to their increased reliability, improved aerodynamics, pilot safety concerns with manned aircraft with regards to force limits caused by velocity/maneuverability, and most importantly only being limited by structural and aerodynamic parameters rather than the pilot’s ability to handle these flight conditions before losing consciousness [1]. This leads to the importance of investigating the fluid-structure interaction limitations of these unmanned aerial systems, and addressing all physical phenomenon that can arise. For instance, the stall effect must still be considered as the physical limits of an airfoil can cause flight disturbances or failure. Furthermore, rotor-powered systems are becoming increasingly utilized due to their unique ability to perform vertical take-off and landing (VTOL) while still having the capability of operating as a fixed-wing aircraft, as well as being cheaper to operate than jet-powered systems [2]. The most popular and well-developed of these vehicles is the tilt-rotor, which alters the rotor configuration to perform vertical take-off and landing or forward flight. This additional degree of freedom, combined with large diameter rotors, causes additional complications and complexities to an already present phenomenon known as whirl flutter.

Whirl flutter is a potentially disastrous aeroelastic instability faced by propeller-driven aircraft that was first witnessed in 1938 but came to prominence again in 1960 [3,4]. Whirl flutter can cause propellers and the surrounding structures to fail and therefore must be accurately predicted. Whirl flutter is the result of a coupling of the structure of the wing, the aerodynamic forces and moments on the wing and the rotor, and the gyroscopic effect of the rotor [5,6]. The rigid nature of the rotors in these VTOL systems introduces additional couplings not seen in traditional helicopters [7]. The onset of whirl flutter is difficult to predict, and the inclusion of the additional degree of freedom from the tilting rotors makes this an even more challenging task. Often, for simplicity, nonlinearities are neglected in the modeling of these systems, but these systems commonly have nonlinear characteristics and behaviors due to the aerodynamics and the inherent structural nonlinearities in the system, such as freeplay caused by construction characteristics, aging, and wear [8]. These nonlinearities make the flutter prediction a more challenging task, and if they are not accurately accounted for, it could cause a catastrophic system failure.

There is a multitude of modern propulsion systems currently in use by aerial vehicles. A broad list of unmanned aerial vehicles using a propeller-driven system from the United States is shown in

Table 1. Broad list of UAVs with either fixed-wing or rotary-wing configuration. #: This column is a list of reference citation numbers used for each aircraft.

Reference Citation #	Application	Model	Manufacturer	Type	Maximum Speed (Km/h)	Wingspan (m)	Length (m)
[9]	Military	Black Hornet	FLIR	Rw	18	0.12	0.168
[10]	Military	Desert Hawk III	Lockheed Martin	Fw	92	1.5
[11]	Military	InstantEye (Gen3)	Physical Sciences	Rw
[12]	Military	Karma	GoPro	Rw
[13]	Military	K-MAX	Kaman	Rw	185	14.7	15.8
[14]	Military	MQ-1 Predator	GA-ASI	Fw	215	16.8	8.22
[15]	Military	MQ-19 Aerosonde	AAI	Fw		11
[16]	Military	MQ-1C Gray Eagle	GA-ASI	Fw	305	17	9
[17]	Military	MQ-25 Stingray	Boeing	Fw
[18]	Military	MQ-4C Triton	Northrop Grumman	Fw	613	40	14.5
[19]	Military	MQ-8 Fire Scout	Northrop Grumman	Rw	250	10.7	12.6
[20]	Military	MQ-9 Reaper	GA-ASI	Fw	445	20	11
[21]	Military	Phoenix 30	UAV Solutions	Rw		0.5	0.5
[22]	Military	RQ-11	AeroVironment	Fw	81	1.4	1
[23]	Military	RQ-12 Wasp	AeroVironment	Fw	83	1	0.8
[24]	Military	RQ-170 Sentinel	Lockheed Martin	Fw
[25]	Military	RQ-180	Northrop Grumman	Fw
[26]	Military	RQ-20 Puma	AeroVironment	Fw	83	2.8	1.4
[27]	Military	RQ-21 Blackjack	Insitu	Fw	167	4.9	2.5
[28]	Military	RQ-23 Tiger Shark	Navmar Applied Sciences Corporation	Fw	148	6.7	4.5
[29]	Military	RQ-4 Global Hawk	Northrop Grumman	Fw		39.9	14.5
[30]	Military	RQ-7 Shadow	AAI	Fw	200	6	3.6
[31]	Military	ScanEagle	Insitu	Fw	150	3.11	1.71
[32]	Military	Silver Fox	BAE, Raytheon	Fw	93	2.4
[33]	Military	Stalker	Lockheed Martin	Fw	72	3.7
[34]	Military	Switchblade	AeroVironment	Fw	157
[35]	Military	T-20 JUMP	Arcturus	Fw	139	5.3	2.7
[36]	Military	T-Hawk	Honeywell	Rw
[37]	Hobby	Mini 2	DJI	Rw	57	0.29	0.25
[38]	Hobby	Matrice 300	DJI	Rw	82.8	0.43	0.42
[39]	Hobby	Mavic 2	DJI	Rw	68	0.18	0.25
[40]	Hobby	Inspire 2	DJI	Rw	94	0.42	0.31
[41]	Hobby	Phantom 4 Pro	DJI	Rw	72	0.35	0.35
[42]	Hobby	Autel EVO II Pro	Autel	Rw	72	0.40	0.40
[43]	Hobby	Brinc Lemur Tactical Drone	Brinc	Rw	80
[44]	Hobby	TRINITY F90+	QUANTUM	Tw	61	2.39	1.20
[45]	Hobby	Sentaero	Censys technologies	Tw	72		1.20
[46]	Hobby	XENO FX	HiTech	Fw	61	1.25	1.00
[47]	Hobby	P-47 Thunderbolt	Durafly	Fw		1.1	0.99
[48]	Hobby	KingTwin 1700	Avios (PNF)	Fw		1.7	1.2
[49]	Hobby	P-51 Warbird	RC Moment	Fw		1.2	1.1
[50]	Hobby	BOREY 20	Uavos	Fw	108	4.4	0.9
[51]	Hobby	ALBATROSS 2.2	Uavos	Fw	275	15.0	6.5
[51]	Hobby	Apus Duo 15	Uavos	Fw	97	15
[51]	Hobby	SAT-i	Uavos	Fw	50	7.30	3.20
[50]	Hobby	SITARIA E	Uavos	Fw	140	5.2	2.8
[52]	Hobby	Albatross UAV	Applied Aeronautics	Fw	129	3.0
[53]	Hobby	PD-2	Ukrspec Systems	Fw/Rw	140	5
[54]	Hobby	EOS C VTOL	Threod	Fw/Rw	122	5.0	1.8
[55]	Hobby	PENGUIN B	UAV Factory	Fw	130	3.3	2.3

, along with the model’s name, manufacturer, propulsion type, maximum speed, wingspan, and overall length. The two types of propeller propulsion systems outlined are rotary-wing (Rw) and fixed-wing (Fw), illustrated in column four of Table 1. Due to the fundamental nature of the mathematical model, the axis of propeller thrust or the size of the aircraft will not affect the results showing when the onset of the whirl flutter phenomenon will occur.

The given data in Table 1 can be categorized by propulsion type, wingspan, length, and maximum flight speed in order to select a planform to generate a mathematical model. These systems rely on motor and propeller propulsion systems, each with unique geometric and dynamic properties. The improvement to motor function and propeller design has made rotor-powered aerial vehicles faster and more maneuverable. These systems, with increased cruise/maximum speeds, have brought new challenges, including a structural phenomenon known as whirl flutter. Therefore, a full mathematical model of structure and aerodynamics is needed to proficiently define vibrational effects on the drone.

Other researchers have begun investigating this topic, but a fully-encompassing study has not yet been carried out. Previous researchers have taken various approaches to this problem. Many have opted to use reduced-order or lumped-parameter models solved analytically, while others have chosen to use computational fluid dynamics (CFD) and multiphysics packages or experiments to investigate this complex topic. Mair et al. [8] carried out an analytical nonlinear stability analysis for a rotor system. This investigation was limited to nonlinearity in the yaw degree of freedom, and nonlinearities are varied for different odd orders of nonlinearity (third-order, fifth-order, etc.). Bifurcation diagrams were presented, and corrected stability boundaries were given, accounting for the additional nonlinearities. In a similar study, Mair et al. [56] presented an improved tilt-rotor model and investigated the stability boundaries and bifurcations when including odd orders of nonlinearity. Both of these analytical studies showed that a correct stability boundary is not found through linear studies, and nonlinear effects must be included.

Rather than using analytical models, other researchers have chosen to use CFD simulations to investigate whirl flutter onset and corresponding stability. Yeo et al. [57] employed two different rotorcraft analysis codes to investigate tilt-rotor whirl flutter. Comprehensive Analytical Model of Rotorcraft Aerodynamics and Dynamics (CAMRAD) II and Rotorcraft Comprehensive Analysis System (RCAS) are both used to investigate whirl flutter and to compare results to one another. The codes compared well to each other, and it was found that blade pitch-flap coupling, rotor lag frequency, rotor rotational speed, and density all have a major impact on whirl flutter onset speed and therefore need to be correctly modeled and selected. Shen et al. [58] also compared two analyses methods for tilt-rotor. In their work, RCAS was compared to Dymore for tilt-rotor whirl flutter prediction. Higgins et al. [59] also employed CFD simulations to investigate whirl flutter onset and stability on the XV-15 and V-22 aircraft. Their solution technique combined the University of Glasgow’s CFD solver (Helicopter Multi-Block) and a NASTRAN structural model. Expanding on the research.

Researchers have not so frequently, especially in recent years, used experimental investigations to study whirl flutter onset and characteristics. Piatak et al. [60] used an experimental study to investigate the effects of pitch-flap coupling and control stiffness on the stability boundaries. In a similar investigation, Wilson et al. [61] proposed a new microwave-based measurement system for vibrations from flutter to aid other investigators in their experimentation.

From the discoveries of tilt-rotor whirl flutter onset, some researchers have turned towards finding a solution in order to control this potentially disastrous phenomenon. Acree et al. [62] investigated rotor designs that could improve the stability margins of these systems using a CAMRAD II model that allows for easy rotor design modifications. Heeg et al. [63] carried out whirl flutter analyses for the proposed X-57 to predict any indication that whirl flutter should be present within its operating condition for different configurations using CAMRAD II and Dymore. Based on their analyses, they found a proposed configuration that should not exhibit any whirl flutter in operation.

This study investigates several parameters governing the onset of whirl flutter to identify critical configurations that can lead to flutter onset. Further research in this work investigates the stability of whirl flutter in the presence of quadratic and cubic structural nonlinearities. In order to better understand the overall system’s response, several nonlinear stiffness models are selected to evaluate the hardening or softening effects of structural stiffnesses in the two degrees of freedom of the system. A bifurcation analysis is then performed to determine the impacts of the structural nonlinearities on the type of instability (supercritical or subcritical) and the amplitude of oscillations. The rest of this paper is organized as follows: The nonlinear aeroelastic modeling is presented in Section 2 and is followed by an investigation of the linear characteristics of the system in Section 3. The nonlinear characteristics of the system are presented in Section 4 and the nonlinear stiffness response with linear characteristics are shown in Section 5. Lastly, conclusions are provided in Section 6.

2. Nonlinear Aeroelastic Modeling

The aeroelastic model consists of a two degrees of freedom system simulating the nacelle and a rigid rotor simulating the propeller [58]. The system is represented by a spring-mass-damper system with a coupling in the yaw and pitch with propeller rotation in a freestream flow environment. A schematic of the system under investigation is presented in Figure 1, illustrating the rotation of the propellers, stiffness, and damping in the yaw and pitch directions of the nacelle motion, and fluid flow normal to the initial state of the system.

The governing equations of motion, given in Equation (1), outline the response of the system in pitch and yaw [8].

$$\left[\begin{array}{cc}{I}_n& 0\\ 0& I_n\end{array}\right]\left\{\begin{array}{c}\ddot{\theta }\\ \ddot{\psi }\end{array}\right\}+\left[\begin{array}{cc}{C}_{\theta }& -I_x\mathsf{\Omega }\\ I_x\mathsf{\Omega }& C_{\psi }\end{array}\right]\left\{\begin{array}{c}\dot{\theta }\\ \dot{\psi }\end{array}\right\}+\left[\begin{array}{cc}{K}_{\theta }& 0\\ 0& K_{\psi }\end{array}\right]\left\{\begin{array}{c}\theta \\ \psi \end{array}\right\}=\left[\begin{array}{c}{F}_{\theta }\\ F_{\psi }\end{array}\right]$$

(1)

where the forcing functions are depicted as aerodynamic reactions based on the geometry of the propellers and rotor angular speed [8].

Quasi-steady aerodynamic forcing for the pitch and yaw motions is considered, as illustrated in Equations (2) and (3). The aerodynamic loads consider the geometric and kinematic responses of the system. For the forcing terms, N_B is the number of blades, $\rho$ denotes the air density, $c_{l,\alpha }$ represents the blade lift slope, R is the rotor radius, $\mathsf{\Omega }$ represents the rotor angular velocity, $c$ is the blade chord length, $a$ is the pivot length to rotor radius ratio, and $V$ is the freestream velocity [64,65]. It should be mentioned that the reference scenario selected for this study has the geometric and dynamic properties shown in

Table 2. Geometric and dynamic parameters of the system [4,8].

Description	Symbol	Value
Rotor radius	R	0.152 m
Rotor angular velocity	$\mathsf{\Omega }$	2291.83 deg-s−1
Freestream velocity	V	6.7 m s−1
Pivot length to rotor radius ratio	$a$	0.25
Rotor moment of inertia	$I_x$	0.000103 kg m2
Nacelle moment of inertia	$I_n$	0.000178 kg m2
Structural pitch damping	$C_{\theta }$	0.05729578 Nm s deg−1
Structural pitch stiffness	$K_{\theta }$	22.9183 Nm deg−1
Structural yaw damping	$C_{\psi }$	0.05729578 Nm s deg−1
Structural yaw stiffness	$K_{\psi }$	22.9183 Nm deg−1
Number of blades	$N_B$	4
Blade chord	c	0.026 m
Blade lift slope	$c_{l,\alpha }$	0.10966 deg−1

[4,66].

$$\begin{gathered} F_{\theta }=\frac{N_B}{2}\left(\frac{1}{2}\rho c_{l,\alpha }{R}^4{\mathsf{\Omega }}^2\right)R\left[-\left(\left(\frac{c}{R}{{\displaystyle \int }}_0^1\frac{{\eta }^4}{\sqrt{{\left(\frac{V}{\mathsf{\Omega }\mathrm{R}}\right)}^2+{\eta }^2}}d\eta \right)+a^2\left(\frac{c}{R}{{\displaystyle \int }}_0^1\frac{{\eta }^2}{\sqrt{{\left(\frac{V}{\mathsf{\Omega }\mathrm{R}}\right)}^2+{\eta }^2}}d\eta \right)\frac{\dot{\theta }}{\mathsf{\Omega }}\right) \\ -\left(\frac{c}{R}{{\displaystyle \int }}_0^1\frac{{\left(\frac{V}{\mathsf{\Omega }\mathrm{R}}\right)}^2{\eta }^2}{\sqrt{{\left(\frac{V}{\mathsf{\Omega }\mathrm{R}}\right)}^2+{\eta }^2}}d\eta \right)\psi +a\mu \left(\frac{c}{R}{{\displaystyle \int }}_0^1\frac{{\eta }^2}{\sqrt{{\left(\frac{V}{\mathsf{\Omega }\mathrm{R}}\right)}^2+{\eta }^2}}d\eta \right)\theta \right] \end{gathered}$$

(2)

$$\begin{gathered} F_{\psi }=\frac{N_B}{2}\left(\frac{1}{2}\rho c_{l,\alpha }{R}^4{\mathsf{\Omega }}^2\right)R\left[-\left(\left(\frac{c}{R}{{\displaystyle \int }}_0^1\frac{{\eta }^4}{\sqrt{{\left(\frac{V}{\mathsf{\Omega }\mathrm{R}}\right)}^2+{\eta }^2}}d\eta \right)+a^2\left(\frac{c}{R}{{\displaystyle \int }}_0^1\frac{{\eta }^2}{\sqrt{{\left(\frac{V}{\mathsf{\Omega }\mathrm{R}}\right)}^2+{\eta }^2}}d\eta \right)\frac{\dot{\psi }}{\mathsf{\Omega }}\right) \\ -\left(\frac{c}{R}{{\displaystyle \int }}_0^1\frac{{\left(\frac{V}{\mathsf{\Omega }\mathrm{R}}\right)}^2{\eta }^2}{\sqrt{{\left(\frac{V}{\mathsf{\Omega }\mathrm{R}}\right)}^2+{\eta }^2}}d\eta \right)\theta +a\mu \left(\frac{c}{R}{{\displaystyle \int }}_0^1\frac{{\eta }^2}{\sqrt{{\left(\frac{V}{\mathsf{\Omega }\mathrm{R}}\right)}^2+{\eta }^2}}d\eta \right)\psi \right] \end{gathered}$$

(3)

In order to develop an effective mathematical model, an unmanned aerial vehicle with aerodynamic and structural properties listed in Table 2 is selected for this investigation. The following geometric and dynamic properties are based on the P-51 and P-47 scale model drones, for which the design conditions for the mathematical model are used. The P-51 and P-47 use a four-blade propeller with dimensions of 9.8 × 6”. The P-47 utilizes a 750 kV brushless motor with a 35A ESC and the P-51 utilizes an 850 kV brushless motor and a 30A ESC, resulting in cruise speeds between 6.6 m/s and 6.9 m/s and rotor angular velocities roughly between 2100 deg/s and 2450 deg/s. These design conditions can also be seen in the scaled model of XC-142A vertical takeoff and landing (VTOL) aircraft [67]. It should be noted that, although this particular mathematical model is based on scaled drones, this type of analysis will work on any rotary-wing manned or unmanned system, as shown in the aeroelastic analysis performed on the AW159 Wildcat Release to Service military document [68]. Ribner [64] presented a deviation of the forces and moments applied to a rotor nacelle system. When deriving the yawed motion based on the simple momentum theory, the second-order $\psi$ is neglected, thereby removing the coupling from the stiffness matrix. A method of averaging, developed by Mair et al [8], was used to select the values for the range provided by Reed III [4].

In order to solve the eigenvalues of the system concerning the equilibrium point, where 0 is a 4 × 4 zero matrix and I is a 4 × 4 identity matrix, the Jacobian matrix is used, as shown in Equation (4), with Equations (5)–(10) being substituted [69].

$$J=\left[\begin{array}{cc}\left[0\right]& \left[I\right]\\ -{\left[\begin{array}{cc}{I}_1& 0\\ 0& I_2\end{array}\right]}^{-1}\left[K\right]& -{\left[\begin{array}{cc}{I}_1& 0\\ 0& I_2\end{array}\right]}^{-1}\left[C\right]\end{array}\right]$$

(4)

$$\begin{gathered} \left[C\right] \\ =\left[\begin{array}{cc}{C}_{\theta }+\frac{c_{K,a}\left(A_3+a^2{A}_{1,i}\right)}{\mathsf{\Omega }}& -I_x\mathsf{\Omega }\\ I_x\mathsf{\Omega }& C_{\psi }+\frac{c_{K,a}\left(A_3+a^2{A}_{1,i}\right)}{\mathsf{\Omega }}\end{array}\right] \end{gathered}$$

(5)

$$\left[K\right]=\left[\begin{array}{cc}{K}_{\theta }-c_{K,a}a{A}_{1,i}& c_{K,a}{A}_2\\ -c_{K,a}{A}_2& K_{\psi }-c_{K,a}a{A}_{1,i}\end{array}\right]$$

(6)

$$c_{K,a}=\frac{1}{2}{N}_{B}R\left(\frac{1}{2}\rho c_{l,\alpha }{R}^4{\mathsf{\Omega }}^2\right)$$

(7)

$$A_{1,i}=\left(\frac{c}{R}\right){{\displaystyle \int }}_0^1\frac{{\left(\frac{U\left(i\right)}{\mathsf{\Omega }\mathrm{R}}\right)}^2}{\sqrt{ {\left(\frac{U\left(i\right)}{\mathsf{\Omega }\mathrm{R}}\right)}^2+n^2 }}dn$$

(8)

$$A_2=\left(\frac{c}{R}\right){{\displaystyle \int }}_0^1\frac{{\left(\frac{U\left(i\right)}{\mathsf{\Omega }\mathrm{R}}\right)}^2{n}^2}{\sqrt{ {\left(\frac{U\left(i\right)}{\mathsf{\Omega }\mathrm{R}}\right)}^2+n^2 }}dn$$

(9)

$$A_3=\left(\frac{c}{R}\right){{\displaystyle \int }}_0^1\frac{n^4}{\sqrt{ {\left(\frac{U\left(i\right)}{\mathsf{\Omega }\mathrm{R}}\right)}^2+n^2 }}dn$$

(10)

After solving the eigenvalues, the undamped natural frequency, $\omega$, and the damping ratio, $\zeta$, can be calculated using Equations (11) and (12), respectively.

$$\omega =\sqrt{Re{\left(\lambda \right)}^2+Im{\left(\lambda \right)}^2}$$

(11)

$$\zeta =\frac{-Re\left(\lambda \right)}{\omega }$$

(12)

Utilizing this mathematical model, several influencing physical variables can be investigated to determine which variable is most significant to initiate whirl flutter effects in the linear system [70]. It should be mentioned that the initial conditions are bound to a thousandth of a radian through each iteration or a reset will occur in an effort to reduce any parametric errors that may accumulate in the recursive nature of the analysis. This work analyzes the modal damping ratio against the ratio of freestream velocity to blade tip velocity for the response of several parameters. The following results are based on two possible whirl flutter modes, either clockwise with the rotor angular velocity (FW—forward) or counter-clockwise, opposite the rotor angular velocity (BW—backward).

3. Linear Characteristics of the System: Onset Speed of Whirl Flutter

Initially, a parametric study for the linear system is performed to determine the main aerodynamic and structural characteristics that affect the onset of instability. A nonlinear analysis is then performed, where nonlinearities are introduced into the nacelle structural properties, and the effects of different combinations of nonlinear functions and locations in the system are discussed. The parameters investigated in this study are the number of blades, length of blades, rotor angular velocity, chord length, structural yaw stiffness, structural pitch stiffness, and nacelle moment of inertia. The mathematical quantities for each of these conditions are independently varied, and the response on the modal damping ratio is plotted. The first variable tested is the number of blades; the resulting modal damping ratio curves are illustrated in Figure 2a, in which the red portion of the plot area is related to the unstable flutter region. The instability region shown in a red tint begins at a modal damping ratio of zero and continues through all negative values. The instability is shown from Equation (12) where the eigenvalues are mathematically represented in combination with the eigenvalue stability theory. As the number of blades is increased from two to six, there is a slight decrease in the U/V_tip of instability crossing. This indicates that as more blades are introduced to aerodynamic systems, whirl flutter is initiated at a slightly lower freestream velocity. Furthermore, it should also be noted that whirl flutter will only occur in the BW angular direction due to the real terms of the eigenvalues [70].

The second study investigates the effects of the blade radius, R, on the onset speed of flutter of the system, as depicted in Figure 2b. For this section, a four-blade system is used, and the length of each blade is varied simultaneously. Based on the results, increasing the length of the blades leads to a decrease in the required freestream velocity to induce flutter. However, there is a notable increase in the modal damping ratio as a result of the increased length, once again, keeping the same I_x, only the BW angular direction leads to flutter.

Next, the effect of altering the rotor’s angular velocity on the system’s instability is studied, as shown in Figure 2c. As expected, by increasing the rotor velocity, the onset of flutter begins at lower freestream velocities. This effect is because the velocity near the tip of the blade is increasing as the rotor velocity increases, which would inevitably make the equality of U/V_tip approach zero. The modal damping ratio is not significantly altered, and BW is the only configuration by which flutter is initiated, as verified by Mair et al. [8].

Concerning the influence of the chord length of the blade on the system’s instability, it is clear from the plotted curves in Figure 2d that a notable variation takes place in the flutter response by increasing the chord length. Indeed, an increase in the chord length results in an increase in the aerodynamic force and, hence, the flutter speed takes place at smaller speeds. Once again, only the BW direction causes flutter instability.

Next, the impacts of the moment of inertia regarding the rotor and nacelle, as well as the effects of the linear pitch and yaw on the onset of flutter, are investigated, as presented in Figure 3. Inspecting these plots, it is clear that if the rotor angular velocity is chosen in the clockwise direction (FW), flutter cannot take place and an increase in the coupled modal damping is obtained. Further, an increase in the nacelle or rotor moment of inertia results in a decrease in the onset speed of instability, as shown in Figure 3a,b, respectively. As for the yaw and pitch stiffnesses, they have a negligible effect on the flutter speed of the system, as depicted in Figure 3c,d, respectively.

Based on this linear investigation on the effects of the system’s parameters on the onset speed of instability, it is shown that the most influential parameters initiating whirl flutter at a lower U/V_tip velocity are rotor angular velocity and length of the blades. On the other hand, the other studied parameters have little effect on the whirl flutter speed. In particular, the linear yaw and pitch stiffnesses have a negligible effect on the onset of instability. It should be mentioned that the clockwise direction (FW) for the rotor angular velocity results in an increase of the total damping in the system and hence no flutter can be present.

4. Bifurcation Analysis: Effects of the Nonlinear Stiffnesses on the System’s Response

There are many sources of nonlinearities in systems subjected to whirl flutter, including material properties, control mechanisms, and fixtures to aerodynamic nonlinearities, such as separated flow and dynamic stall, among others. In this work, the focus is on the possible structural sources of nonlinearities. These nonlinearities can cause either the hardening or the softening of structural behaviors. These later behaviors can introduce several different post- and pre-flutter responses in the bifurcation diagrams of the rotor-based system. In the rest of this study, different nonlinear functions for the restoring moment for pitch and yaw degrees of freedom are considered. The parameters considered for the subsequent analyses can be seen in Table 2.

The first case examined is the cubic restoring moment, depicted in Figure 4a, included symmetrically in both pitch and yaw degrees of freedom. The linear term remains constant as the base value while the cubic term is varied. It follows from the plotted curves in Figure 4a that, as the nonlinearity is increased, there is a hardening effect characterized by an increase in the moment at higher deflections. Mathematically, the cubic term in this defining function creates a more significant moment version angle growth, which rapidly affects the data at values greater than 10°.

The second case of restoring moment considered is the linear-quadratic-cubic case, for which the linear and cubic terms are kept constant and the quadratic coefficient is varied, as depicted in Figure 4b. The introduction of this quadratic term results in an asymmetry in the restoring moment curve, although the restoring moment still presents a hardening effect. The same nonlinear restoring moment is again introduced in both pitch and yaw. The next nonlinear case investigated is the third- and fifth- order case while varying the cubic coefficient. As illustrated in Figure 4c, the first four cubic- and fifth-order cases resemble mainly hardening behaviors. However, the last case, $0.5\psi -20{\psi }^3+500{\psi }^5$, has initial softening between 0° and 5°, after which it begins to exhibit hardening. This case within the cubic- and fifth-order investigation of this system introduces a new parameter to the study by making the softening cubic term dominate the given angles. This cubic softening behavior is investigated to understand its effects on the angles of pitch and yaw on the system.

The final case considered in this work investigates the linear, quadratic, cubic, and fifth-order representation for the yaw and pitch restoring moments. As illustrated in the moment versus angle plot in Figure 4d, the initial deflection curve represents softening behavior that transitions into hardening. It should be noted that the cubic term in this case is negative, which causes the slope of the deflection to be concave downward. It is further understood that the cubic- and fifth-order terms produce inflection points after about 4°, which creates an upward concave curve.

The nonlinear response of this system is shown by the bifurcation diagrams, illustrated in Figure 5, presenting the variation of the maximum angular deflection as a function of the freestream velocity. The presence of the hardening structural nonlinearity in the system introduces limit cycle oscillations to the system in the post-flutter regime. The symmetry of these nonlinear terms in the two degrees of freedom causes an equivalent response in both the pitch and yaw deflections. Considering the linear-cubic case of nonlinearity, Figure 5a shows the Hopf bifurcation occurring at roughly 5.66 m/s, with a supercritical type behavior [71,72,73,74]. As the base value of the cubic term increases, the maximum angular deflection decreases. This response indicates the increase of the hardening effect, which causes the full system to reduce the angular deflection. Additionally, the maximum freestream velocity is limited because there is a possibility of unpredicted stall effects that makes freestream velocities over 10 m/s unfeasible and unrealistic.

As observed from the plotted curves in Figure 5b, the inclusion of the positive quadratic term has a negligible effect on the bifurcation diagram of the system, which can be seen by the lower maximum deflection at the flutter point and end of the simulation. As the nonlinear terms are included in both pitch and yaw, the dynamical response in both degrees of freedom look almost the same. The quadratic nonlinearity in the described case is further examined, as illustrated in Figure 6b. It is clear from the time histories shown in Figure 6b that the pitch and yaw angular deflections are not in agreement. Based on the coupled system proposed, there are inherent asymmetries in the aerodynamic force. However, the coupling between $\theta$ and $\psi$ is not primarily responsible for the described behavior shown in both the time histories and phase portraits in Figure 6b, which will be discussed in-depth later. This behavior is known as the drift effect, which is caused by the presence of the quadratic nonlinearity.

As previously mentioned, the possibility of obtaining a subcritical Hopf bifurcation may exist within linear-cubic-fifth and linear-quadratic-cubic-fifth representations of the restoring moment versus the yaw or pitch angle. This is due to the presence of the softening effect in these representations. Therefore, the bifurcation diagrams’ simulations are run forward and backward in order to detect the subcritical instability if it exists. The results for the pitch deflections in the forward and backward simulations are overlaid with each other and illustrated in Figure 5c,d for both representations. The results clearly show the presence of the subcritical instability, for which there is a hysteresis region when the nonlinear softening effects become more prominent. These results indicate that there are critical values of the negative cubic nonlinearity at which the system changes instability from supercritical to subcritical. It should be mentioned that the positive quadratic nonlinearity, included in the fourth representation (Figure 5c), leads to a reduction in the softening effect nonlinearity, and hence the hysteresis region becomes smaller and smaller when the quadratic nonlinearity is increased.

As briefly mentioned previously, the drift effect is an underlying reason for the unusual response of the system when the quadratic nonlinearity is present in the restoring moment representation, as shown in Figure 5b. The time history of the steady-state response for quadratic and cubic cases previously described is depicted in Figure 6. For each of the five cases, the drift effect is an apparent and distinguished difference between the coupled $\theta$ and $\psi$ deflections, as shown in Figure 6c. In order to visualize the drift effect in case of the linear-quadratic-cubic configuration, the fifth scenario is analyzed, and the results are shown in Figure 6d. Therefore, the quadratic and cubic nonlinear stiffness scenario is reevaluated considering the peak-to-peak measurement. With the inclusion of this shift, the two angular deflections in the pitch and yaw directions follow the same trends and very closely overlap.

Further investigation into the time series of the fifth case is shown in Figure 6b. In the defined scenario, measurements for the angle are made considering the maximum amplitudes. However, as shown in the time history plot, the maximum values of the pitch and yaw angles occur at separate maximum values due to the drift effect and phase difference. In all cases, peak-to-peak amplitude measurements are necessary for accurate results when the quadratic nonlinearity is present in the restoring moment.

The maximum angles for the pitch and yaw at two critical points are illustrated in

Table 3. Points of maximum deflection for all scenarios.

Equation	Max (deg) @ Point after Onset of Flutter (6 m/s) Pitch-θ	Max (deg) @ End of Simulation (10 m/s) Pitch-θ	Max (deg) @ Point after Onset of Flutter (6 m/s) Yaw-ψ	Max (deg) @ End of Simulation (10 m/s) Yaw-ψ
Cubic nonlinearities
$0.5\phi +100{\phi }^3$	1.753	6.720	1.753	6.720
$0.5\phi +200{\phi }^3$	1.239	4.752	1.239	4.752
$0.5\phi +300{\phi }^3$	1.012	3.880	1.012	3.880
$0.5\phi +400{\phi }^3$	0.876	3.360	0.876	3.360
$0.5\phi +500{\phi }^3$	0.784	3.005	0.784	3.005
Quadratic & cubic nonlinearities
$0.5\phi +0.1{\phi }^2+100{\phi }^3$	1.753	6.720	1.753	6.720
$0.5\phi +0.5{\phi }^2+100{\phi }^3$	1.752	6.720	1.752	6.720
$0.5\phi +1.0{\phi }^2+100{\phi }^3$	1.748	6.720	1.748	6.720
$0.5\phi +5.0{\phi }^2+100{\phi }^3$	1.733	6.725	1.735	6.725
$0.5\phi +10{\phi }^2+100{\phi }^3$	1.360	6.870	1.390	6.880
Cubic & fifth-order nonlinearities
$0.5\phi -0.10{\phi }^3+2500{\phi }^5$	4.6525	6.695	4.653	6.695
$0.5\phi -1.0{\phi }^3+2500{\phi }^5$	4.7255	6.755	4.726	6.755
$0.5\phi -5.0{\phi }^3+2500{\phi }^5$	5.065	7.000	5.066	7.000
$0.5\phi -10{\phi }^3+2500{\phi }^5$	5.503	7.305	5.504	7.305
$0.5\phi -20{\phi }^3+2500{\phi }^5$	6.398	7.370	6.399	7.370
Quadratic, cubic, and fifth-order nonlinearities
$0.5\phi +0.05{\phi }^2-20{\phi }^3+500{\phi }^5$	12.389	14.735	12.392	14.735
$0.5\phi +0.10{\phi }^2-20{\phi }^3+500{\phi }^5$	12.386	14.735	12.388	14.735
$0.5\phi +0.20{\phi }^2-20{\phi }^3+500{\phi }^5$	12.197	14.735	12.195	14.725
$0.5\phi +0.50{\phi }^2-20{\phi }^3+500{\phi }^5$	12.188	14.730	12.187	14.710
$0.5\phi +1.00{\phi }^2-20{\phi }^3+500{\phi }^5$	11.256	14.950	11.257	14.845

for the linear-cubic case. The first maximum point shown occurs at a freestream velocity of 6 m/s, considered a point after flutter begins. The second maximum point occurs at 10 m/s, which is selected as the limit of applicability for this mathematical model due to the linear representation of the aerodynamic forces. Every wind speed beyond the 10 m/s threshold is uncertain due to the possibility of the stall effect within the system, which the model does not account for. Each of the four scenarios is outlined in Table 3, with the maximum deflection angle in degrees after the point of flutter, 6 m/s, and at the end of the simulation, 10 m/s, although it should be noted that, as illustrated in Figure 5, not all the bifurcation diagrams are considered to the end of simulation time at 10 m/s. However, this data is still included in the table to illustrate the parametric trends. For instance, in the linear-quadratic-cubic-fifth order configuration, the maximum deflection at the end of the simulation would be around 15°, which, as mentioned before, would introduce the possibility of the stall effect, which has not been accounted for in this mathematic model.

The phase portraits of the pitch and yaw degrees of freedom when considering the four restoring moment representations are illustrated in Figure 7. Clearly, the overall stability of the defined system remains stable within the boundaries shown for each case. It should be noted that the phase portraits are considered at the maximum freestream velocity in each of the four scenarios from Figure 5. The trajectories in phase space suggest an asymmetric behavior, and it is essential to mention that the nonlinearity is symmetric (cubic) and present in pitch and yaw symmetrically. The aerodynamic behavior introduces this asymmetry under propeller rotation and nacelle motion.

The maximum angular deflection for pitch and yaw degrees of freedom based on linear, quadratic, and cubic nonlinear stiffnesses are illustrated in Figure 7b. Based on the bifurcation diagram, the flutter speed occurs at roughly 5.66 m/s, as it has in the previous case due to the fundamental nature of the system with the presence of the supercritical instability. However, unlike the previous scenarios, there is a noticeable dip in the maximum deflection at ~7.6 m/s for case 5. The behavior that occurs at 7.6 m/s may indicate the presence of a secondary bifurcation. To test this, a study of the phase portraits before and immediately after this dip is investigated. Furthermore, it can be seen in Figure 7c that there is a sudden step in the trends of the data at the point of flutter, which could indicate subcritical behavior in the system. As for the fourth representation of the restoring moment, the phase portraits in Figure 7d illustrate the complex behavior of the system. In the phase portrait for the pitch and yaw angles, there are three discernable orbits: a global orbit spanning from −7.5 to 7.5 degrees and two local orbits within.

5. Effects of the System’s Properties and Number of Blades on the Nonlinear Response of the System

A study based on the linear results is considered for the nonlinear first case of the linear-cubic-fifth-order configuration, shown in green in Figure 7, $0.5\phi -0.1{\phi }^3+2500{\phi }^5$. It should be noted that the parameters investigated for this particular case were run with both increasing and decreasing velocities to discover the effects of each parameter on the subcritical or supercritical instability of the system under investigation. The impact of the nacelle moment of inertia is first investigated, as depicted in Figure 8a. It is clear that, as the nacelle moment of inertia increases, a slight variation from a supercritical to subcritical becomes apparent. This is shown by the two extremes for this parameter, I_n equal to 0.0001 kg-m², shown in red, and on the other end, I_n equal to 0.0006 kg-m², shown in cyan. The forward and backward simulations overlap and show supercritical behavior at the point of flutter. However, the forward and backward simulations, shown in green, exhibit subcritical Hopf bifurcation. This bifurcation diagram plot shows a clear trend of increasing subcritical behavior in the system as the nacelle moment of inertia increases, although the transition rate is relatively slow. It should be mentioned that the increase of the nacelle moment of inertia results in an increase in the whirl flutter speed, as demonstrated in the linear analysis Figure 3a.

The second parameter investigated in this study is the rotor moment of inertia. The corresponding bifurcation diagrams are shown in Figure 8b. Inspecting this figure, the increase in the rotor moment of inertia is followed by a decrease in the linear flutter speed, as expected from the linear analysis Figure 3b. Additionally, it is noted that, as the rotor moment of inertia increases, a much faster rate from a supercritical to subcritical instability takes place compared to the nacelle moment of inertia effects. This is illustrated by considering the two boundaries for this parameter, I_x equal to 0.00005 kg-m², shown in red, and I_x equal to 0.0006 kg-m², shown in cyan. The forward and backward simulations in red overlap and show supercritical Hopf bifurcation. However, the cyan forward and backward simulations exhibit subcritical behavior. This plot shows a clear trend of rapidly increasing subcritical behavior in the system as the rotor moment of inertia increases.

The effects of varying the length of the blades (radius) using a four-blade system as the base case are studied in Figure 8c. As demonstrated in the linear analysis Figure 2b, an increase in the blade length is followed by linear flutter speed reduction, which is also observed in the bifurcation diagrams shown in Figure 8c. Further, as the length of the blade increases, the system transitions from a slightly subcritical instability to a supercritical one. It should be noted that the radius shown in Figure 2b of 0.1016 m is unable to run at the given nonlinear parameters of this system. Indeed, the linear flutter speed for this specific radius is higher than 10 m/s. The results of this parameter study suggest a nonlinearly varying flutter speed as the length of the blades increases, as well as the convergent behavior to a supercritical instability.

The plotted bifurcation diagrams shown in Figure 8d show the influence of the number of blades on the system’s instability. This portion of the study investigates the effects from having two to six blades when considering a constant base blade length from Table 2. Clearly, as showed in the linear analysis Figure 2a, the whirl flutter speed increases when the number of blades is decreased. Further, the subcritical nonlinearity is more prominent when the number of blades is decreased, and the flutter speed is higher.

Next, the impacts of the system’s properties and number of blades on the system’s response and instability are investigated when considering the fifth case of the linear-cubic-fifth representation, shown in pink in Figure 5, $0.5\phi -20{\phi }^3+2500{\phi }^5$. It should be mentioned that this scenario showed subcritical Hopf bifurcation in Figure 5. Similar to the plots shown in Figure 8, the same properties are studied when considering forward and backward speeds in order to determine the effects of each parameter on the instability of the given system. Furthermore, the parametric studies for all cases show a variation in the freestream velocity at which flutter occurs. Considering that the case used for this parameter study was initially subcritical at the given base parameters, the effects of each parameter are expected to follow similar trends.

It follows from the plots shown in Figure 9 that the subcritical instability takes place in all configurations, and all parameters investigated do not change the type of instability. However, their effects shown in Figure 8 are present when considering the fourth representation for the restoring moment. Indeed, an increase in the nacelle moment of inertia has a negligible impact on the system’s instability, and all hysteresis regions are similar and only linear effects of the nacelle moment of inertia are present, as depicted in Figure 9a. Concerning the rotor moment of inertia effects on the type of instability and system’s response, it follows from Figure 9b that the linear influence of this parameter is dominant on the flutter speed. In addition, similar to the previous effects shown in Figure 8b, the subcritical behavior becomes more prominent for higher values of the rotor moment of inertia.

The plots in Figure 9c,d illustrate the impacts of the blade length and number of blades on the system’s performance and stability. Inspecting these figures, it is observed that a reduction in the blade length or the blade number result in an increase in the linear flutter speed of the system as well as a more notable subcritical effect on the system’s instability. This trend is particularly important when the number of blades is selected to be two, as shown in Figure 9d.

6. Conclusions

The influence of drones on modern society can be seen around the world, from industry and farming to warfare. By understanding the phenomenon known as whirl flutter, designers and manufactures can improved system lifespans and avoid catastrophic failure. Unmanned ariel systems are taking aviation to new heights, limited only by the structural and aerodynamic parameters rather than the pilot’s ability to handle these flight conditions. This is why, now more than ever, it is important to investigate the structural aeroelastic limitations of these systems and address the phenomenon that can arise. Specifically, whirl flutter, an aeroelastic instability present in propeller-driven aircraft, is a difficult phenomenon to predict. Nonlinearities in the system, both structural and aerodynamic, can contribute to the onset of whirl flutter, complicating its prediction. In this study, a nonlinear reduced-order model was used to investigate the effects of different system’s parameters on whirl flutter and the system’s response and instability to structural nonlinearities. A parametric study was carried out and properties that have a larger impact on the onset of whirl flutter were identified and evaluated. It was found that more blades, shorter blade lengths, increased angular velocity, and increased chord length all induce a lower critical linear flutter velocity, meaning that these are crucial features to optimize and focus on during the designing phase. Other property effects on the system’s linear stability were explored. It was demonstrated that the pitch and yaw stiffnesses have a negligible effect on the linear whirl flutter speed.

Considering one system configuration, the influence of different structural nonlinearities was investigated. Multiple polynomial nonlinearities were introduced into the system, including symmetric cubic hardening nonlinearities for the pitch and yaw degrees of freedom; purely yaw nonlinearity; purely pitch nonlinearity; and a combination of quadratic, cubic, and fifth-order nonlinearities for both degrees of freedom. It was shown that hardening structural nonlinearities introduce limit cycle oscillations in the post-flutter regime. Including a quadratic nonlinearity causes asymmetric oscillations due to the presence of the drift effect. Including softening cubic nonlinearity resulted in the presence of subcritical Hopf bifurcations, which can cause large amplitude deformations before the predicted flutter speed. It was shown that there is a significant effect of the system’s properties and number of blades on the type of instability, depending on the present structural nonlinearities in the system.