<figure><img src="image_1.png" /></figure>

HAL
open science

# Computational methods for whirl flutter analysis

Charles Jacquet, Camille Denoel, Loïc Salles

► To cite this version:

Charles Jacquet, Camille Denoel, Loïc Salles. Computational methods for whirl flutter analysis. 11th International Operational Modal Analysis Conference, Institut National de Recherche en Informatique et en Automatique [Inria], May 2025, Rennes, France. hal-05398750

HAL Id: hal-05398750

https://hal.science/hal-05398750v1

Submitted on 4 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Computational Methods for Whirl Flutter Analysis

Charles Jacquet<fnref n="1" />, Camille Denoël<fnref n="2" />, Loïc Salles<fnref n="2" />,

## ABSTRACT

The rise of Electric Vertical Take-Off and Landing (E-VTOL) vehicles, driven by the demand for fast urban transport, emphasizes the need to understand and mitigate aeroelastic instabilities like whirl-flutter, particularly in designs with distributed electric propulsion systems. Historically, whirl-flutter has caused catastrophic failures, highlighting the need for comprehensive analysis in modern aircraft designs. This work aims to develop a comprehensive computational model using the Finite Element Method to fully capture the dynamics between rotating, expressed in a floating frame of reference, and stationary parts of a structure, enabling a more accurate study of the whirl-flutter phenomenon. This analysis uses the Floquet theory to study the system's stability, particularly in a fixed-rotating frame of reference. The model includes an innovative time-dependent mechanical coupling strategy for the mass, gyroscopic, and centrifugal stiffness structural matrices, therefore fully preserving the dynamics of the structure. The methodology involves the validation of various finite element matrices and components derived from 3D beam elements, followed by the implementation of time-dependent coupling between rotating and stationary components. The developed model is applied to a wing-propeller structure, illustrating its capability to work on complex geometry structures. The results show that the partial coupling between translational degrees of freedom of rotating structures and those of the stationary structure at the hub, is successfully validated. However, it is demonstrated that the Newmark integration scheme does not provide consistent results, highlighting the need for alternative approaches for accurate time integration. This partial success demonstrates the potential of the developed model as a tool for accurately capturing critical dynamics in whirl-flutter analysis, contributing to the design and certification of future urban air mobility vehicles.

Keywords: Whirl Flutter, stability analysis, finite element method, Floquet's theory

<footnote n="1">¹ The von Karman Institute for Fluid Dynamics, Belgium</footnote>

<footnote n="2">² Laboratory Vibration of Turbomachines, University of Liege, Belgium, l.salles@uliege.be</footnote>
# 1. INTRODUCTION

Whirl flutter is an aeroelastic instability that arises in rotor-wing systems due to coupling between aero-dynamic forces, gyroscopic effects, and structural flexibility. This phenomenon has gained importance with the rise of Electric Vertical Take-Off and Landing (E-VTOL) vehicles, where distributed propulsion systems exacerbate stability challenges. Historically, whirl flutter has led to catastrophic failures in propeller-driven aircraft, highlighting the need for accurate computational models.

Previous studies have employed linear time-invariant (LTI) and linear time-periodic (LTP) approaches to analyze whirl flutter. While LTI methods simplify the problem, they fail to capture the time-dependent behavior of rotating systems. LTP methods, particularly those based on Floquet theory, allow for a more comprehensive stability analysis by considering the periodic nature of the problem.

This paper introduces a computational approach using FEM and Floquet theory to model and analyse whirl flutter in rotor-wing systems more accurately. The structure of the paper is as follows: Section 2 details the methodology, including FEM implementation and stability analysis techniques. Section 3 presents numerical results and validation.

# 2. METHODOLOGY

The computational framework for whirl flutter analysis consists of:

## 2.1. Finite Element Modeling

The mathematical derivation is mainly based on the work produced by Vollan and Komzsik [1] with additional information and mathematical development being included.

## 2.2. Full coupling of stationary and rotating parts

Before going into the mathematical derivation of the coupling for the finite element formulation of the problem. It is of interest to understand beforehand, the coupling of particles of mass m. Less complete type of coupling also exists, they might be useful to capture some physical phenomenon. However, there are not of interest in the developed model and therefore only a focus on the full coupling of the particle of mass m is achieved. The coupling of stationary and rotating parts allows translational and rotational displacement on both, stationary and rotating particles. A schematic of the studied coupling is highlighted in Figure 1.

