# Optimal Craig-Bampton mode selection for nonlinear flexible multibody analysis

Océane Topenot, Gaël Chevallier, Gaël Chevallier, Scott Cogan, Christophe Oulerich

► To cite this version:

Océane Topenot, Gaël Chevallier, Gaël Chevallier, Scott Cogan, Christophe Oulerich. Optimal Craig-Bampton mode selection for nonlinear flexible multibody analysis. Vibration, 2025, 8 (4), pp.81. 10.3390/vibration8040081. [hal-05453184]

HAL Id: hal-05453184

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

Submitted on 12 Jan 2026

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

Article

# Optimal Craig-Bampton Mode Selection for Nonlinear Flexible Multibody Analysis

Océane Topenot <fnref n="1" />,<fnref n="2" />*, Gaël Chevallier <fnref n="1" />,<fnref n="3" />,* , Scott Cogan <fnref n="1" /> and Christophe Oulerich <fnref n="2" />

## Abstract

Physics-based simulations are now widely employed in mechanical engineering. Flexible Multibody dynamic Simulations (FMBSs) have proven to be effective in representing the behavior of complex structures with local damping and stiffness nonlinearities. However, due to the broad range of component flexibilities as well as contact behavior between structural elements, time integration analyses can result in high computational burden. The challenge addressed in this article concerns the implementation of an efficient model reduction procedure in order to provide an acceptable tradeoff between calculation time and loss of accuracy in the prediction of system responses and dynamic loads. In most FMBS commercial software, the behavior of linear elastodynamic components is taken into account via imported Craig-Bampton superelements. In this context, dynamic mode selection techniques have been shown to provide a better order reduction than the standard low-frequency truncation. This article provides a review of dynamic mode selection methods that can be found in the literature, followed by a comparison based on simulations of an aircraft engine stator integrated in the full industrial engine model and tested on a speed ramp-up with unbalance.

**Keywords:** Craig-Bampton superelement; multibody simulation; local elastic modes; model order reduction; rotating machinery

## 1. Introduction

Flexible Multibody dynamic Simulations (FMBSs) are a powerful tool for analyzing the behavior of complex machinery under real-world operational conditions, especially in the field of rotating machines where large rotational displacements are involved [1]. The general framework was established several decades ago [2] but is still the topic of recent works [3]: each structural component is first modeled using the Finite Element Method (FEM), then reduced through Component Mode Synthesis (CMS) [4], and then assembled using various types of connections (clamp, pivot, etc.) available in commercial software (ADAMS, SIMPACK, and SIMSCAPE).

In the context of turbomachinery architecture evolution and the growing availability of computational resources, coupled with the expensive nature of experiments, engineers are increasingly leaning towards high-fidelity virtual prototyping. While Flexible Multibody

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

### check for updates

Academic Editor: Evgeny Petrov

Received: 28 October 2025

Revised: 28 November 2025

Accepted: 10 December 2025

Published: 18 December 2025

**Citation:** Topenot, O.; Chevallier, G.; Cogan, S.; Oulerich, C. Optimal Craig-Bampton Mode Selection for Nonlinear Flexible Multibody Analysis. *Vibration 2025*, *8*, 81. https://doi.org/10.3390/vibration8040081

**Copyright:** © 2025 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

<footnote n="1">¹ FEMTO-ST Institute, CNRS, University Marie and Louis Pasteur, 25000 Besançon, France; scott.cogan@univ-fcomte.fr</footnote>

<footnote n="2">² Safran Aircraft Engines, Rond-Point René Ravaud, 77550 Moissy-Cramayel, France; christophe.oulerich@safrangroup.com</footnote>

<footnote n="3">³ Research Center of the French Air Force Academy, French Air Force Academy, 13660 Salon-de-Provence, France</footnote>

<footnote>* Correspondence: oceane.topenot@safrangroup.com (O.T.); gael.chevallier@ecole-air.fr (G.C.)</footnote>
Software (FMBS) solvers exhibit commendable speed given the intricate nature of the models, the rising complexity of engine models leads to a substantial increase in computational burden. Indeed, the structures involve numerous components and connections featuring local damping or nonlinear stiffness. Furthermore, the behavior of these structures must be studied across a broad frequency spectrum. Consequently, the time integration for a single scenario can exceed ten hours, which is a deal-breaker in the design optimization phase.

Although the general framework for FMBS is now fairly standard, several significant scientific obstacles remain, particularly in the context of model reduction and CMS. Firstly, traditionally, CMS assumes linear behavior, but several contributors are seeking to incorporate nonlinearities into the behavior of these components: [5,6], for geometric nonlinearity; [7], for material nonlinearity; or [8], for friction nonlinearities. Secondly, there is great interest in CMS as it allows for parametric evolutions to be captured when dealing with materials that depend on frequency or temperature [9] or when dealing with an optimization process or nonlinear simulation [10]. Finally, research is focusing on linear CMS model improvement. An important way to optimize the reduced model is mode selection, as it reduces the number of variables. Recent approaches are mainly based on energy criteria [11] or truncation error criteria [12]. Recently, Janssen et al. [13] proposed double criteria based on the frequency range of interest and the ability to capture the connection behavior. Our work focuses on this aspect and aims to compare several approaches dedicated to this objective. A key factor influencing computation time is the number of modes considered for modal superposition within each flexible body. Indeed, across all simulations conducted in Section 3.3 using various modes in the modal basis, we observed a quadratic relationship between the number of modes and the time required for time integration.

Traditionally, a Craig-Bampton superelement [4] for each flexible component is imported in FMBS software, consisting in mass and stiffness operators projected onto a basis comprising static modes (Guyan) and fixed-interface modes, also depicted as fixed-interface modes or constrained normal modes. Subsequently, a modal basis is internally recalculated based on the imposed boundary conditions in the model and truncated to the relevant maximal frequency of interest. Depending on the component's geometry and material properties, a significant number of modes may be retained, but only a few actively contribute to the system's response. Indeed, some of them may represent local mode shapes that are either not likely to be excited or contribute negligibly to the loads at the interface.

The local modes in the resulting Reduced-Order Model (ROM) arise from some particular internal Craig-Bampton dynamical modes preserved in the condensation. Various methods exist in the literature, aiming to filter them by means of participation factor computing. Notable methods include the following:

