# Reduction of cyclic symmetric mechanical systems subject to unilateral contact: Application to periodic blade-tip/casing contact interactions

Thibaut Vadcard, Samuel Quaegebeur, Fabrice Thouverez

► To cite this version:

Thibaut Vadcard, Samuel Quaegebeur, Fabrice Thouverez. Reduction of cyclic symmetric mechanical systems subject to unilateral contact: Application to periodic blade-tip/casing contact interactions. Mechanical Systems and Signal Processing, 2025, 237, pp.112974. 10.1016/j.ymssp.2025.112974. hal-05183941

HAL Id: hal-05183941

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

Submitted on 24 Jul 2025

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

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

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

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

Full length article

# Reduction of cyclic symmetric mechanical systems subject to unilateral contact: Application to periodic blade-tip/casing contact interactions

Thibaut Vadcard<sup>1, *</sup>, Samuel Quaegebeur, Fabrice Thouverez

École Centrale de Lyon, CNRS, ENTPE, LTDS, UMR5513, 69130 Écully, France

## ARTICLE INFO

Communicated by J. Baqersad

**Keywords:**

Harmonic balance method

Unilateral contact

Blade-tip/casing interactions

Nonlinear dynamics

Computational mechanics

Cyclic symmetry

Reduction method

## ABSTRACT

This paper proposes a novel reduction technique for industrial bladed disks featuring blade-tip/casing contacts. By analyzing a regularized contact law, a criterion allowing to drastically reduce the size of the problem is derived. It consists in restricting solutions to only a subset of nodal diameters identified in advance. An additional parity rule allowing to reduce the size of problem by half in certain configurations when solving through the harmonic balance method (HBM). These rules are tested on an academic cyclic model featuring unilateral contacts. The reduction is exact both for a regularized law but also for a nonsmooth contact law. The results match perfectly with the reference calculations led on the full system. All the previously mentioned observations remain true for the simulations led on the industrial fan stage ECL5/Catana under rubbing interactions. After including friction forces, the reduction technique is still exact for the industrial model. The previously published traveling wave hypothesis (deducing the full response from a single sector) is shown to perform well for purely traveling excitations. Several loading scenarios demonstrate the applicability of the strategy in numerous configurations. Localized responses are accurately captured by the reduction method in the case of excitation through perturbed traveling waves. This contribution allows to take a major step forward in the numerical modeling of bladed disks featuring blade-tip/casing interactions by taking advantage of cyclic symmetry. The proposed strategy yields an average reduction of 37.6% for the HBM and 25.2% for other methods. For traveling wave excitations, the reduction is rapidly over 90% and goes up as the number of sectors increases.

## 1. Introduction

In order to comply with international regulations that are aiming towards the achievement of net carbon neutrality [1], aircraft engine manufacturers are committed to design continuously more efficient turbomachines. An aeronautical engine is a highly nonlinear and multiphysics system that should be tackled by using advanced modeling techniques. From a structural viewpoint, the nonlinearities play a pivotal role in the dynamics of bladed disks. Notably, the blade-tip/casing contact interface is among the nonlinearities that are the most studied. This interface is mainly addressed for safety concerns about the engine. Indeed, the gaps between the rotating and fixed parts of the engine, here blades and the casing, are decreased in new engines as a way of increasing the aerodynamic efficiency of the turbomachine. This has the impact of yielding structural contacts in nominal operating

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: thibaut.vadcard@ec-lyon.fr (T. Vadcard).</footnote>
conditions [2,3]. These contacts at high speeds are examined with great care by designers since their consequences can threaten the integrity of engine itself, the plane and thus the safety of the passengers.

The realization of full scale experiments of blade-tip/casing contacts in order to evaluate the dynamics response of the engine are oftentimes prohibitively expensive for the manufacturers. In addition, this type of experiments cannot be carried out until the final stage of the design cycle. Such tests cannot serve to actually design the system, and are generally just about safety verifications. Hence, there is a need to setup efficient simulation tools that are able to account for the nonlinear behavior of the engine at early stage of design. This is why many researchers are working towards industry-ready, robust, predictive and accurate numerical strategies in order to characterize these sophisticated and highly nonlinear systems.

Among the strategies used to reduce computational times without hurting the accuracy of the results, the cyclic symmetry [4,5] is widely used in a linear framework. This formulation takes advantage of the angular periodicity of cyclic systems, such as bladed disks. In a perfectly tuned configuration, the dynamics of the bladed disk can be expressed as a sum of decoupled nodal diameters (spatial spectral components). Hence, by analyzing the spectral components present in the loading applied to the system, one can only keep the subset of nodal diameters that are undergoing a loading and thus drastically reduce the size of the system to be solved at no cost in terms of accuracy. However, in a nonlinear framework, the different nodal diameters are not uncoupled anymore and specific strategies should be used. Petrov [6] formulated an hypothesis that states that in the case of a traveling wave excitations, the displacement field and the nonlinear forces have the same shape as the excitation. This so-called traveling wave hypothesis greatly improves the performances of the simulations since the results on a single sector can be used to deduce the response of the full bladed disk [7]. Following this, Petrov also tackled standing wave excitations that can be considered to be traveling wave on a supersector composed of several actual sectors of the bladed disk. Siewert et al. [8] also tackled the problem of friction between shrouds in a cyclic symmetric framework with the traveling wave hypothesis. Quaefebeur [9] went further and derived a criterion allowing to determine the couplings that occur between the different nodal diameters by analyzing the Taylor expansion of the nonlinear forces. This last strategy allows to compute the response of the full system for any excitation without making any approximation on the solution. However, these strategies were designed and tested to be applied only on systems featuring friction interfaces.

Indeed, when considering systems featuring other nonlinearities than friction, such as unilateral contacts or polynomial nonlinearities, no specific criterion exists to determine how the nodal diameters are coupled with each other. Moreover, the traveling wave hypothesis remains untested on cyclic systems with blade-tip/casing contacts: this hypothesis was tested on friction, polynomial nonlinearities and shrouds unilateral contacts. Until now, simulations of blade-tip/casing contacts are carried out either on a single blade [10–16] or on a full bladed disk [17–19]. While simulations on single blades are able to obtain results in a reasonable time, they are unable to reflect the dynamics of the full bladed disk and thus yield an approximation on the vibrations of the system. On the contrary, studies on the full bladed disk do not make any assumption on the solutions but come at an extremely prohibitive computational cost. This contribution aims to tackle this modeling issue and propose a method that allows to model the full bladed disk without having to solve the problem for all unknowns.

The harmonic balance method (HBM) has recently known a growth in popularity when it comes to the characterization of blade-tip/casing interactions [11,20,21]. Its ability to gather information about nonlinear solutions is particularly powerful. Indeed, it is able to state on the stability of periodic orbits [21,22], the connectivity of the solutions [23,24] and it allows to compute nonlinear modal parameters of a system [25,26]. Moreover, it was shown to achieve a similar level of accuracy than the reference explicit time-marching strategies [11] on industrial models at a lower computational cost. For these reasons, it is generally considered that the HBM outperforms the time integration methods for the computation of periodic orbits of a nonlinear system. This article is then centered around the HBM but some of its theoretical elements can still be applied to time-marching simulations.

This paper aims to tackle the challenging subject of reduction of bladed disks in a nonlinear framework. Firstly, analytical derivations are made in order to determine a method of computing a reduction basis on which the solutions of the problem should be projected. This reduction basis is actually composed of a subset of nodal diameters that are expected to be coupled through the unilateral nonlinearity. This analysis relies on the Taylor expansion of a regularized penalty law [20]. Another reduction of the size of the problem is proposed when considering a resolution with the HBM based on the spatio-temporal couplings of nodal diameters with temporal harmonics. Then, the validity of these reductions are assessed on an academic model both for the regularized law that served to derive the reduction rules, but also on a nonsmooth law [27,28] relying on augmented lagrangians. The verification on the latter non-differentiable contact law allows to extend the reduction rule to any contact law that models unilateral laws regardless of the contact formulation that it relies on. Several contact scenarios are studied in order to give an exhaustive view of the reduction capabilities of the strategy. Then, the simulations are taken a step further to characterize a full industrial bladed disk through the cyclic symmetry based reduction strategy, on both contact laws again. The main original contribution of this paper is the new rule that allows to determine in advance the nodal diameters involved in the nonlinear response of a cyclic system subjected to strong unilateral contacts based on the loading's spatial harmonic index for any shape of excitation. It is a generalization of the traveling wave hypothesis proposed by Petrov [6] for all excitation shapes in the case of unilateral contact dominated nonlinear dynamics.

The physical framework relying on the cyclic symmetry is described in Section 2: it includes the formulation of the problem and the derivation of the reduction criterion. Then, the HBM-related resolution strategies are described in Section 3, it also includes the formulation of the supplementary reduction criterion based on the HBM. Then, numerous simulations are led to validate the reduction strategy on a simplified model in cyclic symmetry in Section 4. Finally, all the presented developments are applied to a representative fan blade model in Section 5. It is shown that the procedure is industry-ready and allows to accurately model the dynamics of the full bladed disk at a low computational cost.
## 2. Cyclic symmetry and reduction strategy

### 2.1. Equations of motion in a cyclic symmetric framework

In this section, the equations of motion for a cyclic symmetric mechanical system are derived in terms of spectral generalized coordinates. For each sector $j \in [1, N]$ where N is the number of sectors of the system, the equations of motion read

$$ M_0 \ddot{x}_j + C_0 \dot{x}_j + K_0 x_j = \mathbf{f}_{s,j}^P + \mathbf{f}_{s,j}^n + \mathbf{f}_{nl,j}(x_j(t), \dot{x}_j(t)) + \mathbf{f}_{0,j} + \mathbf{f}_{ex,j}(\omega, t) \quad (1) $$

where $M_0$, $C_0$ and $K_0$ are respectively the mass, damping and stiffness matrix of a single sector of the structure. The full structure is considered to be tuned, meaning that all sectors are identical. $x_j$ is the displacement field of the n degrees of freedom (DoF) of the jth sector. In order to represent the full system, the influence of neighboring sectors is accounted for by the forces applied to the cyclic boundary of the sector $\mathbf{f}_{s,j}^P$ and $\mathbf{f}_{s,j}^n$, related to the forces applied by the previous and the next sector respectively. Note that sector number 1 is adjacent to sectors N and 2, and that sector number N is adjacent to sectors N-1 and 1 so that the cyclic nature is preserved. The system also undergoes nonlinear forces, represented by $\mathbf{f}_{nl,j}$, external pre-loading forces $\mathbf{f}_{0,j}$ and harmonic excitation forces $\mathbf{f}_{ex,j}$. The latter include a dependence in $\omega$ because the present article focuses on the case of a monoharmonic excitation of frequency $\omega$ on the system. All the developments hereafter are valid for both transient and periodic motions but since the harmonic balance method is used to solve the motion equations, only periodic motions will be studied throughout the simulations of this article.

The physical displacement $x_j$ can then be expressed as a function of the real spatial spectral components (nodal diameters) of the cyclic system thanks to the theory of linear representation of cyclic groups of order N [5]. It follows that

$$ x_j = \frac{1}{\sqrt{N}} \hat{x}_0 + \sqrt{\frac{2}{N}} \sum_{h=1}^{\frac{N}{2}-1} \left[ \hat{x}_{h,c} \cos(h(j-1)\alpha) + \hat{x}_{h,s} \sin(h(j-1)\alpha) \right] + \frac{(-1)^{(j-1)}}{\sqrt{N}} \hat{x}_{N/2} \quad \text{if N is even,} \quad (2) $$

$$ x_j = \frac{1}{\sqrt{N}} \hat{x}_0 + \sqrt{\frac{2}{N}} \sum_{h=1}^{\frac{N-1}{2}} \left[ \hat{x}_{h,c} \cos(h(j-1)\alpha) + \hat{x}_{h,s} \sin(h(j-1)\alpha) \right] \quad \text{if N is odd} \quad (3) $$

where $\alpha = \frac{2\pi}{N}$ is the intersector angle and $\hat{x}_{h,c}$ and $\hat{x}_{h,s}$ are the real spectral displacements of order h. With the cyclic Fourier matrix S, physical displacements of all sectors $x = [x_1^T, \dots, x_j^T, \dots, x_N^T]^T$ can be substituted by the spectral components such as $x = S^T \hat{x}$ with $\hat{x} = [\hat{x}_0^T, \dots, \hat{x}_h^T, \dots, \hat{x}_m^T]^T$ with $m = \frac{N}{2}$ if N is even or $m = \frac{N-1}{2}$ if N is odd. All spatial spectral components of order $0 < h < \frac{N}{2}$ contain both cosine and sine contributions such as $\hat{x}_k = [\hat{x}_{h,c}^T, \hat{x}_{h,s}^T]^T$. The notation $\bullet^T$ is the transposition operation.

In this work the cyclic boundaries of the sectors remain linear, thus the different cyclic boundaries are assembled [4,5] leading to a system of equations where all matrices are block-diagonal accounting for the fact that all cyclic indexes are uncoupled. In the nonlinear framework studied, a specific coupling is appearing between the cyclic indexes due to the nonlinear forces applied on the system. This coupling is the main topic of this contribution and is thoroughly addressed throughout the article. In the end, the problem associated to the full system is then expressed such as

$$ \hat{M} \ddot{\hat{x}}(t) + \hat{C} \dot{\hat{x}}(t) + \hat{K} \hat{x}(t) = \hat{f}_{ex}(\omega, t) + \hat{f}_{0} + \hat{f}_{nl}(x(t), \dot{x}(t)) \quad (4) $$

with

$$ \hat{M} = \text{blkdiag} \left( \left\{ \hat{M}^h \right\}_{h=0}^m \right), \quad \hat{C} = \text{blkdiag} \left( \left\{ \hat{C}^h \right\}_{h=0}^m \right), \quad \hat{K} = \text{blkdiag} \left( \left\{ \hat{K}^h \right\}_{h=0}^m \right), \quad (5) $$

$$ \hat{f}_{ex}(\omega, t) = \begin{bmatrix} \hat{f}_{0}^{0} \\ \vdots \\ \hat{f}_{ex}^{h} \\ \vdots \\ \hat{f}_{0}^{m} \end{bmatrix}, \quad \hat{f}_{0} = \begin{bmatrix} \hat{f}_{0}^{0} \\ \vdots \\ \hat{f}_{0}^{h} \\ \vdots \\ \hat{f}_{0}^{m} \end{bmatrix} \quad \text{and} \quad \hat{f}_{nl}(x(t), \dot{x}(t)) = \begin{bmatrix} \hat{f}_{nl}^{0} \\ \vdots \\ \hat{f}_{nl}^{h} \\ \vdots \\ \hat{f}_{nl}^{m} \end{bmatrix}. \quad (6) $$

The notation $\bullet^h$ refers to spatial spectral components of order h, oftentimes referred to as nodal diameters (ND) of order h. It is worth noting that the nonlinear forces expressed in terms of spatial spectral components $\hat{f}_{nl}$ still depends on physical displacements $x(t)$ and velocities $\dot{x}(t)$. Indeed, the evaluation of nonlinear forces, presented further in the article in Section 2.2 should be made on physical coordinates.

This general formulation in cyclic symmetry is actually equivalent of solving the full problem with all sectors and thus yields no computational reduction with respect to the full physical problem. Nevertheless, this strategy enables to perform full computations without defining the whole problem as a single matrix. This particular property is helpful when generating reduced-order models of industrial size models. These models generally cause issues in terms of memory that can be avoided through the use of cyclic symmetry.

Since the nodal diameters are linearly decoupled (i.e. $\hat{M}$, $\hat{C}$ and $\hat{K}$ are block-diagonal), the only link between the different nodal diameters in Eq. (4) is made through the nonlinear forces. Hence, it is possible to analyze the analytical expressions of the nonlinear forces in order to determine the couplings that occur between nodal diameters. By doing so, a subset of nodal diameters $D \subset [0, m]$ is
determined for a specific loading scenario. The latter can be used for a nonlinear calculation instead of the whole range $D^f = \ll 0, m \gg$. Indeed, for certain loading scenarios, some nodal diameters are not involved at all in the dynamics and can thus be withdrawn from the unknown vector without affecting the accuracy of the solutions. This generally yields a significant reduction in computation time and memory usage. Eq. (4) becomes