<figure><img src="image_3.png" /><figcaption>**Figure 1:** Schematic of the full coupled stationary and rotating parts. Reproduced from [1].</figcaption></figure>

In Figure 1, $\sigma$, $\rho$ and $A_\alpha$ represents respectively the translation vector of the stationary part, the trans-
lation vector of the rotating part and transformation matrix of the rotational deformation of the rotating part. The location vector of a rotating particle reads, in its most general form, as:

$$ \bar{\mathbf{r}} = \sigma + \mathbf{B}\mathbf{H}\mathbf{r} + \mathbf{H}\rho + \mathbf{H}\mathbf{A}\alpha + \mathbf{H}\mathbf{r}, \qquad (1) $$

where the {$\bar{\{\}}$} makes reference to the fixed coordinate system. The matrix **B** corresponds to the transformation matrix of the small nodal translations. This matrix reads as:

$$ \mathbf{B} = \begin{bmatrix} 0 & -\bar{\theta} & \bar{\psi} \\ \bar{\theta} & 0 & -\bar{\varphi} \\ -\bar{\psi} & \bar{\varphi} & 0 \end{bmatrix}. \qquad (2) $$

It is also of interest to introduce the transformation matrix **H** which allows to switch from a location vector located in the rotating frame of reference to the fixed coordinate system. This transformation matrix reads:

$$ \mathbf{H} = \begin{bmatrix} \cos \Omega t & -\sin \Omega t & 0 \\ \sin \Omega t & \cos \Omega t & 0 \\ 0 & 0 & 1 \end{bmatrix}. \qquad (3) $$

One also introduces the matrix **A**, which is simply the matrix containing the coordinates of rotating nodes. This matrix reads

$$ \mathbf{A} = \begin{bmatrix} 0 & z' & -y' \\ -z' & 0 & x' \\ y' & -x' & 0 \end{bmatrix}. \qquad (4) $$

### 2.2.1. Mathematical properties of the transformation matrix **H**

It is of interest to introduce some mathematical properties of the transformation matrix **H**.

$$ \mathbf{H}^T \mathbf{H} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} = \mathbf{I}; \qquad (5) $$

$$ \dot{\mathbf{H}} = \Omega \begin{bmatrix} -\sin \Omega t & -\cos \Omega t & 0 \\ \cos \Omega t & -\sin \Omega t & 0 \\ 0 & 0 & 0 \end{bmatrix} = \Omega \bar{\mathbf{H}}; \qquad (6) $$

$$ \ddot{\mathbf{H}} = \Omega^2 \begin{bmatrix} -\cos \Omega t & \sin \Omega t & 0 \\ -\sin \Omega t & -\cos \Omega t & 0 \\ 0 & 0 & 0 \end{bmatrix} = \Omega^2 \bar{\bar{\mathbf{H}}}; \qquad (7) $$

$$ \bar{\mathbf{H}}^T \mathbf{H} = \begin{bmatrix} 0 & 1 & 0 \\ -1 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} = \mathbf{P}; \qquad (8) $$

$$ \mathbf{H}^T \bar{\mathbf{H}} = \begin{bmatrix} 0 & -1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} = \mathbf{P}^T = -\mathbf{P}; \qquad (9) $$

$$ \bar{\mathbf{H}}^T \bar{\mathbf{H}} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{bmatrix} = \mathbf{J} = \bar{\mathbf{H}} \bar{\mathbf{H}}^T \qquad (10) $$
### 2.2.2. Lagrange equation of motion

In order to take into account the rotation of the stationary part, the following transformations allow the rotation of the stationary part to be taken into account:

$$ 
\begin{aligned}
\mathbf{BHr} &= \begin{bmatrix} 0 & -\bar{\theta} & \bar{\psi} \\ \bar{\theta} & 0 & -\bar{\varphi} \\ -\bar{\psi} & \bar{\varphi} & 0 \end{bmatrix} \begin{bmatrix} \cos \Omega t & -\sin \Omega t & 0 \\ \sin \Omega t & \cos \Omega t & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{Bmatrix} x \\ y \\ z \end{Bmatrix} \\
&= \begin{bmatrix} 0 & -\bar{\theta} & \bar{\psi} \\ \bar{\theta} & 0 & -\bar{\varphi} \\ -\bar{\psi} & \bar{\varphi} & 0 \end{bmatrix} \begin{Bmatrix} x \cos \Omega t - y \sin \Omega t \\ x \sin \Omega t + y \cos \Omega t \\ z \end{Bmatrix} \\
&= \begin{Bmatrix} -\bar{\theta}(x \sin \Omega t + y \cos \Omega t) \\ \bar{\theta}(x \cos \Omega t - y \sin \Omega t) \\ -\bar{\psi}(x \cos \Omega t - y \sin \Omega t) + \bar{\varphi}(x \sin \Omega t + y \cos \Omega t) \end{Bmatrix} \\
&= \begin{bmatrix} 0 & 0 & -x \sin \Omega t - y \cos \Omega t \\ 0 & 0 & x \cos \Omega t - y \sin \Omega t \\ x \sin \Omega t + y \cos \Omega t & -x \cos \Omega t + y \sin \Omega t & 0 \end{bmatrix} \begin{Bmatrix} \bar{\varphi} \\ \bar{\psi} \\ \bar{\theta} \end{Bmatrix} \\
&= \mathbf{B}_0 \beta.
\end{aligned} 
\quad (11) $$