* Effective Interface Mass (EIM) and its derivatives VEIM and DEIM proposed by Kammer in [14];

* Participation Factors introduced by Lenoir, Cogan, and Lallement [15] (denoted LPF);

* Optimal Mode Ranking (OMR) [16];

* CMS$\chi$, where CMS stands for Component Mode Synthesis [17];

* Interior Mode Ranking (IMR) [18];

* Energy-Based Ranking (EBR) [11].

It can be noted that some of the filtering methods have been previously compared on simple academic test cases [19,20] highlighting IMR's effectiveness and OMR's poor performance. Additionally, EIM has been successfully applied to the crankshaft and the conrod of Multibody cranktrain model [21]. In the present work, various ranking techniques are compared using an industrial case study: the stator of a next-generation aircraft engine.
To bring it all together, the diagram in Figure 1 outlines the general process of constructing an FMBS model, with a focus on the modified superelement generation. In this process, the fixed-interface modes of the Craig-Bampton superelements are ranked, sorted, and truncated to prioritize the most significant modes in the dynamic response of the structure rather than low-frequency modes. The nonlinear phenomena present in the gearbox do indeed induce high-frequency behavior related to harmonics. Nevertheless, all the methods compared, as well as the reference Craig-Bampton superelement, use the same frequency truncation rule: five times the maximum rotor speed for calculating the fixed interface modes used to construct the superelement and then two and a half times the maximum speed during modal analysis on the resulting superelement as part of the transient calculation. This choice of cutoff frequency is debatable, but it also comes with computational constraints. This study focuses on the frequency range described above, in line with industrial applications and needs.

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Workflow for flexible multibody modeling with CB mode selection. $x_m$ and $x_s$ are respectively master and slave physical DOF, while $q$ represents modal DOF.</figcaption></figure>

In this paper, we aim to compare several of these methods on a realistic industrial case study, an aircraft engine model including a stator, two rotors, and a planetary gear train. The load case is an unbalance response. The evaluation criteria will focus on the computational time savings and the accuracy of forces and displacements at the bearings, which are also the interfaces of the superelements.

In this paper, we begin by recalling the basics of Craig-Bampton condensation and provide a concise overview of the different methodologies found in the literature to optimally select internal dynamical modes. They will be categorized into two groups: the a priori methods, which do not necessitate any prior knowledge of the crucial modes of the structure to be retained, and the a posteriori methods, which require prior knowledge of the important modes of the complete physical system to properly select the appropriate fixed-interface modes.

Subsequently, we assess a subset of these methods using an initial voluntary simplified industrial aircraft engine model, which involves a stator, two rotors, and a linearized planetary gear train. The evaluation aims to identify the most efficient approach based on criteria such as calculation time savings and the accuracy of load predictions across the entire structure.

Finally, we apply the most promising method to a nonlinear flexible multibody system featuring a comprehensive aircraft engine model. The assessment criteria remain centered on computational efficiency and the precision of load predictions throughout the structure.

## 2. Review of Craig-Bampton Mode Selection Techniques

The Craig-Bampton method (CB) stands out as one of the most extensively employed Component Mode Synthesis (CMS) techniques in structural dynamics. Introduced by Donald Craig and Clark Bampton [4], it facilitates the design and the dynamic analysis of structure assemblies.
## 2.1. Theoretical Background

Each component is modeled individually, and the Craig-Bampton method provides a reduced-order model where the interfaces with other components retain physical degrees of freedom.

First, let us recall the conservative dynamic equation:

$$M\ddot{U} + KU = F. \qquad (1)$$

Here, $M$ and $K$ are $n \times n$ matrices representing the mass and stiffness matrices, respectively; $F$ is the vector of external forces; and $U$ denotes the vector of degrees of freedom.

The $n$ degrees of freedom (DOFs) of the finite element structure are categorized into two groups: the “master” DOFs (designated with subscript $m$) and the “slave” DOFs (subscript $s$), such as $n = n_m + n_s$. Master nodes, also referred to as “boundary” nodes, are retained in the FMBS environment for applying forces, modeling connections, or positioning virtual sensors. Conversely, slave nodes, or “interior” nodes, are no longer in use. By reorganizing rows and columns of the finite element matrices, we can express Equation (1) as follows:

$$\begin{bmatrix} M_{mm} & M_{ms} \\ M_{sm} & M_{ss} \end{bmatrix} \begin{Bmatrix} \ddot{u}_m \\ \ddot{u}_s \end{Bmatrix} + \begin{bmatrix} K_{mm} & K_{ms} \\ K_{sm} & K_{ss} \end{bmatrix} \begin{Bmatrix} u_m \\ u_s \end{Bmatrix} = \begin{Bmatrix} R \\ 0 \end{Bmatrix}. \qquad (2)$$

$R$ are the reaction forces resulting from an interface displacement field $u_m$.

The next step involves computing static modes, commonly referred to as Guyan modes, denoted by $\Psi$ and determined by the following:

$$\Psi = -K_{ss}^{-1}K_{sm}. \qquad (3)$$

Then, considering all master nodes constrained, we compute the fixed-interface modes $\Phi$, by solving the eigenvalue problem:

$$K_{ss}\Phi = M_{ss}\Phi\Omega. \qquad (4)$$

$\Omega$ is the diagonal matrix of square eigenvalues. The set of fixed-interface modes $\Phi$ is reduced to a subset associated with natural frequencies lower than a cutoff frequency set by modeling assumptions, yielding $\tilde{\Phi}$; such that the modal properties of the resulting superelement and the full-order model (FOM) in the frequency band of interest for the application remain close. The modal basis used for the coordinates transformation is represented by the matrix $T$:

$$\begin{Bmatrix} u_m \\ u_s \end{Bmatrix} = \underbrace{\begin{bmatrix} I & 0 \\ \Psi & \tilde{\Phi} \end{bmatrix}}_{[T]} \begin{Bmatrix} u_m \\ q \end{Bmatrix} \quad \rightarrow \quad \begin{Bmatrix} u_m = u_m \\ u_s = \Psi u_m + \tilde{\Phi} q \end{Bmatrix}. \qquad (5)$$

$q$ represents the modal degrees of freedom. Finally, the projection of mass and stiffness operators into a reduced space yields the operators of the reduced order model:

$$M^{CB} = T^T MT, \quad K^{CB} = T^T KT, \quad F^{CB} = T^T F, \qquad (6)$$

## 2.2. A Priori Selection Methods

In this section, we will go through the a priori methods for ranking internal modes. Those techniques do not require any prior knowledge on the important modes of the structure.
On the one hand, the Effective Interface Mass (EIM) and Lenoir Participation Factors (LPFs) are fundamentally rooted in the concept of modal masses. The EIM method, however, distinguishes itself by additionally accounting for the influence of mass distribution at the interface.

On the other hand, the Optimal Mode Ranking (OMR) and CMSχ methods are grounded in the evaluation of strain energy within the structure. In essence, both approaches prioritize modes based on the amount of strain energy associated with each mode. Nevertheless, the primary discrepancy between them lies in the treatment of frequency in the ranking process.

For each method, we will express the modal participation factor denoted $S_i$.

### 2.2.1. Lenoir Participation Factors (LPFs)

In 2002, D. Lenoir, S. Cogan and G. Lallement [15] proposed to examine the relationship between the displacement at the interface $u_m$ and the corresponding reaction forces $R$ to define the excitability of the modes of the component with a clamped interface with respect to the inertial forces generated by the interface displacement. To do so, substituting Equation (5) in Equation (2) and premultiplying it by $H^T$ yields the following:

$$ \left( \left[ \begin{array}{cc} \bar{K}_{mm} & 0 \\ 0 & \Omega \end{array} \right] - \omega^2 \left[ \begin{array}{cc} \bar{M}_{mm} & L_{ms} \\ L_{sm} & I \end{array} \right] \right) \left\{ \begin{array}{c} u_m \\ q \end{array} \right\} = \left\{ \begin{array}{c} R \\ 0 \end{array} \right\}, \quad (7) $$

where the eigenmodes of the component clamped at the interface satisfy the following orthonormality relations $\Phi^T M_{ss} \Phi = I$ and $\Phi^T K_{ss} \Phi = \Omega$, and where

$$ \bar{K}_{mm} = K_{mm} + K_{ms} \Psi, $$

$$ \bar{M}_{mm} = M_{mm} + M_{ms} \Psi + \Psi^T M_{sm} + \Psi^T M_{ss} \Psi, $$

$$ L_{ms} = M_{ms} \Phi + \Psi^T M_{ss} \Phi. $$

In this method, the interface mass matrices $M_{mm}$, $M_{ms}$ and $M_{sm}$ are neglected, resulting in

$$ \bar{M}_{mm} = \Psi^T M_{ss} \Psi, $$

$$ L_{ms} = \Psi^T M_{ss} \Phi. $$

Then, we can obtain the relation between the reaction forces and the displacements at the interfaces:

$$ R = \left( \bar{K}_{mm} - \omega^2 \left[ \bar{M}_{mm} + \sum_i \frac{1}{\left(\frac{\omega_i}{\omega}\right)^2 - 1} L_{mi} L_{im} \right] \right) u_m, \quad (8) $$

where $L_{mi} = \Psi^T M_{ss} \Phi_i$ is the $i^{th}$ internal mode participation factor. This relation distinguishes the static component $\bar{M}_{mm}$ and the dynamic component of the reduced interface mass matrix. The product $L_{ms}L_{sm}$ represents the classical modal effective mass matrix.

Finally,

$$ S_i^{LPF} = L_{mi} = \Psi^T M_{ss} \Phi_i. \quad (9) $$

### 2.2.2. Effective Interface Mass (EIM)

Developed by Kammer and Triller in 1996, the Effective Interface Mass (EIM) factor determines the contribution of each constrained normal mode to the dynamic loads at the interface under general displacements of the interior. It is similar to the previous method but does not neglect the interface mass matrices.
The modal participation factor matrix is computed as follows [14]:

$$P = -\Phi^T (M_{ss} \Psi + M_{sm}). \qquad (10)$$

The $i^{th}$ row $P_i$ contains the modal participation factors for the $i^{th}$ fixed-interface mode. They represent multiplication factors for the acceleration inputs from the interface degrees of freedom numbered with the subscript $j$. The larger the $j^{th}$ entry in $P_i$, the more the $i^{th}$ mode will be excited by the $j^{th}$ input. Finally, the EIM factor for internal mode $i$ is obtained through the following equations:

$$\bar{M} = \sum_{i=1}^{n_o} \bar{M}_i = \sum_{i=1}^{n_o} P_i^T P_i = P^T P, \qquad (11)$$

$$S_i^{EIM} = \frac{\mathrm{tr}(\bar{M}_i)}{\mathrm{tr}(\bar{M})}. \qquad (12)$$

The authors extended EIM to consider modal velocity VEIM (13) and modal displacement outputs DEIM (14), weighting low-frequency modes more heavily.

$$S_i^{VEIM} = \frac{\mathrm{tr}(M_{\Omega^{-1}})}{\mathrm{tr}(M_{\Omega^{-1}})}, \qquad (13)$$

$$S_i^{DEIM} = \frac{\mathrm{tr}(M_{\Omega^{-2}})}{\mathrm{tr}(M_{\Omega^{-2}})}, \qquad (14)$$

where $M_{\Omega^{-n}} = P^T \Omega^{-n} P$.

Instead of being ranked by increasing frequency as is the case in the classic Craig-Bampton framework, modes are ranked by decreasing EIM factor. The number of dynamic modes retained in the modal basis is determined by a mass-based criterion. The recommended approach is to retain at least 90 percent of the system's mass. This involves keeping modes with higher EIM values until the cumulative sum reaches 0.9.

### 2.2.3. Optimal Mode Ranking (OMR)

Developed in 2004 by Givoli et al. [16], the participation factor OMR for ranking the constrained normal modes, with eigenvalue $\omega_i$ and eigenvector $\phi_i$, is obtained with the formula below:

$$S_i^{OMR} = \omega_i M_{sm}^T \phi_i \phi_i^T M_{sm} - M_{sm}^T \phi_i \phi_i^T K_{sm} - K_{sm}^T \phi_i \phi_i^T M_{sm} + \frac{1}{\omega_i} K_{sm}^T \phi_i \phi_i^T K_{sm}. \qquad (15)$$

This method had not been tested because it showed mitigated performances on academic cases compared to other methods, probably because not enough weight is given to the frequency of the internal modes as compared to other methods (see Equation (17)).

