# Reduced-order modeling of geometrically nonlinear rotating structures using the direct parametrisation of invariant manifolds

Adrien Martin, Andrea Opreni, Alessandra Vizzaccaro, Marielle Debeurre, Loic Salles, Attilio Frangi, Olivier Thomas, Cyril Touzé

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

## ► To cite this version:

Adrien Martin, Andrea Opreni, Alessandra Vizzaccaro, Marielle Debeurre, Loic Salles, et al. Reduced-order modeling of geometrically nonlinear rotating structures using the direct parametrisation of invariant manifolds. Journal of Theoretical, Computational and Applied Mechanics, 2023, pp.1-31. 10.46298/jtcam.10430. hal-03886793v2

HAL Id: hal-03886793

https://hal.science/hal-03886793v2

Submitted on 25 May 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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
**Identifiers**

DOI 10.46298/jtcam.10430

HAL hal-03886793v2

**History**

Received Dec 6, 2022

Accepted Apr 11, 2023

Published May 25, 2023

**Associate Editor**

Alexander POPP

**Reviewers**

Anonymous

Anonymous

Olivier BAUCHAU

**Open Review**

OAI hal-04087439

**Supplementary Material**

Software

DOI 10.5281/zenodo.7924472

**Licence**

CC BY 4.0

©The Authors

# Reduced-order modeling of geometrically nonlinear rotating structures using the direct parametrisation of invariant manifolds

Adrien MARTIN¹, Andrea OPRENI², Alessandra VIZZACCARO³, Marielle DEBEURRE⁴,
Loïc SALLES⁵,⁶, Attilio FRANGI², Olivier THOMAS⁴, and Cyril TouZÉ¹

¹ IMSIA, CNRS, EDF, CEA, ENSTA Paris, Institut Polytechnique de Paris

² Department of Civil and Environmental Engineering, Politecnico di Milano

³ Department of Engineering Mathematics, University of Bristol

⁴ Arts et Métiers Institute of Technology, LISPEN, HESAM Université

⁵ Department of Aerospace & Mechanical Engineering, University of Liège

⁶ Vibration University Technology Centre, Imperial College London

The direct parametrisation method for invariant manifolds is a nonlinear reduction technique which derives nonlinear mappings and reduced-order dynamics that describe the evolution of dynamical systems along a low-dimensional invariant-based span of the phase space. It can be directly applied to finite element problems. When the development is performed using an arbitrary order asymptotic expansion, it provides an efficient reduced-order modeling strategy for geometrically nonlinear structures. It is here applied to the case of rotating structures featuring centrifugal effect. A rotating cantilever beam with large amplitude vibrations is first selected in order to highlight the main features of the method. Numerical results show that the method provides accurate reduced-order models (ROMs), thus offering remarkable reduction in the computational burden. The hardening/softening transition of the fundamental flexural mode with increasing rotation speed is then investigated in detail and a ROM parametrised with respect to rotation speed and forcing frequencies is detailed. The method is then applied to a twisted plate model representative of a fan blade, showing how the technique can handle more complex structures. Hardening/softening transition is also investigated as well as interpolation of ROMs, highlighting the efficacy of the method.

**Keywords:** rotating structure, invariant manifold, geometric nonlinearity, nonlinear normal modes

## 1 Introduction

Rotating structures are key components in numerous industrial applications ranging from turbo-machines (Rao 1991; Yoo et al. 2001; Legrand et al. 2009; Sinha 2013), helicopter blades (DaSilva and Hodges 1986; Bauchau et al. 2001) to wind turbine (Chaviaropoulos 1996), see also the work of Friswell et al. 2010 for an overview of applications. In recent years, new designs (for example the open rotor concept (Van Zante 2015) in aeronautics or the increasing size of wind turbine) leads to consider even longer blades where large displacements are more easily reached (Hodges and Dowell 1974; Manolas et al. 2015; Rezaei et al. 2018; Delhez et al. 2021; Di Palma et al. 2022), hence exciting geometric nonlinearities. Numerical computation of full-order models with the finite element (FE) method for such problems leads to prohibitive computational times that are generally not compatible with the constraints at the design stage, or asks for specific developments using e.g. parallel computing (Blahoš et al. 2020; Blahoš 2022). Consequently, the development of efficient and accurate reduced-order models (ROM) that are able to tackle the geometric nonlinearity is a key feature.

Reduction methods have a long history, see e.g. (Mignolet et al. 2013; Touzé et al. 2021) for two different review articles on the subject. Numerous methods have been derived within the
framework of Nonlinear Normal Modes (NNMs), first introduced by Rosenberg 1962. In this setting, the definition of an NNM as an invariant manifold (Shaw and Pierre 1991; Shaw and Pierre 1993) allows deriving efficient reduction strategies using either the graph style through the use of the center manifold technique (Pescheck et al. 2002; Jiang et al. 2005), or the normal form approach followed by a truncation on the normal coordinates (Touzé et al. 2004; Touzé and Amabili 2006; Touzé 2014). On the other hand, reduction methods directly applicable to FE structures have also been proposed, using either linear or nonlinear methods, see e.g. (Idelsohn and Cardona 1985; Muravyov and Rizzi 2003; Hollkamp and Gordon 2008).

In recent years, two important steps forward have been realized. On the mathematical point of view, the parametrisation method for invariant manifolds, as first introduced by Cabré et al. 2005 and overviewed with an emphasis on computational applications in the book by Haro et al. 2016, is the key link that allows expressing the different styles of solutions for invariant manifold computations in a unified manner. Rewriting the problem from the invariance equation, it is shown that two main different parametrisations are meaningful: the graph style and the normal form style. Returning to the applications in vibration theory, one understands that Shaw and Pierre method used the graph style, whereas the normal form approach derived in (Touzé et al. 2004; Touzé and Amabili 2006) was relying on the normal form style. Using solutions from the modal space, the parametrisation method has been first applied in vibration theory by Haller and Ponsioen (2016); Ponsioen et al. (2018) who proposed an automated version allowing to overtake the limitations of earlier developments in the field of vibratory systems, which were restricted to low-order asymptotics.

The second important step has been concerned with direct applications of the methods based on invariant manifold theory to finite FE problems. The direct normal form (DNF) has been proposed in (Vizzaccaro et al. 2021b), bypassing the need of expressing the dynamics in the modal space as a first step, which was a needed computation to apply the reduction method proposed in (Touzé et al. 2004; Touzé and Amabili 2006). The main consequence is that one can use nonlinear mappings to go directly from the physical space to an invariant-based span of the phase space. The DNF has been applied to Micro Electro Mechanical Systems (MEMS) featuring internal resonance (Opreni et al. 2021), and to problems with cyclic symmetry in (Quaegebeur et al. 2023). Interestingly, the DNF as developed in (Vizzaccaro et al. 2021b) can be used in a non-intrusive manner. The parametrisation method has also been applied to FE structures in (Jain and Haller 2022) and in (Vizzaccaro et al. 2022; Opreni et al. 2022c), with numerous applications ranging from wing dynamics (Jain and Haller 2022) to MEMS (Vizzaccaro et al. 2022; Opreni et al. 2022c; Opreni et al. 2022b; Di Cristofaro et al. 2022). In particular, the terms due to the presence of an external forcing are fully taken into account in (Opreni et al. 2022c), enlarging the field of applications and showing the robustness of the method in delivering accurate predictions up to large amplitudes with a minimal representation of the ROM. This direct version of the parametrisation method for invariant manifolds, directly applicable to FE structures discretized with 3D elements, has also been proposed as a public open code, MORFE (Opreni et al. 2022a). Elaborating on direct application of the parametrisation method, SSMtool also provides an open public code with beam and plate elements (Jain et al. 2023).

A number of reduction methods for nonlinear dynamical problems featuring geometric nonlinearity have been compared together in a series of papers, see e.g. (Haller and Ponsioen 2017; Vizzaccaro et al. 2021a; Wang and Mignolet 2021; Shen et al. 2021a; Shen et al. 2021b; Shen et al. 2021c; Touzé et al. 2021). More specifically, ad-hoc reduction methods are compared against techniques based on invariant manifold computations, since the latter are known to give the most accurate results from strong mathematical theorems ensuring that the long-term dynamics lies in these reduced subspace. Interestingly, it has been clearly demonstrated that methods, such as the implicit condensation and expansion (ICE) (Hollkamp and Gordon 2008; Frangi and Gobat 2019), or the quadratic manifold with modal derivatives (Jain et al. 2017), are able to deliver accurate results only if a slow-fast assumption between master and slave modes, is at hand.

The aim of this paper is to apply the direct parametrisation for invariant manifold (DPIM), as developed in (Vizzaccaro et al. 2021b; Opreni et al. 2022c) and released in the software MORFE (Opreni et al. 2022a), to rotating structures. The literature on rotating structures shows numerous applications of reduction methods. A special emphasis is generally put on the case of
contact and friction, see e.g. (Laxalde and Thouverez 2009; Petrov 2011; Khalifeh et al. 2017; Joannin et al. 2018; Thorin et al. 2018; Guérin et al. 2018), since this problem is the most stringent one in design, vibration and fatigue analysis. As a consequence, most of the reduction techniques used to tackle such localized nonlinearity fall within the family of substructuring methods, including number of variants ranging from the original Guyan and Craig-Bampton methods to component mode synthesis (Guyan 1965; Craig and Bampton 1968; MacNeal 1971; de Klerk et al. 2008; Allen et al. 2020). On the other hand, geometric nonlinearity is distributed and has been less frequently tackled in the literature concerned with rotating structures. While it is taken into account in (Thomas et al. 2016) with a reduction perspective, both geometric and contact nonlinearities are considered for example in (Delhez et al. 2021). In this last example, the structural geometric nonlinearity is handled with modal derivatives, known to have limitations in the case where the slow/fast assumption is not met (Vizzaccaro et al. 2021a).

The direct parametrisation method for invariant manifold (DPIM) can be adapted to handle the case of rotating structures with geometric nonlinearity. The added terms due to the rotation effects can be separated in two: the centrifugal and the Coriolis effects, appearing through different terms in the equations of motion (Rao 1991; Friswell et al. 2010). Also, depending on the configuration, one of these two terms might not be excited, or one can be dominant on the other one. In order to study them separately, this paper assumes that only the centrifugal effects are present, while the Coriolis force can be safely ignored. The DPIM will be adapted to handle this specific case and the added terms, which also leads to consider a new static position. Two different examples will then be thoroughly analyzed: a rotating beam, and a twisted plate, the geometry of which has been selected in order to resemble a fan blade.

# 2 Problem formulation and reduction method

## 2.1 Equations of motion

A structure rotating with angular speed $\Omega$ and experiencing large amplitude vibrations is considered and described in the rotating frame. The rotation brings about two different effects that can be separately analyzed: the centrifugal effect and the Coriolis force (Rao 1991). The expressions of these forces show that the Coriolis effect can be safely neglected in numerous situations, in particular for structures rotating along an axis that is orthogonal to the main vibrations, which is the case considered here. In order to give a more quantitative assessment, Appendix A details the expression of the Coriolis term and highlights its negligible effect in the numerical cases that are considered in the present study. As a consequence, only the centrifugal force will be considered in the main text. With such assumptions, the general equations of motion can be found for example in (Rao 1991; Géradin and Rixen 2015; Friswell et al. 2010). Starting from the semi-discrete problem, spatially discretised for example with the finite element procedure, the equations of motion write for the unknown displacement $\mathbf{p}(t, \Omega)$, which is an $N$-dimensional vector representing the degrees of freedom, as

$$ \mathbf{M}\ddot{\mathbf{p}}(t, \Omega) + \mathbf{C}\dot{\mathbf{p}}(t, \Omega) + [\mathbf{K} + \mathbf{N}(\Omega)]\mathbf{p}(t, \Omega) + \mathbf{g}(\mathbf{p}(t, \Omega), \mathbf{p}(t, \Omega)) + \mathbf{h}(\mathbf{p}(t, \Omega), \mathbf{p}(t, \Omega), \mathbf{p}(t, \Omega)) = \mathbf{f}^{\text{ext}}(t) + \mathbf{f}^{\Omega}(\Omega). \quad (1) $$

In this equation, $\mathbf{M}$, $\mathbf{C}$ and $\mathbf{K}$ stand respectively for the mass, damping and stiffness matrices, symmetric of size $N \times N$; $\mathbf{N}(\Omega)$ and $\mathbf{f}^{\Omega}(\Omega)$ are respectively the spin softening matrix and the centrifugal load vector, both due to the centrifugal effect brought by rotation; $\mathbf{f}^{\text{ext}}$ is the external load. The matrix $\mathbf{N}(\Omega)$ is symmetric negative semi-definite, $\mathbf{K}$ is positive semi-definite and $\mathbf{M}$ is positive definite. The geometric nonlinearity is distributed and expressed through the nonlinear quadratic and cubic terms $\mathbf{g}$ and $\mathbf{h}$. In the framework of elastic bodies with large-amplitude vibrations, having only quadratic and cubic nonlinearity stems from the equations of motion and is not an approximation, see e.g. (Holzapfel 2000; Touzé et al. 2021).

The space discretisation is assumed to be dealt with a matrix collecting the expansion space functions $\mathbf{H}$. In this paper, the FE method is more specifically addressed such that $\mathbf{H}$ refers to the matrix of interpolation functions associated to the FE procedure (Bathe 1996). The two terms
representing the centrifugal effect can then be explicitly rewritten as (Rao 1991; Géradin and Rixen 2015)

$$ \mathbf{N} = \int_V \rho \mathbf{H}^\mathrm{T} \Omega^2 \mathbf{H} \, dV, \qquad (2a) $$

$$ \mathbf{f}^{\Omega} = \int_V \rho \mathbf{H}^{\mathrm{T}} \Omega^2 \mathbf{y} \, dV \qquad (2b) $$

where $\Omega$ is the skew-symmetric second-order tensor of rotation, while $\mathbf{y}$ represents the position vector of an arbitrary point of the structure in the rotating frame. The geometric nonlinearity arises from the nonlinear relationship between the strain and the displacements, and gives rise to only quadratic and cubic nonlinear terms if using 3D finite elements implementing a full Lagrangian formulation. The nonlinear terms $\mathbf{g}$ and $\mathbf{h}$ can then be expressed in the general form

$$ \mathbf{g}(\mathbf{p}, \mathbf{p}) = \sum_{i,j=1}^{N} \mathbf{g}_{ij} p_i p_j \quad \text{and} \quad \mathbf{h}(\mathbf{p}, \mathbf{p}, \mathbf{p}) = \sum_{i,j,k=1}^{N} \mathbf{h}_{ijk} p_i p_j p_k \qquad (3) $$