As a result, the augmented generalised coordinate vector reads:

$$ g = \begin{Bmatrix} \sigma \\ \beta \\ \rho \\ \alpha \\ \mathbf{r} \end{Bmatrix}, \quad (12) $$

and the governing matrix **M** of the motion becomes:

$$ \mathbf{M} = [ \mathbf{I} \quad \mathbf{B}_0 \quad \mathbf{H} \quad \mathbf{HA} \quad \mathbf{H} ]. \quad (13) $$

As a result, making use of all the previously established relationships, it is possible to derive the Lagrange
equation of motion for the fully coupled case. The equation reads as

$$
\begin{equation}
\begin{aligned}
& \begin{bmatrix}
[M_\sigma] & m[B_0] & m[H] & m[H][A] \\
m[B_0]^T & m[B_0]^T[B_0] & m[B_0]^T[H] & m[B_0]^T[H][A] \\
m[H]^T & m[H]^T[B_0] & [M_\rho] & m[A] \\
m[A]^T[H]^T & m[A]^T[H]^T[B_0] & m[A]^T & [M_\alpha]
\end{bmatrix}
+
\begin{bmatrix}
0 & m[\bar{B}_0] & m[\bar{H}] & m[\bar{H}][A] \\
0 & m[B_0]^T[\bar{B}_0] & m[B_0]^T[\bar{H}] & m[B_0]^T[\bar{H}][A] \\
0 & m[H]^T[\bar{B}_0] & [C_\rho] & m[P]^T[A] \\
0 & m[A]^T[H]^T[\bar{B}_0] & m[A]^T[P]^T & [C_\alpha]
\end{bmatrix}
+
2\Omega
+
\begin{bmatrix}
0 & m[\bar{B}_0] & m[\bar{H}] & m[\bar{H}][A] \\
0 & m[B_0]^T[\bar{B}_0] & m[B_0]^T[\bar{H}] & m[B_0]^T[\bar{H}][A] \\
0 & m[H]^T[\bar{B}_0] & [Z_\rho] & -m[J][A] \\
0 & m[A]^T[H]^T[\bar{B}_0] & -m[A]^T[J] & [Z_\alpha]
\end{bmatrix}
+
\Omega^2
+
\begin{bmatrix}
-m[\bar{H}]\{r\} \\
-m[B_0]^T[\bar{H}]\{r\} \\
\{f_{cp}\} \\
\{f_{c\alpha}\}
\end{bmatrix}
=
\Omega^2
\begin{bmatrix}
\ddot{\sigma} \\
\ddot{\beta} \\
\ddot{\rho} \\
\ddot{\alpha}
\end{bmatrix}
\end{aligned}
\tag{14}
\end{equation}
$$

## 2.3. Numerical implementation of the coupling between stationary and rotating parts

In addition to understanding the phenomena involved, this work's main difficulty is implementing this static/rotational coupling numerically. The aim of this sub-section is therefore to explain the numerical implementation of the coupling and the points that need to be taken into account to reproduce the results obtained.

The whirl-flutter phenomenon is due to a coupling between a rotating object and its attachment point. It is possible, particularly in the case of a wing and propeller system, that there is a static part that is not rotating (the wing for aerospace structure) which has no influence on the coupling but does influence the dynamics of the results obtained by having a cohesion of the deformation modes between the wing and the propeller. The wing is not taken into account in the coupling algorithm. The propeller system can comprise a pylon attached to the wing and a nacelle made up of a rotor and in the case of the study carried out. The propeller is modelled as a pylon (flexible or not) with a hub at its end where the various blades of the propeller are attached to it. The coupling occurs between the blades making up the propeller and the hub. It is important to remember that the dynamics of the blades are not directly linked in this study. In conclusion, coupling occurs only on the propeller system. At the same time, the wing plays a role in the structure expressed in the inertial frame of reference and plays a role in the global dynamics and stability of the structure.