$$ \hat{\mathbf{M}}^{D} \ddot{\hat{\mathbf{x}}}(t) + \hat{\mathbf{C}}^{D} \dot{\hat{\mathbf{x}}}(t) + \hat{\mathbf{K}}^{D} \hat{\mathbf{x}}(t) = \hat{\mathbf{f}}_{\text{ex}}^{D}(\omega, t) + \hat{\mathbf{f}}_{0}^{D} + \hat{\mathbf{f}}_{\text{nl}}^{D}(\mathbf{x}(t), \dot{\mathbf{x}}(t)) \quad (7) $$

with

$$ \hat{\mathbf{M}}^D = \text{blkdiag} \left( \left\{ \hat{\mathbf{M}}^h \right\}_{h \in D} \right), \quad \hat{\mathbf{C}}^D = \text{blkdiag} \left( \left\{ \hat{\mathbf{C}}^h \right\}_{h \in D} \right), \quad \hat{\mathbf{K}}^D = \text{blkdiag} \left( \left\{ \hat{\mathbf{K}}^h \right\}_{h \in D} \right), \quad (8) $$

$$ \hat{\mathbf{f}}_{\text{ex}}^{D} = \left[ \left\{ \hat{\mathbf{f}}_{\text{ex}}^{h} \right\}_{h \in D} \right], \quad \hat{\mathbf{f}}_{0}^{D} = \left[ \left\{ \hat{\mathbf{f}}_{0}^{h} \right\}_{h \in D} \right] \quad \text{and} \quad \hat{\mathbf{f}}_{\text{nl}}^{D} = \left[ \left\{ \hat{\mathbf{f}}_{\text{nl}}^{h} \right\}_{h \in D} \right]. \quad (9) $$

In the general case, the choice of an appropriate set of nodal diameter is a difficult task. However, by analyzing the nonlinearities at stake, it is possible to derive rules that provide ways of determining the right nodal diameters to keep in $D$. For instance, Quagebeur et al. derived a specific rule allowing to determine the couplings between the spatial spectral components of a cyclic structure with friction nonlinearities [9]. This rule analyzes the spectral components of the excitations and the number of sectors of the system in order to determine the subset $D$ of nodal diameters allowing for eventual reduction of the system to be solved. This subset is used by only accounting for the nodal diameter that can be coupled with the excitation as well as with each other by the nonlinearity. In this work, the same type of rule is derived for nonlinearities that model unilateral contacts, it is presented in Section 2.3.

The developments proposed here fall within the framework of the hypotheses of a perfectly tuned bladed disk with no mistuning as a way of taking advantage of the properties of cyclic symmetry. The impact of possible mistuning on the reduction approach is currently beyond the scope of the article and will require future developments.

## 2.2. Evaluation of nonlinear forces in the physical space

It is reminded that the nonlinear forces oftentimes depend on physical displacements and should be evaluated in the physical space. In order to go from the spectral (cyclic) variables to physical (sectors) variables, one should involve the spatial Fourier matrix $S$ and its inverse. As a way of computing the spatial spectral components of the nonlinear forces, the cyclic variables linked to the displacements $\hat{\mathbf{x}}^D(t)$ are firstly transformed to physical variables on all sectors $\mathbf{x}(t)$, then the nonlinear forces are computed in the physical space on each sector $\mathbf{f}_{nl}(t)$ and these physical variables are then transformed back to the cyclic space $\hat{\mathbf{f}}_{nl}^D(t)$.

Firstly, let us assume that $D$ is chosen to be the full set $D = D^f$. Then, for sampling purposes, it is necessary to evaluate the nonlinear forces on all sectors $S^f = \{1, 2, \dots, N\}$ in order to be able to reconstruct the whole spectral decomposition.

However, given Eq. (7), the unknown vectors $\hat{\mathbf{x}}^D(t)$ and $\hat{\mathbf{f}}_{nl}^D(\mathbf{x}(t), \dot{\mathbf{x}}(t))$ are composed of all spectral components that are retained in the subset $D \subset D^f$. Then, if this time $D \neq D^f$ (does not contain every spectral component), the equivalence between physical and cyclic spaces yields that the less nodal diameters are considered, the less sectors are necessary for the computation of the spectral components retained in $D$. Indeed, by suppressing some spectral components of the full basis, the information over all sectors becomes redundant.

As a way of determining the set of sectors $S \subset S^f$, over which the nonlinear forces should be sampled (composed of a number of $j^{\max}$ sectors), it is necessary to analyze the content of $D$. For each non-degenerated spectral component in $D$ (i.e. $h = 0$ and $h = N/2$ if $N$ is even and only $h = 0$ if $N$ is odd), a single sector should be evaluated. For each degenerated spectral components $0 < h < N/2$ in $D$, two sectors should be evaluated. For instance, a structure composed of $N = 24$ sectors for which $D$ is arbitrarily set to $D = \{0, 2, 4, 6, 8, 10, 12\}$, the number of sectors needed for the evaluation of all spectral components of the nonlinear forces is $j^{\max} = 2 \times 1 + 5 \times 2 = 12$ (2 non-degenerated and 5 degenerated spectral components). Then, one should evaluate the nonlinear forces on all sector numbers belonging to $S = \{1, 2, \dots, j^{\max}\}$.

Following these guidelines, it is possible to establish a relation between spectral components restricted to $D$ and the physical information over the sectors of $S$. In order to clarify the notations, complementary sets $\bar{D} = D^f \setminus D$ and $\bar{S} = S^f \setminus S$ are defined. The relation is obtained by partitioning $\hat{\mathbf{x}} = S\mathbf{x}$ over the different sets such as

$$ \begin{bmatrix} \hat{\mathbf{x}}^D \\ \hat{\mathbf{x}}^D \end{bmatrix} = \begin{bmatrix} \mathbf{S}_S^D & \mathbf{S}_S^D \\ \mathbf{S}_S^D & \mathbf{S}_S^D \end{bmatrix} \begin{bmatrix} \mathbf{x}^S \\ \mathbf{x}^S \end{bmatrix}. \quad (10) $$

By definition, restricting the spectral components to a specific set $D$ means that $\hat{\mathbf{x}}^D = \mathbf{0}$. Then, the inverse transformation is expressed as $\mathbf{x} = S^T \hat{\mathbf{x}}$. Since the Fourier matrix $S$ is orthogonal, it yields

$$ \begin{bmatrix} \mathbf{x}^S \\ \mathbf{x}^S \end{bmatrix} = \begin{bmatrix} (\mathbf{S}_S^D)^T & (\mathbf{S}_S^D)^T \\ (\mathbf{S}_S^D)^T & (\mathbf{S}_S^D)^T \end{bmatrix} \begin{bmatrix} \hat{\mathbf{x}}^D \\ \mathbf{0} \end{bmatrix} \quad (11) $$

and it comes that

$$ \mathbf{x}^S = (\mathbf{S}_S^D)^T \hat{\mathbf{x}}^D. \quad (12) $$
Since the set S was built to be equivalent to D in the physical space, $(S_S^D)^T$ is square and invertible, yielding the direct transformation from S to D is

$$ \hat{\mathbf{x}}^D = \left[ (S_S^D)^T \right]^{-1} \mathbf{x}^S. \qquad (13) $$

This shows that restricting the number of spectral components in D yields a double advantage: firstly it decreases the overall number of unknowns of the problem and secondly it reduces the computational cost of evaluating nonlinear forces since it requires their evaluation on less sectors than for the full system. This is illustrated in Fig. 1 later in the article.

## 2.3. Selection of nodal diameters for unilateral contacts

In this section, a complex formulation is employed for the sake of simplifying the developments, however both complex and real cyclic symmetry feature the same properties. In this work, the resolution is carried out with real cyclic symmetry. The complex formulation is only used in this section and it relies on the complex counterpart of the spectral decomposition in Eq. (3).

In the complex formalism, the spectral decomposition of the displacement field of each sector $x_j$ reads

$$ x_j = \frac{1}{\sqrt{N}} \sum_{h=0}^{N-1} (\hat{\mathbf{x}}^h e^{iah(j-1)}) \quad \text{for } j \in [1, N], \qquad (14) $$

with $\hat{\mathbf{x}}^h = \hat{\mathbf{x}}^{-N-h}$ for $h \in [0, m]$. Similarly as in the real formalism, a block-diagonal system can be built for all $h \in D^f$ such as

$$ \hat{\mathbf{M}}^h \hat{\mathbf{x}}^h(t) + \hat{\mathbf{C}}^h \hat{\mathbf{x}}^h(t) + \hat{\mathbf{K}}^h \hat{\mathbf{x}}^h(t) = \hat{\mathbf{f}}^h_{\text{ex}}(\omega, t) + \hat{\mathbf{f}}^h_0 + \hat{\mathbf{f}}^h_{\text{nl}}(\mathbf{x}(t), \hat{\mathbf{x}}(t)). \qquad (15) $$

As stated, when nonlinearities are taken into account the nonlinear forces create interactions between the different nodal diameters and Eq. (15) must be solved for all nodal diameter at once.

In this work, unilateral contact nonlinearities are considered. For theoretical purposes, a regularized penalty law [20] is considered throughout these developments, but it is shown in this article that these results also apply to non-regularized laws. Furthermore, a single nonlinear DoF per sector is assumed, its displacement is noted $x_{nl,j}$. The extension to other DoF is straightforward, based on the developments that follow. The nonlinear forces for sector j are evaluated with the following

$$ f_{nl,j}(t) = -\kappa \frac{(g_j(t) - x_{nl,j}(t))}{2} + \sqrt{\left( \kappa \frac{(g_j(t) - x_{nl,j}(t))}{2} \right)^2 + \gamma^2}, \qquad (16) $$

where $\kappa$ is the contact stiffness and $\gamma$ the regularization parameter. Thermo-mechanical and wear effects are neglected in this study because if a certain symmetry is predicted for the contact interactions, then both the thermal field in the bladed disk and the wear pattern on the abradable layer would be governed by the same symmetry. Hence, since no supplementary couplings are expected, these multi-physics coupling are not accounted for as a way of simplifying the simulations for easier interpretation.

The gap function can be written in terms of cyclic components such as

$$ g_j = \frac{1}{\sqrt{N}} \sum_{h=0}^{N-1} (\hat{g}^h e^{iah(j-1)}) \quad \text{for } j \in [1, N]. \qquad (17) $$

The gap is a known quantity, as such all the components $(\hat{g}^h)_{h \in D^f}$ are also known. As an example, for a perfectly constant gap over one period, only the component $\hat{g}_0$ is different from 0. To simplify the following derivations, the new unknown of the problem is set to

$$ \hat{u}^h = \hat{g}^h - \hat{x}_{nl}^h. \qquad (18) $$

The equation of motion, accounting for the nonlinearities reads

$$ \hat{\mathbf{M}}^h \hat{\mathbf{u}}^h + \hat{\mathbf{C}}^h \hat{\mathbf{u}}^h + \hat{\mathbf{K}}^h \mathbf{u}^h = \hat{\mathbf{f}}^h_{\text{ext}} + \hat{\mathbf{f}}^h_0 + \hat{\mathbf{f}}^h_{\text{nl}} - \hat{\mathbf{K}}^h \mathbf{g}^h \quad \text{for } h \in D^f, \qquad (19) $$

where the force applied to the only nonlinear DoF of the system inside $\hat{\mathbf{f}}^h_{nl}$ is expressed as

$$ \hat{f}_{nl}^h = \frac{1}{\sqrt{N}} \sum_{j=1}^{N} (f_{nl,j} e^{-iah(j-1)}) \quad \text{for } h \in D^f. \qquad (20) $$

Substituting the term $f_{nl,j}$ in Eq. (20) by Eq. (16) (where the physical quantities have been replaced by the cyclic ones) for a given $q \in D^f$ reads

$$ \hat{f}_{nl}^q = \frac{1}{\sqrt{N}} \sum_{j=1}^{N} \left[ \left( -\frac{\kappa}{2} \sum_{h=0}^{N-1} \hat{u}^h e^{ih(j-1)\alpha} + \frac{\kappa}{\gamma} \sqrt{1 + \frac{1}{\gamma^2} \left( \sum_{h=0}^{N-1} \hat{u}^h e^{ih(j-1)\alpha} \right)^2} \right) e^{-iq(j-1)\alpha} \right]. \qquad (21) $$

The second term is expanded as a Taylor series at 0

$$ \sqrt{1 + y^2} = 1 + \sum_{l=0}^{\infty} \beta_l y^{2l}. \qquad (22) $$
This leads to

$$ \hat{f}_{nl}^{q} = \frac{1}{\sqrt{N}} \sum_{j=1}^{N} \left[ \left( -\frac{\kappa}{2} \sum_{h=0}^{N-1} \hat{u}^{h} e^{i(h-q)(j-1)\alpha} + \frac{\kappa}{\gamma} \left( 1 + \sum_{l=0}^{\infty} \left( \frac{\beta_l}{\gamma^2} \left( \sum_{h=0}^{N-1} \hat{u}^{h} e^{i(h-1)\alpha} \right)^{2l} \right) \right) \right) e^{-iq(j-1)\alpha} \right]. \quad (23) $$

Notice that the first term in $\hat{u}_h$ is different from 0 when $h = q$. This means that if a nodal diameter responds, it will provide a component to the nonlinear forces on the same nodal diameter. Moreover the constant term $\frac{\kappa}{\gamma}$ provides a nonlinear force component on its 0th nodal diameter. Finally, the following expression explicits the last term (noted $\hat{f}_{nl,3}^q$) of Eq. (23)

$$ \hat{f}_{nl,3}^{q} = \frac{\kappa}{\sqrt{N\gamma^3}} \sum_{j=1}^{N} \left( \sum_{l=0}^{\infty} \left( \beta_l \left( \sum_{h=0}^{N-1} \hat{u}_h e^{i(h-1)\alpha} \right)^{2l} \right) \right) e^{-iq(j-1)\alpha} \quad (24) $$

$$ = \frac{\kappa}{\sqrt{N\gamma^3}} \sum_{l=0}^{\infty} \sum_{j=1}^{N} \left( \left( \beta_l \left( \sum_{h_1=0}^{N-1} \hat{u}_{h_1} e^{i h_1 (j-1) \alpha} \right) \right) \cdots \left( \sum_{h_{2l}=0}^{N-1} \hat{u}_{h_{2l}} e^{i h_{2l} (j-1) \alpha} \right) e^{i \alpha (j-1) - q} \right). \quad (25) $$

Developing the sums in Eq. (25) gives

$$ \hat{f}_{nl,3}^{q} = \frac{\kappa}{\sqrt{N\gamma^3}} \sum_{l=0}^{\infty} \left( \beta_l \sum_{k \in \mathcal{T}_l} \left( \alpha_k \prod_{p_k=1}^{Q_k} (\hat{u}_{h_{p_k}})^{\beta_{p_k}} \left( \sum_{j=1}^{N} e^{i\alpha(j-1)(\sum_{p_k=1}^{Q_k} h_{p_k})} \right)^{-q} \right) \right). \quad (26) $$

The set $\mathcal{T}_l$ describes all the possible combinations of the sums' multiplications. Each combination $k$ involves one or multiple nodal diameters noted $h_{p_k}$ and their multiplicity $\beta_{p_k}$. The terms $\alpha_k$ denotes the number of occurrence of such a term. Since the degree of the polynomial in Eq. (26) is $\sum_{p_k=1}^{Q_k} \beta_{p_k}$. This degree must be equal to the sum of the degrees of all monomials associated with the $2l$ products of Eq. (25) (from $h_1$ to $h_{2l}$). Thus, the following relation must be verified

$$ \sum_{p_k=1}^{Q_k} \beta_{p_k} = 2l. \quad (27) $$

It corresponds to the sum of the degrees of the different monomials in the combination $k$ considered. The nodal diameter couplings that are created through all terms of Eq. (26) are now analyzed in Section 2.4.

Note that the parameter $\gamma$ does not intervene directly in the relationship Eq. (27) that allows to determine the couplings between nodal diameters. No numerical evaluation is needed to derive the couplings, hence, it is expected that the regularization of the contact law does not influence how nodal diameters interact because it can be as small as possible. In the results section of this article, this particular property of the proposed method is verified by applying the criterion on configurations where the contact is modeled with the non-regularized DLFT-HBM contact law, thus modeling a case where $\gamma = 0$.

## 2.4. Determination of the nodal diameter coupling

This section proposes an algorithm to determine the nodal diameter coupling occurring when contact nonlinearities are taken into account, it relies on the equations set up is Section 2.3. Firstly, $D_q$ denotes the set of nodal diameters that include the nodal diameters corresponding to $\hat{g}$ and the external forces, that is known data. Based on Eq. (19), all these nodal diameters are expected to respond and are thus kept in the basis. Moreover, based on Eq. (23), due to the first term, the 0th nodal diameter also has to be included in the basis, now noted $D_{q,0}$. This basis is being completed through an iterative process with Eq. (24): a new nodal diameter noted $p$ must be included in the basis if