with $\mathbf{g}_{ij}$ and $\mathbf{h}_{ijk}$, vectors of nonlinear coupling coefficients.

Let us now further assume that the rotation speed $\Omega$ is constant. In such case the centrifugal effect leads to consider a new static equilibrium position, and vibrations around this new state. The displacement vector can thus be rewritten as

$$ \mathbf{p}(t, \Omega) = \mathbf{u}_0(\Omega) + \mathbf{u}(t) \qquad (4) $$

where $\mathbf{u}_0$ represents the static equilibrium position, while $\mathbf{u}$ is the new vector of dynamic degrees of freedom (dofs). The deformed configuration with constant rotation speed is found by solving

$$ [\mathbf{K} + \mathbf{N}(\Omega)]\mathbf{u}_0(\Omega) + \mathbf{g}(\mathbf{u}_0(\Omega), \mathbf{u}_0(\Omega)) + \mathbf{h}(\mathbf{u}_0(\Omega), \mathbf{u}_0(\Omega), \mathbf{u}_0(\Omega)) = \mathbf{f}^{\Omega}(\Omega). \qquad (5) $$

Numerically, this nonlinear static problem can be solved using a Newton-Raphson procedure. On the other hand, plugging Equation (4) into Equation (1), and using Equation (5), one obtains that the nonlinear vibrations around the deformed configuration are governed by

$$ \mathbf{M}\ddot{\mathbf{u}}(t) + \mathbf{C}\dot{\mathbf{u}}(t) + \mathbf{K}_t(\Omega) \mathbf{u}(t) + \mathbf{g}(\mathbf{u}(t), \mathbf{u}(t)) + 3\mathbf{h}(\mathbf{u}_0(\Omega), \mathbf{u}(t), \mathbf{u}(t)) \\ + \mathbf{h}(\mathbf{u}(t), \mathbf{u}(t), \mathbf{u}(t)) = \mathbf{f}^{\mathrm{ext}}(t). \qquad (6) $$

In this equation, the nonlinear terms are expanded with the modified tangent stiffness matrix

$$ \mathbf{K}_t(\Omega) = \mathbf{K} + \mathbf{N}(\Omega) + 2\mathbf{g}(\mathbf{u}_0(\Omega), \mathbf{I}) + 3\mathbf{h}(\mathbf{u}_0(\Omega), \mathbf{u}_0(\Omega), \mathbf{I}). \qquad (7) $$

The nonlinear dynamics of the structure rotating with constant speed and excited through an external force, is thus the solution of Equation (6). Direct computations are difficult because of the size of the problem (N the number of dofs), which can be prohibitively large for refined meshes. The aim of the present contribution is to use NNMs as invariant manifolds to drastically reduce the size of the problem to be solved; and more specifically to adapt the DPIM introduced in (Vizzaccaro et al. 2022; Opreni et al. 2022c) and make it available to solve rotational problems with centrifugal effect. As a matter of fact, one can observe that Equation (6) is written under the form used in (Vizzaccaro et al. 2021b; Opreni et al. 2022c). One has just to take care of: (i) computing the rotating position without vibrations using Equation (5), (ii) modifying the tangent stiffness matrix using Equation (7), and (iii) updating the quadratic nonlinear term to take into account the quadratic part coming from cubic nonlinearity in Equation (6). Once that realized, the DPIM method presented in (Vizzaccaro et al. 2021b) for the autonomous part, and in (Opreni et al. 2022c) for the nonautonomous part, can be directly applied to the Equation (6).
## 2.2 Direct Parametrisation of Invariant Manifold

The direct parametrisation of invariant manifolds, applied to nonlinear vibratory systems, allows one to compute the NNMs and the reduced dynamics on an invariant manifold attached to a hyperbolic fixed point of the system, operating directly from the FE dofs, i.e. without the need of first computing the dynamics in the modal basis as a prerequisite. The method is able to handle both damping and forcing in order to propose efficient ROMs for dynamical solutions. The NNMs are tangent to their linear counterparts for small amplitudes, which means that the method use the modal information, and allows to handle larger amplitude by finely capturing the curvatures of the invariant manifold where the dynamics lies. Before recalling the main steps of the derivation from a user's perspective, the linear characteristics needed to start the calculation are briefly summarized, together with the assumptions used about damping and forcing. For a more thorough presentation of the method, the interested reader is referred to (Vizzaccaro et al. 2022; Opreni et al. 2022c).

### 2.2.1 Linear characteristics and forcing

The eigenmodes of the underlying linear conservative autonomous system of Equation (6), are defined as the solutions of

$$[\mathbf{K}_t - \omega_j^2 \mathbf{M}] \boldsymbol{\phi}_j = \mathbf{0}. \qquad (8)$$

The linear normal modes are assumed to be normalized with respect to the mass matrix, leading to the classical relationships

$$\boldsymbol{\phi}_i^\top \mathbf{M} \boldsymbol{\phi}_j = \delta_{ij} \quad \text{and} \quad \boldsymbol{\phi}_i^\top \mathbf{K}_t \boldsymbol{\phi}_j = \delta_{ij} \omega_i^2. \qquad (9)$$

In computational mechanics, it is customary to introduce Rayleigh damping in order to take losses into account (Géradin and Rixen 2015). In this contribution, since the tangent stiffness matrix depends on the rotation speed, see Equation (7), it appears more convenient to select only mass-proportional damping. Note that this is not a limitation of the reduction method which can also handle stiffness-proportional damping (Vizzaccaro et al. 2021b; Vizzaccaro et al. 2022). The damping matrix is thus selected here as

$$\mathbf{C} = \beta \mathbf{M} \qquad (10)$$

where $\beta$ is a constant coefficient, whose value will be tuned to control a given mode with a specific ratio. Indeed one has, for a selected mode labelled $i$: $\beta = 2\zeta_i\omega_i$, where $\zeta_i$ is the nondimensional damping ratio associated to the $i^{th}$ mode. Consequently, in the numerical results shown in the next sections, $\zeta_i$ will be selected for a given mode $i$ (generally the one that is driven by the external forcing), typical values in the range 0.01 to 0.2% will be selected, which in turn will fix the value of $\beta$.

For the damped system, the eigenvalues have now real and imaginary parts and are denoted $\lambda_j = -\zeta_j\omega_j + i\omega_j\sqrt{1 - \zeta_j^2}$ and $\lambda_{j+N} = \bar{\lambda}_j, \forall j = 1, \dots, N$. Note that these correspond to the eigenvalues of the problem written at first-order, leading to pairs of complex conjugate eigenvalues (Géradin and Rixen 2015). They are needed in order to express the reduced dynamics on the invariant manifold. In the rest of the paper, it is also assumed that the damping ratio is small. The important consequence is that the resonance relationships, a key feature of the reduction method (Touzé and Amabili 2006; Touzé 2014), are checked only on the imaginary part of the eigenvalues, in line with the small damping assumption.

In the simulations presented hereafter, the case of an harmonic external forcing with a single excitation frequency, is considered, such that the external load in the equation of motion Equation (6) can be expressed as

$$\hat{\mathbf{f}}_{\text{ext}}(t) = \mathbf{f}_{\text{ext}} \frac{e^{+i\hat{\omega}t} + e^{-i\hat{\omega}t}}{2}, \qquad (11)$$

with $\mathbf{f}_{\text{ext}}$, a vector representing the discretised space dependence of the forcing and $\hat{\omega}$ the excitation frequency. In most of the numerical examples, the forcing frequency will be selected in
the vicinity of an eigenfrequency $\omega_j$, in order to compute the frequency-response curve (FRC) of the system, which collects the maximum amplitude of the response as function of the excitation frequency. In such case, it has also been assumed for simplicity that the spatial dependence of the forcing follows that of the $j^1$ eigenmode shape directly excited. Consequently the vector $\mathbf{f}^{\text{ext}}$ can be expressed in this case as

$$ \mathbf{f}^{\text{ext}} = \alpha \mathbf{M} \phi_j, \qquad (12) $$

with $\alpha$, the amplitude of the forcing. Note that this assumption is not a limitation of the method, and numerical results with a pointwise forcing are reported in Appendix C.

### 2.2.2 Nonlinear mappings and reduced dynamics

The DPIM relies on the introduction of two nonlinear mappings, respectively between the displacement field $\mathbf{u}$ and its associated velocity field $\mathbf{v} = \dot{\mathbf{u}}$, in relation with a reduced coordinate describing the dynamics on the associated invariant manifold. This reduced coordinate $\mathbf{z}$ is denoted as the normal coordinate. It has the size of $2n$, with $n$ the number of master modes selected to derive the ROM. In practical applications, $n$ is assumed to be very small. For nonlinear dynamics with a unimodal Duffing-like behaviour, $n = 1$. On the other hand, when more complex dynamics arise with the excitation of internal resonance and, for instance, the appearance of quasi-periodic solutions, then $n$ needs to be increased in order to take into account all the internally resonant modes.

The two nonlinear mappings are expressed by first separating the contributions coming from the autonomous system (without forcing) to that of the non-autonomous one. This method is used to compute two separate expansions and treat individually the resonance conditions appearing in each contribution. Consequently, the two mappings read

$$ \mathbf{u} = \psi(\mathbf{z}) + \hat{\psi}(\mathbf{z}, \hat{\omega}, t) \qquad (13a) $$

$$ \mathbf{v} = v(\mathbf{z}) + \hat{v}(\mathbf{z}, \hat{\omega}, t) \qquad (13b) $$

where $\psi$ and $v$ refers to the autonomous part of the expansion, while $\hat{\psi}$ and $\hat{v}$ to the non-autonomous part. For an unforced problem, only the autonomous terms are to be taken into account and the searched invariant manifold are independent of time. On the other hand, taking the forcing into account leads to a time dependence of the invariant manifolds, which can show small deformations during a period of the forcing in addition to the main oscillatory forced motion. To take this property into account, $\hat{\psi}$ and $\hat{v}$ are written with an explicit dependence upon time $t$ and excitation frequency $\hat{\omega}$.

The reduced dynamics, which is an unknown at this stage of the development, is also expanded following the same separation, and is thus assumed to read

$$ \dot{\mathbf{z}} = \mathbf{f}(\mathbf{z}) + \hat{\mathbf{f}}(\mathbf{z}, \hat{\omega}, t) \qquad (14) $$

where $\mathbf{f}(\mathbf{z})$ refers to the autonomous part (without forcing), while $\hat{\mathbf{f}}(\mathbf{z}, \hat{\omega}, t)$ represents the non-autonomous part.

Each of the six terms of the previous expansions are expressed as high-order polynomials of $\mathbf{z}$. Using Equations (13) and (14), and inserting them in the equation of motion rewritten in first-order form, one can eliminate time to obtain the so-called *invariance equation* (Haro et al. 2016; Touzé et al. 2021; Vizzaccaro et al. 2022; Opreni et al. 2022c). Collecting terms of like powers of $\mathbf{z}$, one is then able to write the homological equations, for each order of the expansion. The linear solution is easy to derive and shows that the linear part of the mappings is colinear with the master eigenmodes, while the linear terms of the reduced dynamics retrieve the usual linear modal oscillator equations. In order to simplify the calculations associated with the non-autonomous part, it is assumed that the dependence as function of the excitation frequency $\hat{\omega}$ is small. Consequently, the mappings are computed for a single excitation frequency, instead of recomputing the whole procedure for all excitation frequency, as first proposed for example in (Jiang et al. 2005). This assumption is important to simplify the computational burden and gives access to efficient ROMs able to compute FRC (Opreni et al. 2022c). Finally, the time dependence
of the non-autonomous terms, see Equations (13) and (14), disappears when considering harmonic forcing, since the invariance equation for the non-autonomous part is written for each harmonic of the forcing term (Opreni et al. 2022c). More specifically, one thus rewrite the non-autonomous part of the mappings expressed in Equation (13) as

$$ \hat{\psi}(\mathbf{z}, \hat{\omega}, t) = \hat{\psi}_{+}(\mathbf{z}) e^{i\hat{\omega}t} + \hat{\psi}_{-}(\mathbf{z}) e^{-i\hat{\omega}t}, \quad \hat{v}(\mathbf{z}, \hat{\omega}, t) = \hat{v}_{+}(\mathbf{z}) e^{i\hat{\omega}t} + \hat{v}_{-}(\mathbf{z}) e^{-i\hat{\omega}t}. \quad (15) $$

Since a single excitation frequency is used in this contribution, only two invariance equations, related to the associated spectral component $\pm i\hat{\omega}$, need to be solved for the non-autonomous part. With all these assumptions and results, one can write the general solution for the nonlinear mapping on the displacement as:

$$ \psi(\mathbf{z}) = \Phi\mathbf{z} + \sum_{p=2}^{o_a} [\psi(\mathbf{z})]_p, \quad \hat{\psi}(\mathbf{z}, \hat{\omega}, t) = \sum_{p=0}^{o_{na}} [\hat{\psi}(\mathbf{z}, \hat{\omega}, t)]_p, \quad (16) $$

where $\Phi$ is the matrix gathering the master modes only, which has thus a dimension $N \times 2n$, with $n$ the number of selected master modes, sorted as

$$ \Phi = [\phi_1 \ \phi_2 \ \cdots \ \phi_n \ \phi_1 \ \phi_2 \ \cdots \ \phi_n]. \quad (17) $$

In these expressions, the notation $[\bullet]_p$ refers to a generic polynomial expression of order $p$; $o_a$ and $o_{na}$ to the maximal order of expansion retained respectively for the autonomous and the non-autonomous part. In the rest of the paper, the notation $O(o_a, o_{na})$ will be systematically used to refer to the selected orders in the numerical examples. One can observe that the displacement part of the nonlinear mapping is tangent to the master eigenmodes, the nonlinear terms starting from order 2. On the other hand, the non-autonomous development starts at order zero. As discussed in (Opreni et al. 2022c), the zero order for the non-autonomous part is consistent with earlier approximations in reduced-order modeling strategies, see e.g. (Touzé and Amabili 2006), where the modal forcing is directly appended to the right-hand side of the reduced dynamics. Higher orders allows approximating the deformations of the time-dependent invariant manifold (also referred to as the *whisker*), resulting in more accurate results for larger amplitudes, as well as taking into account different forcing conditions including parametric excitation (Opreni et al. 2022c; Ponsioen et al. 2019).