To visualize the latter, a representation of the coupling for the mass structural matrix that undergoes a coupling between rotating and stationary components is shown in Figure 2. Analysis of the Figure reveals the coupling between the rotor hub and the blades. As well as the coupling between the blades and the rotor. All this shows that the coupling only affects the propeller system and not the stationary wing.

### 2.3.1. Mass matrix coupling

The coupling algorithm for the mass matrix is described here below. First of all, the mass matrix of the blade, i.e. rotating structure coupling in Figure 2, is achieved. For the mass matrix, no coupling takes
place on this structural submatrix. Therefore, the elementary mass matrix for this submatrix is expressed as

$$ M_{el} = \rho A \int_{0}^{l} \rho A N^{T} N dx. $$

To assemble the different elementary matrices, expressed in the element frame, must before being assembled, be expressed in the rotating frame of reference. This is done thanks to the transformation matrix $T_e$. The transformation from the element frame of reference to the rotating frame of reference is done by

$$ M_{eS} = T_{e}^{T} M_{el} T_{e}. $$

The assembly concerns all degrees of freedom of the element.

Concerning the coupling between the Hub and the blades rotating, the coupling matrix that must be taken into account when looking at the coupling between the blades and the rotating hub translational degree of freedom is the matrix $H$. As the different elementary matrices are expressed and computed in the element frame of reference, the coupling matrix $H$, which is expressed in the rotating frame of reference must transform to the element frame of reference thanks to the transformation matrix $R_e$. This transformation is expressed as

$$ H_{el} = R_e H R_e^T. $$

Therefore, the elementary mass matrix, when tackling the coupling between the hub and the blades rotating reads

$$ M_{el} = \rho A \int_{0}^{l} N^{T} R_e H R_e^T N dx. $$

To assemble in the rotating frame of reference, the elementary mass matrix $M_{el}$ must be transformed from the elementary to the rotating frame of reference. This is done similarly as before

$$ M_{eS} = T_{e}^{T} M_{el} T_{e}. $$

The assembly is made by adding to the translational degrees of freedom of the hub (located on the rows) and on the degrees of freedom of the rotating elements (the blades), the elementary mass matrix expressed in the rotating frame, the first three rows and all the columns of the first three rows of $M_{eS}$.

Similarly, the coupling between the blades and the hub is achieved using this time the coupling matrix $H^T$. As what was done previously, the coupling matrix which is expressed in the rotating frame of reference must be transformed in the local frame of reference. Therefore, the elementary mass matrix reads:

$$ M_{el} = \rho A \int_{0}^{l} N^{T} R_e H^{T} R_e^T N dx. $$

To assemble in the rotating frame of reference, the elementary mass matrix $M_{el}$ must be transformed from the elementary to the rotating frame of reference. This is similarly done as before

$$ M_{eS} = T_{e}^{T} M_{el} T_{e}. $$

The assembly is made by adding on the degrees of freedom of the rotating elements (the blades) and the translational degrees of freedom of the hub (located on the rows), the elementary mass matrix expressed in the rotating frame of reference; the first three columns and all the rows of $M_{eS}$.

Once all the elements requiring coupling have been traversed, the various structural sub-matrices are assembled. Before that, the coupling matrices and the matrix corresponding to the rotating structure (in our case, the blades) must have the rows and columns corresponding to the hub's degrees of freedom being removed. In fact, in the rotating frame of reference, the blades are considered to be clamped to the rotating hub.
<figure><img src="image_4.png" /><figcaption>**Figure 2:** Representation of a structural matrix undergoing time coupling between the rotating structure and the hub. The different rows and columns (not represented) represents the different degrees of freedom of the complete structure.</figcaption></figure>

### 3. NUMERICAL RESULTS AND DISCUSSION

To validate the proposed method, a rotor-wing system with known stability characteristics is analyzed. The key findings are as follows:

### 4. SIMPLE VALIDATION MODEL

In the first instance, the aim is to validate the coupled finite element model established previously partially. The coupling is tested on a ground resonance model established by Hammond et al [2]. The studied model is represented in Figure 3.

<figure><img src="image_5.png" /><figcaption>**Figure 3:** Validation case: Ground resonance model. Reproduced from Hammond et al [2].</figcaption></figure>