$$ \sum_{i=1}^{2l} h_i - p \equiv 0[N] \quad \text{with} \quad h_i \in D_{q,0} \quad \text{and} \quad l \in \mathbb{N}. \quad (28) $$

As $l$ goes up to infinity, such an iterative process is cumbersome. However, a few mathematical considerations allow to perform such procedure in an efficient way. First, it is assumed that the nodal diameter basis only contains the diameters 0 and $h_{ex}$ (the nodal diameter associated with the external loading). The 0th nodal diameter is not involved in the iterative process as it does not create further coupling. Then, Eq. (28) becomes

$$ 2lh_{ex} \equiv p[N], \quad l \in \mathbb{N}. \quad (29) $$

Following this, there exists an integer $K$ such that the remainder of the euclidian division of $2h_{ex}$ by $l$ (noted $p$), for $l \ge K$, is the same as the one of the euclidian division of $2h_{ex}$ by $l$ with $l \in [1, K]$. Then, it comes that

$$ 2(l+K)h_{ex} \equiv p[N]. \quad (30) $$

As $2lh_{ex} \equiv p[N]$, one obtains the following

$$ 2Kh_{ex} \equiv 0[N] \quad (31) $$
Therefore, K is the first integer for which $2Kh_{ex}$ is a multiple of N. The determination of K is then made thanks to the implementation of

$$K = \frac{\operatorname{lcm}(2h_{ex}, N)}{2h_{ex}} \qquad (32)$$

with $\operatorname{lcm}(\cdot, \cdot)$ returns the lowest common multiple of the two arguments. Note that Eq. (32) cannot be used when $h_{ex} = 0$, but in that case $K = 1$ and only the nodal diameter 0 is involved in the response.

Based on this result, the updated set $D_{q,0}$ must contains 0, $h_{ex}$, and all even multiples of $h_{ex}$ up to $2K$. However, these new nodal diameters can get coupled once again with $h_{ex}$ and thus the set $D_{q,0}$ must also contain all the uneven multiple of $h_{ex}$ up to $2K$. However thanks to Eq. (31), the remainder of the euclidian division of $(2K - l)h_{ex}$ by N is equal to the remainder of the euclidian division of $l h_{ex}$ by N for $l \in [1, K]$. As a consequence, in order to determine the nodal diameter coupling when the initial set $D_{q,0}$ contains the nodal diameter $h_{ex}$, one must first determine K through Eq. (31) and then include all the multiples of $h_{ex}$ up to K only. Then, the final set D is built following these three steps

1. the list $v_{K,h_{ex}}$ of all multiples of $h_{ex}$ up to K is firstly built such as

$$v_{K,h_{ex}} = h_{ex} \times \underbrace{[0, 1, 2, \dots, q, \dots, K]}_{v_K} = [0, h_{ex}, 2h_{ex}, \dots, qh_{ex}, \dots, Kh_{ex}], \qquad (33)$$

2. then, the remainder of each element of $v_{K,h_{ex}}$ by the euclidian division by N, noted $r_q \forall q \in v_K$, is stored in the vector $v_R = [r_0, r_1, r_2, \dots, r_q, \dots, r_K]$,

3. finally, the nodal diameters associated with each $r_q$ are computed and stored in $D = [d_0, d_1, \dots, d_q, \dots, d_K]$ with $d_q = r_q$ if $r_q < m$ and $d_q = N - r_q$ if $r_q > m \forall q \in v_K$

Examples of application of this strategy are found throughout this article: each time a computation in cyclic symmetry is conducted with this criterion, the value of K and the final basis used is given.

In the case where the initial set $D_{q,0}$ contains several nodal diameters, two cases arise. First, if the included nodal diameters share a common divider strictly greater than 1, then the previous procedure can be employed to obtain a reduction by using $h_{ex} = \gcd(h_{ex}^1, h_{ex}^2)$, where $\gcd(\cdot, \cdot)$ returns the greater common divider of the two arguments. The other case corresponds to the situation where the nodal diameters are prime with one another $(\gcd(h_{ex}^1, h_{ex}^2) = 1)$. In this case, unfortunately all nodal diameters respond and no reduction can be achieved.

An additional case where no reduction can be achieved by this criterion is the configuration where N is a prime number. Then, no engine order (except the 0th) are subject to any reduction and full calculations should be made. This property is a consequence of the fact that $K = N$ when N is prime and it yields that $D = D^t$.

### 3. Numerical framework: harmonic balance method with nodal diameters selection

The resolution of Eq. (7) can be carried out indifferently with time domain or frequency domain methods. The determination of the set D does not make any assumption on the shape of the solution, thus the strategy presented in Section 2.4 is valid for all numerical strategies. In this work, the retained solving method for the application of the diameter selection strategy method is the HBM. Firstly the general formulation of the HBM is applied to Eq. (4), then a particular form of the HBM is used to take advantage of the selection of certain nodal diameters. A new criterion is also proposed in addition to the one in Section 2.4. It takes advantage of the coupling between nodal diameters and temporal harmonic to halve the size of the system. Finally, specific developments concerning the response to traveling wave excitations are carried out.

#### 3.1. General harmonic balance method

The HBM is used to solve Eq. (4) for periodic solutions and allows fast solving and flexibility. It relies on the assumption that, under an excitation of angular frequency $\omega$, the time-dependent unknown $\hat{x}(t)$ is written as a truncated Fourier series of order $N_h$ and fundamental angular frequency $\omega$

$$\hat{x}(t) \approx \frac{1}{2}\hat{a}_0 + \sum_{k=1}^{N_h} (\hat{a}_k \cos(k\omega t) + \hat{b}_k \sin(k\omega t)) \qquad (34)$$

This decomposition allows to perform a Fourier-Galerkin projection on Eq. (4) in order to obtain an algebraic system of equations related to cyclic coordinates in the frequency domain

$$Z(\omega)\hat{\mathbf{x}} - \hat{\mathbf{f}}_{nl}(\hat{\mathbf{x}}) - \hat{\mathbf{f}}_{ex} - \hat{\mathbf{f}}_0 = \mathbf{0}. \qquad (35)$$

The new unknown of the problem is then the so-called multiharmonic vector $\hat{\mathbf{x}}$ containing the real Fourier decomposition of the spectral components of the displacement field $\mathbf{x}(t)$, written as in Eq. (36a). The vectors associated with external and nonlinear forces
are also organized as multiharmonic vector which expressions are given in Eqs. (36b), (36c) and (36d)

$$ \hat{\mathbf{x}} = \left[ \frac{1}{2} \hat{\mathbf{a}}_0^\mathrm{T} \hat{\mathbf{a}}_1^\mathrm{T} \hat{\mathbf{b}}_1^\mathrm{T} \dots \hat{\mathbf{a}}_{N_h}^\mathrm{T} \hat{\mathbf{b}}_{N_h}^\mathrm{T} \right]^\mathrm{T}, \qquad (36a) $$

$$ \hat{\mathbf{f}}_{nl} = \left[ \frac{1}{2} (\hat{\mathbf{a}}_0^{nl})^\mathrm{T} (\hat{\mathbf{a}}_1^{nl})^\mathrm{T} (\hat{\mathbf{b}}_1^{nl})^\mathrm{T} \dots (\hat{\mathbf{a}}_{N_h}^{nl})^\mathrm{T} (\hat{\mathbf{b}}_{N_h}^{nl})^\mathrm{T} \right]^\mathrm{T}, \qquad (36b) $$

$$ \hat{\mathbf{f}}_{\mathrm{ex}} = [\mathbf{0} \ (\hat{\mathbf{a}}_1^{\mathrm{ex}})^\mathrm{T} \ (\hat{\mathbf{b}}_1^{\mathrm{ex}})^\mathrm{T} \ \mathbf{0} \ \dots \ \mathbf{0} \ \mathbf{0}]^\mathrm{T}, \qquad (36c) $$

$$ \hat{\mathbf{f}}_0 = [(\hat{\mathbf{a}}_0^{\mathrm{T}} \ \mathbf{0} \ \dots \ \mathbf{0} \ \mathbf{0})]^\mathrm{T}. \qquad (36d) $$

$Z$ is a block-diagonal matrix accounting for the linear behavior of the system and is composed of ($N_h + 1$) blocks, defined such as

$$ Z_0 = \hat{\mathbf{K}} \quad \text{and} \quad Z_k(\omega) = \begin{bmatrix} \hat{\mathbf{K}} - (k\omega)^2 \hat{\mathbf{M}} & k\omega \hat{\mathbf{C}} \\ -k\omega \hat{\mathbf{C}} & \hat{\mathbf{K}} - (k\omega)^2 \hat{\mathbf{M}} \end{bmatrix} \quad \forall k \in [1, N_h]. \qquad (37) $$

## 3.2. Modifications induced by the selection of specific nodal diameters

However, the HBM allows a specific methodology that consists in selecting a subset of all spectral components, or nodal diameters, to reduce the size of the system that is actually solved. While in the time domain equation the restriction of the unknown to a set of nodal diameters $D$ affects the whole equation in the same way, it is not necessarily the case for HBM. Indeed, since the blocks $Z_k$ of the dynamic stiffness matrix are independent, it is possible to use different subsets $D^k$ specific to each temporal harmonic $k$ as a way of taking advantage of the spatio-temporal couplings that occur in cyclic structures. Then, these sets $D^k$ are used to build specific constitutive matrices $\hat{\mathbf{M}}_k$, $\hat{\mathbf{C}}_k$ and $\hat{\mathbf{K}}_k$ and corresponding loading and unknown vectors of each $k$th harmonic. Then, Eq. (37) is rewritten such as

$$ Z_0 = \hat{\mathbf{K}}_0 \quad \text{and} \quad Z_k(\omega) = \begin{bmatrix} \hat{\mathbf{K}}_k - (k\omega)^2 \hat{\mathbf{M}}_k & k\omega \hat{\mathbf{C}}_k \\ -k\omega \hat{\mathbf{C}}_k & \hat{\mathbf{K}}_k - (k\omega)^2 \hat{\mathbf{M}}_k \end{bmatrix} \quad \forall k \in [1, N_h]. \qquad (38) $$

with

$$ \hat{\mathbf{M}}_k = \mathrm{blkdiag} \left( \left\{ \hat{\mathbf{M}}^h \right\}_{h \in D^k} \right), \quad \hat{\mathbf{C}}_k = \mathrm{blkdiag} \left( \left\{ \hat{\mathbf{C}}^h \right\}_{h \in D^k} \right), \quad \hat{\mathbf{K}}_k = \mathrm{blkdiag} \left( \left\{ \hat{\mathbf{K}}^h \right\}_{h \in D^k} \right) \qquad (39) $$

The different loading vectors of Eq. (6) should also be rewritten as

$$ \hat{\mathbf{f}}_{\mathrm{ex},k} = \left[ \left\{ \hat{\mathbf{f}}_{\mathrm{ex}}^h \right\}_{h \in D^k} \right] \quad \text{and} \quad \hat{\mathbf{f}}_{0,k} = \left[ \left\{ \hat{\mathbf{f}}_0^h \right\}_{h \in D^k} \right] \qquad (40) $$

depending on the temporal harmonic number $k$ they are related to. Following these notations, the unknown vector is also reduced and defined as $\hat{\mathbf{f}}_k = [\hat{\mathbf{x}}_k^h]_{h \in D^k}$ and $\hat{\mathbf{f}}_{nl,k} = [\hat{\mathbf{f}}_{nl}^h]_{h \in D^k}$.

If all subsets $D^k$ for $k \in [0, N_h]$ are $D^k = D^f$, then this strategy is equivalent to solving the full system with all sectors.

In the literature, Petrov showed that for cyclic symmetric systems undergoing nonlinearities and traveling wave excitations, solution displacement fields can be searched as traveling waves of the same shape as the excitation [6]. This particular hypothesis allows to select only a single index in each of the $D^k$, significantly reducing the computational costs without accuracy loss. Knowing the direction of propagation of the traveling wave also allows to divide the size of the system by two.

At this point, it is possible to perform a condensation of the linear DoF by taking advantage of the algebraic nature of the HBM equation and the fact that contact nonlinearities are localized [28]. Then, the system is restricted to nonlinear DoF and the computations are significantly sped up. The reduced system of equations formulated on sets of nodal diameters $D^k$ finally reads

$$ Z_{k,r} \hat{\mathbf{x}}_{k,r} - \hat{\mathbf{f}}_{nl,k,r} - \hat{\mathbf{f}}_{\mathrm{ex},k,r} - \hat{\mathbf{f}}_{0,k,r} = \mathbf{0} \quad \forall k \in [0, N_h] \qquad (41) $$

that is aggregated into a single equation

$$ Z_r \hat{\mathbf{x}}_r - \hat{\mathbf{f}}_{nl,r} - \hat{\mathbf{f}}_{\mathrm{ex},r} - \hat{\mathbf{f}}_{0,r} = \mathbf{0}. \qquad (42) $$

This equation is solved by means of a Newton-Raphson nonlinear iterative solver. The jacobian matrix of Eq. (42) should thus be computed at each iteration.

The proposed reduction technique is not consuming any supplementary time since it consists in assuming some unknowns of the full problem (unaccounted nodal diameters) to be zero. Two main computational reductions are identified for the HBM as a consequence of the proposed method at two stages of the nonlinear solving process:

* The size of the full problem is reduced, yielding savings on the application of the condensation of linear DoF [28] and during the Newton-Raphson solving process. This scales like the cube of the system size $O(n^3)$ because it involves matrix inversions and multiplications.

* The number of nonlinear nodes of the interface is reduced because the number of sectors in $S$ is lower than the full number of sectors. It yields a linear saving $O(n)$ since the nonlinear forces need to be evaluated on less nodes.
Both computational gains are related to two stages of the solving process that are carried out once per iteration (for an arc-length continuation), then according to the partition between linear and nonlinear DoF (contact interface) in the system, the reduction is expected to scale between $\mathcal{O}(n)$ and $\mathcal{O}(n^3)$. The size of the problem $n$ is itself proportional to $\frac{j^{\max}}{N}$ with $j^{\max}$ the number of sectors for AFT-CP, so reductions in terms of time saving are expected to scale between $\frac{j^{\max}}{N}$ and $(\frac{j^{\max}}{N})^3$.

### 3.3. Parity criterion for unilateral nonlinearities

For unilateral nonlinearities such as the one studied in Sections 2.3 and 2.4, the procedure of selection of diameters harmonic-wise can help further reducing of the size of the problem. The additional criterion and its implementation are provided in this section, however the underlying mathematical developments are left in Appendix A.

The additional reduction consists in taking advantage of how the different harmonics of the nonlinear forces couple themselves. Indeed, when certain conditions described hereafter are met, it is observed that all diameters among the subset $D$ determined in Section 2.4 do not respond on all harmonics. A distinction comes from the parity of the harmonic number: a subset $D^e$ responds only on even harmonics while a subset $D^o$ responds only on odd harmonics. The sets $D^e$ and $D^o$ respect the relations

$$D^e \cup D^o = D \text{ and } D^e \cap D^o = \emptyset. \qquad (43)$$

This means that all diameters in $D$ are needed, but only a subset $D^e$ or $D^o$ can be selected according to the parity of the harmonic.

In order to know if this parity criterion can be used, one should examine the value of $Kh_{ex}$ where the definition of $K$ is given in Eq. (32) and $h_{ex}$ is the nodal diameter of the excitation. Two cases arise:

if $Kh_{ex} \equiv 0$ [$N$]: no reduction is expected and all harmonics should use the whole $D$ set such as

$$D^k = D \quad \forall k \in [0, N_h], \qquad (44)$$

if $Kh_{ex} \equiv \frac{N}{2}$ [$N$]: the reduction according to the parity of the harmonic number can be used and it comes that

$$D^k = D^e \text{ if } k \text{ is even,} \qquad (45)$$

$$D^k = D^o \text{ if } k \text{ is odd,} \qquad (46)$$

The sets $D^e$ and $D^o$ are built such as $D^e$ contains all $d_q$ associated to even values of $q$ in $v_K$ while $D^o$ is composed of all $d_q$ associated to odd values of $q$ in $v_K$.

Overall, this parity criterion allows to further divide the size of the system by a factor 2 when the condition $Kh_{ex} \neq 0$ [$N$] is respected. It is noted that the demonstration of Appendix A shows that this condition can only occur when $N$ is even.