### 2.2.4. CMS "khi"

Developed in 2007 by Liao et al. [17] as a declination of OMR, the participation factor CMS$\chi$ for ranking the constrained normal modes $\phi_i$ is obtained with the formula below:

$$S_i^{CMS\chi} = \frac{1}{\omega_i} (M_{sm} - M_{ss} \Psi)^T \phi_i \phi_i^T (M_{sm} - M_{ss} \Psi). \qquad (16)$$

A relation between some participation factors is demonstrated in [20] between Kammer's methods, OMR, and CMS$\chi$:
$$S_i^{\text{EIM}} = \omega_i^2 S_i^{\text{VEIM}} = \omega_i^4 S_i^{\text{DEIM}} = \omega_i^2 S_i^{\text{CMS}\chi} = \omega_i^{-2} S_i^{\text{OMR}}, \quad (17)$$

where $\omega_i$ is the frequency of $\phi_i$. It is shown that CMS$\chi$ is equivalent to VEIM while being more computationally expensive. As a consequence, it will not be tested.

## 2.3. A Posteriori Selection Methods

In this section, we will go through the a posteriori methods for ranking internal modes. In contrast to a priori methods, they require prior knowledge of the important modes of the structure or excitation in order to choose the appropriate fixed-interface modes to keep for the condensation, ensuring good representability of the physical model by the superelement for the application case. But, ultimately, we want to perform the selection without initial calculations. For the a posteriori methods, we therefore opted for a definition of force that contains all possible components. This approach is clearly disadvantageous for a posteriori methods but corresponds to the reality of the intended use.

### 2.3.1. Energy-Based Ranking (EBR)

The mode selection criterion proposed in [11] ranks the interior modes by using scalar coefficients representing the contribution of each interior mode to the mean kinetic and potential elastic energy stored by the complete system in a period of the external force. These coefficients provide a measure of the importance of each interior mode to the forced response of the full-order system and account explicitly for both the frequencies and the spatial distribution of the force.

The system is supposed to be excited on the master DOFs by a set of periodic external nodal forces $f_e$ consisting of a sum of $n_f$ harmonic components $f_k$:

$$f(t) = \left\{ \begin{matrix} f_e(t) \\ 0 \end{matrix} \right\} = \sum_{k=1}^{n_f} f_k(t) = \sum_{k=1}^{n_f} \left\{ \begin{matrix} f_{1,k} \cos(\omega_k t + \alpha_{1,k}) \\ \vdots \\ f_{n_{m,k}} \cos(\omega_k t + \alpha_{n_{m,k}}) \\ 0_{n_{m+1,k}} \\ \vdots \\ 0_{n,k} \end{matrix} \right\}, \quad (18)$$

where $f_k$ is the $k^{th}$ harmonic component of the periodic force and $f_{i,k}$, $\omega_{i,k}$, and $\alpha_{i,k}$ are respectively, the amplitude, the angular frequency, and the relative phase of the $k^{th}$ harmonic component exciting the $i^{th}$ DOF.

The contribution of the $i^{th}$ interior mode to the system mean energy can be evaluated through the scalar coefficients:

$$S_i^{\text{EBR}} = \sum_{k=1}^{n_f} \left| \frac{1}{2} (\omega_i^2 + \omega_k^2) (G_{\eta_{i,k}} f_k)^2 + \omega_k^2 |f_k^T G_{e_k}^T| D_{\beta_{e,k} \gamma_{i,k}} M_{ms} \phi_i |G_{\eta_{i,k}} f_k| \right|, \quad (19)$$

with $G_{e_k}$ and $G_{\eta_k}$ arising from the partition of the contributions of $u_m$ and $u_s$ in

$$[K^{\text{CB}} - \omega_k^2 M^{\text{CB}}]^{-1} = \begin{bmatrix} G_{e_k} \\ G_{\eta_k} \end{bmatrix}, \quad (20)$$

and with $D_{\beta_{e,k}\gamma_{i,k}} = \text{diag}(\cos(\beta_{i,k} - \gamma_{e,k}))$, $i = 1, \dots, s$, $e = 1, \dots, m$, where $\beta_{i,k}$ is the relative phase of the response of the $i^{th}$ DOF to the $k^{th}$ harmonic component and $\gamma_{j,k}$ represents the relative phase of the response of the $j^{th}$ interior modal coordinate to the $k^{th}$ harmonic component.
This method requires knowing the external force to be applied before generating the reduced-order model.

### 2.3.2. Interior Mode Ranking (IMR)

The mode selection criterion proposed in [18], referred to as the Interior Mode Ranking (IMR) method, aims to rank the constrained normal modes according to the contribution they provide to the dynamics of one or more selected vibration modes of the full system with actual boundary conditions. The eigenvector matrix of the full system is denoted V and has the size $n \times n$.

The participation coefficients of the interior modes i in the dynamics of the $\zeta^{th}$ normal mode of the full order system is as follows:

$$S_{\zeta,i}^{IMR} = \Gamma_{2\zeta,i} \left( \frac{1}{\left( \frac{\omega_i}{\omega_{\zeta}} \right)^2 - 1} \right) m^{\mathrm{T}}(j\omega_f) \phi_i \quad i = 1, \dots, n_s, \qquad (21)$$

where $\phi_i$ is the $i^{th}$ fixed interface normal mode shape, $\Gamma_{2\zeta,i}$ arises from the partitioning of the transformation matrix $V^T M H$ into $\Gamma_1$ (size $n \times n_m$) and $\Gamma_2$ (size $n \times n_s$) that splits the contributions of both $u_m$ and $u_s$, and the frequency of the fictitious force $\omega_f$ is assumed equal to the angular frequency of the mode of interest $\omega_{\zeta}$ plus a small frequency shift and $m(j\omega_f)$ is defined by the following:

$$m(j\omega_f) = \left\{ \frac{\omega_{\zeta}^2 M_{sm} + \omega_f^2 M_{ss} \Psi}{\omega_f^2 (\omega_{\zeta}^2 - \omega_f^2)} v_{\zeta,m} + \frac{M_{ss}}{\omega_f^2} v_{\zeta,s} \right\}, \qquad (22)$$

where $v_{\zeta,m}$ and $v_{\zeta,s}$ arise from the partitioning on master and slave DOFs of the modal shape of the $\zeta^{th}$ mode of the full system.