The nonlinear mappings for the velocity fields are also expanded in a similar manner as

$$ v(\mathbf{z}) = \Phi\Lambda\mathbf{z} + \sum_{p=2}^{o_a} [v(\mathbf{z})]_p \quad \text{and} \quad \hat{v}(\mathbf{z}, \hat{\omega}, t) = \sum_{p=0}^{o_{na}} [\hat{v}(\mathbf{z}, \hat{\omega}, t)]_p \quad (18) $$

where $\Lambda$ is the matrix gathering the complex conjugate eigenvalues of the sole master modes

$$ \Lambda = \mathrm{diag}[\lambda_1 \ \cdots \ \lambda_n \ \bar{\lambda}_1 \ \cdots \ \bar{\lambda}_n]. \quad (19) $$

Finally, the reduced dynamics is also expanded following

$$ \dot{\mathbf{z}} = \Lambda\mathbf{z} + \sum_{p=2}^{o_a} [\mathbf{f}(\mathbf{z})]_p + \sum_{p=0}^{o_{na}} [\hat{\mathbf{f}}(\mathbf{z}, \hat{\omega}, t)]_p. \quad (20) $$

One can note here that the linear part of the reduced dynamics follows the usual linear modal equation. Nonlinear terms account for the nonlinear oscillatory behaviour expressing the dynamics on the invariant manifold.

### 2.2.3 Styles of parametrisation

The linear solutions for the unknown mappings and reduced dynamics have been extracted from the general solutions in the previous section, in order to highlight the important property of identity-tangency, underlining that the NNMs as invariant manifold are the natural continuation of the linear modes to the nonlinear regime.
The solutions of the homological equations of arbitrary order *k*, starting from the quadratic terms with *k* = 2, is made difficult due to various reasons. First, the homological equations are under-determined, such that the solution is not unique. Indeed, as shown theoretically for example in (Haro et al. 2016), there are infinitely many solutions. Each of the possible solutions is a specific style of parametrisation. Nevertheless, among all these solutions, two very specific styles appear more naturally in the developments.

The first solution is denoted as the *graph style*, since it assumes an identity relationship between the master coordinate and the normal coordinate, meaning that the slave modal coordinates are expressed as functions of the master ones. In vibration theory, this style of solution is the one used by Shaw and Pierre in their first developments on NNMs (Shaw and Pierre 1991; Shaw and Pierre 1993; Touzé et al. 2021). One limitation of this graph style is its inability capture folding points of the invariant manifold, see e.g. (Haller and Ponsioen 2016; Vizzaccaro et al. 2022) for illustration of this fact.

The second style is the *normal form style* already used in (Touzé et al. 2004; Touzé and Amabili 2006). In this case, there is a nonlinear relationship between the master modal coordinate and the normal coordinate, allowing the parametrisation to go over folding points. As also explained in (Vizzaccaro et al. 2022), different normal form styles can be computed: the complex normal form style, or the real normal form style, that can be interesting in vibration problems since it allows to get more easily to real equations. In the present paper, only the complex normal form (CNF) will be used in the numerical examples. Interestingly, the best choice of parametrisation for a given problem is still an open problem, the reader is referred for example to (Stoychev and Römer 2022) for dedicated examples where different styles are compared.

The second main problem in deriving arbitrary order solutions of the homological equations is to compute direct solutions, in order to bypass the expression in the modal coordinates. In (Vizzaccaro et al. 2022; Opreni et al. 2022c), this is achieved by using a bordering technique to solve for the homological equations. On the other hand, a norm-minimizing solution is adopted in (Jain and Haller 2022) to cope with the same problem.

Once the homological equations are solved at an arbitrary order, one has then access to all the coefficients in Equations (16), (18) and (20), meaning that the nonlinear mappings and the reduced dynamics, both for autonomous and non-autonomous terms, are available at high order. In the sequel, the selected order for the computations will be denoted as $O(o_a, o_{na})$.

## 2.3 Numerical techniques

The aim of the forthcoming numerical results is to compare the predictions given by the ROM to a reference solution obtained from the full-order model (FOM) simulation. The comparisons will be made on the frequency-response curves, obtained by exciting the rotating structure in the vicinity of one of its eigenmode by a harmonic force. In these cases, the spatial dependence of the forcing is assumed to be that of the excited eigenmode. Being orthogonal to the other linear modes, this assumption ensures that no direct energy is transferred to the other modes by the forcing.

The FRCs of the full-order model are computed via a continuation method directly implemented on all the dofs of the FE discretisation. The numerical method uses a harmonic balance technique to express the time-domain problem as an expansion onto the harmonics of the periodic solution. The continuation is then realized through an arclength method coupled to an alternating frequency/time (AFT) procedure (Cameron and Griffin 1989), and this solution will be referred to as HBFEM for Harmonic Balance on the Finite Element Method. As a general rule, five harmonics will be used in this HBFEM procedure as it provides a converged result with sufficient accuracy, given the range of amplitudes tested. The computational time associated to such simulations is generally very important (the order of one day). Consequently the computation of the stability of the solution branches is not done for the HBFEM solutions in order to reduce the burden. The ROMs are derived using the direct parametrisation implemented in MORFE (Opreni et al. 2022a). Autonomous problems are used to compute the backbone curves, while non-autonomous ROMs allows computing FRCs for direct comparison with FOM. The solution branches of periodic orbits are also computed thanks to a continuation method for the ROM.
## 2.4 Problem description

A one-meter long rotating cantilever beam, shown in Figure 1, is considered. The geometric dimensions are reported in Table 1. Its length is large as compared to the thickness and the width (respectively 2 cm and 3 cm) in order to excite easily geometric nonlinearity. The material is assumed to follow a linear elastic constitutive law, characterized by its Young modulus $E = 104$ GPa and Poisson ratio $\nu = 0.3$. The density is $\rho = 4400$ kg/m$^3$, those values being typical of titanium alloy. The axis of rotation is set to a distance $R = 10$ cm of the clamped end of the beam. The beam is discretised with 27-nodes hexahedral 3D finite elements, which are full quadratic Lagrangian elements. Thirty elements are used along the main length, while four elements are considered in the cross-section. This results in a total number of 1525 nodes and 4575 dofs.

<figure><img src="image_3.png" /><figcaption>**Figure 1** Schematic representation of the cantilever rotating beam.</figcaption></figure>

<table><caption>**Table 1** Dimensions, material properties and FE discretisation of the beam</caption><thead><tr><th>L [m]</th><th>h [m]</th><th>b [m]</th><th>R [m]</th><th>ρ [kg/m³]</th><th>E [N/m²]</th><th>ν [-]</th><th>nb. elem. [x × y × z]</th></tr></thead><tbody><tr><td>1</td><td>0.02</td><td>0.03</td><td>0.1</td><td>4400</td><td>104 × 10<sup>9</sup></td><td>0.3</td><td>30 × 2 × 2</td></tr></tbody></table>