The fact that this supplementary reduction only applies to cyclic systems with an even number of sectors actually lies in the fact that a certain symmetry can be achieved in the response (allowing to consider only half of the sectors) over the bladed disk when the non-degenerated nodal diameter $\frac{N}{2}$ is included in the basis. Since systems with odd numbers of sectors do not have this nodal diameter in their basis, this symmetry can never be achieved and thus no further reduction can be obtained.

Due to the necessity of splitting the diametral basis among different harmonics, this method only applies to an HBM solving process. Throughout the rest of the article, this criterion is used for the simulations when it is relevant and applicable. As a way of visualizing the reduction procedure, an illustration of the concept of selecting a subset of sectors $S$ with and without the parity criterion is proposed in Fig. 1 where the sectors colored in blue represent the set $S$. The case with a number of sectors $N = 8$ and $h_{ex} = 2$ chosen. In this configuration, $K = 2$ and the parity criterion holds because $Kh_{ex} = 4 \equiv N/2$. First of all, without considering the parity criterion the reduction method predicts that it is needed to sample the nonlinear forces over sectors 1 to 4. However, the fact that the parity criterion is met actually means that the information over sectors 3 and 4 (■) is redundant with the information over sectors 1 and 2 (■) which allows to account only for the latter in the reduced problem.

### 3.4. Evaluation of nonlinear forces: alternating Frequency/Time-Cyclic/Physical

This section details the strategy retained for the evaluation of the nonlinear forces throughout the HBM resolution. Two different contact laws are used in the article and both rely on the alternating frequency/time procedure (AFT) [29]. Firstly, the regularized penalty law (referred to as Regularized-HBM (R-HBM)) that was used for the analytical developments in Section 2.3 is employed to compute nonlinear responses. Secondly, a stiffer contact law relying on augmented lagrangian [27,28], denominated Dynamic Lagrangian Frequency/Time (DLFT-HBM), is also used. Both contact laws constitute the state-of-the-art in terms of computing unilateral contact dynamics through HBM [20,28], especially for blade-tip/casing contacts [11,20,30]. The DLFT-HBM strategy is, unlike the R-HBM, a nonsmooth strategy which makes it more numerically sensitive but also more accurate. The use of this second method to evaluate the contact forces allows to put to the test the criterions derived in Sections 2.4 and 3.3. Then, this article aims to demonstrates that the criterions not only work for a regularized contact law, but that they can also be extended to the nonsmooth and non-differentiable contact laws that approximate unilateral contact.

In order to compute the nonlinear forces in a HBM framework, it is oftentimes necessary to resort to the AFT scheme [29]. In this work, a supplementary step is required since the HBM unknown is composed of spectral components. Indeed, the cyclic components
<figure><img src="image_3.png" /><figcaption>Fig. 1. Illustration of the choice of S for N = 8 and h<sub>ex</sub> = 2 with and without the parity criterion.</figcaption></figure>

should be transformed to the physical space for the computation of nonlinear forces and back to the cyclic space for the evaluation of the residual at each iteration. This new alternation type of is denominated Alternating Frequency/Time-Cyclic/Physical (AFT-CP).

The transformation between the frequency and time domain is done by using the real inverse discrete Fourier transform (IDFT) matrix $F^{-1}$ so that $\{\bullet(t_i)\}_{i=1...N_t} = F^{-1}\tilde{\sigma}$. Whereas, the transformation from the time domain to the frequency domain is made by means of the real discrete Fourier transform (DFT) matrix $F$ and the relation between both domains is $\tilde{\sigma} = F\{\bullet(t_i)\}_{i=1...N_t}$. A time period is sampled into $N_t$ (user-defined parameter) instants. It should be chosen as a compromise between computation time and accuracy of the procedure: increasing $N_t$ favors convergence and prevents aliasing problems but increases computational costs [31].

In this work, the link between the cyclic and physical spaces is made in the frequency domain. It relies on $\tilde{G}$ and $\tilde{\overline{G}}$ that are built out of sub-matrices of the spatial Fourier transform matrix $S$. The construction of these matrices requires to determine in advance the sets $S^k$ (defined by the integers $j_k^{\max}$) associated with each $D^k$ following the rule presented in Section 2.2. In practice, since the evaluation of nonlinear forces in the time domain on a set of sectors requires all harmonics for all sectors, the set of sectors over which nonlinear forces are computed is actually $S = (\bigcup_{k=0}^{N_h} S^k)$. Then, the matrices $\tilde{G}$ and $\tilde{\overline{G}}$ are defined as

$$ \tilde{G} = \text{blkdiag}(\tilde{G}^k) \otimes \mathbf{I}_{n_{nl}} \quad \text{and} \quad \tilde{\overline{G}} = \text{blkdiag}(\overline{\tilde{G}}^k) \otimes \mathbf{I}_{n_{nl}} \quad \text{with } k \in [0, N_h] \quad (47) $$

where $n_{nl}$ is the number of nonlinear DOF of a single sector. Hereafter, the notations are consistent with what is defined in Section 2.2. The inverse transformation (from cyclic to physical variables) is straightforward and

$$ \tilde{\overline{\mathbf{G}}}^k = \begin{bmatrix} (\mathbf{S}_S^{D^k})^\top & \mathbf{0} \\ \mathbf{0} & (\mathbf{S}_S^{D^k})^\top \end{bmatrix} \quad \text{if } k > 0 \text{ or } \tilde{\mathbf{G}}^k = (\mathbf{S}_S^{D^k})^\top \quad \text{if } k = 0. \quad (48) $$

However, for the forward transformation, since $S$ contains a larger number of sectors than each individual harmonic set $S^k$, some of the information contained over sectors of $S$ is redundant and should be truncated to switch from physical to cyclic variables. This truncation from $S$ to $S^k$ is made by the matrix $\mathbf{R}_{S^k}$ defined as

$$ \mathbf{R}_{S^k} = [\mathbf{I}_k^{\max} \quad \mathbf{0}] \quad \text{if } S^k \neq S \quad \text{and} \quad \mathbf{R}_{S^k} = \mathbf{I}_k^{\max} \quad \text{if } S^k = S. \quad (49) $$

Following this definition, it comes that the direct transformation $\tilde{G}$ is built out of blocks $\tilde{\mathbf{G}}^k$ such as

$$ \tilde{\mathbf{G}}^k = \begin{bmatrix} \left[ (\mathbf{S}_{S^k}^{D^k})^\top \right]^{-1} \mathbf{R}_{S^k} & \mathbf{0} \\ \mathbf{0} & \left[ (\mathbf{S}_{S^k}^{D^k})^\top \right]^{-1} \mathbf{R}_{S^k} \end{bmatrix} \quad \text{if } k > 0 \text{ or } \tilde{\mathbf{G}}^k = \left[ (\mathbf{S}_{S^k}^{D^k})^\top \right]^{-1} \mathbf{R}_{S^k} \quad \text{if } k = 0. \quad (50) $$

Note that the block $\left[ (\mathbf{S}_{S^k}^{D^k})^\top \right]$ is invertible since sets $S^k$ and $D^k$ were chosen so that this sub-matrix of $S$ is square. The distinction with respect to $k$ comes from the presence of cosine and sine terms in real HBM for non-zero harmonic numbers $k$.
Now that all matrices are defined for the alternation between the different spaces, the AFT-CP methodology is presented hereafter. Say that an estimate of the solution is known in the shape of a multiharmonic vector corresponding to spectral components $\tilde{\mathbf{x}}$, the overall alternating strategy is composed of five steps:

1. the vector $\hat{\mathbf{x}}$ is transformed to the physical space in the frequency domain by means of $\tilde{\mathbf{x}} = \tilde{\mathbf{G}}\hat{\mathbf{x}}$. At this stage, $\tilde{\mathbf{x}}$ contains all real Fourier coefficients of the nonlinear DoF of the system on sectors $j \in [1, j^{\max}]$ where $j^{\max} = \max_{k \in [0, N_h]} j_k^{\max}$,

2. the multiharmonic vector of displacements is transformed to the time domain by using the IDFT matrix such that $\mathbf{x}(t) = \mathbf{F}^{-1}\tilde{\mathbf{x}}$ and $\dot{\mathbf{x}}(t) = \mathbf{F}^{-1}\nabla\tilde{\mathbf{x}}$ where $\nabla$ is the frequency-domain derivation operator,

3. nonlinear forces are evaluated in the time domain on physical sectors according to the expression of the nonlinearity $\mathbf{f}_{nl}(\mathbf{x}(t), \dot{\mathbf{x}}(t))$,

4. the time-dependent nonlinear forces are transformed back to the frequency domain by using the DFT matrix such as $\tilde{\mathbf{f}}_{nl} = \mathbf{F}^{-1}\mathbf{f}_{nl}(t)$,

5. the frequency domain physical nonlinear forces are then used to compute the spectral components following $\hat{\mathbf{f}}_{nl} = \tilde{\mathbf{G}}\tilde{\mathbf{f}}_{nl}$.

This type of alternation actually nests the classical AFT from steps 2 to 4. Steps 1 and 5 are added in order to perform the back and forth transformations between cyclic and physical domains. Note that, in this exploratory work the solutions are computed regardless of their stability. Stability analyses of periodic solutions for blade-tip/casing contacts were recently explored [22], but the implementation of these methods exceeds the scope of this paper.

### 3.5. Special case of traveling wave excitations

In this section, the particular case of traveling wave excitation is dealt with. For this type of loading, further reductions are possible. This reduction method relies on the hypothesis that the solution displacement field of the system has of the same spatial decomposition than the excitation. This idea was initially proposed by Petrov [6]. In this work, Petrov's method is assessed on the industrial application featuring blade-tip/casing interactions and serves as a basis of comparison for the more generic method proposed here.

In the following, a traveling wave excitation of order $h_{ex}$ of the form

$$ \mathbf{f}_{ex,j}^{TW}(t) = f_{ex} \cos \left( t + \frac{\alpha(j-1)h_{ex}}{\omega} \right) \quad \forall j \in [1, N] \qquad (51) $$

is applied to a cyclic system. The excitation is defined in the fixed frame such as it can model stator wakes. Then, assuming the bladed disk is rotating in the same direction as blades are numbered, the fixed frame excitation propagates backwards in terms of blade numbering: 1, N, N - 1, ...

In the case of such an excitation, the subsets $D^k$ should be chosen by reading the zig-zag diagram, such as

$$ D^k = \{ \mathcal{Z}(kh_{ex}) \} \qquad (52) $$

where $\mathcal{Z}(h)$ is the function that stands for the zig-zag diagram. In the particular case of traveling waves, the spectral indexes in $D^k$ are signed, allowing to know the specific propagation direction of the traveling waves.

If there exists a spectral index $-m \le q < 0$ such that $(q-h) \equiv 0[N]$ then $\mathcal{Z}(h) = q$, otherwise $\mathcal{Z}(h) = q$ where $q$ respects $(q+h) \equiv 0[N]$. Both conditions are exclusive and only a single spectral component $-m \le q \le m$ can be obtained by this method.

For the case of the traveling wave hypothesis, the excitation should be a pure traveling wave hypothesis, otherwise the strategy is not applicable. For all cases where the external loading contains a disturbance or a standing wave component, the nodal diameter coupling strategy proposed in Section 2.4 should be used to ensure that all necessary cyclic components are accounted for. This sensitivity is assessed on the industrial model for a disturbed traveling excitation shape in Section 5.4.

## 4. Results on a simplified bladed disk model

This section presents the results obtained on a simplified model composed of a cyclic chain of oscillators that undergo unilateral contacts. This model is approximating a bladed disk that experiences blade-tip/casing contacts. The demonstration of the nodal diameters selection method proposed in this article is made on several blade numbers N as well as several loading scenarios for both the regularized and the nonsmooth contact law. In all cases, the arc-length continuation procedure [29] is employed to compute the forced responses that feature numerous turning points.

### 4.1. Presentation of the model

The model of simplified bladed disk is adapted from the literature [32]. Each sector is composed of 5 masses linked with stiffnesses. A single sector of this model is illustrated in Fig. 2. The links between the different sectors are made through the stiffnesses between $m_5^j$ and $m_4^{j+1}$ where j denotes the sector number. All the stiffnesses and masses of a sector are given in Table 1.

The excitation force on each sector is applied to the mass $m_1$ that models the tip of the blade. Modal damping is used such that all modes of the system are damped with the same modal damping ratio of $\xi = 1 \cdot 10^{-3}$. For this model, no pre-loading is applied to
<figure><img src="image_4.png" /><figcaption>Fig. 2. Simplified model of a bladed disk (sector *j*).</figcaption></figure>

Source: Adapted from Ref. [32].

<table><caption>**Table 1**
Structural parameters of a sector.</caption><thead><tr><th></th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th></tr></thead><tbody><tr><td>m (kg)</td><td>0.2</td><td>0.3</td><td>0.4</td><td>1.2</td><td>1.2</td></tr><tr><td>k (10<sup>6</sup>N m<sup>-1</sup>)</td><td>2</td><td>1</td><td>40</td><td>0.6</td><td>50</td></tr></tbody></table>

the system: $\mathbf{f}_0 = \mathbf{0}$. The number of sectors of the system is parameterized by the variable $N$ allowing to study the dynamics of the system with different number of sectors. This allows to validate the rule derived in Section 2.4 for multiple values of $N$.

For the HBM simulations, the number of harmonics is set to $N_h = 10$ and the number of time instants is $N_t = 10000$. For R-HBM simulations, $\kappa = 1 \cdot 10^8$ N m$^{-1}$ and $\gamma = 2.5$ N. For DLFT-HBM simulations, a wall of stiffness $\kappa = 1 \cdot 10^8$ N m$^{-1}$ is added and $\epsilon = 1 \cdot 10^9$ N m$^{-1}$. Both contact laws are then simulating the dynamics of the model in Fig. 2 despite the difference in contact treatment.

## 4.2. Contact scenario

The nonlinear DoF are all the masses $m'_1$ with $j \in [1, N]$ that are impacting a wall at a given distance $g_{0,j}(t)$. Unless otherwise mentioned, the gap function is constant along the period with $g_{0,j}(t) = c_j = 1 \cdot 10^{-3}$ m for all sectors, making the casing purely circular like displayed in Fig. 3(a). In some contact scenarios, the casing is chosen to be a rigid diametral shape, such that

$$g_{0,j}(t) = c_j + p_j \cos \left( h_c \Omega \left( t - j \frac{2\pi}{N\Omega} \right) \right) \quad (53)$$

where $c_j$ is the mean value of the gap function of node $j$, $p_j$ is the amplitude of oscillation of the gap function around $c_j$, $\Omega$ is the rotation speed (synchronous with the excitation frequency $\omega$) and $h_c$ is the nodal diameter index of the casing. An example of diametral casing is displayed in Fig. 3(b) as an illustration for $h_c = 5$. In order to initiate contact interactions, the system is also excited for several engine orders by standing waves excitations. The definition of a traveling wave excitation of order $h_{ex}$ is Eq. (51) and a standing wave excitation is the sum of two traveling wave excitations of the same order traveling in the opposite directions, such as

$$\mathbf{f}_{\text{ex},j}^{\text{SW}}(t) = f_{\text{ex}} \left[ \cos \left( t + \frac{\alpha(j-1)h_{\text{ex}}}{\omega} \right) + \cos \left( t - \frac{\alpha(j-1)h_{\text{ex}}}{\omega} \right) \right] \quad \forall j \in [1, N]. \quad (54)$$

## 4.3. Even number of sectors: $N = 24$

The first application chosen for the presented methodologies is the case study of Fig. 2 with $N = 24$ sectors. Since $N$ is even, this case is suitable for the presentation of the parity criterion for unilateral nonlinearities in HBM (see Section 3.3). Following this
<figure><img src="image_5.png" /><figcaption>Fig. 3. Casing shapes and contact scenario.</figcaption></figure>

<table><caption>**Table 2**
Reduction for the case of the simplified system with $N = 24$ sectors for 10 harmonics (full size: $N \times n_{nl} \times (2N_h + 1) = 504$ unknowns).</caption><thead><tr><th>Engine order</th><th>0</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th><th>11</th><th>12</th></tr></thead><tbody><tr><th>Value of K</th><td>1</td><td>12</td><td>6</td><td>4</td><td>3</td><td>12</td><td>2</td><td>12</td><td>3</td><td>4</td><td>6</td><td>12</td><td>1</td></tr><tr><th>Reduced size (DoF)</th><td>21</td><td>252</td><td>126</td><td>84</td><td>63</td><td>252</td><td>42</td><td>252</td><td>63</td><td>84</td><td>126</td><td>252</td><td>21</td></tr><tr><th>Sectors for AFT-CP</th><td>1</td><td>12</td><td>6</td><td>4</td><td>3</td><td>12</td><td>2</td><td>12</td><td>3</td><td>4</td><td>6</td><td>12</td><td>1</td></tr><tr><th>Reduction (%)</th><td>95.8</td><td>50.0</td><td>75.0</td><td>83.3</td><td>87.5</td><td>50.0</td><td>91.7</td><td>50.0</td><td>87.5</td><td>83.3</td><td>75.0</td><td>50.0</td><td>95.8</td></tr></tbody></table>