This method requires knowing the modes of the full system that will be excited.

## 3. Improved Performances of Flexible Multibody Analysis

### 3.1. Presentation of the Industrial Application Cases

To reduce the duration of flexible multibody simulations, we have identified and reviewed several methods to improve the construction of the Craig-Bampton superelement. For comparative analysis, we will now apply these methods to a stator component v1 introduced into a deliberately simplified flexible multibody model featuring two rotors and a linearized planetary gearbox (Model A in Figure 2). The Finite Element (FE) model of this particular stator component exhibits numerous local elastic modes, facilitating the distinction of the efficiency of various methods.

Subsequently, after identifying the most promising method, we apply it to another stator component v2, consisting of an improved version of stator v1, completed with a nacelle and a pylon. Stator v2 is integrated into a complete turbomachine model consisting of three rotors and a fully detailed nonlinear gear train (Model B in Figure 2). The FE model of stator v2 is of better quality than v1 and exhibits less local elastic modes. Concerning the planetary gear train, each planet, the planet carrier, the sun gear, and the ring have flexible bodies. Detailed gear force elements are employed, accounting for the three-dimensional motion of components and meshing effects. Internal mode filtering is exclusively applied to the statoric parts because the rotors are modeled as beams, where all nodes are considered as master nodes, allowing no fixed-interface modes. Regarding the gear train, due to cyclic symmetry, we assume that there are no local deformations except for the teeth deflection, which we aim to preserve.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Industrial cases for Craig–Bampton internal mode selection techniques comparison (A) and final real-world application (B). The orange and green arrows indicate, respectively, the applied speed and imposed torque on the shaft. The black comb symbol shows that the system is rigidly linked to the ground.</figcaption></figure>

To classify the methods, the test case involves transient analysis of an unbalance response. A mass is attached to the low-pressure turbine, and the angular velocity of the rotors is linearly imposed from 0 to the maximum speed, while a resistive torque is applied in accordance with real-world operating conditions. The features under comparison include the evolution of loads in the bearings and suspensions of the engine with angular velocity.

## 3.2. *Practical Implementation and Results*

The full physical operators of the components are initially extracted as op4 files from the MSC Nastran (v2018) Finite Element model using a dedicated DMAP script. These files are processed using MATLAB (v2019b) routines to obtain them as MATLAB files. We then conduct a reordering of rows and columns to gather master nodes on the one hand and the slave nodes on the other hand. We also remove non-independent nodes from the model (attached with MPC—Multiple Point Constraint—or RBE—Rigid Body Element—for example).

Both stator versions have 13 master nodes, which means 78 master DOFs, resulting in the same amount of static modes, which are computed according to (3). Constrained normal modes are computed, ranked, and sorted by decreasing participation factors for each method. Considering that the stator (v1 or v2) features over 150,000 degrees of freedom, and due to limited computational resources, the computation of all fixed-interface modes is hardly feasible. Consequently, our approach involves calculating these modes up to five times the maximum frequency of interest. Modes beyond this limit are unlikely to exhibit significant participation factors.

Finally, the ranked and sorted modes are truncated based on criteria defined in the literature or through trial and error to evaluate performance in the industrial application case. Full mass, stiffness, and damping operators are projected onto a modified Craig–Bampton basis, retaining the most significant modes according to each ranking method. The reduced-order matrices are written in a binary file for compatibility with the multibody software and imported into the whole-engine model.
As an illustrative example, let us focus on the application of Effective Interface Mass (EIM). Figure 3 compares the cumulative sum of EIM values for modes sorted by increasing frequency, as in the classic Craig–Bampton method, against modes sorted by decreasing EIM. This comparison reveals that many low-frequency modes in the structure exert minimal influence on the system’s response at interface and can be safely excluded from the projection modal basis.

<figure><img src="image_6.png" /><figcaption>Figure 3. Cumulative sum of effective interface mass participation factor.</figcaption></figure>

Kammer and Triller recommend truncation to ensure a minimum of 90% effective interface mass. In a conservative approach, maintaining 95% of Effective Interface Mass (EIM) would necessitate selecting approximately 1400 modes when sorted by increasing frequency. However, when ordered by decreasing EIM participation factor, only 500 modes would be enough. In the latter case, the resulting superelements would have fewer modal DOFs, with the number of physical DOFs remaining the same.

Figure 4 shows the Modal Assurance Criterion (MAC; see (23)) matrix between the modes of the EIM superelement and the full physical model. To make the figure easier to read, the matrix is reduced to the first 13 modes. They are representative of the entire matrix. While the low-frequency modes of a classic Craig–Bampton superelement correspond to those of the full-order model and would show a diagonal with ones, here, certain modes (highlighted in red boxes) present in the FOM are absent in the EIM superelement. These modes correspond to local elastic modes, while global modes are well-represented (MAC > 0.95 attesting to the resemblance of modal shapes (green boxes)).

$$MAC(\{\phi_s\}, \{\phi_r\}) = \frac{|\{\phi_r\}^T \{\phi_s\}|^2}{(\{\phi_r\}^T \{\phi_s\})(\{\phi_s\}^T \{\phi_r\})} \quad (23)$$

We have successfully reduced the number of modes of the superelement in the relevant frequency range from 494 to 340. This reduction significantly cuts down the computational time required for the multibody simulation, all while maintaining a very good accuracy of the results. The excluded modes, having negligible contributions to the system’s response, do not compromise the overall outcomes.

Furthermore, it is worth discussing the selection criteria for EIM. In Figure 5, only MAC values greater than 0.95 are plotted, comparing the mode shapes of resulting EIM superelements with selection criteria 90%, 95% and 610 modes to those from the full Finite Element (FE) model modal analysis. The short dashed horizontal lines indicate the number of degrees of freedom of the substructure, while the long dashed horizontal lines represent
the number of internal modes. The vertical line indicates the number of modes to be used
in the physical model, for modal superposition in the frequency range of interest.

<figure><img src="image_7.png" /><figcaption>Figure 4. Modal assurance criterion matrix of the stator superelement obtained with EIM ranking vs. the full physical model.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 5. Impact of the chosen criterion for EIM in terms of internal modes and degrees of freedom.</figcaption></figure>