Figure 2 reports some of the linear characteristics of the selected rotating beam. The eigenfrequency dependence on the rotation speed (Campbell diagram) is shown in Figure 2(a), where the modes have been labelled as follows: F for flap-wise bending i.e. bending vibration along the thinner direction (y direction in Figure 1; E for edge-wise bending, i.e. bending in thicker direction (z direction in Figure 1); while T stands for torsion. As expected, the eigenfrequencies along edge-wise bending are larger than along flap-wise, and all bending frequencies are increasing with rotation speed as a result of the centrifugal stiffening. On the other hand, the eigenfrequency of the first torsional mode is almost not affected by rotation.

Figures 2(c) and 2(d) show the dependence of the eigenmode shapes as function of the rotation speed, for the first two bending modes in flap-wise direction (1F and 2F), and the torsion mode. As expected, the torsion mode shape is almost independent of the rotation speed. On the other hand, the bending modes show an important dependence, in line with previously reported results, see e.g. (Thomas et al. 2016). Finally, Figure 2(b) displays the longitudinal displacement of the rotating cantilever for three different rotation speeds, illustrating how the centrifugal effect elongates the beam.

## 2.5 Nonlinear vibrations of the fundamental bending mode and hardening/softening behaviour

In this section, a detailed study of the nonlinear vibrations of the first flap-wise bending mode, denoted as 1F, is reported. The convergence of the ROM is first illustrated. Then, interpolation of the coefficients of the reduced dynamics is done in order to study with more accuracy the hardening/softening transition with rotation speed. Finally, it is shown that a light interpolation with very few points is needed to reproduce with sufficient accuracy the behaviour of the FOM on a large range of rotation speeds.

### 2.5.1 Validation of the ROM

The nonlinear vibrations of mode 1F is known to be of the hardening type without rotation (Nayfeh and Pai 2004; Touzé and Thomas 2004), while the centrifugal effect leads to a transition from hardening to softening behaviour (Thomas et al. 2016). Figure 3 shows a comparison of the FOM
<figure><img src="image_4.png" /><figcaption>**Figure 2** Linear characteristics of the rotating cantilever beam illustrating the dependence with rotation speed Ω. (a) Campbell diagram. Color code for the eigenfrequencies: [—] 1F; [—] 1E; [—] 2F; [—] 2E; [—] 3F; [—] 3E; [—] 4F; [—] 1T; [—] 4E. (b) Static position versus Ω. Colour code: [—] 0 rpm; [—] 2500 rpm; [—] 5000 rpm. (c)-(d) Mode shape dependence upon rotation speeds for modes 1F and 1T, same colour code.</figcaption></figure>

and the ROM for five different rotation speeds, which is given in revolutions per minutes (rpm): 0 rpm, 500, 1000, 1500 and 2000 rpm. For these simulations, the damping has been set to 0.05 % on the 1F mode, and the forcing amplitude is chosen so that the linear transverse amplitude of the beam is 60 cm, in order to reach important vibration amplitudes where the geometric nonlinearity is excited. As shown in Figure 3(a), the transverse displacement is up to half the length of the beam, while the in-plane displacement reported in Figure 3(b) is up to 8 cm. Note that, since the system is stiffening due to rotation, the forcing amplitude factor α and the damping ratio needs to be increased with the rotation speed, in order to maintain a given level of vibration amplitude and the same amount of losses. The corresponding values of the forcing α, for each rotation speed, are reported in the caption of Figure 3. On the other hand, the damping ratio is increased such that it is always equal to 0.05 % when the rotation speed is varied.

Figure 3 shows a perfect agreement between the FOM simulation and the ROM, even for these large values of the vibration amplitude. The autonomous ROM with an order 7 in the asymptotic expansion has been used to draw out the backbone curves for each of the rotation speeds considered, while the FRCs have been computed with an order O(7, 6), meaning that the non-autonomous part is expanded up to order 6. This first result shows the accuracy of the reduction method, which is able to recover, with a single master mode and thus a single oscillator equation, the nonlinear behaviour of the FOM and the transition from hardening to softening. The computing time to obtain the complete FRC of the FOM is an heavy calculation and lasts approximately 1.5 day. On the other hand, the computation of the FRC with the ROM consisting of a single master mode is almost immediate and is obtained within 30 seconds. Computations have been realized with a standard desktop computer equipped with Intel® Xeon® Silver 4110 processor and 32 GB RAM.

Figure 4 reports the convergence of the non-autonomous part of the ROM, illustrating the importance of the time dependence of the forced invariant manifold on the reduced dynamics.
<figure><img src="image_5.png" /><figcaption>**Figure 3** FRC of the rotating cantilever beam. Comparison between the full order model (HBFEM) and the reduced order model for nonlinear vibrations along mode 1F. Three different curves are reported for each rotation speed: HBFEM solution, 5 harmonics, with continuous lines [---]; backbone curve computed from the autonomous ROM, order $O(7)$, continuation with 7 harmonics, dotted lines [-----]; FRC computed with non-autonomous DPIM $O(7, 6)$, dashed lines, 7 harmonics [---]. Colour code for rotation speeds and corresponding forcing amplitude factors: [---] 0 rpm, $\alpha = 4.77$; [---] 500 rpm, $\alpha = 5.26$; [---] 1000 rpm, $\alpha = 6.67$; [---] 1500 rpm, $\alpha = 8.89$; [---] 2000 rpm, $\alpha = 11.79$.</figcaption></figure>

Indeed, previous numerical studies generally used a zero-order approximation for the terms brought by the forcing, as shown for example in (Touzé and Amabili 2006; Jain and Haller 2022), and the importance of considering higher-order terms for the forcing has been recently illustrated in (Opreni et al. 2022c). Note that, for the sake of brevity, the convergence of the autonomous part is not addressed in this Section, since numerical results underlined that, for the range of amplitudes reached (up to 0.6 times the length of the beam), a parametrisation up to order 5 gives a good result with very small discrepancy as compared to the FOM, while an order 7 allows obtaining a converged result for the ROM; these numerical outcomes being completely in line with previous studies on a cantilever reported in (Vizzaccaro et al. 2022) for the autonomous problem, and in (Opreni et al. 2022c) for the forced system.

Figure 4 shows the convergence by increasing the order of the non-autonomous part, while keeping the autonomous part at order 7, for two different rotation speeds (beam at rest, 0 rpm, and beam rotating at 2000 rpm), and for two different vibration amplitudes. A high amplitude,

<figure><img src="image_6.png" /><figcaption>**Figure 4** Convergence of the non-autonomous part of DPIM for the 1F mode. (a) FRC of the transverse displacement for two rotation speeds and different order of the DPIM formulation for the non-autonomous terms. (b) Zoom on the large amplitudes. (c) Zoom on the moderate amplitudes. Color code: HBFEM with 5 harmonics with dashed line [---]; DPIM $O(7, 0)$ to DPIM $O(7, 6)$, 7 harmonics in the continuation, with continuous lines from softer to darker colours: [---] DPIM $O(7, 0)$; [---] DPIM $O(7, 2)$; [---] DPIM $O(7, 4)$; [---] DPIM $O(7, 6)$. Two forcing amplitudes at 0 rpm: [---] 0 rpm, $\alpha = 4.77$; [---] 0 rpm, $\alpha = 2.38$. Two forcing amplitudes at 2000 rpm: [---] 2000 rpm, $\alpha = 11.79$; [---] 2000 rpm, $\alpha = 5.89$.</figcaption></figure>

corresponding to the case highlighted in Figure 3, is contrasted to a smaller amplitude, resulting in transverse vibrations of 1/4 the length of the beam. Again, the values of the forcing amplitude are rotation dependent because of the global stiffening of the structure, they have thus been adjusted as function of the rotation speed in order to keep a vibration amplitude of the same magnitude. The selected values are reported in the caption. For the small vibration amplitude, the
DPIM ROM with order $O(7, 0)$ slightly overpredicts the maximum vibration amplitude, and this discrepancy increases importantly for the high amplitude tested, such that the ROM $O(7, 0)$ is not reliable anymore. On the other hand, increasing the order to 2 for the non-autonomous part is sufficient to recover a good approximation, while a fully converged result is obtained for the large amplitude case with $O(7, 4)$. One can see that the conclusions about this convergence are independent from the rotation speed, underlining that once the vibration amplitudes are reaching half the length of the cantilever beam, then an order 0 on the non-autonomous part is clearly not sufficient to produce accurate ROM predictions, and an order 2 on the non-autonomous part is at least needed.

For completeness, an additional comparison of the ROM with a full order model is provided in Appendix C, where the case of a much thinner beam is considered. In this case, the reference results are taken from previously published results given in (Thomas et al. 2016), where a different model is used to make the computation, since 1D beam elements implementing the Timoshenko assumptions, are considered. The results in this case needs to take into account the difference in the models since the calculations are here reported with 3D FE.

### 2.5.2 Accurate interpolation of the autonomous ROM

In this section, an accurate interpolation of the coefficients of the ROM is detailed, in order to investigate more precisely the transition from the hardening to the softening behaviour. Since this transition is governed by the backbone curve of the mode, only the autonomous ROM is considered here, which means that damping and forcing are discarded in this section. They will be reintroduced in Section 2.5.4. The ROM is here obtained with a single master mode assumption. Moreover, the complex normal form style (CNF) is used. With a single master mode, the normal coordinate has two entries, that one can denote as $z = [z_1 \ z_2]^T$. Since $z_2$ is the complex conjugate of $z_1$, $z_2 = \bar{z}_1$, the reduced dynamics can be formally written only for the coordinate $z_1$ as (Vizzaccaro et al. 2022; Opreni et al. 2022c):

$$ \dot{z}_1 = i\omega_1 z_1 + if^{(3)} z_1^2 \bar{z}_1 + if^{(5)} z_1^3 \bar{z}_1 + \dots + if^{(o_a)} z_1^{\frac{o_a+1}{2}} \bar{z}_1^{\frac{o_a-1}{2}} + O(z_1^{o_a+1}). \quad (21) $$

In this expression, $o_a$ is the order retained for the approximation of the autonomous part, and $f^{(p)}$ refers to the order-$p$ coefficient computed by the parametrisation method with CNF. One can note in particular that only odd orders are appearing, since they correspond to the only trivially resonant monomials to be considered in such a normal form approach (Touzé 2014). The coefficients are purely imaginary (which has been underlined by factorizing i), and the equation for $z_2 = \bar{z}_1$ is simply the complex conjugate of Equation (21).

The coefficients of the ROM, recalled in Equation (21), will be interpolated from order 1 (corresponding to the radian eigenfrequency) up to order 11, with respect to the rotation speed $\Omega \in [0, 1500]$rpm. Note that such interpolation has already been proposed in previous studies, see e.g. (Sternchüss and Balmès 2006; Khalifeh et al. 2017) where the tangent stiffness matrix is directly interpolated in order to draw out models for arbitrary values of rotation speed in a selected interval. In (Sternchüss and Balmès 2006; Khalifeh et al. 2017), the interpolation was realized with polynomial containing only even powers of the rotation speed, as a consequence of the expressions of the spin softening matrix and the centrifugal effects, see Equations (2a) and (2b), that shows an explicit dependence on $\Omega^2$. Here the perspective is a bit different since the coefficients of the ROM are directly interpolated. In the expressions of these coefficients, see e.g. (Touzé et al. 2004; Touzé and Amabili 2006; Vizzaccaro et al. 2022; Opreni et al. 2022c) where some analytical details are given for a few of them, complex algebraic manipulations appear both at the numerator and the denominators, such that it appears difficult to predict beforehand if all these coefficients might only depend on even powers of the rotation speed. Consequently this assumption is not retained here for the analysis. One can also note that the dependence of the eigenfrequencies with $\Omega$ is a result from the eigenvalue problem, where the tangent stiffness matrix is given by Equation (7). Because of the appearance of the nonlinear quadratic and cubic operators in Equation (7), there is no clear evidence that the eigenfrequencies should only depend on even powers of $\Omega$.
<table><caption>**Table 2** Examples of Chebyshev nodes in the range [0, 1500].</caption><tr><td>1 point</td><td>750</td></tr><tr><td>2 points</td><td>219.67 1280.33</td></tr><tr><td>3 points</td><td>100.48 750 1399.52</td></tr><tr><td>4 points</td><td>57.09 462.99 1037.01 1442.91</td></tr><tr><td>5 points</td><td>36.71 309.16 750 1190.84 1463.29</td></tr><tr><td>6 points</td><td>25.56 219.67 555.89 944.11 1280.33 1474.44</td></tr></table>

The interpolations are computed using Lagrange polynomials of order *l*, computed on sets of *l* + 1 Chebyshev points. Appendix B recalls the formula needed to find the Chebyshev points in an arbitrary interval [*a*, *b*], as well as the expressions of the Lagrange polynomials. For the sake of illustration, the Chebyshev nodes for the selected interval of rotation speeds [0, 1500]rpm, are given in Table 2.

Figure 5 shows the variation of the coefficients of the ROM from order 3 to order 11. A reference behaviour for each of these coefficients has been computed with a high order polynomial by using 50 Chebyshev points. Then the low-order polynomials of order 2 to 5 are shown for each of them. Qualitatively, it appears that an interpolation with a degree 2 is rather insufficient to ensure a very accurate result. On the other hand, a polynomial with degree 3 seems to offer a sufficiently accurate interpolation of the ROM on the selected interval.

<figure><img src="image_7.png" /><figcaption>**Figure 5** Interpolation of the coefficients of the autonomous ROM. (a)-(e) coefficients $f^{(p)}$ of orders $p = 3$ to 11. (f) Root mean square error for the coefficients of the ROM as function of the number of Chebyshev points used to compute the polynomial interpolation. Colour code for (a)-(e): [●] reference interpolation with 50 points (polynomial of degree 49); [—] interpolation, polynomial of degree 2; [—] degree 3; [—] degree 4. Colour code for (f): [●●] coef. order 1 (eigenfrequency); [●●●] order 3; [●●●●] order 5; [●●●●●] order 7; [●●●●●●] order 9; [●●●●●●●] order 11.</figcaption></figure>

In order to give a more quantitative description of the convergence of the interpolated coefficients, Figure 5(f) shows the Root-Mean-Square Error (RMSE) between the reference (computed using 50 interpolation points), and the lower degree interpolations, for all the coefficients from order 1 (eigenfrequency) to order 11. Because the type of nonlinearity is changing from hardening to softening, the coefficients of the ROM are vanishing at certain values of the rotation speed. Consequently it is not possible to use a classical relative error since this would lead to infinite values at the vanishing points. The RMSE allows coping with such problems. It is defined as (Atkinson 1989)

$$E(z) = \frac{1}{\sqrt{b-a}} \sqrt{\sum_{i=1}^{n_p} [f(z_i) - r(z_i)]^2} = \frac{\|\mathbf{f} - \mathbf{r}\|_2}{\sqrt{b-a}}, \quad (22)$$

where **f** is the function to interpolate; **r** is the computed interpolation of **f**. Both functions **f**
and **r** are defined on the interval [a, b], and $n_p$ is the number of points used to discretise the interval [a, b]. One limitation of this definition for the RMSE is the fact that the amplitude of the data is not taken into account, which could lead to very disperse curves if comparing different interpolations, which is the case here since many coefficients are investigated. In order to obtain a normalized criterion usable for all the coefficients, the following definition of the RMSE has been finally used:

$$E_n(z) = \frac{1}{\max \mathbf{f} - \min \mathbf{f}} \sqrt{\sum_{i=1}^{n_p} [f(z_i) - r(z_i)]^2} = \frac{\|\mathbf{f} - \mathbf{r}\|_2}{\max \mathbf{f} - \min \mathbf{f}} \quad (23)$$

Figure 5(f) confirms the convergence trend of the interpolation. One can observe in particular that the important gain of using a polynomial of degree 3 as compared to degree 2, which was qualitatively observed from the behaviour of the different coefficients, is confirmed quantitatively. Secondly, the dependence of the coefficients with only even powers of the rotation speed, does not appear clearly in the trend, underlining that using arbitrary degree of polynomial is indeed a good choice. Finally, the figure shows that a large number of points might be needed to obtain a high accuracy. This question will be further investigated in Section 2.5.4, but with a different perspective. Indeed, the question will be to investigate the minimum number of interpolation points needed to obtain a sufficient accuracy of the ROM on a given interval of rotation speed. Furthermore, the non-autonomous part of the ROM will also be considered in Section 2.5.4, thus enlarging the scope of interpolating the ROM to the amplitude of the forcing, such that ROMs could be available to compute FRC on a given interval of rotation speed and forcing amplitude.

### 2.5.3 Detailed study of the hardening/softening transition

The accurately interpolated ROMs are here used in order to investigate more closely the transition from hardening to softening behaviour of the 1F mode with increasing rotation speed. As underlined in the previous section, each coefficient $f^{(p)}$ from $p = 3$ has a vanishing point on the selected interval [0, 1500]rpm. This behaviour is more closely reported in Figures 6(a) and 6(b), where the pair of opposites coefficients corresponding to each of the two equations of the reduced dynamics (i.e. $z_1$ and $z_2$ which are simply complex conjugates, leading to opposite coefficients), are shown for each order. The critical rotation speed $\Omega_c^{0p}$ is defined as the value for which each of these coefficients $f^{(p)}$ vanishes. Their exact values can be computed thanks to the interpolated models, and are reported in Table 3 for the different orders from $p = 3$ to $p = 11$.

<figure><img src="image_8.png" /><figcaption>**Figure 6** Hardening/softening behaviour of the rotating beam. (a) Behaviour of the nonlinear coefficient $f^{(p)}$, for $p = 3$ to 11, as function of the rotation speed. (b) Zoom on the behaviour around 500 rpm. (c) Backbone curves for increasing values of $\Omega$. Colour code for (a)-(b): [—] $f^{(3)}$; [—] $f^{(5)}$; [—] $f^{(7)}$; [—] $f^{(9)}$; [—] $f^{(11)}$. Colour code for (c): six different rotation speeds: [—] $\Omega = 400$ rpm; [—] $\Omega = 468.32$ rpm; [—] $\Omega = 578.23$ rpm; [—] $\Omega = 661.98$ rpm; [—] $\Omega = 852.57$ rpm; [—] $\Omega = 1010$ rpm. For each rotation speeds, three backbones obtained by truncating the ROM at different orders are shown: [—] order 3; [—] order 7; [—] order 11.</figcaption></figure>

The most important critical rotation speed needed to predict the change from hardening to softening behaviour is $\Omega_c^{03}$. Indeed the cubic term in Equation (21) drives the first-order
perturbation of the backbone and thus dictates the type of nonlinearity (Touzé et al. 2004). This coefficient has been used in past studies to predict for instance the type of nonlinearity of shallow spherical-cap shells as function of the geometry (Touzé and Thomas 2006), or for imperfect plates (Touzé et al. 2008). According to this criterion, the change of hardening/softening behaviour for this cantilever rotating beam occurs at $\Omega_c^{03} = 1008.63$ rpm, in line with the FRCs shown in Section 2.5.1.

<table><caption>**Table 3** Critical rotation speeds $\Omega_c^{0p}$ of the various orders.</caption><thead><tr><th>Ω<sub>c</sub><sup>03</sup> [rpm]</th><th>Ω<sub>c</sub><sup>05</sup> [rpm]</th><th>Ω<sub>c</sub><sup>07</sup> [rpm]</th><th>Ω<sub>c</sub><sup>09</sup> [rpm]</th><th>Ω<sub>c</sub><sup>11</sup> [rpm]</th></tr></thead><tbody><tr><td>1008.63</td><td>696.52</td><td>627.44</td><td>529.02</td><td>407.46</td></tr></tbody></table>

Higher order developments allows one to obtain a more accurate picture of this transition, since all orders in the asymptotic expansion of the backbone are at hand. Interestingly, the critical rotation speed (value at which the coefficient changes sign) for the next orders appears at smaller and smaller values. As stated in Table 3, the order 5 coefficient vanishes at 696 rpm, order 7 at 627 rpm, etc. This means that the transition starts for smaller values than the critical one, but this effect is also related to larger amplitudes of vibration. Indeed, the first coefficients to vanish are the highest ones, and then propagates to the lowest. Since the highest coefficients gets more important at large vibration amplitude values on the backbone, this means that even though the behaviour is hardening for small vibration amplitude, the backbone starts to bend and change curvature at very large amplitudes.

This phenomenon is illustrated in Figure 6(c), showing the evolution of the backbone curves for different rotation speeds that have been selected in order to fall between the critical rotation speeds gathered in Table 3. Also, in order to highlight more clearly the effects of the different orders in the computation of the backbone, three different solutions are shown for each rotation speed, namely $O(3)$, $O(7)$ and $O(11)$. For rotation speeds smaller than $\Omega_c^{03} = 1008.63$ rpm, the behaviour is indeed hardening but one can observe the effect of the higher orders on the backbones that start perturb it from very large amplitudes by changing the curvature. Note that this phenomena appears at extremely large amplitudes: here the backbones are plotted up to vibration amplitudes reaching the length of the beam. At $\Omega = 852$ rpm, only the cubic term is positive while all others are negative. Consequently, even if the behaviour is first hardening it turns to softening type, meaning that there is a region before $\Omega_c^{03}$ where the behaviour can be denoted as “mixed”. Finally for $\Omega = 1010$ rpm, the cubic order term is almost zero such that the $O(3)$ backbone is quasi linear. Taking into account the higher orders shows clearly the softening curvature.

### 2.5.4 Minimal interpolation for predictive modeling of the non-autonomous ROM

In this section, the interpolation method is further investigated and extended to the non-autonomous problem, hence taking into account the external forcing. The aim is to illustrate the minimal interpolation needed in order to deliver accurate ROMs that are able to compute the FRCs of the 1F mode for a range of rotation speed as well as a range of forcing amplitudes.

The method proposed in (Opreni et al. 2022c) in order to take into account the time dependence of the invariant manifold, is to insert a first-order development of the non-autonomous part in the invariance equation, following previous studies where the same idea was proposed (Haro and de la Llave 2006; Breunung and Haller 2018). As a main consequence of this choice, the added terms in the reduced dynamics take the form of a product of the forcing term (at power 1) with monomials (of arbitrary even orders) of the normal coordinates. No terms of higher orders in the forcing (e.g. forcing at power two or three) can be derived because of the starting assumption. In the context of interpolating the ROMs, one can take advantage of this property to easily extend to the case of forcing. Since only linear terms with the forcing appears in the reduced dynamics, then only a single forcing amplitude is needed to interpolate linearly all the terms involving the forcing. In practice, the coefficients of the non-autonomous ROM are thus computed for a unitary value of the forcing amplitude factor $\alpha$, i.e. for $\alpha = 1$. Then, any value of the forcing can then be used to compute the FRC of the interpolated ROM, one has just to multiply the coefficients of the non-autonomous terms, by the selected $\alpha$ value. Of course, one needs just to take care of not
using too large values of $\alpha$ in order to stay in the range where the ROM is valid (i.e. small strains to excite only geometric nonlinearities).

Based on the results shown in the previous section for the autonomous problem, it is awaited that an interpolated ROM computed with only 3 points (degree two interpolation) will lead to small errors, while using 4 points (degree 3) should be sufficient. This is illustrated in Figure 7, where the interpolated ROM is used to compute the FRC of the rotating cantilever at rotation speeds that are different from the Chebyshev nodes. To that purpose, $\Omega = 500$ and 1000 rpm have been selected as falling out of the Chebyshev points reported in Table 2 for 3, 4 and 5 points. In each case, the FRCs computed by the interpolated ROMs are compared to a reference FOM solution.

<figure><img src="image_9.png" /><figcaption>**Figure 7** Comparisons of FRCs obtained from interpolated ROMs with FOM solution. In (c), backbones obtained by normalizing the frequency axis in order to highlight the difference in terms of nonlinearity. Colour code for (a)-(b): Reference HBFEM solution computed with 5 harmonics [-----]. Interpolated ROMs computed with 3 calculation points: [—] DPIM $O(7, 6)$ for the FRC, [---] backbone $O(7)$. Interpolated ROMs computed with 4 points: [—] DPIM $O(7, 6)$ for the FRC, [---] backbone $O(7)$. Colour code for (c): [—] backbone $O(7)$, ROM interpolated with 3 points; [—] backbone $O(7)$, 4 points; [---] backbone $O(7)$, 5 points.</figcaption></figure>

With only 3 points to interpolate the coefficients of the ROM (such that each variations are approximated by a polynomial of degree 2), one can observe in Figure 7 that the eigenfrequency is not well recovered, such that a shift of the FRC is obtained in both cases for $\Omega = 500$ and 1000 rpm. This shift is however rather small (the relative error in the eigenfrequency being less than 1 %), and can be recovered by a single linear calculation if one wants to feed the ROM with the exact value. Interestingly, the small errors in the nonlinear coefficients of the ROM are almost not visible on the FRCs which clearly depicts the correct curvature. This is more quantitatively assessed in Figure 7(c) where the normalized backbones are reported for each case of interpolated ROMs. As a general result, one can observe that using 4 points to interpolate the ROM gives perfect results, while using 3 points leads to very small discrepancies that are largely acceptable. It is also very important to underline that the non-autonomous part is interpolated using a single forcing amplitude, from which all the coefficients can be computed.

As a conclusion, very efficient interpolations can be computed thanks to the parsimonious representation of reduced dynamics offered by the complex normal form (CNF). Indeed, the DPIM with CNF offers a reduced dynamics that has the minimal number of monomials and coefficients to handle, since reducing the dynamics to the sole resonant complex monomials. Numerical results underline that the interpolated ROMs can extend easily to the non-autonomous case. FRCs can thus be computed with accuracy on a large range of rotation speeds and forcing amplitudes. Note that these properties can be efficiently used for design purposes. Here the rotation speed can be seen as a varying parameter but other cases can be studied, for example parametrized shape defects can be inserted, see e.g. (Marconi et al. 2020), or varying material properties, etc.

## 2.6 Nonlinear vibrations of the second bending mode

In this section, the nonlinear vibrations of the second bending mode are investigated and the reduction method is applied by selecting mode 2F as master. The forcing frequency is in the vicinity of the second flapwise bending mode, and it spatial dependence is that of the corresponding mode shape, as in the previous case. The damping has been set to 0.05 %. Since
the eigenfrequency is larger, the problem is stiffer and it is more difficult to reach very large amplitudes of vibration. Even by importantly increasing the values of the forcing amplitude $\alpha$ as compared to the previous section, the range of maximum vibration amplitudes investigated for mode 2F is about 6 cm at the tip.

Figure 8 compares the FRC produced by the full-order model as compared to that computed by the ROM consisting of a single master mode. As expected from previous studies, the nonlinear behaviour of the second mode is of the softening type without rotation (Nayfeh and Pai 2004; Touzé and Thomas 2004), and does not show important variations with the rotation speed (Thomas et al. 2016). In Figure 8, three different values of rotation speed have been selected, namely $\Omega = 0$, 2500 and 5000 rpm. One can observe a mild variation of the softening behaviour which is slightly decreasing.

<figure><img src="image_10.png" /><figcaption>**Figure 8** FRC of the rotating cantilever beam. Comparison between the full order model (HBFEM) and the reduced order model for nonlinear vibrations along mode 2F. Three curves are reported for each rotation speed: HBFEM solution, 5 harmonics [---]; backbone curve computed from the autonomous ROM, order $O(7)$, continuation with 7 harmonics [-----]; FRC computed with non-autonomous DPIM $O(7, 6)$ and 7 harmonics [---]. Colour code for the rotation speeds and corresponding forcing amplitude factors: [---] 0 rpm, $\alpha = 18.69$; [---] 2500 rpm, $\alpha = 36.70$; [---] 5000 rpm, $\alpha = 87.86$. [---] Backward continuation with decreasing excitation frequencies, 5000 rpm, $\alpha = 87.86$.</figcaption></figure>

The general trend of the behaviour is very well reproduced by the ROM. In this case, a solution with order $O(7, 6)$ is reported, where a converged result is obtained, in line with the previous section dealing with mode 1F. However, the full order model shows an important bifurcated branch which is probably due to an internal resonance. Interestingly, this bifurcated branch is present for the three tested rotating speeds. At 0 rpm, it occurs at large amplitudes, around 6.6 cm, not shown in the present figure. Then, by increasing $\Omega$, this bifurcated branch move along the amplitude, but is still present, for example it occurs at an amplitude around 4 cm for $\Omega = 2500$ and 5000 rpm. A distinctive feature is also given by the fact that for the 5000 rpm case, two bifurcated branches are found with the HBFEM solution, depending on the choice of increasing or decreasing the frequency sweep during the continuation, see Figure 8 where the forward continuation led to the branch reported in orange, while the red solution branch is found by decreasing the frequency. This underlines a complex bifurcation scenario, difficult to analyze.

By analyzing the modal amplitudes along the bifurcated branch, it has been found an important participation of the first four flapwise bending modes. Different ROMs have been built by including up to the first four bending modes as master. However, none of them have been able to recover the bifurcated branch shown by the FOM. Instead, the computed FRCs displayed by the tested ROMS with up to four master modes showed the exact same results as those presented in Figure 8, and are thus not added to the plot. This underlines that the resonance relationship leading to the bifurcated branch is probably driven by a higher-order internal resonance that has not been found and might also involve other high frequency modes. It is also likely that such a resonance might appear on the nonlinear frequencies rather than on the linear ones that are here tested. And detecting internal resonances on the nonlinear frequencies is a difficult task since one needs to compute the backbones of each mode to get the exact dependence upon amplitude, and test the resonance relationships for different amplitude values.
Even though the ROM has not been able to recover this complex bifurcated branch, such phenomena are known to be difficult to analyze and reproduce. Also they appear for very specific values in the parameter space, whereas the most likely scenario is to recover the unimodal behaviour reproduced by the ROM, which indeed operates a filtering of the complex, high-order internal resonances. On the other hand, one can underline that the ROM allows a perfect prediction of the amplitudes of the HBFEM model for all other excitation frequencies.

## 3 Second test case: a twisted plate

In this section, the reduction method is applied to a more complex structure in order to demonstrate the applicability to a more industrial test case. To that purpose, a rotating twisted clamped-free plate, resembling a fan blade, is used. The twisting creates important quadratic nonlinearities that makes the problem much more difficult to handle than the cantilever beam because of important shell effects. The nonlinear vibrations of the fundamental flexural mode are investigated.

### 3.1 Presentation and linear behaviour

The selected structure, with a simplified geometry as compared to a real fan blade, is shown in Figure 9(a). The dimensions have been selected by taking inspirations from existing industrial blades. New aircraft engines have a growing input diameter, such that the fan diameter also increases. As example, we can cite the CFM leap engine which fan blades are ~ 0.75 m long; while next generation engine with the open rotor project have fan blades which are up to 1.2 m long. As a compromise between these two cases, a 1 m-long plate has been selected, while the width has been taken equal to 20 cm. A uniform thickness of 1.5 cm has also been selected. A global twist angle of 60° has been imposed between the foot and the tip, with a quadratic twisting dependence upon the length, in order to have a more important twist at the tip. The boundary conditions are clamped at the foot and free on other three sides, and the axis of rotation has been shifted from a distance of 30 cm from the foot. The axis of rotation is also chosen aligned with the foot of the blade. Material properties (Young modulus, Poisson ratio and density) have the same values as for the beam, in order to represent an isotropic blade made of Titanium alloy (see Table 1). This is again a simplification as compared to real blades that are generally made of composite materials. However, the aim of the present study is to highlight an application of the reduction method to a more complex geometry. Additional effects like anisotropy due to the composite material of such blades is beyond the scope of the present developments.

The blade is meshed with 120 27-nodes hexaedric elements. The studied model is thus composed of 1365 nodes and 4095 dofs.

Some of the linear characteristics of the blade are reported in Figure 9. The mode shapes of the first three flexural modes are shown together with the first torsional mode, for the case without rotation, in Figure 9(c). The centrifugal effect leads to a new static position of the blade which is illustrated in Figure 9(b) for Ω = 3000 rpm. Interestingly for this geometry, the plate elongates but has a twisted angle that increases with Ω. The moment produced by centrifugal loading can be divided into two components: one related to geometric characteristics that twist the blade and one component related to the effect of axial and torsion forces that untwist the blade. In the case of a perfect twisted flat plate, the centre of pressure and gravity lie on the same line and consequently the second untwist component of the moment due to centrifugal loading is null (Golinval 2002). In summary, the geometric property of the flat plate explains the twisted behaviour observed with the increase in rotational speed.

The Campbell diagram is shown in Figure 9(d). In line with the behaviour found with the cantilever beam, the eigenfrequencies of the flexural modes are increasing with the rotation speed, due to the global stiffening effect. The quantitative values of the radian eigenfrequencies of mode 1F are reported in Table 4 for different rotation speeds; these values will be used to normalize the frequency axis in the next section. For the twisted plate, the first torsional mode shows a dependence with Ω, in contrast with the observation reported for the beam. A crossing between the second flexural mode and the torsional mode occurs at Ω ≈ 850 rpm, which should lead to a 1 : 1 internally resonant dynamics. It is worth noticing that this Campbell diagram is for
<figure><img src="image_11.png" /><figcaption>(a) Geometry</figcaption></figure>

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

1F

2F

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

1T

1E

<figure><img src="image_14.png" /><figcaption>(b) Static equilibrium at $\Omega = 3000$ rpm</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(d) Campbell diagram</figcaption></figure>

**Figure 9** Linear characteristics of the twisted plate representative of a fan blade. In (b), static equilibrium position at $\Omega = 3000$ rpm (with colours) compared to the undeformed geometry (grey); amplitudes magnified by a factor 10 for readability. (d) colour code: [—] 1F; [—] 2F; [—] 1T; [—] 1E; [—] 3F.

a single blade. Once mounted on a disk to form a full fan disk, the global behaviour is more complicated and out of the scope of the present study. Nevertheless, one can underline that when mounted in order to form a full fan disk, cyclic symmetry can be applied to reduce the model, see e.g. (Quaegbeur et al. 2023).

<table><caption>**Table 4** Eigenfrequency of first flexural mode 1F for various rotation speeds.</caption><tr><td>Ω [rpm]</td><td>0</td><td>1000</td><td>2000</td><td>3000</td></tr><tr><td>ω<sub>1F</sub> [rad/s]</td><td>60.77</td><td>100.51</td><td>166.04</td><td>234.37</td></tr></table>

## 3.2 Nonlinear vibrations of the fundamental bending mode

Figure 10 shows the frequency-response curves obtained for the first bending mode, and for four different rotation speeds: 0, 1000, 2000 and 3000 rpm. The results obtained by the FOM are contrasted to those obtained by assuming a single master mode in the reduction procedure. The damping value has here been set to 0.05 % on the first bending mode and applied through the mass matrix. The FRC for HBFEM model has been obtained using five harmonics. The non-autonomous ROMs have been computed with an order $O(5, 4)$, convergence studies revealing that the result is converged with such an order, given the amplitudes tested.

Without rotation, the blade shows a hardening behaviour. A flat untwisted plate is known to display a hardening behaviour. However, the imposed twist creates important quadratic nonlinearities that might change the behaviour. However, for the twist angle selected (60°), the behaviour is still of the hardening type. As for the beam, the centrifugal effect leads to a transition from hardening to softening behaviour, occurring between 1000 and 2000 rpm for this FE discretisation.

A perfect match is obtained between the ROM and the FOM, for all the rotation speed tested, and up to the vibration amplitude shown here which is already large for such a blade (20 cm vibration at the tip). As for the beam, and because of the global stiffening of the structure, it has been selected to increase the value of the forcing $\alpha$ in order to maintain the peak vibration amplitude to the same level, around 0.2 m. As reported in Figure 10, $\alpha$ is thus varying from 1.98
at 0 rpm, to 31.5 at 3000 rpm. For this last value, the HBFEM solution depicts a small bifurcated branch that is not recovered by the single master mode analysis of the ROM. As for the mode 2F of the beam, this bifurcated branch might involve a high-order internal resonance that is difficult to qualify properly. Nevertheless, it occurs on a very limited range of the parameter values. Different ROMs including up to four master modes have been tested in order to see if a simple model could capture this bifurcated branch. However, none of them has been able to recover what seems to be a high-order internal resonance.

<figure><img src="image_16.png" /><figcaption>**Figure 10** FRC of the twisted plate (fan blade). Comparison between the full-order model (HBFEM) and the reduced-order model for nonlinear vibrations along mode 1F for four rotational velocities: [—] 0 rpm, $\alpha = 1.98$; [—] 1000 rpm, $\alpha = 5.57$; [—] 2000 rpm, $\alpha = 15.61$; [—] 3000 rpm, $\alpha = 31.5$. Reference HBFEM solution [—] and DPIM $O(5, 4)$ ROM solution [---].</figcaption></figure>

### 3.3 Effect of the discretisation

In this section, the effect of increasing the number of dofs of the FE model is investigated. This point has not been questioned before since the goal was to demonstrate that the ROM is able to recover the behaviour of the FOM, starting from the same discretisation. Since it has been demonstrated in the previous section that a perfect match between ROM and FOM is at hand, here only the ROM is studied, and the convergence of the solution with respect to the discretisation is investigated. All the computations reported in this section has been realized with a standard laptop having 16 GB of RAM, thus underlining that involved computations with model size comparable to industrial needs can be performed with a light computational environment thanks to the proposed reduction method. This is not the case for a FOM simulation that needs important memory access.

To that purpose, different size meshes are given in Table 5, ranging from a coarse one with 4095 dofs to a refine one with 86 751 dofs, pertaining to reasonable mesh that can be handled with 16 GB RAM. Table 5 also shows the convergence of the first bending mode 1F radian eigenfrequency (third row), while the last row shows the difference (in percent) between the most accurate calculation (mesh 40 × 3 × 25 with 86 751 dofs taken as reference), and the current configuration. In order to take into account rotation, the convergence is here shown for a case where the rotation speed is equal to 3000 rpm.

<table><caption>**Table 5** Convergence of the first flexural mode eigenfrequency $\omega_{1F}$ with respect to mesh refinement.</caption><thead><tr><td>mesh number</td><td>FE discretisation</td><td>size [dof]</td><td>ω<sub>1F</sub> at Ω = 3000 rpm [rad/s]</td><td>ε [%]</td></tr></thead><tbody><tr><td>1</td><td>10 × 2 × 6</td><td>4095</td><td>234.37</td><td>1.25</td></tr><tr><td>2</td><td>12 × 2 × 8</td><td>6375</td><td>233.41</td><td>0.84</td></tr><tr><td>3</td><td>15 × 2 × 10</td><td>9765</td><td>232.62</td><td>0.5</td></tr><tr><td>4</td><td>20 × 2 × 15</td><td>19 065</td><td>232.02</td><td>0.24</td></tr><tr><td>5</td><td>30 × 2 × 20</td><td>37 515</td><td>231.62</td><td>0.065</td></tr><tr><td>6</td><td>40 × 3 × 25</td><td>86 751</td><td>231.48</td><td>/</td></tr></tbody></table>

Figure 11(a) reports the computational times needed to construct the ROM on a laptop equipped with 16 GB RAM, for an order $O(5, 4)$ of the reduction. The computing time for the 37 515 dofs mesh is still reasonable with less than 7 minutes to construct the ROM. Note also that this calculation corresponds to the offline stage, needed to construct the ROM from the mesh, and it has to be done only once. Then to compute FRCs, corresponding to the online stage, the
difference between the FOM simulation and the ROM is impressive, order of the day versus less
than one minute.

Figure 11(b) shows the convergence of the FRCs with respect to mesh refinement. The two last points correspond to meshes 4 and 5 in Table 5, and underlines that convergence of the FRC is reached. The convergence study has been made for four different rotation speeds, but for the sake of clarity, only two are shown in Figure 11(b), namely 0 and 1000 rpm. The convergent behaviour for 2000 and 3000 rpm followed the exact same trend. As a final result, the FRC for these four rotation speed of a spatially converged solution (mesh 5, 37 515 dofs) is reported in Figure 11(c). In all the computations reported here, DPIM with order O(5, 4) is considered. Interestingly, the global behaviour is a bit different from the one reported in the previous section with mesh 1 (4095 dofs), underlining that the spatial discretisation of the previous model was not converged. Here the behaviour is globally softer, and the hardening/softening transition with increasing rotation speeds occurs before 1000 rpm.

<figure><img src="image_17.png" /><figcaption>**Figure 11** Computing time and convergence analysis with respect to mesh refinement for the fan blade. (a) Computing time to construct the ROM (offline stage) DPIM *O*(5, 4) with respect to the number of dofs. (b) FRC convergence with increasing mesh refinement, for two different rotation speeds: [---] Ω = 0 rpm, α = 1.972; [---] Ω = 1000 rpm, α = 5.56. Colour code for mesh refinements: [---] mesh 1; [---] mesh 2; [---] mesh 3; [---] mesh 4; [---] mesh 5. (c) FRC obtained with the DPIM *O*(5, 4) ROM, mesh 5, four rotation speeds: [---] Ω = 0 rpm, α = 1.972; [---] Ω = 1000 rpm, α = 5.56; [---] Ω = 2000 rpm, α = 15.53; [---] Ω = 3000 rpm, α = 31.2.</figcaption></figure>

3.4 Interpolation of the ROM

To conclude this analysis, the interpolation procedure used for the beam is here employed, in order to find the minimal interpolation needed to achieve accurate results on a given range of rotation speeds and forcing frequencies. As highlighted for the beam case, the peculiarity of the treatment of the non-autonomous problem, see (Opreni et al. 2022c), allows us to take benefit of the linear dependence of the non-autonomous terms with the forcing amplitude. Consequently, a single simulation, with a unitary forcing amplitude $\alpha = 1$, is needed and can be used to deduce the coefficients of the non-autonomous terms for any forcing amplitude. The interpolation of the autonomous terms is conducted as in the previous case, using Lagrange interpolation on Chebyshev points, for a given, selected interval of $\Omega$. Note that more adverse conditions as compared to the beam are selected here. The plate is twisted, hence creating quadratic nonlinearities, and the increase of the twist with the rotation speed leads to enhance the importance of the quadratic terms with rotation. Also, the reduction method using complex normal form (CNF) uses only resonant monomials and hence only odd power terms in the reduced dynamics; meaning that all the even terms are computed and reported to the next orders in the algebraic automated manipulations.

Figure 12 shows the results obtained, with the idea of finding the minimal number of interpolation points in order to recover accurate predictions. The interpolation interval is [0, 3000] rpm. Four selected rotation speeds, falling out of the Chebyshev nodes, have been taken in order to test the predictions given by the interpolated ROMs out of the learning data set: $\Omega = 0, 1000, 2000$ and $3000$ rpm. For each of these cases, a reference solution is obtained by calculating directly the ROM solution at the given rotation speed, and this result is compared to the one given by interpolated ROMs. The computations have been made using the finer mesh
with 37 515 dofs. Interpolated ROMs with small number of calculation points are tested against
this reference, from 4 points (interpolation with Lagrange polynomial of degree 3) up to 7 points
(degree 6). The reference ROM is in red in Figure 12, and the lighter interpolations are in blue,
with light blue for the smallest degree (3) to dark blue for the highest one (6).

<figure><img src="image_18.png" /><figcaption>**Figure 12** FRC of interpolated ROMs with increasing degrees compared to reference solution, for the first flexural mode 1F of the fan blade, for four rotation speeds. (a) $\Omega = 0$ rpm, $\alpha = 1.97$; (b) $\Omega = 1000$ rpm, $\alpha = 5.56$; (c) $\Omega = 2000$ rpm, $\alpha = 15.53$; (d) $\Omega = 3000$ rpm, $\alpha = 31.2$. [—] Reference solution DPIM $O(7, 6)$ computed at the exact rotation speed. Interpolated ROMs DPIM $O(7, 6)$ with increasing number of learning points: [—] 4 points (degree 3); [—] 5 points; [—] 6 points; [—] 7 points.</figcaption></figure>

Since the case is more difficult than the beam, the convergence of the interpolated ROMs
is longer to obtain. This is the consequence of a more complex behaviour of the coefficients
of the ROM with Ω, reflecting also the more difficult calculations needed to obtain the higher
orders and eliminating the even power terms. Nevertheless, convergence is achieved and one can
observe that 7 points of computation (interpolation with a polynomial of order 6) is enough to
reproduce correctly the FRCs of the system. Again, it is also observed that the convergence of the
eigenfrequency is key and drive the analysis of the convergence, while the nonlinear behaviour is
correctly estimated with a lower order interpolation. This result indicates that the method can be
advantageously used and be parameterized to obtain a converged collection of ROMs with a
varying parameter.

**4 Conclusion**

In this contribution, the direct parametrisation method for invariant manifold has been applied for model-order reduction of geometrically nonlinear rotating structures featuring centrifugal effect. The key features of the method is to be directly applicable to FE discretisation, and to propose arbitrary order expansions, hence allowing one to obtain converged results up to large amplitude vibrations. Applicability of the method to two different cases, a cantilever beam and a twisted plate representative of a fan blade, has been illustrated. In each cases, it has been shown that the method is able to perfectly retrieve the solution curve provided by a full-order solution. In terms of computational gains, the effect is dramatic since FOM simulations typically last one day, versus less than one minute for the ROM.

ROMs with one master mode have a very limited number of coefficients, hence making easy
to interpolate their values such that one is able to have reduced models for a large range of
continuous values of the rotation speed. The interpolation has been applied in order to investigate more closely the hardening/softening behaviour transition of the first bending mode of the cantilever beam. It has been highlighted that higher-order terms cancels first, while the cubic term, which dictates the type of nonlinearity, vanishes last, for increasing values of Ω. The minimal number of points needed to provide interpolated ROMs with sufficient accuracy, has then been investigated. An interesting feature of the reduction method for the non-autonomous term is to provide only linear terms as function of the forcing. Therefore, a single, unitary amplitude forcing simulation, is needed to interpolate for any forcing value. For the coefficients of the autonomous system, it has been found that 4 points are sufficient in order to interpolate the cantilever beam on a range of rotation speed. For the twisted plate, the problem at hand is much more difficult, with quadratic nonlinearities, that are also increasing with rotation since the twisting of the position at rest increases with Ω. For that structure, 6 to 7 interpolation points are needed to provide highly accurate results for any rotational speed.

All the results presented in this article highlights the excellent accuracy and the very important gains in computing time that are offered by the reduction method. It underlines that using nonlinear normal modes as invariant manifold of the phase space is an efficient tool for fast and accurate reduced-order modeling of geometrically nonlinear structures.

The next targeted development is to incorporate the Coriolis force in the reduction process and implement it in the software MORFE. No theoretical limitation hinders this development, and results have already been reported in the literature where the parametrisation method is used with Coriolis effect in (Li et al. 2022).

For applications to engineering structures like helicopter or turbine blades, the proposed reduction method can thus be directly applied to deal with the structural nonlinearities. However, in such problems, nonlinearities arising from the aerodynamic loads need also to be taken into account. While this is yet out of the scope of the present study, one can see that numerous different options are available to extend the results presented here. A first option would be to consider heuristic models for the aerodynamics effect, that can be appended to the equations of the ROM. The opposite option would be to consider the fully coupled fluid-structure interaction problem and introduce the parametrisation method in such framework, which asks for very important developments.

## A Effect of the Coriolis force

This appendix is devoted to the justification of neglecting the Coriolis force in the case under study in the present paper, by comparing its effect to the centrifugal force. The equations of motion for a rotating structure considering the Coriolis force leads to change Equation (1) for

$$ M\ddot{\mathbf{p}} + [\mathbf{C} + \mathbf{G}]\dot{\mathbf{p}} + [\mathbf{K} + \mathbf{N}]\mathbf{p} + \mathbf{g}(\mathbf{p}, \mathbf{p}) + \mathbf{h}(\mathbf{p}, \mathbf{p}, \mathbf{p}) = \mathbf{f}^{\text{ext}} + \mathbf{f}^{\Omega}, \quad (A.1) $$

where the added term is the Coriolis matrix G which depends on the rotation speed. Its expression in a finite element approach is given by

$$ \mathbf{G} = 2 \int_V \rho \mathbf{H}^\top \Omega \mathbf{H} dV, \quad (A.2) $$

where Ω is the rotation matrix. The rotation matrix is by construction skewsymmetric and thus gives also this property to the Coriolis matrix. This expression can be compared to Equation (2a) for evaluating the relative effects of Coriolis term with respect to centrifugal ones. In the specific case of the beam studied in Section 2.4, the rotation matrix can be expressed as

$$ \Omega = \begin{bmatrix} 0 & -\Omega_z & 0 \\ \Omega_z & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}, \quad (A.3) $$