criterion, the expected reduction in terms of total size and number of sectors to evaluate inside the AFT-CP procedure are presented in Table 2.

For this case, almost all engine orders benefit from the 50% reduction allowed by the criterion on the parity of Section 3.3. Only $h_{ex} = 8$ does not fall into this category: $K = 3$ and $Kh_{ex} \equiv 0$ [N], thus, no reduction is achieved on the parity of the harmonics. Overall, the mean over all engine orders of the expected reduction is 75%. This number is significantly impacted by the parity criterion and would be of only 51.3% without it. Note that reductions presented in Table 2 do not take into account the type of excitation that is applied to the system. Indeed, by using the traveling wave hypothesis formulated in Section 3.5, the achieved reduction is higher. As a reminder, when considering the traveling wave hypothesis, the size of the system (and the AFT-CP) is brought down to the size of a single sector, yielding a reduction of 95.8% for any engine orders. Note that the value of K is systematically equal to the numbers of sectors needed for the calculation. This means that K has a physical meaning related to a relation of symmetry along the blades of the cyclic system.

For this case study, as a way of demonstrating all the reduction capacities of the HBM method for different loading scenarios, several types of calculations are made. Computations are carried out for standing wave excitations in Section 4.3.1 and for a mixed contact scenario involving a traveling wave excitation and a casing deformation in Section 4.3.2. For all cases, the result obtained with the proposed reduction is always compared to the full calculations (with all sectors) that serve as reference. The latter is directly built by assembling all sectors without the use of cyclic symmetry. Note that a perfect match was observed between full calculations and the traveling wave hypothesis reduction, but the results are not shown for the sake of brevity.

### 4.3.1. Standing wave excitation

In this section, simulations are carried out on the nonlinear dynamics of the system under a standing wave excitation. This type of excitation allows to observe how the nonlinearity couples the different nodal diameters together. Hence, it is possible to verify that the analytical developments made in Sections 2.4 and 3.3 are relevant and actually are exact reductions of the problem. Standing wave excitations are applied to the system for $h_{ex} = 1$ with $f_{ex} = 2$ N, and for $h_{ex} = 4$ with $f_{ex} = 1$ N.

Following the rules of determination of the nodal diameters to account for in each harmonic, the basis chosen are

$$ \begin{cases} D^k = \{0, 2, 4, 6, 8, 10, 12\} & \text{if } k \text{ is even} \\ D^k = \{1, 3, 5, 7, 9, 11\} & \text{if } k \text{ is odd} \end{cases} \quad \forall k \in [0, N_h] \qquad (55) $$

for $h_{ex} = 1$ ($K = 12 \Rightarrow Kh_{ex} = 12 \equiv \frac{N}{2}$ [N]) and

$$ \begin{cases} D^k = \{0, 8\} & \text{if } k \text{ is even} \\ D^k = \{4, 12\} & \text{if } k \text{ is odd} \end{cases} \quad \forall k \in [0, N_h] \qquad (56) $$
<figure><img src="image_6.png" /><figcaption>Fig. 4. Standing wave excitation on the simplified bladed disk with $N = 24$ sectors for various engine orders, linear (……), R-HBM: cyclic symmetry (—), full (---), DLFT-HBM: cyclic symmetry (—), full (---).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 5. Cyclic components over all harmonics at the resonance point of the full calculation in Fig. 4(b), standing wave with $h_{ex} = 1$. Nodal diameter: 0 (■), 1 (■), 2 (■), 3 (■), 4 (■), 5 (■), 6 (■), 7 (■), 8 (■), 9 (■), 10 (■), 11 (■), 12 (■).</figcaption></figure>

for $h_{ex} = 4$ ($K = 3 \Rightarrow Kh_{ex} = 12 \equiv \frac{N}{2} [N]$). The results of the different HBM computations are then displayed for all sectors in Fig. 4 for both contact laws and values of $h_{ex}$.

Once again, a perfect match is observed between the full calculations and the reduced calculations with the proposed approach. The standing wave nature of the response of the system can be seen to the presence of curves of different amplitude (linear and nonlinear), related to the level of forcing that a given sector is actually experiencing due to the standing wave excitation. Similarly as for the case of traveling wave excitations, the cyclic contributions to the displacement field at the resonance for the full DLFT calculation is displayed in Fig. 5 for $h_{ex} = 1$ and in Fig. 6 for $h_{ex} = 4$.
<figure><img src="image_8.png" /><figcaption>Fig. 6. Cyclic components over all harmonics at the resonance point of the full calculation in Fig. 4(d), standing wave with $h_{ex} = 4$. Nodal diameter: 0 (■), 4 (■), 8 (■), 12 (■).</figcaption></figure>

For the case of $h_{ex} = 1$ of Fig. 5, the basis was correctly identified a priori in Eq. (55). This case demonstrates that even if no multiplicity is observed between $h_{ex}$ and $N$, the HBM can take advantage of the criterion on the parity of the harmonics to still reduce the computational burden of 50%. This can be observed in Fig. 5 where all ND are involved in the whole solution, however one half only responds on even harmonics and the other half on odd harmonics. The sampling strategy of AFT-CP is then able to perform the computations on only half of the sectors, yielding a significant gain in terms of computational time (50%). This type of strategy would not be possible for time-integration strategies since all cyclic components would have been needed to compute the response.

Similarly, the cyclic contributions responding in Fig. 6 for the case $h_{ex} = 4$ were correctly identified in Eq. (56) and the split among even and odd harmonics is also highlighted. For this couple ($N, h_{ex}$) a reduction of 87.5% is reached, and the results are strictly equivalent to the full model. This section has demonstrated the relevance of the method used to determine the nodal diameter couplings as a way of reducing the computational effort.

For both Figs. 6 and 5, one can note that the response is purely stationary, meaning that forward and backward contributions are strictly equal. Hence, when such responses emerge (standing wave excitations) in the rest of the article, only the first graph is displayed for the sake of brevity.

### 4.3.2. Mixed contact scenario

Finally, the system's response is analyzed under a loading scenario that involves both the excitation through a traveling wave and also a nodal diameter shape of the surrounding casing, as presented in Fig. 3(b). When considering them independently, the excitation and the casing both induce a traveling wave response in cyclic system. However, if $h_{ex} \neq h_c$ a coupling between the two types of loading happens and a partly standing wave response is expected to emerge. This scenario could not have been studied without a specific criterion in cyclic symmetry. Indeed, Petrov's method [6] is restricted to a single traveling wave loading.

The chosen configuration is a case where the basis associated to $h_{ex}$ and $h_c$ are different. Indeed, trivial cases occur when one of the two basis is included in the other ($D_{ex}^k \subset D_c^k$ or $D_c^k \subset D_{ex}^k$). In the non trivial case, one should consider the basis that is formed by the greatest common divisor (GCD) of $h_{ex}$ and $h_c$. In order to demonstrate this property, a case with $h_{ex} = 6$ and $h_c = 8$ is studied. The basis for each of the two are respectively

$$ \begin{cases} D^k = \{0, 12\} & \text{if } k \text{ is even} \\ D^k = \{6\} & \text{if } k \text{ is odd} \end{cases} \quad \forall k \in [0, N_h] \qquad (57) $$

for the excitation and

$$ D^k = \{0, 8\} \quad \forall k \in [0, N_h] \qquad (58) $$

for the casing, no parity criterion exists for an excitation over 8ND. In order to cover all the possible couplings that can happen, one should choose the basis that would emerge from a 2ND excitation, such as

$$ D^k = \{0, 2, 4, 6, 8, 10, 12\} \quad \forall k \in [0, N_h]. \qquad (59) $$

It should be noted that this basis not only contains the nodal diameters in Eqs. (57) and (58) but also supplementary nodal diameters that are expected to appear through the coupling of both basis (4ND and 10ND). Also, the parity criterion has not been used for this expression. Then, the basis Eq. (59) is used since it includes all couplings generated by Eqs. (57) and (58). The results of the simulations are displayed in Fig. 7 for both DLFT-HBM and R-HBM on two sectors.

Firstly, the forced responses obtained in cyclic symmetry match the full calculations for both methods. The responses of DLFT-HBM and R-HBM are not exactly identical because of the difference in contact treatment. This configuration being a highly nonlinear very sensitive case study, it highlights the discrepancies. The standing wave nature of the displacement field can be seen by the fact that all sectors do not experience the same response. This fact can be seen more clearly in Fig. 8 where the cyclic components are displayed.
<figure><img src="image_9.png" /><figcaption>Fig. 7. Traveling wave excitation on the simplified bladed disk with $N = 24$ sectors for $h_{ex} = 6$ and a casing deformed with $h_c = 8$, linear (......), R-HBM: cyclic symmetry (—), full (---), DLFT-HBM: cyclic symmetry (—), full (---).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8. Cyclic components over all harmonics at the resonance point of the full calculation in Fig. 7, traveling wave with $h_{ex} = 6$ and casing with $h_c = 8$. Nodal diameter: 0 (●), 2 (■), 4 (■), 6 (■), 8 (■), 10 (■), 12 (■).</figcaption></figure>

The decomposition of the solution at the resonance in terms of cyclic components reveals that the basis Eq. (59), emerging from the coupling of the excitation and the casing, was correctly identified a priori. Indeed, the only nodal diameters that are involved in the response are in the basis. Moreover, as discussed beforehand, a parity criterion is observed on the response: nodal diameters {0, 4, 8, 12} respond on even harmonics while nodal diameters {2, 6, 10} intervene only on odd harmonics. Detecting the coupling of the casing and the excitation in terms of parity would have then allowed to divide the computational effort by two in this particular case. Moreover, it is observed that the couplings of Eqs. (57) and (58) make appear 4ND and 10ND contributions in the simulations even though it is in neither of the two basis, justifying that the basis should be built according to the GCD of $h_{ex}$ and $h_c$ in order to cover all possible couplings.
<figure><img src="image_11.png" /><figcaption>**Fig. 9.** ECL5/Catana fan stage model, cyclic symmetry boundary (■), clamped (■), axial blockage (■), nonlinear nodes (□), observation node for NFRC (□), excitation node (□).</figcaption></figure>

The test case studied in depth here was featuring *N* = 24 sectors and was prone to large reductions in all types of loading scenarios. Indeed, *N* being even ensures that all even values of *h*<sub>ex</sub> will provide a significant reduction of around at least 50%.

This section compared full calculations with the reduced methodology proposed in this article. The results are in perfect
agreement and thus these numerical simulations give further credibility to the theoretical demonstration. Notice that the numerical
simulations showed that the reduction strategy also works for non-regularized contact law. The arbitrary case of *N* = 24 was dealt
with here, but the proposed methodology can be applied for a wide number of cases. In order to cover all configurations, an odd
number of sectors of *N* = 45 is also dealt with in Appendix B. The latter is dealt with for the sake of completeness as the reduction
rates are quite lower for systems with odd numbers of sectors.

**5. Results on an industrial model**

Now that the reduction capabilities of the proposed strategy was extensively demonstrated on academic examples, the study is
taken further to tackle a more realistic industrial configuration. The same range of loading scenarios as for the academic case is
applied to demonstrate the ability of the method to be employed on large-scale systems involving blade-tip/casing contacts.

**5.1. Presentation of the model**

The model chosen to demonstrate the capabilities of the strategy is the open test case fan stage ECL5/Catana designed at Ecole
Centrale de Lyon [33,34] to be a step towards more representative open numerical models of ultrahigh bypass ratio architectures.
The stage is composed of 16 blades built out of unidirectional carbon fibers and epoxy composite plies. The disk is made of titanium.
This model was recently the subject of experimental [35] and numerical [36] concerning structural dynamics. The full finite-element
mesh is composed of 87,310 nodes for a total of 19,642 elements. It is displayed in Fig. 9(a) for information. Both the material data
and the boundary conditions are those given in Ref. [33]. The system is considered in cyclic symmetry and the cyclic boundaries,
as long as other boundary conditions, are displayed in Fig. 9(b). In this work, the contacts between the blades and the disk are
linearized and no relative motion is possible. The Coriolis effect is neglected throughout this study since blade-dominant modes
are studied.

In order to decrease the computational costs of simulating the dynamics of this system, a Craig–Bampton reduced-order modeling
technique [37] is applied to the system for each nodal diameter. A total of 11 boundary nodes are retained in the reduced-order
model: $n_{nl} = 9$ boundary nodes are kept at the tip of the blade at the intrados side for contact treatment purposes, one node is kept
at the blade-tip at mid chord in order to apply the excitation and one observation node is retained in order to monitor the bending
amplitude away from the nonlinear interface. The contact nodes are evenly spaced along the blade-tip. The displacement of the
observation node in the circumferential direction is systematically used throughout the article in order to plot the forced response
curves. 20 internal modes are retained for the Craig–Bampton procedure. Since the system is formulated using the cyclic symmetry,
the number of unknowns associated to non-degenerated nodal diameters (OND and 8ND) is 53 and 106 for all other diameters.
Modal damping is used to build the damping matrix. The first five modes are damped with a ratio $\xi_{1-5} = 1 \cdot 10^{-3}$ and the remaining
modes with a ratio of $\xi_{5+} = 5 \cdot 10^{-3}$. The mode of interest in this study is the first bending mode.

Given the number of sectors of the ECL5 fan stage of *N* = 16, the expected reductions with the proposed methodology are
provided in Table 3. Since *N* is even, it is subject to the parity criterion, and parity reduction is predicted along all nodal diameters.
<caption>**Table 3**</caption>

Reduction for the case of the industrial system with N = 16 sectors for 10 harmonics (full size: N × 3n<sub>nl</sub> × (2N<sub>h</sub> + 1) = 9072 unknowns).

<table><thead><tr><th>Engine order</th><th>0</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th></tr></thead><tbody><tr><th>Value of K</th><td>1</td><td>8</td><td>4</td><td>8</td><td>2</td><td>8</td><td>4</td><td>8</td><td>1</td></tr><tr><th>Reduced size (DoF)</th><td>567</td><td>4536</td><td>2268</td><td>4536</td><td>1134</td><td>4536</td><td>2268</td><td>4536</td><td>567</td></tr><tr><th>Sectors for AFT-CP</th><td>1</td><td>8</td><td>4</td><td>8</td><td>2</td><td>8</td><td>4</td><td>8</td><td>1</td></tr><tr><th>Reduction (%)</th><td>93.8</td><td>50.0</td><td>75.0</td><td>50.0</td><td>87.5</td><td>50.0</td><td>75.0</td><td>50.0</td><td>93.8</td></tr></tbody></table>

## 5.2. Contact scenario

In this article, the dynamics of the casing at the contact interface is neglected. The casing is considered to be a rigid mathematical profile. The casing is located at a distance $c_j = 4 \cdot 10^{-4}$ m of each boundary node after accounting for the pre-loading forces $f_0$. The latter are centrifugal forces that deform the blade into its working conditions. The system is rotating at an angular speed of $\Omega = 1,151.9$ rad s$^{-1} = 11,000$ rpm. These characteristics are chosen to reflect the same order of magnitude that was experimentally measured on the PHARE#2 test rig [38]. However, blade-tip/casing contacts were not observed during the experimental work. Thus, the blade-tip contact is numerically triggered by decreasing the nominal gap and applying a sufficient level of excitation to the system.

The excitation is applied to a node at the tip of the blade on each of the 16 blades in the circumferential direction in order to study the response of bending modes. The excitation is considered to be synchronous with the rotation of the blade ($\omega = \Omega$). The contact scenario considered for this study is similar as the one that was considered on the academic model in Section 4.2.

When contact interactions occur, dry friction forces are applied to the blade. Given the high relative speeds at stake in such interactions, sliding is assumed at all times. This hypothesis was validated several times by correlating numerical studies with experimental data [39,40]. The friction forces are expressed on the jth node at time $t_j$, such as

$$f_{nl,j}^{\theta}(t_j) = -\mu \frac{v_j^{\theta}(t_j) + \rho_j \Omega}{\sqrt{(v_j^{\theta}(t_j) + \rho_j \Omega)^2 + (v_j^z(t_j))^2}} f_{nl,j}^{N}(t_j) \quad (60)$$