As expected, the less restrictive the selection criteria, the more internal modes are kept,
the better the correspondence between the FOM and the ROM within the studied frequency
band. However, this also leads to reduced potential for decreasing computation time.
Therefore, the ultimate decision entails finding a balance between the desired prediction
accuracy and the need for saving computation time.

3.3. *Performance Benchmarking of Selection Techniques*

For the purpose of comparison, a reference model was used to evaluate both compu-
tational time savings and the loss of precision in bearing load estimations. This reference
model employed the Craig-Bampton method with a truncation frequency set to five times
the maximum frequency of interest, resulting in 1650 dynamic modes within the superele-
ments, in addition to the 78 static modes. This approach yielded a total of 494 modes
considered after integration into the MBS model and truncation to the frequency of interest.
The time integration for this reference model required approximately 20 h.

The process described in the previous subsection is replicated for most a priori meth-
ods, employing various truncation criteria. The a posteriori method IRM was omitted from
testing, as it is deemed unsuitable for our case. To determine the appropriate modes from
the full physical model, one would need to conduct an expensive simulation and invest
effort in analyzing modal participations. Given the presence of several hundred modes, identifying the relevant ones becomes a challenging task. The other a posteriori method, referred to as EBR, was tested by selecting arbitrary excitation forces. This approach was necessitated by the lack of prior knowledge about the locations of loads without running the simulation. Consequently, we are not confident in the method's relevance in this case, as the excitation forces are chosen without prior information.

For each mode selection method and chosen selection criterion, the number of retained dynamic modes for condensation is specified in Table 1, along with the difference between the resulting number of modes used for modal superposition in the multibody simulation and in the reference case. The differences between the methods rely in the choice of internal modes and its intricate link to the resulting number of modes in modal basis of the ROM. For instance, all methods were tested while retaining the first 292 internal modes sorted by decreasing participation factors. Notably, the number of resulting statoric modes used for time integration varies. These discrepancies can be further investigated in Figure 6, where MAC > 0.98 are plotted between FOM and ROM for the various methods, highlighting their distinct selectivity tendencies. We can, for example, confirm the statement that VEIM assigns more significance to low-frequency modes than EIM. Another observation is that all methods filtering low-frequency modes seem to target approximately the same FOM modes for elimination.

<table><caption>**Table 1.** Summary table of computation time savings and precision on load estimation in bearings among various internal mode selection methods for a benchmark on model A. NIDM: Number of fixed-interface modes. Diff: Difference in nb of statoric modes compared to ref. Time: Normalized computation time (ref. = 1).</caption><thead><tr><th rowspan="2">Method</th><th rowspan="2">N</th><th rowspan="2">Diff</th><th rowspan="2">Mass (%)</th><th rowspan="2">Time</th><th colspan="3">Max. Rel. Error on Loads at Bearing (%)</th><th colspan="3">Max. Rel. Error on Displacement at Bearings (%)</th><th rowspan="2">Score (/20)</th></tr><tr><th>Worst</th><th>Best</th><th>Mean</th><th>Worst</th><th>Best</th><th>Mean</th></tr></thead><tbody><tr><td rowspan="2">CB</td><td>610</td><td>0</td><td></td><td>1.01</td><td>4.2</td><td>0.2</td><td>1.1</td><td>2</td><td>1</td><td>1</td><td>8.9</td></tr><tr><td>292</td><td>-187</td><td></td><td>0.45</td><td>26.4</td><td>3.5</td><td>10.2</td><td>11</td><td>3</td><td>7</td><td>5.5</td></tr><tr><td rowspan="3">LPF</td><td>610</td><td>-112</td><td>98.3</td><td>0.66</td><td>4.4</td><td>0.2</td><td>1.5</td><td>3</td><td>1</td><td>2</td><td>11.9</td></tr><tr><td>327</td><td>-231</td><td>95.0</td><td>0.37</td><td>5.3</td><td>1.1</td><td>2.6</td><td>2</td><td>1</td><td>2</td><td>13.7</td></tr><tr><td>292</td><td>-257</td><td>94.2</td><td>0.30</td><td>5.8</td><td>1.3</td><td>3.3</td><td>4</td><td>1</td><td>3</td><td>13.7</td></tr><tr><td rowspan="3">EIM</td><td>610</td><td>-105</td><td>96.7</td><td>0.63</td><td>3.9</td><td>0.3</td><td>1.2</td><td>2</td><td>1</td><td>1</td><td>12.5</td></tr><tr><td>485</td><td>-154</td><td>95.0</td><td>0.51</td><td>4.4</td><td>0.5</td><td>1.7</td><td>3</td><td>1</td><td>2</td><td>13.2</td></tr><tr><td>292</td><td>-245</td><td>90.0</td><td>0.33</td><td>13.0</td><td>1.9</td><td>5.0</td><td>5</td><td>2</td><td>3</td><td>11.7</td></tr><tr><td rowspan="2">VEIM</td><td>292</td><td>-198</td><td>99.7</td><td>0.51</td><td>5.4</td><td>0.5</td><td>1.8</td><td>2</td><td>1</td><td>2</td><td>13.1</td></tr><tr><td>100</td><td>-375</td><td>95.0</td><td>0.12</td><td>45.9</td><td>9.4</td><td>24.9</td><td>11</td><td>5</td><td>8</td><td>8.8</td></tr><tr><td rowspan="2">EBR</td><td>338</td><td>-167</td><td>99.0</td><td>0.58</td><td>9.0</td><td>1.6</td><td>4.0</td><td>3</td><td>2</td><td>2</td><td>10.2</td></tr><tr><td>292</td><td>-209</td><td>98.4</td><td>0.45</td><td>22.0</td><td>2.1</td><td>9.3</td><td>9</td><td>2</td><td>6</td><td>6.2</td></tr></tbody></table>

The percentage of time saved compared to the reference model is also provided. The relation between time savings and the amount of statoric modes in the frequency range for analysis is highlighted in Figure 7. Indeed, the calculation time depends linearly on the number of degrees of freedom *N* for solving second-order differential systems, but the algebraic-differential equations associated with mechanical connections introduce a “complexity” that leads to solution times proportional to the cube of the number of degrees of freedom. Recent work [22] shows that it is possible to reduce the time required to solve constraint equations without, however, achieving proportionality to *N*. In our case, proportionality to *N*² has been observed empirically.
Furthermore, we compute the maximum relative error $\epsilon$ observed during the complete ramp-up phase using Equation (24), where $x$ represents the observed value and $x_{ref}$ denotes the reference value:

$$ \epsilon = \max_{simu} \left( \frac{abs(x - x_{ref})}{\max(x_{ref}) - \min(x_{ref})} \right). \qquad (24) $$

Out of the eight load features, $\epsilon$ is given in the table for the worst load feature, as well as the best and mean out of the eight.

Finally, a score $S = S_1 + S_2$ out of 20 is calculated based on the accuracy and time saving according to the following rules:

* if $\epsilon_{mean} < 10\%$, $S_1 = 10 - \epsilon_{mean}$ else, $S_1 = 0$

* 0/10 if $t < t_{ref}$, $S_2 = 10 \frac{t}{t_{ref}}$ else $S_2 = 0$

The best marks are obtained for the EIM and LPF methods. The latter is applied in the next section on the statoric parts of the next generation of civil aircraft engines CFM RISE for industrial validation. The EIM method is indeed equivalent to the LPF method, which does not neglect the mass terms at the interface. The similar results obtained using these two methods show that the mass at the interface has little effect in the case of the stators studied. This conclusion can be generalized to any structure with distributed mass and relatively small interfaces.

<figure><img src="image_9.png" /><figcaption>Figure 6. Modal Assurance Criterion (MAC) between the Full Order Model (FOM) and various Reduced Order Models (ROMs) using different mode selection techniques. Only MAC values greater than 0.98, which attest to a good correspondence between modal shapes, are displayed, allowing for a comparison of the selectivity of the different methods.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 7. Influence of the number of modes in the modal basis on the simulation duration.</figcaption></figure>
### 3.4. Validation

The most promising methods for fixed-interface mode selection had been identified thanks to a benchmark based on the intentionally simplified yet relevant model A (cf. Figure 2).

Lenoir's method with selection criterion set to 95% of the mass is applied to model B, which represents the most detailed FMBS model of RISE to date, including notably the nonlinear effects of the planetary gear box. The gain in computational time was 17.6% (22.25 h instead of 27 h), which is particularly noteworthy considering the minimal additional effort required. Indeed, it only requires adding a single line of code to call a DMAP script implementing the LPF method to the Nastran file used for superelement construction. Moreover, it is worth noting that the performed improvement pertains to only a small portion of the DOFs of the entire system, exhibiting linear behavior. As for the accuracy of prediction, it was excellent for the application discrepancies inferior to 1%, observed on all the load features while performing unbalance response on the low-pressure turbine. The explanation lies in the fact that much of the high-frequency behavior is captured by static modes (Guyan). Since the nonlinearities are not internal to the superelement, the latter filters the excitations that reach it. Guyan modes ensure the transmission of forces from one bearing to another, regardless of the choice of frequency truncation. The results obtained on model B show, through their fidelity to the reference simulation, that there is no need to question the selection method when nonlinearity occurs outside the reduced model. Indeed, the dominant frequency content in the response (in terms of amplitude) is clearly present in the range.

## 4. Conclusions and Future Work

Flexible Multibody Simulations are efficient for predicting loads in complex nonlinear structures in operation. However, despite the ongoing advancements in computational resources, FMBS can impose a substantial computational burden. Through our investigation, we have observed that the computational time required for temporal simulation of a simplified aircraft engine model exhibits a quadratic relationship with the number of retained modes of the stator's flexible body for modal superposition. Consequently, reducing the number of superelement modes within the frequency range of interest is essential to reduce computational time, particularly when the structure contains local elastic modes that have minimal significance in the overall system response.

In this study, we have adapted the classical Craig-Bampton reduction methods by incorporating various existing techniques from the literature, which offer a refined selection of dynamic modes. To assess the efficiency of these methods, we conducted a comparative analysis based on results obtained from a speed ramp-up with unbalanced response in a realistic aircraft engine model. This evaluation primarily considered time-saving benefits and the accuracy of load predictions.

Our findings highlight two standout methods: the Effective Interface Mass method, as proposed by Kammer [14], and a method proposed by Lenoir and Cogan in [15]. These methods have demonstrated the ability to achieve significant computational time reductions, up to 65 percent, while introducing a maximal relative error of 4 percent in load predictions for all bearings during the ramp-up phase, which is satisfactory for the application.

In the scope of this work and with regard to the time available, the methods OMR and CMSχ have not been employed due to their recognized lower efficiency on academic test scenarios or similarities with other tested methods. Likewise, the exploration of the IMR method has been omitted due to its restricted relevance to the specific application as it necessitates prior knowledge of the modes to retain from the full-order system in order to select the most suitable fixed-interface modes for Craig-Bampton's condensation.
The list of tested methods is not exhaustive and only concerns improvements of the Craig-Bampton reduction method, which is most widely used in commercial multibody software. Nevertheless, there are alternative approaches, such as moment-matching by projection on Krylov subspaces or SVD-based reduction techniques to take into account the time-varying boundary conditions in the reduction step [23–25]. Those are worth including in this benchmark. However, due to considerations regarding integration feasibility with existing industrial methods, time constraints, and the already satisfactory results obtained, these methods were not explored as part of this study.

In the future, it is certainly the use of superelements constructed from data, so-called “surrogate models” [26], that could make it possible to select fewer variables to simulate, further reduce calculation times, and improve accuracy.

**Author Contributions:** Conceptualization, S.C., G.C. and C.O.; software, O.T.; validation, O.T. and C.O.; formal analysis, O.T.; investigation, O.T.; writing—original draft preparation, O.T.; writing—review and editing, S.C., G.C. and C.O.; supervision, S.C., G.C. and C.O. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Safran Aircraft Engines as part of a CIFRE agreement with the FEMTO-ST Institute.

**Data Availability Statement:** The datasets presented in this article are not available due to industrial confidentiality.

**Conflicts of Interest:** The authors declare no conflicts of interest.

## Abbreviations

The following abbreviations are used in this manuscript:

FMBS Flexible Multibody Simulation

FEM Finite Element Method

MAC Modal Assurance Criterion

DOF Degree Of Freedom

CMS Component Mode Synthesis

CB Craig-Bampton

EIM Effective Interface Mass

LPF Lenoir Participation Factor