where only the rotation along z axis has been considered.
In order to better assess the magnitude of the Coriolis force, it appears easier to use its more classical expression given by

$$ \mathbf{f}_C = -2m\omega \times \mathbf{v} = -2m \begin{bmatrix} 0 \\ 0 \\ \omega_z \end{bmatrix} \times \begin{bmatrix} v_x \\ v_y \\ v_z \end{bmatrix} = 2m \begin{bmatrix} \omega_z v_y \\ -\omega_z v_x \\ 0 \end{bmatrix}. \quad (A.4) $$

It appears clearly that this force has non-vanishing components only in the x and y directions, following the axis representation given in Figure 1. The y direction is that of the flexural motions of the beam and may have the most important effect on the vibrations of the structure. However, this term is proportional to $v_x$ the longitudinal vibration velocity, which is negligible since axial vibrations are very small in the present context of nonlinear flexural vibration. In the same line, the effect of the term in x direction, even though proportional to $v_y$ which mat be large, is awaited to be very small since acting on the stiffer direction of the structure. Consequently the effect of this term is globally awaited to be very small.

In order to illustrate this assumption, we compute the eigenfrequencies of the rotating beam with and without Coriolis effect. The result is given in Figure A.1 for rotation speed between 0 and 5000 rpm. We can clearly observe that the impact of Coriolis effect on these eigenfrequencies is negligible. In order to give a quantitative comparison, the maximum difference observed between the two cases considered (with and without Coriolis force) for the first nine eigenfrequencies shown in Figure A.1 is 0.8 %. The Coriolis effect can thus safely be neglected in the general Equation (1).