for the circumferential direction and

$$f_{nl,j}^{z}(t_j) = -\mu \frac{v_j^{z}(t_j)}{\sqrt{(v_j^{\theta}(t_j) + \rho_j \Omega)^2 + (v_j^z(t_j))^2}} f_{nl,j}^{N}(t_j) \quad (61)$$

for the axial direction. $v_j^\theta$ and $v_j^z$ are the circumferential and axial components of the vibrational velocity, $\mu = 0.15$ is the dry friction coefficient and $\rho_j$ are the radial distances between each of boundary node and the rotation axis. Note that, since R-HBM uses a regularized contact law, the residual normal force associated with the regularization induces a residual friction force as well.

For all contact laws, since permanent sliding is assumed, the tangent force expressed in the friction direction is just the product of the normal contact forces $f_{nl}^N$ and the friction coefficient $\mu$. In Eqs. (60) and (61), a projection of the tangent force in the cylindrical frame is performed depending on the rotation speed and the vibrational velocities of the nodes. Even though the expressions are nonlinear with respect to $v_j^z$ and $v_j^\theta$, the friction direction is highly driven by the rotation speed of the whole system. Indeed, for the given application $\rho_j \Omega \gg v_j^z$ and $\rho_j \Omega \gg v_j^\theta$. Under this hypothesis, it comes that $f_{nl,j}^\theta \approx \mu f_{nl,j}^N$ and $f_{nl,j}^z \approx 0$ and thus the friction forces being proportional to normal forces it does not create any supplementary couplings of nodal diameters. Indeed, adding friction only changes the direction of the resultant contact force from $f_{nl} = [1 \ 0 \ 0]$ to $f_{nl}^N$ to $f_{nl} = [1 \ -\mu \ 0]$ $f_{nl}^N$ which does not modify the nonlinear expression of the term that creates the couplings $f_{nl}^N$, but only the direction of the contact force. Following this reasoning, friction forces are assumed to create no further couplings than what was already identified in Section 2.4 and the same procedure is used for the industrial model than what was used in Section 4. This assumption is verified numerically in the next section through the comparison between the reduced model and the full system. It is important to note that it is unrealistic that the contact reaches a sticking state with the high relative velocities at stake and thus the coupling criterion should not encompass Coulomb's Law in order to maximize the reduction capabilities of the method. The application studied here is vibro-impact dominated and one should refer to Ref. [9] to tackle friction dominated dynamics applications. The latter would only add unneeded complexity for blade-tip/casing contacts.

For this study, $N_h = 10$ harmonics and $N_t = 2000$ instants are used for the HBM. For DLFT-HBM, the value $\epsilon = 1 \cdot 10^7$ N m$^{-1}$ is used. The friction coefficient is set to $\mu = 0.15$. In this section, only the simulations run with DLFT-HBM are shown. The simulations were also run with the R-HBM and provided near identical results. They are not presented here for the sake of brevity.

This contact scenario is not necessarily representing a real life occurrence of blade-tip/casing contact interactions, however the overly stiff setup (rigid casing) is used to back the methodology on the most challenging configuration possible regarding contact. Moreover, the use of a simple harmonic force on the blade to initiate contact is chosen in order not to introduce unneeded complexity for the validation of the coupling criterion. Indeed, it is reminded that the proposed reduction methodologies of Sections 2.4 and 3.3 are not restricted to this contact scenario, as shown by the results on the academic example in Section 4, and that the strategy could still be applied to any synchronous loading (casing deformation or any external excitation).
<figure><img src="image_12.png" /><figcaption>Fig. 10. Traveling wave excitation on the Catana industrial blade model for various engine orders, linear (……), DLFT-HBM: cyclic symmetry (—), full (---).</figcaption></figure>

## 5.3. Traveling wave excitations

The first case study for this industrial model is the configuration where the casing is purely circular and the system is excited by a traveling wave excitation. In order to demonstrate the applicability of the traveling wave hypothesis formulated in Section 3.5, the cases of $h_{ex} = 1$ and $h_{ex} = 4$ (odd and even) are dealt with. Like for the studies in Section 4, the results are obtained through DLFT-HBM for the cyclic symmetric model as well as the full model. For these values of $h_{ex}$, the forcing amplitude is $f_{ex} = 2$ N and the results are displayed in Fig. 10.

First of all, the same response is obtained through both methodologies, thus consolidating the relevance of both methodologies to tackle blade-tip/casing contact problems. Moreover, as for the academic case in Section 4, there is an exact match between the full computations and the reduced calculations made with the traveling wave hypothesis formulation. The latter yields a 93.75% reduction, allowing to perform calculations in a few minutes, where the full calculations usually last for between 5 to 15 h, depending on the choice of continuation parameters, contact treatment algorithm and other HBM parameters. Indeed, the scaling of the reduction rate with the computational time is not linear: numerous operations made during the solving process scale in $O(n \log n)$ or even $O(n^2)$ and higher order. Then, a reduction of the size of the system of a certain percentage, say $R_{\%}^{size}$, can be translated into an effective reduction in terms of time $R_{\%}^{time}$ that is significantly higher. Such that generally speaking $R_{\%}^{time} > R_{\%}^{size}$, depending on the implementation. While this effect is not systematically highlighted on small systems, it underlines that relatively small reductions in terms of size can yield significant speed-ups of the overall execution of a calculation.

For instance, for the case $h_{ex} = 4$, the full calculation took 7 h and 14 min while the cyclic symmetric one only lasted for 1 min a 31 s, making the reduction of computing time of over $R_{\%}^{time} = 99.6\%$. This is almost 6 points greater than the expected size reduction of $R_{\%}^{size} = 93.75\%$. Obviously, this reduction is highly dependent on implementation and hardware. Throughout the article, only effective size reductions are presented but it must be kept in mind that the announced reductions are most likely underestimated.

## 5.4. Partly standing wave excitation

Now that the capabilities of the overall strategy is demonstrated on the case of traveling wave excitations, a new scenario is considered. In turbomachinery, the order of the traveling excitations is generally determined by the number of blades in the immediately adjacent stages of the studied stage. However, it is possible for a stage to be excited by stages further than its immediate neighbors. In this section, the loading scenario that is studied is one where the blade is under a backward excitation for 90% (adjacent stage) and under a forward excitation for 10% (further stage) for the same diametral index. This means that the excitation is a sum of a backward traveling wave (80%) and of a standing wave (20%). For this case, the engine order considered is $h_{ex} = 4$ with $f_{ex} = 1$ N and the excitation is such as

$$f_{ex,j}(t) = f_{ex} \left[ 0.9 \cos \left( t + \frac{\alpha(j-1)h_{ex}}{\omega} \right) + 0.1 \cos \left( t - \frac{\alpha(j-1)h_{ex}}{\omega} \right) \right] \quad \forall j \in [1, N]. \qquad (62)$$

Since the excitation includes stationary components, the traveling wave hypothesis is not suitable for this case. Hence, it is necessary to resort to the criterion of Section 2.4. It is then predicted that the following basis is needed for such a problem:

$$\begin{cases} D^k = \{0, 8\} & \text{if } k \text{ is even} \\ D^k = \{4\} & \text{if } k \text{ is odd} \end{cases} \quad \forall k \in [0, N_h] \qquad (63)$$

Indeed, this case is subject to the parity criterion since $K = 2$ and $Kh_{ex} = 8 \equiv \frac{N}{2} [N]$. Using this basis, the results of the simulations are displayed in Fig. 11 for DLFT-HBM. Firstly, it is observed that the match is once again excellent between the full and reduced
<figure><img src="image_13.png" /><figcaption>Fig. 11. Partly standing wave excitation on the Catana industrial blade model for $h_{ex} = 4$, linear (……), DLFT-HBM: cyclic symmetry (—), full (---), equivalent backward traveling wave excitation (—).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 12. Cyclic components over all harmonics at the resonance point of the full calculation in Fig. 11 (sector 1), 90% forward and 10% backward traveling wave with $h_{ex} = 4$. Nodal diameter: 0 (■), 4 (■), 8 (■).</figcaption></figure>

calculations even if the contact law is nonsmooth. The standing characteristics of the response are clearly identified by comparing the responses of sectors 1 (Fig. 11(a)) and 2 (Fig. 11(b)). An interesting behavior even occurs at the resonance of sector 1. Indeed, the cyclic decomposition given in Fig. 12 shows that a purely standing wave displacement field emerges from the excitation that is only 20% stationary. Actually, this behavior is due to the fact that the dynamic motion of the system reaches a point where only one in two blades experiences contact along a period at this resonance. On the frequency response curves of Fig. 11, all even numbered blades (sector 2 in Fig. 11(b)) experience a linear behavior while odd sectors feature resonances around 293 Hz. The associated branch for the even sectors also experience a turning point at the same frequency. These points can be visualized in the zooms of Fig. 11. One can see that it is superimposed with the classical linear response of the system.

This type of response underlines the need of investigating the dynamics of the system using general nodal diameter coupling criteria such as the one presented in this article. Indeed, even if the system is under an almost purely backward traveling wave, the small standing components introduced for this scenario creates new couplings with the nonlinearity and makes a pure standing wave motion emerge. Hence, it highlights that no part of the excitation should be neglected in a nonlinear framework because small contributions can still induce specific dynamics in the system. Moreover, this demonstrates that the traveling wave hypothesis is not robust to small disturbances as shown by the superimposed equivalent backward traveling wave excitation in Fig. 11 that predicts wrong vibrational amplitudes on this case.

Besides this observation, the diametral basis chosen in Eq. (63) was correctly identified beforehand and matches perfectly with the cyclic components emerging in the full response in Fig. 12. The predicted parity criterion is also respected like expected. This shows that on a fairly complex configuration, the rules derived in Sections 2.4 and 3.3 are still suitable in partly standing wave excitations scenarios on large scale industrial models.

In order to demonstrate the usefulness of the proposed strategy once again, a similar example is studied with $h_{ex} = 7$. The percentages of forward and backward excitations are kept the same and $f_{ex} = 0.3$ N. $h_{ex} = 7$ is actually prime with $N = 16$ and no reduction is expected in terms of overall diameters in the response. However, the parity criterion yields that only half of the diameters are considered on odd and even harmonics. This means that even if $h_{ex}$ and $N$ are prime with one another, a reduction of 50% is reached (as shown in Table 3). The basis considered for these computation is thus

$$ \begin{cases} D^k = \{0, 2, 4, 6, 8\} & \text{if } k \text{ is even} \\ D^k = \{1, 3, 5, 7\} & \text{if } k \text{ is odd} \end{cases} \quad \forall k \in [0, N_h] \qquad (64) $$
<figure><img src="image_15.png" /><figcaption>Fig. 13. Partly standing wave excitation on the Catana industrial blade model for $h_{ex} = 7$, linear (......), DLFT-HBM: cyclic symmetry (—), full (---), equivalent backward traveling wave excitation (—).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 14. Cyclic components over all harmonics at the resonance point of the full calculation in Fig. 13 (sector 1), 90% forward and 10% backward traveling wave with $h_{ex} = 7$. Nodal diameter: 0 (□), 1 (■), 2 (●), 3 (▲), 4 (■), 5 (▼), 6 (▼), 7 (■), 8 (▲).</figcaption></figure>

The responses obtained for this configuration are displayed in Fig. 13. The response of the system to a pure backward traveling excitation of amplitude $f_{ex} = 0.3$ N is also displayed for comparison purposes. For this case, the reduction method also produces an exact approximation of the solution obtained on the full model. This further supports the reduction capabilities of the method exposed in this article for the reduction of blade-tip/casing problems using cyclic symmetry. The forced responses displayed in Fig. 13 highlight that the purely traveling wave excitation response is very different from the response on the case where 10% of forward traveling excitation is introduced. This shows the level of sensitivity that blade-tip/casing problems experience when it comes to the loading scenario. Actually, when looking at the cyclic components of the response at the resonance of Fig. 13(a) in Fig. 14, one can see that the response is actually a standing wave (both forward components in Fig. 14(a) and backwards components in Fig. 14(b) are equal). For instance, at this point only two sectors out of the $N = 16$ sectors actually experience contact along a period of vibration. This means that even if the excitation is mainly composed of a backward traveling wave, the small forward traveling perturbation yields a localized nonlinear response on a few sectors. These localizations are actually symptoms of the appearance of standing wave displacement fields. One should also note that the parity criterion (that yields all the reduction properties for this example) is strictly respected as predicted by the analysis of the loading that allowed to find the basis in Eq. (64). At this stage, one can note that, as predicted by the analysis of the nonlinear forces applied to the blade in sliding friction, no further couplings are created between nodal diameters by the addition of friction. Indeed, the nodal diameters determined by the criterion of Section 2.4 is valid whether sliding friction is accounted or not. The studied application prevents the apparition of sticking phases during friction due to the high velocities generated by the rotation of the system.

As a way of providing more insights about the phenomenon of localization observed in Fig. 14, the cyclic components of the first harmonic are displayed in Fig. 15 for the whole response for the full calculation carried out with DLFT-HBM in Fig. 13. The results displayed in Fig. 15 show that in the linear domain, i.e. at low and high frequencies, the backward 7ND component is significantly higher than the forward 7ND component and this reflects the excitation applied to the system. However, once the nonlinearity is activated, all other nodal diameters have both their cyclic components increased of the same amount, meaning that these nodal diameters are actually responding on a standing wave. The behavior is quite different with 7ND (associated to the excitation $h_{ex} = 7$) since the backward contribution suddenly drops and the forward contribution raises rapidly. This continues until the resonance where both contributions are almost equal, as already highlighted in Fig. 14. This demonstrates how the response of the system goes from an almost pure traveling wave displacement field in the linear domain (90%) to a standing wave response at its nonlinear resonance.
<figure><img src="image_17.png" /><figcaption>Fig. 15. Cyclic components of the first harmonic for the whole response of the full calculation in Fig. 13(a), 90% forward and 10% backward traveling wave with $h_{ex}$ = 7. Nodal diameter: 1 ( ), 3 ( ), 5 ( ), 7 ( ), backward components ( - - ), forward components ( — ).</figcaption></figure>

Overall, the reduction capabilities of the strategy was here demonstrated on the industrial test case for various loading scenarios: traveling wave excitations and partly standing wave excitations. A localization behavior was shown to emerge in the cases where the traveling excitation was not pure. Though, this behavior was correctly captured by the reduction strategy. The emergence of such solutions vibrational solutions was not previously highlighted and this contributions demonstrates the importance of modeling accurately the excitation in order not to simplify the response of the bladed disk. Indeed, this shows some limitations of Petrov's traveling wave hypothesis [6] in the case of a disturbed excitation. Moreover, as shown in the results presented in this section, a perfect match is observed between reduced and full calculations with the nonsmooth DLFT-HBM contact law. This allows to conclude that the reduction criterion is not sensitive to the regularized parameter that helped performing the analytical derivations in Section 2.3.

## 6. Conclusion

This paper demonstrates the reduction capabilities of the cyclic symmetry in order to tackle vibro-impact applications, as much as small size systems than on industrial models. Firstly, the nodal diameter coupling criterion is derived by analyzing a regularized contact law. Its implementation is quite straightforward and allows to drastically decrease the size of the system to be solved. The criterion relies on the analysis of the Taylor expansion of the contact law. The reduction is dependent on both the number of sectors of the system and the nodal diameters that are excited by the external loadings.

A method to further reduce the computing times is also proposed by analyzing how the different harmonics of the harmonic balance method interact with the nodal diameters. It results in a criterion that is able to reduce even more the size of the problem by a factor 2 for certain combinations of numbers of sectors and nodal diameters of the excitation.

These two original contributions are then put to the test on numerous forced response calculations. In order to cover all possible configurations, academic models with both even and odd numbers of sectors are studied through different loading scenarios: standing wave excitations and loadings along multiple nodal diameters (traveling wave excitation and casing deformation). The two criterion used to build the basis of diameters that serve to reduce the size of the problem are shown to produce results that match exactly with the full reference calculations. Indeed, the reduction is exact and yields no loss in precision whatsoever.

These observations of the response are also true when analyzing the results obtained on the state-of-the-art fan stage ECL5/Catana under blade-tip/casing contacts. In fact, even though the friction forces are added in the modeling of blade-tip/casing interactions, the criterion holds true and the results on the full model are matching exactly with the reduced ones, even on the most complex loading scenarios. The reduction strategy also allows to access the intricate dynamics of responses involving strong couplings between several nodal diameters even when the excitation is a quasi-traveling wave. It was shown the traveling wave hypothesis still performs well when the excitation is a pure traveling wave. However, the new criterion enables to study the robustness of the traveling wave hypothesis to small perturbations and even showcase a limitation of Petrov's method [6].