The different parameters used to perform the validation are available in Table 1.
<table><caption>**Table 1:** Ground resonance model parameters (Hammond et al [2]).</caption><thead><tr><th>Parameter</th><th>Value</th><th>Units</th></tr></thead><tbody><tr><td>Number of blade N</td><td>4</td><td>-</td></tr><tr><td>Rotor radius R</td><td>5.64</td><td>m</td></tr><tr><td>Operational rotor speed Ω<sub>0</sub></td><td>31.42</td><td>rads<sup>-1</sup></td></tr><tr><td>Blade mass m<sub>b</sub></td><td>94.9</td><td>kg</td></tr><tr><td>Blade mass moment S<sub>b</sub></td><td>281.1</td><td>kgm</td></tr><tr><td>Blade mass moment of inertia I<sub>b</sub></td><td>1084.7</td><td>kgm<sup>2</sup></td></tr><tr><td>Lag hinge offset e</td><td>0.3048</td><td>m</td></tr><tr><td>Lag spring k<sub>b</sub></td><td>0</td><td>Nmrad<sup>-1</sup></td></tr><tr><td>Longitudinal hub mass M<sub>x</sub></td><td>8026.6</td><td>kg</td></tr><tr><td>Lateral hub mass M<sub>y</sub></td><td>3283.6</td><td>kg</td></tr><tr><td>Number of blades N<sub>b</sub></td><td>4</td><td>-</td></tr><tr><td>Longitudinal hub spring k<sub>x</sub></td><td>1240481.08</td><td>N/m</td></tr><tr><td>Lateral hub spring k<sub>y</sub></td><td>1240481.08</td><td>N/m</td></tr></tbody></table>

## 4.1. Rigid model

A rigid model has already been developed by former student Sviatoslav Tezikov [3]. In his work, a rigid model is studied which does not take into account the deformation of the helicopter blades and therefore their potential impact on the stability of the system under study. The equations governing the motions of such a system have been derived by Hammond [2] and read:

$$ (m_x + N m_b) \ddot{x}_h + c_x \dot{x}_h + k_x x_h = S_b \sum_{i=1}^{N} \left[ (\ddot{\zeta}_i - \Omega^2 \zeta_i) \sin \psi_i + 2\Omega \dot{\zeta}_i \cos \psi_i \right]; $$

$$ (m_y + N m_b) \ddot{y}_h + c_y \dot{y}_h + k_y y_h = -S_b \sum_{i=1}^{N} \left[ (\ddot{\zeta}_i - \Omega^2 \zeta_i) \cos \psi_i - 2\Omega \dot{\zeta}_i \sin \psi_i \right]. $$

This system of equations consist of N+2 differential equations which are coupled. The periodic coefficient arise because the blades equations are developed in a rotating frame of reference whereas the hub equations are developed in a fixed frame of reference. An alternative consists of writing the hub equations of motion in the rotating frame of reference. The equations thus read:

$$ \ddot{x} + \eta_h \dot{x} + (\omega_h^2 - \Omega^2) \bar{x} - 2\Omega \dot{\bar{y}} - \Omega \eta_h \bar{y} = \nu_h^2 \sum_{j=1}^{N} \left[ (\zeta_j - \Omega^2 \zeta_j) \sin \frac{2\pi}{N}(j-1) + 2\Omega \zeta_j \cos \frac{2\pi}{N}(j-1) \right]; $$

$$ \ddot{y} + \eta_h \dot{y} + (\omega_h^2 - \Omega^2) \bar{y} + 2\Omega \dot{\bar{x}} + \Omega \eta_h \bar{x} = -\nu_h^2 \sum_{j=1}^{N} \left[ (\zeta_j - \Omega^2 \zeta_j) \cos \frac{2\pi}{N}(j-1) - 2\Omega \zeta_j \sin \frac{2\pi}{N}(j-1) \right], $$

where the bar notation refers to the global frame of reference and where the following parameters have been introduced

$$ \nu_h^2 = \frac{S_b}{m_x + N m_b}; $$

$$ \omega_h^2 = \frac{k_x}{m_x + N m_b}; $$

$$ \eta_h = \frac{c_x}{m_x + N m_b}. $$
## 4.2. Validation of the different numerical scheme

To validate the various numerical scheme, the various presented numerical schemes are tested and compared to results available in the literature. The different results obtained with the various numerical schemes are highlighted in Figure 4, Figure 5, Figure 6 and Figure 7.