<figure><img src="image_19.png" /><figcaption>Figure A.1 Campbell diagram of the beam. [—] without Coriolis effect; [•] with Coriolis effect. Color codes for the eigenfrequencies: [—] 1F; [—] 1E; [—] 2F; [—] 2E; [—] 3F; [—] 3E; [—] 4F; [—] 1T; [—] 4E.</figcaption></figure>

## B Interpolation: Lagrange polynomials and Chebyshev points

The interpolations are computed thanks to Lagrange polynomials, see e.g. (Atkinson 1989) for a detailed presentation. Here only the main formula are recalled for the sake of completeness. Considering $n + 1$ points $(x_0, y_0), \dots, (x_n, y_n)$, the only polynomial with minimal degree $n$, passing exactly through the given points, is the Lagrange polynomial, whose expression reads

$$ \mathbf{p}(x) = \sum_{j=0}^{n} y_j \left[ \prod_{i=0, i \neq j}^{n} \frac{x - x_i}{x_j - x_i} \right]. \quad (B.1) $$

When dealing with interpolation problems, it is customary to use the Chebyshev nodes since they minimize the interpolation error (Atkinson 1989). Chebyshev points are defined on a normalized interval $[-1, 1]$ as

$$ x_i = \cos\left(\frac{2i + 1}{2n + 2}\pi\right) \quad \text{for } i = 0, 1, \dots, n. \quad (B.2) $$