In both cases, the analysis of the cyclic components of the response validate a *posteriori* the choice of the different basis of nodal diameters for the reduced calculations. As a whole, the analytical derivation about the nodal diameter couplings were consolidated by numerous examples and the criterion appears to be industry-ready for the simulation of complex blade-tip/casing interactions on full bladed disks at a low computational cost. Moreover, it appears that these derivation made on the regularized contact law are still valid when considering a truly nonsmooth contact law. This means that the reductions still hold when the contact laws are non differentiable.
<figure><img src="image_18.png" /><figcaption>Fig. 16. Expected reduction as a function of N until N = 150, even number of sectors with the parity rule (■), odd number of sectors (□), even number of sectors (○).</figcaption></figure>

As a way of concluding this study, a wide view of the reduction rates achieved by the proposed criterion is given hereafter. In the case of a traveling wave excitation, the traveling wave hypothesis systematically yields a drastic reduction of a rate

$$R_{\%}^{\text{TW}} = 100 \left( 1 - \frac{1}{N} \right), \qquad (65)$$

that is over 90% as soon as the system is composed of more than N = 10 blades.

For the more general case where the loading is not necessarily a traveling wave, the mean reductions achieved along all non-zero engine orders for numbers of sectors up to 150 are displayed in Fig. 16. First of all, the reduction strategy is unfortunately not suitable when the number of sectors of the system is a prime number: no reduction can be achieved. Moreover, it can be seen that even numbers of sectors are much more subject to reduction than odd numbers of sectors. Considering a resolution through harmonic balance method that allows to use the parity rule, the mean reduction of even numbers of sector reaches 62.4% while non-prime odd numbers of sectors are reduced of 23.4% on average and it drops down to 12.5% when considering the prime numbers. Though, the overall mean reduction is 37.6% which is quite large when considering the fact that the reduction is exact. If a resolution through other methods than the harmonic balance method is preferred, the mean reduction among even number of sectors drops down to 37.7% and the overall reduction falls to an average of 25.2%.

Overall, this work has participated in making a major step forward in the efficient simulation of the nonlinear vibrations of bladed disk featuring blade-tip/casing contacts. Indeed, such simulations were conducted on a single blade or on the full bladed disk to date, making them unrealistic or extremely cumbersome. Hence, the development of a reduction strategy such as the one presented here enables to carry out increasingly realistic calculations at an affordable computational burden.

Among the perspectives of this work, the stability analysis of the response of these cyclic systems featuring vibro-impact is one of the main avenues of research. Indeed, this study has focused on solving the equation of motion through the harmonic balance method to validate the reduction procedure in terms of nodal diameters. It has shown that reduction is exact since the same results are obtained with the cyclic symmetric formulation as well as with the full system. However, this study does not give any information about the stability of the presented solutions and a stability analysis would be required to state on the nature of these solutions. Moreover, a cross-verification of the stability of these nonlinear solutions was not carried out in this article because it exceeded its scope. Nevertheless, it would be necessary in order to validate the physical nature of these solutions. The characterization of eventual branches in terms of cyclic components would also allow to understand more finely the dynamics response of such cyclic systems. The use of the reduction strategy in conjunction with time-marching strategies to tackle transient dynamics of cyclic systems is also among the following steps of this work. Indeed, the nodal diameter criterion is not relying on any periodicity hypothesis and is expected to be also relevant for solution methods related to non-periodic vibrations.

In order to generalize the reduction strategy, it would also be necessary to investigate how it is possible to extend the nodal diameters couplings to systems with a broken cyclic symmetry : structural mistuning, variation in contact conditions, geometric mistuning. Dealing with this engineering challenge would be the next step towards the industrialization of this strategy.

The potential applications of this work go beyond the world of aeronautical turbomachinery: any nonlinear system dominated by unilateral contact with a cyclic symmetry and subject to an external loading along known spatial harmonic indices could benefit from this type of reduction. Among these systems, one can cite epicyclic gearings, automotive turbocompressors, rotodynamic pumps found in rocket engines or even turbomolecular pump.
CRediT authorship contribution statement

Thibaut Vadcard: Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Conceptualization. Samuel Quaegebeur: Writing – review & editing, Methodology, Conceptualization. Fabrice Thouverez: Writing – review & editing, Supervision, Funding acquisition.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

The authors are grateful for the financial support of the ANR (project ANR-22-CPJ2-0061-01) and the Horizon Europe HE-ART project (Hybrid Electric propulsion system for regional AiRcraft, grant agreement ID: 101102013). The authors thank LMFA for providing the finite element mesh of the ECL5 fan stage. Funded by the European Union. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the Clean Aviation Joint Undertaking. Neither the European Union nor the granting authority can be held responsible for them.

Appendix A. Developments associated with the parity criterion for the selection of nodal diameters in the case of unilateral nonlinearities

When using the HBM, the choice of the sets $D^k$ can lead to further reduction by evaluating which harmonic will respond on which nodal diameter. This section aims to demonstrate the existence of such a relation for unilateral nonlinearities. For simplicity, only a single nodal diameter will be considered in the initial basis $D_{q,0}$. Extension to multiple initial nodal diameters is then straightforward.

Consider $h_{ex}$, the initial nodal diameter associated to K the smallest integer which respects the following

$$2Kh_{ex} \equiv 0 [N]. \qquad (A.1)$$

This definition is derived in Section 2.4. As a reminder, to determine the full nodal diameter basis, one has to proceed with the following steps

* Create a vector $v_K = [0, 1, 2, ..., q, ..., K]$ (a scalar of this vector is noted $q$)

* Multiply it by $h_{ex}$ leading to a vector $v_{K,h_{ex}} = [0, h_{ex}, 2h_{ex}, ..., qh_{ex}, ..., Kh_{ex}]$.

* for each component $qh_{ex}$, its reminder through the Euclidian division by $N$ is evaluated, noted $r_q$. Then, if $0 < r_q < m$ with $m$ the maximum number of nodal diameters for the given number of sectors, the excited diameter is $d_q = r_q$. Otherwise, if $m < r_q$ then the excited diameter is $d_q = N - r_q$. A final vector is obtained with the list of the interacting nodal diameter $D = [d_0, d_1, d_2, ..., d_q, ..., d_K]$ (where $d_0 = 0$, and $d_1 = h_{ex}$).

Two scenarios exist:

* if $Kh_{ex} \equiv 0 [N]$, then no reduction is obtained with the HBM,

* if $Kh_{ex} \equiv \frac{N}{2} [N]$, then the nodal diameters $d_q$ with $q$ even only respond on even harmonics, and the nodal diameters $d_q$ with $q$ odd only respond on odd harmonics.

This second scenario allows to divide the number of unknowns by 2. The following sections aim to demonstrate this property: each scenario is dealt with separately in Appendices A.1 and A.2 respectively.

A.1. First scenario

For the first scenario, one assumes that

$$Kh_{ex} \equiv 0 [N]. \qquad (A.2)$$

The number $K$ is then necessarily an odd number. Indeed, if $K$ were to be even, it could be factorized by an integer $K_2$ such that $K = 2K_2$. Substituting in Eq. (A.2) gives

$$2K_2h_{ex} \equiv 0 [N]. \qquad (A.3)$$

This breaks the definition of $K$, i.e. $K$ is defined as the smallest integer verifying Eq. (A.1). Therefore, $K_2$ cannot exist and $K$ is necessarily an odd number in the following. Then, it comes that

$$v_K = [0, 1, \dots, \frac{K-1}{2}, \frac{K+1}{2}, \dots, K-1, K]. \qquad (A.4)$$
Following the expression of $v_K$ in this case, a symmetry then occurs in the interacting nodal diameters $D$ between both ends of the list: $d_K = d_0 = 0$, then $d_{K-1} = d_1 = h_{ex}$, etc. As a consequence, each nodal diameter appears twice in $D$ and it is such that

$$D = \left[ 0, h_{ex}, \dots, d_{\frac{K-1}{2}}, d_{\frac{K-1}{2}}, \dots, h_{ex}, 0 \right]. \qquad (A.5)$$

Notice that both indices 1 and $K-1$ (respectively an even and an odd number) are associated to the same nodal diameter $h_{ex}$.

Consider now the HBM in an iterative procedure as described in Section 2.4. Initially, only the nodal diameter $h_{ex}$ constitutes the basis and responds only on the first harmonic $e^{\pm i\omega t}$. Then through Eq. (24), the second diameter $2h_{ex}$ responds on its second harmonic $e^{\pm 2i\omega t}$, and so on. Notice that $K-1$ is even and thus the nodal diameter $r_{K-1} = h_{ex}$ will also have an even harmonic component. Due to the nodal diameter interaction, the harmonics gets coupled and eventually all harmonics respond on all nodal diameters.

## A.2. Second scenario

For this second scenario, one assumes that

$$Kh_{ex} \equiv r[N], \quad r \neq 0. \qquad (A.6)$$

Since $2Kh_{ex} \equiv 0$ [$N$], Eq. (A.6) implies that $2Kh_{ex}$ is proportional to $N$ with an odd number. Moreover, since $2Kh_{ex}$ is an even number then $N$ must be even, and $K$ must be odd. Therefore, it comes that

$$\begin{aligned} 2Kh_{ex} &= (p+1)N \\ Kh_{ex} &= \frac{N}{2} + 2pN, \end{aligned} \qquad (A.7)$$

hence $r = \frac{N}{2}$.

Now, assume that there exists an integer $1 < q < K$ such that

$$qh_{ex} \equiv h_{ex}[N]. \qquad (A.8)$$

This relation is equivalent to

$$(q \pm 1) h_{ex} \equiv 0 \quad [N]. \qquad (A.9)$$

To ensure that $K$ is the smallest integer verifying Eq. (A.1), the relation $q \pm 1 > K$ must be respected. Therefore, either of $q > K+1$ or $q > K-1$ is true. Neither of the two relations is not compatible with $q < K$. As a consequence, there exists no integer $1 < q < K$ such that Eq. (A.8) holds true. Therefore, in this scenario, each nodal diameter that is contained in $D$ is only present for a single value $q$ unlike for the first scenario in Appendix A.1. Hence it comes that

$$v_K = [0, 1, \dots, q, \dots, K] \qquad (A.10)$$

and

$$D = \left[ 0, h_{ex}, \dots, d_q, \dots, \frac{N}{2} \right]. \qquad (A.11)$$

To determine which harmonic will respond on which nodal diameter, the HBM can be employed with an iterative process. Without loss of generality, the excited diameter $d_1 = h_{ex}$ is assumed to be excited along its first harmonic $e^{i\omega t}$. Through the first iteration (determination of the coupling with only $d_1$), the nodal diameter $d_q$ will systematically have a component along the $q$th harmonic $e^{\pm i\omega q t}$. After this first iteration, the entire set of interacting nodal diameters is determined. Due to the aforementioned property, after this first iteration, each nodal diameter is only associated to a single harmonic component. Notice also that if $q$ is even (resp. odd) then $d_q$ is associated with an even (resp. odd) harmonic.

The next iterations of the procedure to determine the harmonic couplings in HBM consist in identifying how the nodal diameter will combine with each other and hence how their temporal harmonic content will be combined. To determine this, Eq. (26) is recalled

$$\hat{f}_{nl,3}^{q} = \frac{\kappa}{\sqrt{N}\gamma^3} \sum_{l=0}^{\infty} \left( \beta_l \sum_{k \in \mathcal{F}_l} \left( \alpha_k \prod_{p_k=1}^{Q_k} (\hat{u}_{h_{p_k}})^{\beta_{p_k}} \left( \sum_{j=1}^{N} e^{i\alpha_{j-1}} \left( \sum_{p_k=1}^{Q_k} h_{p_k} \beta_{p_k} - q \right) \right) \right) \right) \qquad (A.12)$$

The second iteration concerns $l=1$, hence $\sum_{p_k=1}^{Q_k} \beta_{p_k} = 2$. In this case, the interaction between only two nodal diameters is studied. The nodal diameter $d_q$ can be obtained through the combination of two other nodal diameter $d_a$ and $d_b$ if $a \pm b \equiv q \mod N$. Two scenarios arise

* If $q$ is even, then necessarily $a$ and $b$ must both be either odd or even. In either case, the combination of their temporal harmonic $e^{i\omega a t}$ and $e^{i\omega b t}$ will be an even harmonic $e^{i\omega(a\pm b)t}$ associated to $d_q$.

* If $q$ is odd, then necessarily $a$ and $b$ one is even and the other odd, leading to an odd harmonic $e^{i\omega(a\pm b)t}$ associated to $d_q$.
<caption>**Table B.4**</caption>

Reduction for the case of the simplified system with $N = 45$ sectors for 10 harmonics (full size: $N \times n_{nl} \times (2N_h + 1) = 945$ unknowns).

<table><thead><tr><th>Engine order</th><th>0</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th><th>11</th></tr></thead><tbody><tr><th>Value of K</th><td>1</td><td>45</td><td>45</td><td>15</td><td>45</td><td>9</td><td>15</td><td>45</td><td>45</td><td>5</td><td>9</td><td>45</td></tr><tr><th>Reduced size (DoF)</th><td>21</td><td>945</td><td>945</td><td>315</td><td>945</td><td>189</td><td>315</td><td>945</td><td>945</td><td>105</td><td>189</td><td>945</td></tr><tr><th>Sectors for AFT-CP</th><td>1</td><td>45</td><td>45</td><td>15</td><td>45</td><td>9</td><td>15</td><td>45</td><td>45</td><td>5</td><td>9</td><td>45</td></tr><tr><th>Reduction (%)</th><td>97.8</td><td>0.0</td><td>0.0</td><td>66.7</td><td>0.0</td><td>80.0</td><td>66.7</td><td>0.0</td><td>0.0</td><td>88.9</td><td>80.0</td><td>0.0</td></tr><tr><th>Engine order</th><th>12</th><th>13</th><th>14</th><th>15</th><th>16</th><th>17</th><th>18</th><th>19</th><th>20</th><th>21</th><th>22</th></tr><tr><th>Value of K</th><td>15</td><td>45</td><td>45</td><td>3</td><td>45</td><td>45</td><td>5</td><td>45</td><td>9</td><td>15</td><td>45</td></tr><tr><th>Reduced size (DoF)</th><td>315</td><td>945</td><td>945</td><td>63</td><td>945</td><td>945</td><td>105</td><td>945</td><td>189</td><td>315</td><td>945</td></tr><tr><th>Sectors for AFT-CP</th><td>15</td><td>45</td><td>45</td><td>3</td><td>45</td><td>45</td><td>5</td><td>45</td><td>9</td><td>15</td><td>45</td></tr><tr><th>Reduction (%)</th><td>66.7</td><td>0.0</td><td>0.0</td><td>93.3</td><td>0.0</td><td>0.0</td><td>88.9</td><td>0.0</td><td>80.0</td><td>66.7</td><td>0.0</td></tr></tbody></table>

Following these elements, the second iteration shows that a nodal diameter $d_q$ can be generated through another interaction (involving two nodal diameter) but the parity of its associated harmonic does not change.

Through an iterative process, the $(k+1)$th coupling relation $(l=k)$ is conditioned by $\sum_{p_k=1}^{Q_k} \beta_{p_k} = 2k$. Here, the interaction between $2k$ nodal diameters is studied. A similar result is obtained: nodal diameter $d_q$ can be obtained through this interaction but the parity property must still hold. This is caused by the fact that $2k$ is an even number.

In the end, if $q$ is even, only even harmonics will respond on the nodal diameter $d_q$ and conversely if $q$ is odd.

## Appendix B. Academic case with an odd number of sectors: $N = 45$

The case study of this appendix is the same model as in Section 4.3 but with $N = 45$ sectors. It is dealt with for completeness: the reductions between even and odd numbers of sectors are significantly different. The fact that $N$ is odd allows to showcase the reduction properties of the proposed method. This time, no calculation is able to take advantage of the parity criterion on the harmonics as mentioned in Section 3.3. With this number of sectors, the expected reductions are displayed in Table B.4 for all engine orders. It can be seen that the reductions are lower that for the previous case in Section 4.3: the mean reduction only reaches 38.1% (mean reduction of 75% for $N = 24$). Moreover, Table B.4 reveals that 12 values of $h_{ex}$ are not subject to any reduction at all.