OMR Optimal Mode Ranking

IMR Interior Mode Ranking

EBR Energy-Based Ranking

FOM Full-Order Model

ROM Reduced-Order Model

DMAP Direct Matrix Abstraction Program

## References

1. Pappalardo, C.M.; Del Giudice, M.; Oliva, E.B.; Stieven, L.; Naddeo, A. Computer-aided design, multibody dynamic modeling, and motion control analysis of a quadcopter system for delivery applications. *Machines* **2023**, *11*, 464. [CrossRef]

2. Shabana, A.A. Flexible multibody dynamics: Review of past and recent developments. *Multibody Syst. Dyn.* **1997**, *1*, 189–222. [CrossRef]

3. Wang, G.; Niu, Z.; Feng, Y. Improved Craig-Bampton Method Implemented into Durability Analysis of Flexible Multibody Systems. *Actuators* **2023**, *12*, 65. [CrossRef]

4. Craig, R.R., Jr.; Bampton, M.C. Coupling of substructures for dynamic analyses. *AIAA J.* **1968**, *6*, 1313–1319. [CrossRef]

5. Mashayekhi, F.; Zucca, S. Modal derivatives for efficient vibration prediction of geometrically nonlinear structures with friction contact. *Appl. Sci.* **2024**, *14*, 3936. [CrossRef]

6. Bui, T.A.; Park, J.; Kim, J.S. A review of combining component mode synthesis and model order reductions for geometrically nonlinear analysis. *J. Mech. Sci. Technol.* **2024**, *38*, 4699–4711. [CrossRef]
7. Penas, R.; Balmes, E.; Gaudin, A. A unified non-linear system model view of hyperelasticity, viscoelasticity and hysteresis exhibited by rubber. *Mech. Syst. Signal Process.* **2022**, *170*, 108793. [CrossRef]

8. Festjens, H.; Chevallier, G.; Dion, J. Nonlinear model order reduction of jointed structures for dynamic analysis. *J. Sound Vib.* **2014**, *333*, 2100–2113. [CrossRef]

9. Festjens, H.; Gaël, C.; Franck, R.; Jean-Luc, D.; Remy, L. Effectiveness of multilayer viscoelastic insulators to prevent occurrences of brake squeal: A numerical study. *Appl. Acoust.* **2012**, *73*, 1121–1128. [CrossRef]

10. Vlachas, K.; Garland, A.; Quinn, D.D.; Chatzi, E. Parametric reduced-order modeling for component-oriented treatment and localized nonlinear feature inclusion. *Nonlinear Dyn.* **2024**, *112*, 3399–3420. [CrossRef]

11. Palomba, I.; Richiedei, D.; Trevisani, A. Energy-based optimal ranking of the interior modes for reduced-order models under periodic excitation. *Shock Vib.* **2015**, *2015*, 348106. [CrossRef]

12. Yuan, J.; Salles, L.; El Haddad, F.; Wong, C. An adaptive component mode synthesis method for dynamic analysis of jointed structure with contact friction interfaces. *Comput. Struct.* **2020**, *229*, 106177. [CrossRef]

13. Janssen, L.A.; Fey, R.H.; Besselink, B.; van de Wouw, N. Mode selection for component mode synthesis with guaranteed assembly accuracy. *J. Sound Vib.* **2024**, *589*, 118596. [CrossRef]

14. Kammer, D.C.; Triller, M.J. Selection of Component Modes for Craig-Bampton Substructure Representations. *J. Vib. Acoust.* **1996**, *118*, 264–270. [CrossRef]

15. Lenoir, D.; Cogan, S. Analyse d'incertitudes structurales complexes par sous-structuration. In Proceedings of the XIIIème Colloque GAMI Vibrations, Chocs et Bruits, Lyon, France, 12-14 June 2002.

16. Givoli, D.; Barbone, P.E.; Patlashenko, I. Which are the important modes of a subsystem? *Int. J. Numer. Methods Eng.* **2004**, *59*, 1657–1678. [CrossRef]

17. Liao, B.S.; Bai, Z.; Gao, W. The important modes of subsystems: A moment-matching approach. *Int. J. Numer. Methods Eng.* **2007**, *70*, 1581–1597. [CrossRef]

18. Palomba, I.; Richiedei, D.; Trevisani, A. Mode selection for reduced order modeling of mechanical systems excited at resonance. *Int. J. Mech. Sci.* **2016**, *114*, 268–276. [CrossRef]

19. Palomba, I.; Richiedei, D.; Trevisani, A. A ranking method for the selection of the interior modes of reduced order resonant system models. In *Proceedings of the Engineering Systems Design and Analysis; American Society of Mechanical Engineer: New York, NY, USA, 2014; Volume 45844*, p. V002T07A028.

20. Kim, S.M.; Kim, J.G.; Chae, S.W.; Park, K. Evaluating mode selection methods for component mode synthesis. *AIAA J.* **2016**, *54*, 2852–2863. [CrossRef]

21. Ricci, S.; Troncossi, M.; Rivola, A. Modal Selection Through Effective Interface Mass with Application to Flexible Multibody Cranktrain Dynamics. *J. Comput. Nonlinear Dyn.* **2014**, *9*, 011002. [CrossRef]

22. Featherstone, R. *Rigid Body Dynamics Algorithms*; Springer: Berlin/Heidelberg, Germany, 2008.

23. Fehr, J.; Eberhard, P. Simulation process of flexible multibody systems with non-modal model order reduction techniques. *Multibody Syst. Dyn.* **2011**, *25*, 313–334. [CrossRef]

24. Holzwarth, P.; Eberhard, P. SVD-based improvements for component mode synthesis in elastic multibody systems. *Eur. J. Mech.-A/Solids* **2015**, *49*, 408–418. [CrossRef]

25. Frie, L.; Dieterich, O.; Eberhard, P. Model Order Reduction for Elastic Multibody Systems with Fast Rotating Flexible Bodies. In *Proceedings of the ECCOMAS Thematic Conference on Multibody Dynamics, Jyväskylä, Finland, 12–15 December 2021; Budapest University of Technology and Economics: Budapest, Hungary, 2021; pp. 432–443*.

26. Boschetti, G.; Sinico, T. Designing digital twins of robots using simscape multibody. *Robotics* **2024**, *13*, 62. [CrossRef]

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.