In order to translate the location of these nodes to an arbitrary interval $[a, b]$, one can use

$$ x_i = \frac{a+b}{2} + \frac{b-a}{2} \cos\left(\frac{2i+1}{2n+2}\pi\right). \quad (B.3) $$
## C Comparison with existing results on a thin beam

In this section, a further comparison of the ROM with full-order simulation results is proposed. Existing results already published in (Thomas et al. 2016) are selected to compare against the ROM predictions. Moreover, since the FRCs in (Thomas et al. 2016) considered only direct time integration, the outcome of the model recently proposed in (Debeurre et al. 2023), which couples continuation methods to the Timoshenko beam model, are also used, in order to draw out a complete comparison. The considered structure is a simple straight beam of dimensions 1 m × 10 cm × 5 mm, made of titanium. The damping model is mass-proportional, with a damping ratio of 0.5 % on the studied mode and the beam is forced by a sine point load applied at the tip of the beam, with constant transverse direction in the global frame. In (Thomas et al. 2016), the FE discretisation is made thanks to 20 Timoshenko elements with linear shape functions and reduced integration, to avoid shear locking. Since time integration was used in (Thomas et al. 2016), only some points of the FRC were obtained, by waiting for the steady state. To summarize, the outcomes of two full-order models, both using Timoshenko beam elements, are here compared against ROMs derived thanks to the high-order DPIM, discretised with 3D elements. A special attention is thus devoted to the different kinematical assumptions used in the respective models.

The beam selected in (Thomas et al. 2016) has a rather large width (1/10 of the length), such that plate effects might not be completely ignored. To track the effects of the width and fairly compare the outcomes of the Timoshenko beam models to the 3D simulations, two different beams are tested with the 3D FE discretisation. The first one has the exact same dimension as the one used in (Thomas et al. 2016). It is meshed with 30 × 2 × 8 27-nodes hexaedral FE (leading to 5185 nodes and 15 555 dofs), and is referred to as 3DP (the acronym standing for 3D model with plate-like behaviour due to the important width). A second beam, with the same dimensions except for the width which has been decreased, is also considered. This second model has for final dimensions: 1 m × 1 cm × 5 mm, and is discretised with a 60 × 2 × 2 H27 mesh, equivalent to 9075 dofs. This second 3D model is called 3DB for 3D structure with beam-like behaviour. In order to keep other quantities equal (forcing, damping, bending rigidity) between the models, the nondimensionalisation rules, as given for example in (Debeurre et al. 2023), have been used.

Table C.1 compares the eigenfrequencies computed with the different models: Timoshenko beam model from (Thomas et al. 2016), and the two 3D models with large and small width: 3DP and 3DB. One can see in particular that the converged mesh for the 3DP model does not coincide

<table><caption>**Table C.1** Eigenfrequencies of the first two bending modes 1F and 2F. Comparison between the Timoshenko beam model from (Thomas et al. 2016), the 3DP model (3D elements, plate like behaviour due to important width: 1 m × 10 cm × 5 mm) and the 3DB model (3D elements, small width: 1 m × 1 cm × 5 mm).</caption><thead><tr><th></th><th>0 rpm</th><th>500 rpm</th><th>1000 rpm</th><th>2000 rpm</th><th>3000 rpm</th></tr></thead><tbody><tr><td>f<sub>1F</sub> [Hz] (Thomas et al. 2016)</td><td>3.927</td><td>10.288</td><td>19.068</td><td>36.876</td><td>54.720</td></tr><tr><td>f<sub>1F</sub> [Hz] 3DP</td><td>3.961</td><td>10.322</td><td>19.114</td><td>36.95</td><td>54.826</td></tr><tr><td>ε<sub>1F</sub> [%]</td><td>0.875</td><td>0.325</td><td>0.244</td><td>0.201</td><td>0.194</td></tr><tr><td>f<sub>1F</sub> [Hz] 3DB</td><td>3.932</td><td>10.294</td><td>19.078</td><td>36.895</td><td>54.751</td></tr><tr><td>ε<sub>1F</sub> [%]</td><td>0.127</td><td>0.058</td><td>0.052</td><td>0.052</td><td>0.057</td></tr><tr><td>f<sub>2F</sub> [Hz] (Thomas et al. 2016)</td><td>24.644</td><td>33.414</td><td>51.229</td><td>92.268</td><td>134.822</td></tr><tr><td>f<sub>2F</sub> [Hz] 3DP</td><td>24.837</td><td>33.575</td><td>51.371</td><td>92.425</td><td>135.021</td></tr><tr><td>ε<sub>2F</sub> [%]</td><td>0.78</td><td>0.481</td><td>0.277</td><td>0.17</td><td>0.148</td></tr><tr><td>f<sub>2F</sub> [Hz] 3DB</td><td>24.642</td><td>33.416</td><td>51.234</td><td>92.276</td><td>134.835</td></tr><tr><td>ε<sub>2F</sub> [%]</td><td>0.008</td><td>0.00599</td><td>0.0976</td><td>0.00867</td><td>0.0096</td></tr></tbody></table>

perfectly with the beam model, the 3DP model being slightly stiffer. This is also confirmed on the FRCs shown in Figure C.2, where the 3DP model shows an enhanced hardening behaviour as compared to the Timoshenko beam model.

Interestingly, this difference is more pronounced at Ω = 0 rpm, and then decreases if the rotation speed Ω increases (see Table C.1). The same trend was observed for the nonlinear FRC, the two models being closer as Ω increases. Those FRCs results are not reported here for a sake of conciseness. In Figure C.2, the results of (Thomas et al. 2016) are also compared to computations
obtained with the same Timoshenko beam FE discretisation, where a continuation method with harmonic balance and asymptotic-numerical method is appended to the calculations in order to obtain the complete solution branches of the FRCs (Debeurre et al. 2023). In this case, the match

<figure><img src="image_20.png" /><figcaption>**Figure C.2** FRC of the cantilever beam studied in (Thomas et al. 2016). Comparison between the ROM obtained with DPIM and numerical results on the full-order model obtained from (Debeurre et al. 2023) (continuation method) and (Thomas et al. 2016) (direct time integration): [-] results from (Thomas et al. 2016), time integration; [-----] Timoshenko model with continuation from (Debeurre et al. 2023); [---] DPIM $O(5, 4)$ with large width, model 3DP displaying plate-like effect; [—] DPIM $O(5, 4)$ with small width, model 3DB.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure C.3** Frequency-response curves for the rotating beam, comparison between the predictions given by the DPIM applied to 3DB model, and the results from (Thomas et al. 2016). (a)-(f): first bending mode 1F at $\Omega = 0, 50, 100, 500, 1000$ and 3000 rpm, respectively. (g)-(i): second bending mode 2F at $\Omega = 0, 1000$ and 3000 rpm, respectively. Colour code: [—] DPIM $O(5, 4)$ on the 3DP structure; [\*] direct time integration from (Thomas et al. 2016).</figcaption></figure>

with the time integration is perfect for the computed points, meaning that the two Timoshenko models solved out with different numerical techniques (continuation vs. direct time integration),
provide the exact same result. One can note, as usual with time integration, that due to the shrinking of the basin of attraction at the top of the FRCs, the points obtained from the numerical time integration reported in (Thomas et al. 2016) for the second mode, failed to reach the top of the curve.

Now comparing the results of the second 3DB model, with smaller width, one can observe first in Table C.1 a much better coincidence with the Timoshenko beam model, the differences in the eigenfrequencies being severely decreased. Comparing the FRCs in Figure C.2, the match is perfect between this 3DB model and the computations reported in (Thomas et al. 2016). In fact, the linear and nonlinear dynamics of geometrically exact straight beams do not depend on their width, as seen with dimensionless quantities in (Debeurre et al. 2023). However, a beam model (associated to plane stress assumption in the width direction) tends to a plate model (with plain strain assumptions), slightly stiffer, as the width increases, the beam model becoming slightly erroneous for “plate like” beams. As a consequence, the following comparison between the test case of (Thomas et al. 2016) and the ROMs obtained with the DPIM will be done with the 3DB geometry, which has a smaller width.

The comparison of the DPIM (computed with the 3DB “beam like” geometry) with the results of (Thomas et al. 2016) are shown in Figure C.3 for varying rotating speeds. It can be observed that the match between the two models is perfect in term of bending of the FRC due to the geometrical nonlinearities associated to the centrifugal force. Some discrepancies can be noticed at the top of the FRC associated to jumps, which were not reached by the time-integrations of (Thomas et al. 2016). Those comparisons perfectly validate the four models and associated computational methods (1D beam elements in time integration and continuation, 3D elements with continuation (HBFEM solution) and in DPIM).

## References

Allen, M., D. Rixen, M. van der Seijs, P. Tiso, T. Abrahamsson, and R. Mayes (2020). *Substructuring in Engineering Dynamics*. Vol. 594. CISM International Centre for Mechanical Sciences. Springer. [DOI], [OA].

Atkinson, K. (1989). *An Introduction to Numerical Analysis*. 2nd ed. Wiley. ISBN: 9780471624899.

Bathe, K. J. (1996). *Finite Element Procedures*. Prentice-Hall.

Bauchau, O., C. Bottasso, and Y. Nikishkov (2001). Modeling rotorcraft dynamics with finite element multibody procedures. *Mathematical and Computer Modelling* 33(10-11):1113-1137. [DOI].

Blahoš, J. (2022). Parallel Harmonic Balance Method for Analysis of Nonlinear Mechanical Systems. PhD thesis. Imperial College London. [HDL].

Blahoš, J., A. Vizzaccaro, L. Salles, and F. El Haddad (2020). Parallel Harmonic Balance Method for analysis of nonlinear dynamical systems. *Volume 11: Structures and Dynamics: Structural Mechanics, Vibration, and Damping; Supercritical CO2* (London, UK, Sept. 21, 2020–Sept. 25, 2020). GT2020-15392. ASME. [DOI], [OA].

Breunung, T. and G. Haller (2018). Explicit backbone curves from spectral submanifolds of forced-damped nonlinear mechanical systems. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences* 474(2213):20180083. [DOI], [OA].

Cabré, X., E. Fontich, and R. de la Llave (2005). The parameterization method for invariant manifolds III: overview and applications. *Journal of Differential Equations* 218(2):444-515. [DOI], [OA].

Cameron, T. M. and J. H. Griffin (1989). An Alternating Frequency/Time domain method for calculating the steady-state response of nonlinear dynamic systems. *Journal of Applied Mechanics* 56(1):149-154. [DOI], [OA].

Chaviaropoulos, P. (1996). Development of a state-of-the art aeroelastic simulator for horizontal axis wind turbines: Part 1: Structural aspects. *Wind Engineering* 20(6):405-421.

Craig, R. and M. Bampton (1968). Coupling of substructures for dynamic analyses. *AIAA Journal* 6(7):1313-1319. [DOI], [HAL].

DaSilva, M. C. and D. Hodges (1986). Nonlinear flexure and torsion of rotating beams, with application to helicopter rotor blades - I. Formulation. *Vertica* 10 (2):151-169. [HAL].
De Klerk, D., D. Rixen, and S. Voormeeren (2008). General framework for dynamic substructuring: History, review and classification of techniques. *AIAA Journal* 46(5):1169-1181. [DOI].

Debeurre, M., A. Grolet, B. Cochelin, and O. Thomas (2023). Finite element computation of nonlinear modes and frequency response of geometrically exact beam structures. *Journal of Sound and Vibration* 548:117534. [DOI], [HAL].

Delhez, E., F. Nyssen, J.-C. Golinval, and A. Batailly (2021). Reduced-order modeling of blades with geometric nonlinearities and contact interactions. *Journal of Sound and Vibration* 500:116037. [DOI], [HAL].

Di Cristofaro, D., A. Opreni, M. Cremonesi, R. Carminati, and A. Frangi (2022). An Arbitrary Lagrangian Eulerian approach for estimating energy dissipation in micromirrors. *Actuators* 11(10):298. [DOI], [OA].

Di Palma, N., B. Chouvion, and F. Thouverez (2022). Parametric study on internal resonances for a simplified nonlinear blade model. *International Journal of Non-Linear Mechanics* 141:103941. [DOI], [OA].

Frangi, A. and G. Gobat (2019). Reduced-order modelling of the non-linear stiffness in MEMS resonators. *International Journal of Non-Linear Mechanics* 116:211-218. [DOI].

Friswell, M., J. Penny, S. Garvey, and A. Lees (2010). *Dynamics of Rotating Machines*. Cambridge University Press. [DOI].

Golinval, J.-C. (2002). *Lecture notes in Mechanical Design of Turbomachinery*. University of Liège.

Guyan, R. J. (1965). Reduction of stiffness and mass matrices. *AIAA Journal* 3(2):380-380. [DOI], [HAL].

Guérin, N., A. Thorin, F. Thouverez, M. Legrand, and P. Almeida (2018). Thermomechanical model reduction for efficient simulations of rotor-stator contact interaction. *Journal of Engineering for Gas Turbines and Power* 141(2). [DOI], [HAL].

Géradin, M. and D. Rixen (2015). *Mechanical Vibrations: Theory and Applications to Structural Dynamics*. 3rd ed. Wiley. ISBN: 9781118900192.

Haller, G. and S. Ponsioen (2016). Nonlinear normal modes and spectral submanifolds: existence, uniqueness and use in model reduction. *Nonlinear Dynamics* 86(3):1493-1534. [DOI], [ARXIV].

Haller, G. and S. Ponsioen (2017). Exact model reduction by a slow-fast decomposition of nonlinear mechanical systems. *Nonlinear Dynamics* 90(1):617-647. [DOI], [ARXIV].