For the simulations, the approach is similar as in Section 4.3: firstly the response of the system under standing wave excitations is analyzed to validate the choice of the nodal diameter basis and then a loading scenario involving both a traveling wave excitation and a casing deformation is studied. Note that a perfect match was observed between full calculations and the traveling wave hypothesis reduction, but the results are not shown for the sake of brevity.

### B.1. Standing wave excitations

In this section, more complex excitations are applied to the system: standing wave excitations. In order to demonstrate the reduction capabilities of the methodology, two values of $h_{ex}$ where a reduction is predicted are picked from Table B.4. Simulations are led for $h_{ex} = 9$ with $f_{ex} = 1$ N and $h_{ex} = 15$ with $f_{ex} = 1$ N. The two basis associated with these values of $h_{ex}$ are

$$D^k|_{k \in \llbracket 0, 10 \rrbracket} = \{0, 9, 18\} \text{ for } h_{ex} = 9 \quad \text{and} \quad D^k|_{k \in \llbracket 0, 10 \rrbracket} = \{0, 15\} \text{ for } h_{ex} = 15. \qquad (B.1)$$

For $h_{ex} = 9$, the procedure of Section 2.4 yields $K = 5$ and $K = 3$ for $h_{ex} = 15$. With these basis, the results are displayed in Fig. B.17 for both R-HBM and DLFT-HBM.

Again, the reduction criterion on the nodal diameters produces an exactly matching response for both the regularized and the nonsmooth contact laws. As a way of verifying the choice of the basis of Eq. (B.1), the cyclic contributions are displayed in Fig. B.18(a) for $h_{ex} = 9$ and in Fig. B.18(b) for $h_{ex} = 15$. For this case, since the response is purely stationary (forward contributions are strictly equal to backward contributions), no distinction is made on the direction of propagation and only one of the two graphs is displayed for each $h_{ex}$. For both cases, the basis identified in Eq. (B.1) are exactly matching with what the response of the full model features. This shows that the coupling rule derived in Section 2.4 is also relevant for cases with an odd number of sectors. Contrarily with the case in Section 4.3.1, there is no parity rule among the cyclic components and the whole basis responds on all harmonics. This property is also correctly predicted by the parity rule of Section 3.3: the criterion is only possible for even number of sectors, which is not the case here.

### B.2. Mixed contact scenario

Finally, the academic test case is characterized through an even more complex contact scenario: it is excited through a traveling wave and the casing it interacts with is deformed with a diametral shape. The system is excited along a traveling wave excitation on the 12th engine order ($h_{ex} = 12$) with an amplitude of $f_{ex} = 2$ N and the casing is deformed with a 3ND shape ($h_c = 3$) with $p_j = 0.05c_j$. Both excitations are traveling waves, however both orders are expected to couple themselves and make a displacement
<figure><img src="image_19.png" /><figcaption>Fig. B.17. Standing wave excitation on the simplified bladed disk with N = 45 sectors for various engine orders, linear (……), R-HBM: cyclic symmetry (—), full (—-), DLFT-HBM: cyclic symmetry (—), full (—-).</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. B.18. Cyclic components over all harmonics at the resonance point of the full calculation in Fig. B.17(b) and in Fig. B.17(d), standing wave excitations. Nodal diameter: 0 (■), 9 (■), 15 (■), 18 (■).</figcaption></figure>

field of standing wave shape emerge. In order to build the basis of the nodal diameters that intervene in the response, one should take a decomposition that is related to the GCD of $h_c$ and $h_{ex}$, 3 in this case, and the basis is

$$D^k|_{k \in \{0,1,0\}} = \{0, 3, 6, 9, 12, 15, 18, 21\}. \qquad (B.2)$$

Even if the GCD of $h_{ex}$ and $h_c$ is equal to $h_{ex}$, the nodal diameter basis is extended. Indeed, if no deformation of the casing were present, a simple traveling wave basis could have been chosen. The computations are then led for both DLFT-HBM and R-HBM and displayed in Fig. B.19 for two sectors (1 and 4).

In this scenario and with the chosen diametral basis, there is again a perfect match between the computations led in cyclic symmetry and the full calculations. This agreement for such of a sophisticated contact scenario demonstrates the usefulness of the proposed strategy. It offers reduction capabilities in numerous configurations. For the present case, a reduction of 66.7% is achieved. As a way of demonstrating that the diametral basis is correctly identified, the cyclic contributions are displayed in Fig. B.20. Contrarily with standing wave excitations, the displacement field is not purely stationary: the system experiences a vibrational motion composed of multiple traveling waves of different orders as well as stationary components that cannot be easily predicted by analyzing the loadings (excitation and casing) that are applied. Nevertheless, the nodal diameters responding in Fig. B.20
<figure><img src="image_21.png" /><figcaption>Fig. B.19. Traveling wave excitation on the simplified bladed disk with $N = 24$ sectors for $h_{ex} = 12$ and a casing deformed with $h_c = 3$, linear (· · · · ·), R-HBM: cyclic symmetry (—), full (---), DLFT-HBM: cyclic symmetry (—), full (---).</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Fig. B.20. Cyclic components over all harmonics at the resonance point of the full calculation in Fig. B.19, traveling wave with $h_{ex} = 12$ and a casing deformed with $h_c = 3$. Nodal diameter: 0 (■), 3 (■), 6 (■), 9 (■), 12 (■), 15 (■), 18 (■), 21 (■).</figcaption></figure>

correspond to the basis identified in Eq. (B.2) justifying once again the relevance of the nodal diameter coupling criterion for unilateral nonlinearities.

Data availability

The authors do not have permission to share data.
References

[1] International Energy Agency, Net zero by 2050, a roadmap for the global energy sector, 2021, URL https://www.iea.org/reports/net-zero-by-2050.

[2] A. Millecamps, J.-F. Brunel, P. Dufrénoy, F. Garcin, M. Nucci, Influence of thermal effects during blade-casing contact experiments, in: Proceedings of the ASME IDETC-CIE 2009, 2009, pp. 855-862, http://dx.doi.org/10.1115/DETC2009-86842.

[3] A. Millecamps, A. Batailly, M. Legrand, F. Garcin, Snecma's viewpoint on the numerical and experimental simulation of blade-tip/casing unilateral contacts, in: Proceedings of the ASME Turbo Expo 2015, 2015, http://dx.doi.org/10.1115/GT2015-42682.

[4] R. MacNeal, R. Harder, J. Mason, NASTRAN cyclic symmetry capability, NASTRAN: Users' Exp. (1973) https://ntrs.nasa.gov/api/citations/19740006491/downloads/19740006491.pdf.

[5] D.L. Thomas, Dynamics of rotationally periodic structures, Internat. J. Numer. Methods Engrg. 14 (1) (1979-01) 81-102, http://dx.doi.org/10.1002/nme.1620140107.

[6] E.P. Petrov, A method for use of cyclic symmetry properties in analysis of nonlinear multiharmonic vibrations of bladed disks, J. Turbomach. 126 (1) (2004-01) 175-183, http://dx.doi.org/10.1115/1.1644558.

[7] E. Petrov, Analysis of sensitivity and robustness of forced response for nonlinear dynamic structures, Mech. Syst. Signal Process. 23 (1) (2009-01) 68-86, http://dx.doi.org/10.1016/j.ymssp.2008.03.008.

[8] C. Siewert, L. Panning, J. Wallaschek, C. Richter, Multiharmonic forced response analysis of a turbine blading coupled by nonlinear contact forces, J. Eng. Gas Turbines Power 132 (8) (2010-05) http://dx.doi.org/10.1115/1.4000266.

[9] S. Quaegebeur, B. Chouvion, F. Thouerez, Model reduction of nonlinear cyclic structures based on their cyclic symmetric properties, Mech. Syst. Signal Process. 145 (2020-11) 106970, http://dx.doi.org/10.1016/j.ymssp.2020.106970.

[10] Y. Colaïtis, A. Batailly, Development of a harmonic balance method-based numerical strategy for blade-tip/casing interactions: Application to NASA rotor 37, J. Eng. Gas Turbines Power 143 (11) (2021) http://dx.doi.org/10.1115/1.4051967.

[11] T. Vadcard, Y. Colaïtis, A. Batailly, F. Thouerez, Assessment of two harmonic balance method-based numerical strategies for blade-tip/casing interactions: Application to nasa rotor67, J. Eng. Gas Turbines Power 144 (12) (2022) http://dx.doi.org/10.1115/1.4055416.

[12] M. Legrand, C. Pierre, Numerical investigation of abradable coating wear through plastic constitutive law: Application to aircraft engines, in: Proceedings of the ASME IDETC-CIE 2009, 2009, pp. 907-916, http://dx.doi.org/10.1115/DETC2009-87669.

[13] R.J. Williams, Simulation of blade casing interaction phenomena in gas turbines resulting from heavy tip rubs using an implicit time marching method, in: Proceedings of the ASME 2011 Turbo Expo, 2011, pp. 1007-1016, http://dx.doi.org/10.1115/GT2011-45495.

[14] K.E. Turner, C. Padova, M. Dunn, Airfoil deflection characteristics during rub events, J. Turbomach. 134 (1) (2010) 11, http://dx.doi.org/10.1115/1.4003257.

[15] S.K. Sinha, Non-linear dynamic response of a rotating radial Timoshenko beam with periodic pulse loading at the free-end, Int. J. Non Linear Mech. 40 (1) (2005) 113-149, http://dx.doi.org/10.1016/j.ijnonlinmec.2004.05.019.

[16] E. Piollet, F. Nyssen, A. Batailly, Blade/casing rubbing interactions in aircraft engines: Numerical benchmark and design guidelines based on NASA rotor 37, J. Sound Vib. 460 (2019) 114878, http://dx.doi.org/10.1016/j.jsv.2019.114878.

[17] A. Batailly, M. Legrand, P. Cartraud, C. Pierre, Assessment of reduced models for the detection of modal interaction through rotor stator contacts, J. Sound Vib. 329 (26) (2010-12) 5546-5562, http://dx.doi.org/10.1016/j.jsv.2010.07.018.

[18] E. Delhez, F. Nyssen, J.-C. Golinval, A. Batailly, Numerical study of bladed structures with geometric and contact nonlinearities, J. Sound Vib. 544 (2023) 117382, http://dx.doi.org/10.1016/j.jsv.2022.117382.

[19] F. Nyssen, Y. Colaïtis, A. Batailly, Numerical investigation of a mistuned industrial bladed disk dynamics with structural contacts using time and frequency methods, J. Sound Vib. 535 (2022-09) 117077, http://dx.doi.org/10.1016/j.jsv.2022.117077.

[20] Y. Colaïtis, A. Batailly, The harmonic balance method with arc-length continuation in blade-tip/casing contact problems, J. Sound Vib. 502 (2021) 116070, http://dx.doi.org/10.1016/j.jsv.2021.116070.

[21] E.P. Petrov, Stability analysis of multiharmonic nonlinear vibrations for large models of gas turbine engine structures with friction and gaps, J. Eng. Gas Turbines Power 139 (2) (2016) http://dx.doi.org/10.1115/1.4034353.

[22] Y. Colaïtis, A. Batailly, Stability analysis of periodic solutions computed for blade-tip/casing contact problems, J. Sound Vib. 538 (2022) 117219, http://dx.doi.org/10.1016/j.jsv.2022.117219.

[23] L. Salles, B. Staples, N. Hoffmann, C. Schwingshackl, Continuation techniques for analysis of whole aeroengine dynamics with imperfect bifurcations and isolated solutions, Nonlinear Dynam. 86 (3) (2016) 1897-1911, http://dx.doi.org/10.1007/s11071-016-3003-y.

[24] T. Vadcard, F. Thouerez, A. Batailly, On the detection of nonlinear normal mode-related isolated branches of periodic solutions for high-dimensional nonlinear mechanical systems with frictionless contact interfaces, Comput. Methods Appl. Mech. Engrg. 419 (2024) 116641, http://dx.doi.org/10.1016/j.cma.2023.116641.

[25] M. Krack, Nonlinear modal analysis of nonconservative systems: Extension of the periodic motion concept, Comput. Struct. 154 (2015) 59-71, http://dx.doi.org/10.1016/j.compstruc.2015.03.008.

[26] D. Laxalde, F. Thouerez, Complex non-linear modal analysis for mechanical systems: Application to turbomachinery bladings with friction interfaces, J. Sound Vib. 322 (4-5) (2009) 1009-1025, http://dx.doi.org/10.1016/j.jsv.2008.11.044.

[27] S. Nacivet, C. Pierre, F. Thouerez, J. Lezequel, A dynamic Lagrangian frequency-time method for the vibration of dry-friction-damped systems, J. Sound Vib. 265 (1) (2003) 201-219, http://dx.doi.org/10.1016/S0022-460X(02)01447-5.

[28] T. Vadcard, A. Batailly, F. Thouerez, On harmonic balance method-based Lagrangian contact formulations for vibro-impact problems, J. Sound Vib. 531 (2022) 116950, http://dx.doi.org/10.1016/j.jsv.2022.116950.

[29] G. Von Groll, D. Ewins, The harmonic balance method with arc-length continuation in rotor/stator contact problems, J. Sound Vib. 241 (2) (2001) 223-233, http://dx.doi.org/10.1006/jsvi.2000.3298.

[30] T. Vadcard, F. Thouerez, A. Batailly, Computation of isolated periodic solutions for forced response blade-tip/casing contact problems, J. Eng. Gas Turbines Power 146 (4) (2024) 041011, http://dx.doi.org/10.1115/1.4063704.

[31] L. Woiwode, N.N. Balaji, J. Kappauf, F. Tubita, L. Guillot, C. Vergez, B. Cochelin, A. Grolet, M. Krack, Comparison of two algorithms for harmonic balance and path continuation, Mech. Syst. Sig. Process. 136 (2020) 106503, http://dx.doi.org/10.1016/j.ymssp.2019.106503.

[32] C. Joannin, B. Chouvion, F. Thouerez, J.-P. Ousty, M. Mbaye, A nonlinear component mode synthesis method for the computation of steady-state vibrations in non-conservative systems, Mech. Syst. Sig. Process. 83 (2017) 75-92, http://dx.doi.org/10.1016/j.ymssp.2016.05.044.

[33] V. Pagès, P. Duquesne, S. Aubert, L. Blanc, P. Ferrand, X. Ottavy, C. Brandstetter, UHBR open-test-case fan ECL5/CATANA, Int. J. Turbomach. Propuls. Power 7 (2) (2022) http://dx.doi.org/10.3390/ijtpp7020017.

[34] C. Brandstetter, V. Pages, P. Duquesne, X. Ottavy, P. Ferrand, S. Aubert, L. Blanc, UHBR open-test-case fan ECL5/catana, part 1: Geometry and aerodynamic performance, in: 14th European Conference on Turbomachinery Fluid Dynamics and Thermodynamics, in: ETC14, European Turbomachinery Society, 2021, http://dx.doi.org/10.29008/etc2021-626.

[35] K. Billon, L. Sanchez, G. Bouvard, C. Gibert, L. Blanc, F. Thouerez, Experiment on tuned open-test-case fan ECL5/catana: structural characterization under vacuum conditions, in: Proceedings of the ASME Turbo Expo 2024, 2024, http://dx.doi.org/10.1115/GT2024-127539.
[36] T. Vadcard, S. Quaegebeur, F. Thouverez, Numerical analysis of the coupling mechanisms between blade/disk friction and blade-tip/casing contacts on the state-of-the-art ECL5/catana fan stage, Int. J. Non-Linear Mech. 170 (2025-03) 104970, http://dx.doi.org/10.1016/j.ijnonlinmec.2024.104970.

[37] R.R. Craig, M.C.C. Bampton, Coupling of substructures for dynamic analyses., AIAA J. 6 (7) (1968) 1313–1319, http://dx.doi.org/10.2514/3.4741.

[38] A.P. Schneider, A.-L. Fiquet, B. Paoletti, X. Ottavy, C. Brandstetter, Experiments on tuned UHBR open-test-case fan ECL5/CATANA: Performance and aerodynamics, J. Turbomach. 146 (8) (2024) http://dx.doi.org/10.1115/1.4064231.

[39] A. Batailly, Q. Agrapart, A. Millecamps, J.-F. Brunel, Experimental and numerical simulation of a rotor/stator interaction event localized on a single blade within an industrial high-pressure compressor, J. Sound Vib. 375 (2016) 308–331, http://dx.doi.org/10.1016/j.jsv.2016.03.016.

[40] I. Favretti, L. Roux, A. Batailly, Validation of the numerical simulation of rotor/stator interactions in aircraft engine low-pressure compressors, J. Eng. Gas Turbines Power 147 (5) (2024-12) http://dx.doi.org/10.1115/1.4066789.