<figure><img src="image_6.png" /><figcaption>**Figure 4:** Relative modal damping ratio obtained using Newmark integration scheme using rigid model.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5:** Relative modal damping ratio obtained using Runge-Kutta integration scheme using rigid model.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6:** Relative modal damping ratio obtained using Newmark V2 integration scheme using rigid model.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7:** Relative modal damping ratio obtained using α generalised integration scheme using rigid model.</figcaption></figure>

The different numerical schemes employed provided results in agreement with results presented by Ignacio et al in [4]. The different results are obtained neglecting the external damping on the hub $c_x$ and $c_y$ and the blade damping $c_b$ but the coupling term corresponding to the gyroscopic term are taken into account. However, one can observe that the results obtained using the Runger-kutta integration scheme are quite different from the ones obtained using the Newmark integration scheme. In fact, the zone of instability which corresponds to a positive modal damping is larger than the zone of instability obtained using the Newmark integration scheme. In addition, Newmark's integration scheme detects two zones of instability separated each separated by a stability zone. These results seem to be more in accordance with results obtained in the literature. The difference in the instability zone is probably coming from the particularity of the different integration schemes. In fact, the Runge-Kutta integration scheme employs a variable time step to limit the relative error between two integration steps whereas the Newmark integration scheme employs a constant time step. The results obtained using alpha generalised numerical scheme are a bit higher than the one obtained using Newmark numerical scheme but this is due to the fact that the spectral radius is taken into in the computation of the value of the parameter $\beta$ and $\gamma$ of the scheme, which makes them different from the one of the Newmark scheme.

Although Floquet theory and the analysis of Floquet exponents have enabled us to conclude whether the
rigid model under study is stable or not, it is interesting to study the behavior of the dynamic system under study at a speed of rotation concluded by Floquet theory to be stable or unstable. To do so, The Newmark integration scheme is employed using as an initial condition an initial displacement, over one blade, of 0.001 m and no initial velocity. Figure 8 and Figure 10 each represent the movement in the x and y degrees of freedom plane of the hub. It can be seen that when the rigid system is rotated at a speed of 3 rad/s, which is a stable speed according to Floquet's theory, the movement remains bounded around its equilibrium position, which is confirmed by analyzing the movement over ten periods of revolutions (Figure 9), whereas when this system is rotated at a speed of rotation said to be unstable by Floquet's theory, it can be seen that the hub's movement is not bounded and could lead to its destruction if operated at this speed of rotation. This critical behavior is visible when looking at Figure 11, when the movement of the hub is studied over ten periods of revolution of the propeller.

<figure><img src="image_10.png" /><figcaption>**Figure 8:** Hub displacement over one period of revolution when rotor rotating at $\Omega = 3$ rad/s.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9:** Hub displacement over ten periods of revolution when rotor rotating at $\Omega = 3$ rad/s.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10:** Hub displacement over one period of revolution when rotor rotating at $\Omega = 30$ rad/s.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 11:** Hub displacement over ten periods of revolution when rotor rotating at $\Omega = 30$ rad/s.</figcaption></figure>

# 5. CONCLUSIONS

This study introduces an advanced computational framework for whirl flutter analysis, leveraging FEM and Floquet theory. The results demonstrate improved accuracy in stability prediction, particularly through the use of time-dependent coupling. Key conclusions include:

* Time-dependent coupling enhances accuracy in whirl flutter analysis.

* Floquet theory effectively identifies instability regions in rotor-wing systems.
* Runge-Kutta integration provides superior results but at a higher computational cost.

Future work will focus on optimizing computational efficiency and integrating aerodynamic forces for more comprehensive analysis. Experimental validation will also be pursued to further confirm the effectiveness of the proposed approach.

## REFERENCES

[1] Arne Vollan and Louis Komzsik. *Computational Techniques of Rotor Dynamics with the Finite Element Method*. 2017. ISBN 13:978-1-4398-4770-1.

[2] C E Hammond, Langley Directorate, U S Armv, Air Moblltv, and R8d Laboratorv. An application of floquet theory to prediction of mechanical Instability. *Journal of the American Helicopter Society*, 19:14–23, 1973.

[3] Sviatoslav Tezikov. Computational stability analysis of whirl flutter with application to civil urban mobility. 2023.

[4] José A. Ignácio da Silva, Douglas D. Bueno, and Gustavo L.C.M. de Abreu. On the controllers’ design to stabilize ground resonance helicopter. *JVC/Journal of Vibration and Control*, 25:2894–2909, 12 2019. ISSN 17412986. doi: 10.1177/1077546319873797.