Haro, À., M. Canadell, J.-L. Figueras, A. Luque, and J. M. Mondelo (2016). *The Parameterization Method for Invariant Manifolds. From Rigorous Results to Effective Computations*. Applied Mathematical Sciences. Springer International Publishing. [DOI].

Haro, À. and R. de la Llave (2006). A parameterization method for the computation of invariant tori and their whiskers in quasi-periodic maps: rigorous results. *Journal of Differential Equations* 228(2):530-579. [DOI].

Hodges, D. and E. H. Dowell (1974). *Nonlinear equations of motion for the elastic bending and torsion of twisted nonuniform rotor blades*. Tech. rep. NASA technical note. [OA].

Hollkamp, J. and R. Gordon (2008). Reduced-order models for non-linear response prediction: Implicit condensation and expansion. *Journal of Sound and Vibration* 318(4-5):1139-1153. [DOI].

Holzapfel, G. (2000). *Nonlinear Solid Mechanics: A Continuum Approach for Engineering*. Wiley. ISBN: 9780471823193.

Idelsohn, S. and A. Cardona (1985). A reduction method for nonlinear structural dynamic analysis. *Computer Methods in Applied Mechanics and Engineering* 49(3):253-279. [DOI].

Jain, S., M. Li, and G. Haller (2023). *SSM tool 2.3: Computation of invariant manifolds in high-dimensional mechanics problems*. URL: https://github.com/jain-shobhit/SSMTool.

Jain, S. and G. Haller (2022). How to compute invariant manifolds and their reduced dynamics in high-dimensional finite-element models? *Nonlinear Dynamics* 107(2):1417-1450. [DOI], [OA].

Jain, S., P. Tiso, J. Rutzmoser, and D. Rixen (2017). A quadratic manifold for model order reduction of nonlinear structural dynamics. *Computers & Structures* 188:80-94. [DOI], [ARXIV].

Jiang, D., C. Pierre, and S. Shaw (2005). Nonlinear normal modes for vibratory systems under harmonic excitation. *Journal of Sound and Vibration* 288(4-5):791-812. [DOI].

Joannin, C., F. Thouverez, and B. Chouvion (2018). Reduced-order modelling using nonlinear modes and triple nonlinear modal synthesis. *Computers & Structures* 203:18-33. [DOI], [HAL].
Khalifeh, E., E. Piollet, A. Millecamps, and A. Batailly (2017). Non-linear modeling of centrifugal stiffening effects for accurate bladed component reduced-order models. *Volume 7B: Structures and Dynamics* (Charlotte, United States, June 26, 2017–June 30, 2017). ASME. [DOI], [HAL].

Laxalde, D. and F. Thouverez (2009). Complex non-linear modal analysis for mechanical systems: Application to turbomachinery bladings with friction interfaces. *Journal of Sound and Vibration* 322(4-5):1009–1025. [DOI], [HAL].

Legrand, M., C. Pierre, P. Cartraud, and J.-P. Lombard (2009). Two-dimensional modeling of an aircraft engine structural bladed disk-casing modal interaction. *Journal of Sound and Vibration* 319(1-2):366–391. [DOI], [HAL].

Li, M., S. Jain, and G. Haller (2022). Nonlinear analysis of forced mechanical systems with internal resonance using spectral submanifolds – Part I: Periodic response and forced response curve. *Nonlinear Dynamics* 110(2):1005–1043. [DOI], [OA].

MacNeal, R. (1971). A hybrid method of component mode synthesis. *Computers & Structures* 1(4):581–601. [DOI], [HAL].

Manolas, D. I., V. A. Riziotis, and S. G. Voutsinas (2015). Assessing the importance of geometric nonlinear effects in the prediction of wind turbine blade loads. *Journal of Computational and Nonlinear Dynamics* 10(4). [DOI].

Marconi, J., P. Tiso, and F. Braghin (2020). A nonlinear reduced-order model with parameterized shape defects. *Computer Methods in Applied Mechanics and Engineering* 360:112785. [DOI], [OA].

Mignolet, M., A. Przekop, S. Rizzi, and M. Spottswood (2013). A review of indirect/non-intrusive reduced-order modeling of nonlinear geometric structures. *Journal of Sound and Vibration* 332(10):2437–2460. [DOI].

Muravyov, A. and S. Rizzi (2003). Determination of nonlinear stiffness with application to random vibration of geometrically nonlinear structures. *Computers & Structures* 81(15):1513–1523. [DOI].

Nayfeh, A. and F. Pai (2004). *Linear and Nonlinear Structural Mechanics*. Wiley. [DOI], [OA].

Opreni, A., A. Vizzaccaro, A. Martin, A. Frangi, and C. Touzé (2022a). *MORFEinvariantManifold.jl: Nonlinear model order reduction of large dimensional Finite Element models using the direct parametrisation method for invariant manifolds*. [URL].

Opreni, A., M. Furlan, A. Bursuc, N. Boni, G. Mendicino, R. Carminati, and A. Frangi (2022b). One-to-one internal resonance in a symmetric MEMS micromirror. *Applied Physics Letters* 121(17):173501. [DOI].

Opreni, A., A. Vizzaccaro, A. Frangi, and C. Touzé (2021). Model order reduction based on direct normal form: application to large finite element MEMS structures featuring internal resonance. *Nonlinear Dynamics* 105:1237–1272. [DOI], [OA].

Opreni, A., A. Vizzaccaro, C. Touzé, and A. Frangi (2022c). High order direct parametrisation of invariant manifolds for model order reduction of finite element structures: application to generic forcing terms and parametrically excited systems. *Nonlinear Dynamics*. [DOI], [HAL].

Pescheck, E., C. Pierre, and S. Shaw (2002). A new Galerkin-based approach for accurate non-linear normal modes through invariant manifolds. *Journal of Sound and Vibration* 249(5):971–993. [DOI], [HAL].

Petrov, E. (2011). A high-accuracy model reduction for analysis of nonlinear vibrations in structures with contact interfaces. *Journal of Engineering for Gas Turbines and Power* 133(10). [DOI].

Ponsioen, S., T. Pedergnana, and G. Haller (2019). Analytic prediction of isolated forced response curves from spectral submanifolds. *Nonlinear Dynamics* 98(4):2755–2773. [DOI], [ARXIV].

Ponsioen, S., T. Pedergnana, and G. Haller (2018). Automated computation of autonomous spectral submanifolds for nonlinear modal analysis. *Journal of Sound and Vibration* 420:269–295. [DOI], [ARXIV].

Quaegebeur, S., B. Chouvion, N. Di Palma, and F. Thouverez (2023). Model reduction of a cyclic symmetric structure exhibiting geometric nonlinearity with a normal form approach. *European Journal of Mechanics - A/Solids* 97:104822. [DOI], [OA].

Rao, S. J. (1991). *Turbomachine Blade Vibration*. New Age International. ISBN: 9780470217641.

Rezaei, M., H. Zohoor, and H. Haddadpour (2018). Aeroelastic modeling and dynamic analysis of
a wind turbine rotor by considering geometric nonlinearities. *Journal of Sound and Vibration* 432:653-679. [DOI].

Rosenberg, R. M. (1962). The normal modes of nonlinear *n*-degree-of-freedom systems. *Journal of Applied Mechanics* 29(1):7-14. [DOI], [HAL].

Shaw, S. and C. Pierre (1991). Non-linear normal modes and invariant manifolds. *Journal of Sound and Vibration* 150(1):170-173. [DOI], [HAL].

Shaw, S. and C. Pierre (1993). normal modes for non-linear vibratory systems. *Journal of Sound and Vibration* 164(1):85-124. [DOI], [OA].

Shen, Y., N. Béreux, A. Frangi, and C. Touzé (2021a). Reduced-order models for geometrically nonlinear structures: Assessment of implicit condensation in comparison with invariant manifold approach. *European Journal of Mechanics - A/Solids* 86:104165. [DOI], [OA].

Shen, Y., N. Kesmia, C. Touzé, A. Vizzaccaro, L. Salles, and O. Thomas (2021b). Predicting the type of nonlinearity of shallow spherical shells: Comparison of direct normal form with modal derivatives. *Second International Nonlinear Dynamics Conference (NODYCON 2021)* (Rome (Virtual), Italy, Feb. 16, 2021-Feb. 19, 2021). Springer, pp 361-371. [DOI], [HAL].

Shen, Y., A. Vizzaccaro, N. Kesmia, T. Yu, L. Salles, O. Thomas, and C. Touzé (2021c). Comparison of reduction methods for finite element geometrically nonlinear beam structures. *Vibration* 4(1):175-204. [DOI], [OA].

Sinha, S. (2013). Rotordynamic analysis of asymmetric turbofan rotor due to fan blade-loss event with contact-impact rub loads. *Journal of Sound and Vibration* 332(9):2253-2283. [DOI], [HAL].

Sternchüss, A. and E. Balmès (2006). On the reduction of quasi-cyclic disk models with variable rotation speeds. *International Conference on Noise and Vibration Engineering (ISMA2006)* (Leuven, Belgium, Sept. 18, 2006-Sept. 20, 2006). [HAL].

Stoychev, A. and U. J. Römer (2022). Failing parametrizations: what can go wrong when approximating spectral submanifolds. *Nonlinear Dynamics* 111:5963-6000. [DOI], [OA].

Thomas, O., A. Sénéchal, and J.-F. Deü (2016). Hardening/softening behavior and reduced-order modeling of nonlinear vibrations of rotating cantilever beams. *Nonlinear Dynamics* 86(2):1293-1318. [DOI], [HAL].

Thorin, A., N. Guérin, M. Legrand, F. Thouverez, and P. Almeida (2018). Nonsmooth thermoelastic simulations of blade-casing contact interactions. *Journal of Engineering for Gas Turbines and Power* 141(2). [DOI], [HAL].

Touzé, C. and O. Thomas (2004). Reduced-order modeling for a cantilever beam subjected to harmonic forcing. *EUROMECH Colloquium 457: Nonlinear modes of vibrating systems* (Fréjus, France, June 7, 2004-June 9, 2004), pp 165-168. [HAL].

Touzé, C. and M. Amabili (2006). Non-linear normal modes for damped geometrically non-linear systems: Application to reduced-order modeling of harmonically forced structures. *Journal of Sound and Vibration* 298(4-5):958-981. [DOI], [HAL].

Touzé, C. and O. Thomas (2006). Non-linear behaviour of free-edge shallow spherical shells: effect of the geometry. *International Journal of Non-Linear Mechanics* 41(5):678-692. [DOI], [HAL].

Touzé, C., O. Thomas, and A. Chaigne (2004). Hardening/softening behaviour in non-linear oscillations of structural systems using non-linear normal modes. *Journal of Sound and Vibration* 273(1-2):77-101. [DOI], [HAL].

Touzé, C. (2014). Normal form theory and nonlinear normal modes: Theoretical settings and applications. *Modal Analysis of Nonlinear Mechanical Systems*. Ed. by G. Kerschen. Vol. 555. CISM courses and lectures. Springer, pp 75-160. [DOI], [HAL].

Touzé, C., C. Camier, G. Favraud, and O. Thomas (2008). Effect of imperfections and damping on the type of nonlinearity of circular plates and shallow spherical shells. *Mathematical Problems in Engineering* 2008:1-19. [DOI], [OA].

Touzé, C., A. Vizzaccaro, and O. Thomas (2021). Model order reduction methods for geometrically nonlinear structures: A review of nonlinear techniques. *Nonlinear Dynamics* 105(2):1141-1190. [DOI], [OA].

Van Zante, D. (2015). Progress in open rotor research: A U.S. perspective. *Aircraft Engine; Fans and Blowers; Marine* (Montréal, Canada, June 15, 2015-June 19, 2015). Vol. 1. ASME. [DOI], [HDL].

Vizzaccaro, A., A. Opreni, L. Salles, A. Frangi, and C. Touzé (2022). High order direct parametrisation of invariant manifolds for model order reduction of finite element structures: application to
large amplitude vibrations and uncovering of a folding point. *Nonlinear Dynamics* 110(1):525-571. [DOI], [OA].

Vizzaccaro, A., L. Salles, and C. Touzé (2021a). Comparison of nonlinear mappings for reduced-order modeling of vibrating structures: normal form theory and quadratic manifold method with modal derivatives. *Nonlinear Dynamics* 103(4):3335-3370. [DOI], [OA].

Vizzaccaro, A., Y. Shen, L. Salles, J. Blahoš, and C. Touzé (2021b). Direct computation of nonlinear mapping via normal form for reduced-order models of finite element nonlinear structures. *Computer Methods in Applied Mechanics and Engineering* 384:113957. [DOI], [ARXIV].

Wang, X. and M. Mignolet (2021). Discussion on “A non-intrusive model-order reduction of geometrically nonlinear structural dynamics using modal derivatives”. *Mechanical Systems and Signal Processing* 159:107638. [DOI].

Yoo, H. H., J. Hun Park, and J. Park (2001). Vibration analysis of rotating pre-twisted blades. *Computers & Structures* 79(19):1811-1819. [DOI].

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the authors—the copyright holder. To view a copy of this license, visit creativecommons.org/licenses/by/4.0.

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

**Authors' contributions** AM carried out the main part of study, defined the examples, performed the numerical simulations and drafted the manuscript. AO and AV developed the methodology and built the main parts of the Julia code implementing the reduction method, and AO developed the first version of the HBFEM code which has been updated for rotation in collaboration with AM. MD performed all the simulations shown in Appendix C related to the Timoshenko beam model with continuation. LS supervised the work, discussed applications to blades, and helped in designing and understanding the twisted plate model. AF supervised the work and help in the development of the methodology. OT helped in all discussions related to the comparisons with the thin beam example and wrote Appendix C. CT supervised the work, carried out most of the writing and developed the methodology. All authors read and approved the final manuscript.

**Supplementary Material** Data points necessary to reproduce the figures as well as a Software Heritage Identifier pointing to the version of the code used in the paper along with input parameters (meshes) are available at the [permalink](https://doi.org/10.5281/zenodo.7924472).

**Acknowledgements** The funding from AID (Agence de l'Innovation de Défense), project REMODEL, contract number 2020 65 0057 ENSTA, is gratefully acknowledged.

**Ethics approval and consent to participate** Not applicable.

**Consent for publication** Not applicable.

**Competing interests** The authors declare that they have no competing interests.

**Journal's Note** JTCAM remains neutral with regard to the content of the publication and institutional affiliations.