# 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

Thibaut Vadcard, Samuel Quaegebeur, Fabrice Thouverez

## Abstract

For the first time, the interaction between two contact nonlinearities, blade-tip/casing contact and blade/disk friction, is investigated simultaneously. This contribution addresses the associated complex coupling mechanisms on the state-of-the-art industrial fan stage ECL5/Catana open test case. Both contact interfaces are accounted for in nonlinear calculations using the harmonic balance method. First, a static analysis is performed to investigate the coupling mechanism between the friction coefficient and the static blade-tip gap. Then, the dynamics response of the blade is analyzed in forced response for two reference configurations: (1) with blade-tip/casing contacts only and (2) with blade/disk friction only. Particular attention is paid to the contact kinematics of both interfaces and the vibrational motion of the blade. The forced responses are analyzed using different local indicators of the contact interfaces. As the friction coefficient is usually uncertain, the robustness of the numerical strategy is then demonstrated by performing a large number of simulations on different root friction coefficients. In the end, even more computations are performed, as a part of sensitivity studies on first order parameters, in order to map out the different relevant quantities of both interfaces: number of contacting nodes at the tip, tip contact resultant force, energy dissipation. These large scale mappings allow to identify trends and inter-dependences between the kinematics of each interface. The originality of this paper lies in the thorough analysis of this seldom studied nonlinear interaction, that is crucially important for the design of safe and efficient aircraft engines. This work demonstrates that accounting for the blade/disk friction interface in numerical simulations strongly mitigates the dynamics of the blade-tip/casing contacts and can even prevent it for a specific range of friction coefficients.

## Keywords

harmonic balance method, friction damping, blade-tip/casing interactions, nonlinear dynamics, computational mechanics
# Analyse numérique des mécanismes de couplages entre le frottement aube/disque et le contact aube/carter sur un étage de soufflante de dernière génération ECL5/Catana

Thibaut Vadcard, Samuel Quaegebeur, Fabrice Thouverez

## Résumé

Pour la première fois, l'interaction entre deux non-linéarités de contact, le contact aube/carter et le frottement aube/disque, est étudiée simultanément. Cette contribution aborde les mécanismes de couplage complexes associés sur l'étage de soufflante industriel de dernière génération ECL5/Catana, un cas test ouvert. Les deux interfaces de contact sont prises en compte dans les calculs non linéaires à l'aide de la méthode de l'équilibrage harmonique. Tout d'abord, une analyse statique est réalisée pour étudier le mécanisme de couplage entre le coefficient de frottement en pied de pale et le jeu statique entre la tête de pale et le carter. Ensuite, la réponse dynamique de la pale est analysée en réponse forcée pour deux configurations de référence : (1) avec des contacts entre la tête de pale et le carter uniquement et (2) avec un frottement entre la pale et le disque uniquement. Une attention particulière est accordée à la cinématique de contact des deux interfaces et au mouvement vibratoire de la pale. Les réponses forcées sont analysées à l'aide de différents indicateurs locaux des interfaces de contact. Le coefficient de frottement étant généralement incertain, la robustesse de la stratégie numérique est ensuite démontrée en effectuant un grand nombre de simulations avec différents coefficients de frottement au pied de pale. Enfin, d'autres calculs sont effectués, dans le cadre d'études de sensibilité sur les paramètres de premier ordre, afin de cartographier les différentes quantités pertinentes des deux interfaces : nombre de points de contact au sommet d'aube, force résultante de contact au sommet d'aube et dissipation d'énergie. Ces cartographies à grande échelle permettent d'identifier les tendances et les interdépendances entre les cinématiques de chaque interface. L'originalité de cet article réside dans l'analyse approfondie de cette interaction non linéaire rarement étudiée, qui est d'une importance cruciale pour la conception de moteurs aéronautiques sûrs et efficaces. Ce travail démontre que la prise en compte de l'interface de frottement aube/disque dans les simulations numériques atténue fortement la sévérité des contacts entre l'extrémité de l'aube et le carter et peut même l'empêcher pour une gamme spécifique de coefficients de frottement.

## Mots-clés

Méthode d'équilibrage harmonique, amortissement frottant, contact aube/carter, dynamique non linéaire, mécanique numérique

1 - École Centrale de Lyon, Laboratoire de Tribologie et Dynamique des Systèmes, UMR CNRS 5513, 36 avenue Guy de Collongue, Écully, 69134, France
# 1 Introduction

As a way of complying with international regulations aiming to achieve net carbon neutrality [1], aircraft engine manufacturers are working towards the development of increasingly efficient turbomachines. These systems are inherently nonlinear and should be studied with advanced modeling techniques. Nonlinearities play a central role in the dynamics of bladed disks. The blade/disk friction interface and the blade-tip/casing interface are among the most studied. The former is closely linked to the level of damping present in the system: quantifying the additional nonlinear damping brought by friction interfaces is of a strong interest for designers. Meanwhile, the latter should be addressed for safety purposes: the reduction of blade-tip/casing gaps yields a gain in aerodynamic efficiency but makes structural contacts more likely in nominal operating conditions.

Full scale experiments allowing to evaluate the dynamics of bladed disks under friction and blade-tip/casing contacts are too much of a financial burden for the manufacturers. Moreover, setting up such experiments can only come at the validation stage and can not serve as a priori verifications. Thus, accounting for nonlinear phenomena in early stages of design through numerical simulations is now mandatory. This raises the need for the development of industry-ready, robust, predictive and accurate numerical strategies for the characterization of these complex systems.

The dynamics of friction interfaces is widely studied through the Harmonic Balance Method (HBM) [2] in the literature [3, 4, 5, 6, 7, 8, 9, 10, 11]. Frequency domain strategies are well-suited for the analysis of industrial friction problems since the structural dynamics problem usually involves a great number of unknowns that prevents the use of direct time integration strategies. Numerical strategies usually differ from each other by the way the contact law is modeled inside the alternating frequency/time procedure [12]. One the one hand, penalty methods [4, 10, 11] yield a simpler implementation and faster computation times, but also come with modeling defects such as residual penetration and the impossibility of reaching a purely sticking state. On the other hand, the dynamic lagrangian frequency/time (DLFT) method [3, 7, 8, 9] (relying on augmented lagrangians), models more accurately the nonsmooth Coulomb Law and yields more predictive results. Naturally, this accuracy comes with more complex numerical developments, and oftentimes higher computation costs.

The study of the blade-tip/casing interface has recently known a renewed interest [13, 14, 15, 16, 17, 18] through the advent of the HBM as one of the main methodological tools of analysis of nonlinear systems. Indeed, it was shown that the HBM is able to achieve a similar level of accuracy than reference direct time integration algorithms [13, 17] on rubbing applications for fewer computing costs. Moreover, harmonic methods even outperform transient simulations [19, 20] since they allow to gather more information about the dynamic state of the system such as the connectivity of the different solutions [16, 17, 21], their stability [14, 18] or nonlinear modal parameters [22, 23, 24].

When characterizing industrial nonlinear systems, researchers oftentimes focus on a single nonlinearity (geometric [25], friction [26], blade-tip/casing contacts [27], material [28]). This choice is multifactorial and is made because numerous challenges arise when considering multiple nonlinearities or interfaces, namely implementation, interpretation, numerical conditioning and computing costs. However, while models featuring a single nonlinearity allow to perform relatively simple influence analyses, it only reflects a part of the dynamics response of the system. Indeed, strong couplings between the dynamics induced by nonlinearities are expected. This is why the next natural step in the study of nonlinear systems is to consider several nonlinearities in order to enhance the accuracy of the numerical models. As a first step towards this goal, Delhez et al. [29] have highlighted a complex coupling between geometric nonlinearities on blade-tip/casing contacts through transient simulations on a compressor blade. Zucca et al. [30] recently performed calculations on an academic model featuring a friction contact as well as geometrical nonlinearities. Some authors have performed calculation with multiple interfaces of the same nature: friction at shrouds and under platform dampers [31] or under platform dampers and blade attachments [32].

This paper aims to tackle the numerous challenges above-mentioned about the realization and interpretation of these ambitious simulations. This works studies the interdependence of the blade/disk friction interface and blade-tip/casing contacts for a state-of-the-art fan stage. The study includes resolutions through a frequency domain method that finds constant rotational speed solutions for a single sector without wear effects on the casing. The influence of the static state of the blade/disk interface on the blade-tip displacement is firstly underlined. Then both contact interfaces are analyzed independently in order to serve as references for the coupled case, where both contact interfaces are accounted for simultaneously. The responses are analyzed through local indicators at both
inter interfaces, as well as a global indicators. Then, numerous forced responses are computed to perform a sensitivity analysis on both the blade/disk friction coefficient and the level of excitation. Finally, a connection between the blade/disk interface dynamics and the number of contacting nodes at the blade-tip is highlighted.

The numerical strategy used throughout this article is described in Sec. 2. The industrial ECL5/Catana fan stage test case is presented in Sec. 3. A study of the static state of the system is performed in Sec. 4. The dynamics of the blade under a single nonlinearity, firstly blade-tip/casing contacts only and then blade/disk friction only, is analyzed in Sec. 5 to serve as references. Then, the coupled configuration, with both nonlinearities activated at once, is compared in depth to the reference cases in Sec. 6. In the end, a large scale sensitivity analysis over the blade/disk friction coefficient and the excitation level is conducted and commented in Sec. 7.

## 2 Numerical framework

In this work, the time and space normalized equation of motion for a n-DoF (degrees of freedom) nonlinear mechanical system reads

$$ \frac{\alpha}{\beta^2} \mathbf{M} \ddot{\mathbf{x}}(t) + \frac{\alpha}{\beta} \mathbf{C} \dot{\mathbf{x}}(t) + \alpha \mathbf{K} \mathbf{x}(t) = \mathbf{f}_{\text{ex}}(\omega, t) + \mathbf{f}_{0} + \mathbf{f}_{\text{nl}}(\mathbf{x}(t), \dot{\mathbf{x}}(t)) \quad (1) $$

where $\mathbf{x}$ is the normalized unknown displacement field, $\mathbf{M}$, $\mathbf{C}$ and $\mathbf{K}$ are the mass, damping and stiffness matrices respectively. $\mathbf{f}_{\text{ex}}$ is the vector of harmonic excitation forces of fundamental angular frequency $\omega$, $\mathbf{f}_{0}$ are static pre-loading forces accounting for centrifugal effects. Finally, $\mathbf{f}_{\text{nl}}$ is the vector of nonlinear forces that represents contact interactions in this work. Two types of contact are modeled here: blade/disk friction damping and blade-tip/casing unilateral contacts. In this work, the former is modeled by a Coulomb law [33] and the latter respects the Hertz-Moreau-Signorini conditions [33].

Equation Eqs. (1) is normalized in time and space by applying two linear changes of variables: (1) on the time variable so that the normalized time reads $t = \tau/\beta$ with $\tau$ being the physical time, and (2) on the displacement variable so that the normalized displacement is $\mathbf{x} = \mathbf{u}/\alpha$ with $\mathbf{u}$ the physical displacement. This normalized problem yields better numerical conditioning and hence the solver achieves better performances.

In order to solve Eq. (1), the HBM [2] is retained. In this work, only periodic solutions are sought, making the HBM particularly suited for the study of the presented problem. This numerical strategy and all useful related methodological tools are briefly recalled in Sec. 2.1 to put the presented work in its mathematical context.

### 2.1 Harmonic Balance Method

The HBM makes the assumption that, under a monoharmonic excitation $\mathbf{f}_{\text{ex}}(\omega, t) = \mathbf{a}_{1}^{\text{ex}} \cos(\omega t) + \mathbf{b}_{1}^{\text{ex}} \sin(\omega t)$, the periodic solutions of Eq. (1) can be written as Fourier series with $N_h$ harmonics, such that

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

with $\mathbf{a}_k$ and $\mathbf{b}_k$ being the real coefficients of the Fourier decomposition of $\mathbf{x}$. Carrying-out a Fourier-Galerkin projection of Eq. (1) on the Fourier basis $\mathcal{B}^{N_h} = [t \mapsto \frac{1}{2}, t \mapsto \cos(\omega t), t \mapsto \sin(\omega t), \dots, t \mapsto \cos(N_h\omega t), t \mapsto \sin(N_h\omega t)]$ yields an algebraic system of $n(2N_h + 1)$ nonlinear equations, usually concatenated in the following vectorial system

$$ \mathcal{H}(\tilde{\mathbf{x}}, \omega) = \mathbf{Z}(\omega)\tilde{\mathbf{x}} - \tilde{\mathbf{f}}_{\text{nl}}(\tilde{\mathbf{x}}) - \tilde{\mathbf{f}}_{\text{ex}} - \tilde{\mathbf{f}}_{0} = \mathbf{0}. \quad (3) $$

The new unknown of the problem is the so-called multiharmonic vector $\tilde{\mathbf{x}}$, composed of all Fourier coefficients of Eq. (2). Multiharmonic vectors $\tilde{\mathbf{f}}_{\text{nl}}$, $\tilde{\mathbf{f}}_{\text{ex}}$ and $\tilde{\mathbf{f}}_{0}$ of Eq. (1) consist of the Fourier coefficients of $\mathbf{f}_{\text{nl}}$, $\mathbf{f}_{\text{ex}}$ and $\mathbf{f}_{0}$ respectively. All multiharmonic vectors are structured as follows
$$ \tilde{\mathbf{x}} = \left[ \frac{1}{2} \mathbf{a}_0^\mathrm{T} \mathbf{a}_1^\mathrm{T} \mathbf{b}_1^\mathrm{T} \dots \mathbf{a}_{N_h}^\mathrm{T} \mathbf{b}_{N_h}^\mathrm{T} \right]^\mathrm{T} \quad (4a) $$

$$ \tilde{\mathbf{f}}_{nl} = \left[ \frac{1}{2} (\mathbf{a}_0^{nl})^\mathrm{T} (\mathbf{a}_1^{nl})^\mathrm{T} (\mathbf{b}_1^{nl})^\mathrm{T} \dots (\mathbf{a}_{N_h}^{nl})^\mathrm{T} (\mathbf{b}_{N_h}^{nl})^\mathrm{T} \right]^\mathrm{T} \quad (4b) $$

$$ \tilde{\mathbf{f}}_{\mathrm{ex}} = [\mathbf{0} \ (\mathbf{a}_1^{\mathrm{ex}})^{\mathrm{T}} \ (\mathbf{b}_1^{\mathrm{ex}})^{\mathrm{T}} \ \mathbf{0} \ \dots \ \mathbf{0} \ \mathbf{0}]^{\mathrm{T}} \quad (4c) $$

$$ \tilde{\mathbf{f}}_0 = [(\mathbf{a}_0^{\mathrm{0}})^{\mathrm{T}} \ \mathbf{0} \ \dots \ \mathbf{0} \ \mathbf{0}]^{\mathrm{T}} \quad (4d) $$

where $\bullet^\mathrm{T}$ is the transposition operation. The block-diagonal matrix $\mathbf{Z}$ is the dynamic stiffness matrix and accounts for the linear behavior of the system. Its $(N_h + 1)$ blocks are defined such as

$$ \mathbf{Z}_0 = \alpha \mathbf{K} \quad \text{and} \quad \mathbf{Z}_k(\omega) = \alpha \begin{bmatrix} \mathbf{K} - \left(\frac{k\omega}{\beta}\right)^2 \mathbf{M} & \frac{k\omega}{\beta} \mathbf{C} \\ -\frac{k\omega}{\beta} \mathbf{C} & \mathbf{K} - \left(\frac{k\omega}{\beta}\right)^2 \mathbf{M} \end{bmatrix} \quad \forall k \in [1, N_h]. \quad (5) $$

## 2.2 Hybrid condensation procedure

As a way of reducing the computational cost of solving Eq. (3), two exact condensation procedures are employed on the nonlinear system of equations.

The first condensation relies on the property that the linear DoF, not subject to any nonlinear forces, of the system can be reduced, yielding a significant reduction in size of the actual system to be solved. This procedure is fairly common when accounting for localized nonlinearities such as contact. Thus, it is presented in A. After this step, the unknown of the problem is now $\tilde{\mathbf{x}}_{nl}$ where only nonlinear DoF remain and the system reads

$$ \mathbf{Z}_{nl} \tilde{\mathbf{x}}_{nl} = \tilde{\mathbf{f}}_{\mathrm{ex},nl} + \tilde{\mathbf{f}}_{0,nl} + \tilde{\mathbf{f}}_{nl,nl}. \quad (6) $$

In this work, two types of contact interfaces are accounted for: (1) a friction interface at the blade root and (2) a rotor/stator interaction (RSI) interface. The friction interface is modeled through the DLFT procedure [3] and thus requires a condensation into relative coordinates, whereas the RSI interface should remain in absolute coordinates. As a way of setting-up the system in the right form, a so-called hybrid condensation step should be performed. Friction related DoF and RSI related DoF are respectively superscripted $\bullet$^fri and $\bullet$^R/S$. The relative condensation relies on a pairing of DoF from either solids involved in the friction interaction, here from the blade and the disk, respectively identified by subscripts $\bullet$^b and $\bullet$^d$. Following these notations, the vector containing all nonlinear DoF reads $\tilde{\mathbf{x}}_{nl} = [(\tilde{\mathbf{x}}_b^{\mathrm{fri}})^{\mathrm{T}}, (\tilde{\mathbf{x}}_d^{\mathrm{fri}})^{\mathrm{T}}, (\tilde{\mathbf{x}}^{\mathrm{R/S}})^{\mathrm{T}}]^\mathrm{T}$. The condensation procedure relies on the following relation, linking $\tilde{\mathbf{x}}_{nl}$ to the new unknown $\tilde{\mathbf{x}}_r$

$$ \tilde{\mathbf{x}}_r = \begin{bmatrix} \tilde{\mathbf{x}}_{\mathrm{rel}}^{\mathrm{fri}} \\ \tilde{\mathbf{x}}^{\mathrm{R/S}} \end{bmatrix} = \begin{bmatrix} \mathbf{I}_{n^{\mathrm{fri}}} & -\mathbf{I}_{n^{\mathrm{fri}}} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{I}_{n^{\mathrm{R/S}}} \end{bmatrix} \begin{bmatrix} \tilde{\mathbf{x}}_b^{\mathrm{fri}} \\ \tilde{\mathbf{x}}_d^{\mathrm{fri}} \\ \tilde{\mathbf{x}}^{\mathrm{R/S}} \end{bmatrix} = \mathbf{P}^{\mathrm{T}} \tilde{\mathbf{x}}_{nl} \quad (7) $$

where $n^{\mathrm{fri}}$ is the number of DoF of the friction interface, $n^{\mathrm{R/S}}$ the number of DoF of the RSI interface and $\mathbf{I}_\bullet$ the identity matrix of size $\bullet$. Due to the reciprocity of the nonlinear forces on both sides of the friction interface, the matrix $\mathbf{P}$ also links the different nonlinear forces such as

$$ \tilde{\mathbf{f}}_{nl,nl} = \begin{bmatrix} \tilde{\mathbf{f}}_{nl,b}^{\mathrm{fri}} \\ \tilde{\mathbf{f}}_{nl,d}^{\mathrm{fri}} \\ \tilde{\mathbf{f}}_{nl}^{\mathrm{R/S}} \end{bmatrix} = \mathbf{P} \begin{bmatrix} \tilde{\mathbf{f}}_{nl}^{\mathrm{fri}} \\ \tilde{\mathbf{f}}_{nl}^{\mathrm{R/S}} \end{bmatrix} = \mathbf{P} \tilde{\mathbf{f}}_{nl,r} \quad (8) $$

where $\tilde{\mathbf{f}}_{nl}^{\mathrm{fri}} = \tilde{\mathbf{f}}_{nl,b}^{\mathrm{fri}} = -\tilde{\mathbf{f}}_{nl,d}^{\mathrm{fri}}$. By multiplying Eq. (6) by $\mathbf{Z}_{nl}^{-1}$ and substituting $\tilde{\mathbf{f}}_{nl,nl}$ by Eq. (8), it leads to

$$ \tilde{\mathbf{x}}_{nl} = \mathbf{Z}_{nl}^{-1} (\tilde{\mathbf{f}}_{\mathrm{ex},nl} + \tilde{\mathbf{f}}_{0,nl} + \mathbf{P}\tilde{\mathbf{f}}_{nl,r}). \quad (9) $$
The reduced unknown $\tilde{\mathbf{x}}_r$ then appears when Eq. (9) is multiplied by $\mathbf{P}^\top$

$$ \tilde{\mathbf{x}}_r = \mathbf{P}^\top \mathbf{Z}_{nl}^{-1} (\tilde{\mathbf{f}}_{\text{ex,nl}} + \tilde{\mathbf{f}}_{0,nl}) + \mathbf{P}^\top \mathbf{Z}_{nl}^{-1} \mathbf{P} \tilde{\mathbf{f}}_{\text{nl,r}} \quad (10) $$

Multiplying Eq. (10) by $\mathbf{Z}_r = (\mathbf{P}^\top \mathbf{Z}_{nl}^{-1} \mathbf{P})^{-1}$ allows to retrieve the classical form of the frequency domain equation of motion

$$ \mathbf{Z}_r \tilde{\mathbf{x}}_r = \tilde{\mathbf{f}}_r + \tilde{\mathbf{f}}_{nl,r}, \quad (11) $$

with $\tilde{\mathbf{f}}_r = \mathbf{Z}_r \mathbf{P}^\top \mathbf{Z}_{nl}^{-1} (\tilde{\mathbf{f}}_{\text{ex,nl}} + \tilde{\mathbf{f}}_{0,nl})$. For the sake of reducing computational costs, this condensation is applied independently for each harmonic since $\mathbf{Z}_{nl}$ is block-diagonal. However, in the case where one of both solids involved in the friction interface is subject to rigid body modes (a blade for example), a specific procedure should be employed for the static problem (harmonic 0). The associated developments are given in B.

## 2.3 Evaluation of nonlinear forces

When using the HBM, the evaluation of the nonlinear forces in the frequency domain generally requires the use of the alternating frequency/time (AFT) strategy [12]. The latter consists in estimating the nonlinear forces $\mathbf{f}_{nl}(t)$ in the time domain from a given frequency-domain displacement $\tilde{\mathbf{x}}$ and then transforming the time-dependent nonlinear forces into the frequency domain $\tilde{\mathbf{f}}_{nl}$. The step of going from the frequency domain to the time domain is done by using the inverse discrete Fourier transform (IDFT) matrix $\tilde{\mathcal{F}}$. The transformation from the time domain to the frequency domain is done by using the discrete Fourier transform (DFT) matrix $\mathcal{F}$. These matrices are used to transform time-discrete vectors into multiharmonic vectors $\tilde{\mathbf{\bullet}} = \mathcal{F}\{\bullet(t_i)\}_{i=1...N_t}$ and vice-versa $\{\bullet(t_i)\}_{i=1...N_t} = \tilde{\mathcal{F}}\tilde{\mathbf{\bullet}}$.

The number of instants for the discretization of a period in the time domain, noted $N_t$, is a user-defined parameter. In this work, $N_t$ is chosen to be different according to the type of interface, as a compromise between convergence and computational efficiency. However, its value does not strongly affect the solutions of the problem since the residual is still computed in the frequency domain. Usually, having $N_t$ as large as possible favors the convergence of the nonlinear solver and prevents aliasing problems [34].

The strategies to evaluate the nonlinear forces on the friction and RSI interface both rely on the previously published dynamic Lagrangian frequency/time [3] (DLFT) algorithm and were thoroughly characterized separately on both friction [32, 35] and vibro-impact applications [17, 36].

Since the unknown vector $\tilde{\mathbf{x}}_r$ is composed of two interfaces, Eq. (11) is also partitioned such as

$$ \begin{bmatrix} \mathbf{Z}_r^{\text{fri}} \\ \mathbf{Z}_r^{\text{R/S}} \end{bmatrix} \begin{bmatrix} \tilde{\mathbf{x}}_{\text{rel}}^{\text{fri}} \\ \tilde{\mathbf{x}}^{\text{R/S}} \end{bmatrix} = \begin{bmatrix} \tilde{\mathbf{f}}_r^{\text{fri}} \\ \tilde{\mathbf{f}}_r^{\text{R/S}} \end{bmatrix} + \begin{bmatrix} \tilde{\mathbf{f}}_{\text{nl,r}}^{\text{fri}} \\ \tilde{\mathbf{f}}_{\text{nl,r}}^{\text{R/S}} \end{bmatrix}. \quad (12) $$

The procedure for the evaluation of $\tilde{\mathbf{f}}_{nl,r}^{\text{R/S}}$ and $\tilde{\mathbf{f}}_{nl,r}^{\text{fri}}$ are explicit in Sec. 2.3.1 and Sec. 2.3.2 respectively. For the sake of concision, only the formulas related to the forces themselves are recalled here: more details can be found on these strategies in the previously published work [17, 37], in particular for the computation of Jacobian matrices. $\mathbf{Z}_r^{\text{fri}}$ and $\mathbf{Z}_r^{\text{R/S}}$ are rectangular and stand for all the linear couplings (internal forces) between the two interfaces.

### 2.3.1 Blade-tip/casing contact

The strategy for the evaluation of blade-tip/casing contact forces relies on the previously published DLFT [3, 36]. It relies on a prediction-correction scheme detailed hereafter. The second line of Eq. (12) is partitioned according to the three directions of the cylindrical coordinate system such that

$$ \begin{bmatrix} \mathbf{Z}_r^{\text{R/S,N}} \\ \mathbf{Z}_r^{\text{R/S},\theta} \\ \mathbf{Z}_r^{\text{R/S},z} \end{bmatrix} \tilde{\mathbf{x}}_r = \begin{bmatrix} \tilde{\mathbf{f}}_r^{\text{R/S,N}} \\ \tilde{\mathbf{f}}_r^{\text{R/S},\theta} \\ \tilde{\mathbf{f}}_r^{\text{R/S},z} \end{bmatrix} + \begin{bmatrix} \tilde{\mathbf{f}}_{nl,r}^{\text{R/S,N}} \\ \tilde{\mathbf{f}}_{nl,r}^{\text{R/S},\theta} \\ \tilde{\mathbf{f}}_{nl,r}^{\text{R/S},z} \end{bmatrix}. \quad (13) $$
where circumferential and axial directions are superscripted ●θ and ●z respectively. The radial direction, noted by the superscript ●N, is the only direction handled with augmented lagrangian. Thus, the prediction-superscripted ●P-only concerns normal contact forces. This prediction is made in the frequency domain following

$$ \tilde{\mathbf{f}}_{\text{nl,r}}^{\text{R/S,N,p}} = \mathbf{Z}_{\text{r}}^{\text{R/S,N}} \tilde{\mathbf{x}}_{\text{r}} - \tilde{\mathbf{f}}_{\text{r}}^{\text{R/S,N}} - \varepsilon^{\text{R/S}} \tilde{\mathbf{g}}_{\text{x}}. \quad (14) $$

where $\tilde{\mathbf{g}}_{\text{x}}$ is a multiharmonic vector standing for the normal gap function $g_j(t) = g_{0,j} - x_j^{\text{R/S,N}}(t)$ of each contact node j and $\varepsilon^{\text{R/S}}$ a weighting parameter responsible for the respect of the unilateral contact conditions. Note that this prediction depends on the dynamics of the other interface thanks to the term $\mathbf{Z}_{\text{r}}^{\text{R/S,N}} \tilde{\mathbf{x}}_{\text{r}}$ since $\tilde{\mathbf{x}}_{\text{r}} = [\tilde{\mathbf{x}}_{\text{rel}}^{\text{fri}} \tilde{\mathbf{x}}^{\text{R/S}}]^{\top}$. The predicted contact force of Eq. (14) is in the frequency domain and encompasses all contact nodes. Once this prediction is made, contact forces $f_{\text{nl},j}^{\text{R/S,N}}(t)$ of each node j are handled separately in the time domain. At a given instant of the time discretization $t_i$, two cases may arise based on the predicted state of contact:

**Separation** if $f_{\text{nl},j}^{\text{R/S,N,p}}(t_i) > 0$, the jth node is separated at time $t_i$ and:

$$ f_{\text{nl},j}^{\text{R/S,N}}(t_i) = 0, \quad (15a) $$

$$ f_{\text{nl},j}^{\text{R/S},\theta}(t_i) = 0, \quad (15b) $$

$$ f_{\text{nl},j}^{\text{R/S},z}(t_i) = 0. \quad (15c) $$

**Contact** if $f_{\text{nl},j}^{\text{R/S,N,p}}(t_i) < 0$, the jth node is in contact at time $t_i$ and:

$$ f_{\text{nl},j}^{\text{R/S,N}}(t_i) = f_{\text{nl},j}^{\text{R/S,N,p}}(t_i), \quad (16a) $$

$$ f_{\text{nl},j}^{\text{R/S},\theta}(t_i) = -\mu^{\text{R/S}} \frac{v_j^{\theta}(t_i) + \rho_j \Omega}{\sqrt{(v_j^{\theta}(t_i) + \rho_j \Omega)^2 + v_j^{z}(t_i)^2}} |f_{\text{nl},j}^{\text{R/S,N}}(t_i)|, \quad (16b) $$

$$ f_{\text{nl},j}^{\text{R/S},z}(t_i) = -\mu^{\text{R/S}} \frac{v_j^{z}(t_i)}{\sqrt{(v_j^{\theta}(t_i) + \rho_j \Omega)^2 + v_j^{z}(t_i)^2}} |f_{\text{nl},j}^{\text{R/S,N}}(t_i)|. \quad (16c) $$

$\mu^{\text{R/S}}$ is the friction coefficient of the RSI interface, $v_j^{\theta}$ and $v_j^{z}$ are circumferential and axial vibrational velocities respectively, $\rho_j$ is the radial distance between the jth contact node and the rotational axis and $\Omega$ the rotational speed. For this type of interface, the tangential friction forces are computed with the hypothesis that all contacting nodes are in a sliding state, due to the high relative velocities at stake.

After looping on all nodes and all instants, the RSI contact forces are fully evaluated and are then concatenated and transformed back to the frequency domain, yielding the term $\tilde{\mathbf{f}}_{\text{nl,r}}^{\text{R/S}}$ of Eq. (12).

### 2.3.2 Blade/disk friction

Similarly, the strategy for the evaluation of friction contact forces relies on the DLFT algorithm [3]. The contact model used here allows to deal with contact constraints as a node-to-node approach, but other methods relying on microslip contact models are also possible [38] for the description of the contact interface. However, the evaluation of these contact models exceeds the scope of this paper. Here, contrary to RSI, the computation of contact forces in all three directions are handled by the prediction-correction scheme, so that $\tilde{\mathbf{f}}_{\text{nl,r}}^{\text{fri,p}} = \tilde{\mathbf{f}}_{\text{nl,r}}^{\text{fri,c}} + \tilde{\mathbf{f}}_{\text{nl,r}}^{\text{fri,p}}$ with $\tilde{\mathbf{f}}_{\text{nl,r}}^{\text{fri,p}}$ being the prediction and $\tilde{\mathbf{f}}_{\text{nl,r}}^{\text{fri,c}}$ the correction. The contact forces in all directions are predicted for all contact nodes following

$$ \tilde{\mathbf{f}}_{\text{nl,r}}^{\text{fri,p}} = \mathbf{Z}_{\text{r}}^{\text{fri}} \tilde{\mathbf{x}}_{\text{r}} - \tilde{\mathbf{f}}_{\text{r}}^{\text{fri}} - \varepsilon^{\text{fri}} \tilde{\mathbf{x}}_{\text{rel}}^{\text{fri}}. \quad (17) $$

Similarly as for the blade-tip/casing contact forces, the dynamics of the other interface is accounted through to the term $\mathbf{Z}_{\text{r}}^{\text{fri}} \tilde{\mathbf{x}}_{\text{r}}$ because $\tilde{\mathbf{x}}_{\text{r}} = [\tilde{\mathbf{x}}_{\text{rel}}^{\text{fri}} \tilde{\mathbf{x}}^{\text{R/S}}]^{\top}$. For the case of friction, the predictions of tangential forces need to be
re-estimated at each instant based on the value of the correction computed at the last instant, so that the adjusted prediction-superscripted ●p²-reads

$$ \left\{ \begin{array}{l} \mathbf{f}_{nl,j}^{\text{fri,T,p2}}(t_i) = \mathbf{f}_{nl,j}^{\text{fri,T,p}}(t_i) - \mathbf{f}_{nl,j}^{\text{fri,T,c}}(t_{i-1}), \\ \mathbf{f}_{nl,j}^{\text{fri,N,p2}}(t_i) = \mathbf{f}_{nl,j}^{\text{fri,N,p}}(t_i), \end{array} \right. \qquad (18a) $$

$$ \left\{ \begin{array}{l} \mathbf{f}_{nl,j}^{\text{fri,N,p2}}(t_i) = \mathbf{f}_{nl,j}^{\text{fri,N,p}}(t_i), \\ \mathbf{f}_{nl,j}^{\text{fri,T,p2}}(t_i) = \mathbf{f}_{nl,j}^{\text{fri,T,p}}(t_i), \end{array} \right. \qquad (18b) $$

making the assumption that the contact element is in a sticking state. For a given node j at time t_i, three cases then arise according to the values of $f_{nl,j}^{\text{fri,N,p2}}(t_i)$ and $\mathbf{f}_{nl,j}^{\text{fri,T,p2}}(t_i)$:

**Separation** if $f_{nl,j}^{\text{fri,N,p2}}(t_i) > 0$, the jth node is separated at time $t_i$ and:

$$ \left\{ \begin{array}{l} f_{nl,j}^{\text{fri,N}}(t_i) = 0, \\ \mathbf{f}_{nl,j}^{\text{fri,T}}(t_i) = \mathbf{0}. \end{array} \right. \qquad (19a) $$

$$ \left\{ \begin{array}{l} f_{nl,j}^{\text{fri,N}}(t_i) = 0, \\ \mathbf{f}_{nl,j}^{\text{fri,T}}(t_i) = \mathbf{0}. \end{array} \right. \qquad (19b) $$

The corrections are time $t_i$ then read $\mathbf{f}_{nl,j}^{\text{fri,T,c}}(t_i) = -\mathbf{f}_{nl,j}^{\text{fri,T,p}}(t_i)$.

**Sticking** if $f_{nl,j}^{\text{fri,N,p2}}(t_i) < 0$ and $\|\mathbf{f}_{nl,j}^{\text{fri,T,p2}}(t_i)\| < \mu|f_{nl,j}^{\text{fri,N,p2}}(t_i)|$, the jth node is sticking at time $t_i$ and:

$$ \left\{ \begin{array}{l} f_{nl,j}^{\text{fri,N}}(t_i) = f_{nl,j}^{\text{fri,N,p2}}(t_i), \\ \mathbf{f}_{nl,j}^{\text{fri,T}}(t_i) = \mathbf{f}_{nl,j}^{\text{fri,T,p2}}(t_i). \end{array} \right. \qquad (20a) $$

$$ \left\{ \begin{array}{l} f_{nl,j}^{\text{fri,N}}(t_i) = f_{nl,j}^{\text{fri,N,p2}}(t_i), \\ \mathbf{f}_{nl,j}^{\text{fri,T}}(t_i) = \mathbf{f}_{nl,j}^{\text{fri,T,p2}}(t_i). \end{array} \right. \qquad (20b) $$

In this case, the correction at time $t_i$ is obtained through $\mathbf{f}_{nl,j}^{\text{fri,T,c}}(t_i) = \mathbf{f}_{nl,j}^{\text{fri,T,c}}(t_{i-1})$.

**Sliding** if $f_{nl,j}^{\text{fri,N,p2}}(t_i) < 0$ and $\|\mathbf{f}_{nl,j}^{\text{fri,T,p2}}(t_i)\| \ge \mu|f_{nl,j}^{\text{fri,N,p2}}(t_i)|$, the jth node is sliding at time $t_i$ and:

$$ \left\{ \begin{array}{l} f_{nl,j}^{\text{fri,N}}(t_i) = f_{nl,j}^{\text{fri,N,p2}}(t_i), \\ \mathbf{f}_{nl,j}^{\text{fri,T}}(t_i) = -\mu|f_{nl,j}^{\text{fri,N,p2}}(t_i)| \frac{\mathbf{f}_{nl,j}^{\text{fri,T,p2}}(t_i)}{\|\mathbf{f}_{nl,j}^{\text{fri,T,p2}}(t_i)\|} \end{array} \right. \qquad (21a) $$

$$ \left\{ \begin{array}{l} f_{nl,j}^{\text{fri,T}}(t_i) = -\mu|f_{nl,j}^{\text{fri,N,p2}}(t_i)| \frac{\mathbf{f}_{nl,j}^{\text{fri,T,p2}}(t_i)}{\|\mathbf{f}_{nl,j}^{\text{fri,T,p2}}(t_i)\|}. \end{array} \right. \qquad (21b) $$

The value of the correction at time $t_i$ is

$$ \mathbf{f}_{nl,j}^{\text{fri,T,c}}(t_i) = \mathbf{f}_{nl,j}^{\text{fri,T,c}}(t_{i-1}) + \mathbf{f}_{nl,j}^{\text{fri,T,p2}}(t_i) \left( 1 - \mu \frac{|f_{nl,j}^{\text{fri,N,p2}}(t_i)|}{\|\mathbf{f}_{nl,j}^{\text{fri,T,p2}}(t_i)\|} \right). \qquad (22) $$

$\mu^{\text{fri}}$ is the friction coefficient of the blade/disk friction interface. The value of the correction at the first instant is initialized such as $\mathbf{f}_{nl,j}^{\text{fri,T,c}} = \mathbf{0}$. The correction is then updated at each instant, and the hypothesis made at the first time vanishes after looping over two full periods. After looping on all nodes and all instants, the estimated contact forces are concatenated and transformed to the frequency domain, fully defining the term $\tilde{\mathbf{f}}_{nl,r}^{\text{fri}}$ of Eq. (12). The underlying mathematical background exceeds the scope of this article and can be found in the literature [3].

## 2.4 Continuation procedure

In this work, due to the highly nonlinear nature of the simulated physics, the nonlinear frequency response curves (NFRC) may exhibit turning points and are thus computed by using an arc-length parameterization [39] in forced response. Similarly, nonlinear modal analyses of the system, that are conducted thanks to complex nonlinear modes [22], are computed with an arc-length parameterization. Complex nonlinear modes involve a nonlinear damping ratio $\zeta$ that is a new unknown of the problem.

In order to perform such continuations, the computation of analytical gradients of the nonlinear forces with respect to both the excitation pulsation $\omega$ and the nonlinear damping $\zeta$ are used. For the sake of concision, their expressions are not recalled in this article but can be found in the literature [36, 37].
# 3 Industrial application: ECL5 Catana blade model

## 3.1 Finite-elements model

The industrial model retained for this study is the open-test-case fan stage ECL5/Catana designed at Ecole Centrale de Lyon [40] to be representative of near-future fan blades concepts in ultrahigh bypass ratio architectures. The stage is composed of two parts: (1) 16 blades made out of unidirectional carbon fibers and epoxy composite plies and (2) a titanium disk. The blade geometry was extensively studied in both experimental and numerical conditions from aerodynamic [41], aeroelastic [42], acoustic [43] and structural dynamics [44] perspectives.

The full finite-element mesh is composed of a total of 19,642 elements and 87,310 nodes and is displayed in Fig. 1a. The material data associated with the blade composite is given in Ref. [40]. The boundary conditions applied to the disk are the same as in previously published studies [40] and displayed in Fig. 1b. In this work, a single blade attached to a portion of the disk is studied to focus on the interaction of both contact interfaces. Analyzing the full structure is beyond the scope of this paper. The disk cyclic boundaries are thus clamped.

<figure><img src="image_1.png" /><figcaption>**Figure 1.** ECL5/Catana fan stage model, clamped (■), axial blockage (■).</figcaption></figure>

For the sake of reducing computational costs, a reduced-order model is computed by means of the Craig-Bampton procedure [45] on each substructure. For the disk, 54 boundary nodes are retained at the root alongside 20 internal modes. The reduced-order model of the blade includes (a) $n^{\text{fri}} = 54$ nonlinear boundary nodes at the root for friction management purposes, (b) $n^{\text{R/S}} = 9$ nonlinear boundary nodes at the tip of the blade for RSI management purposes, (c) 1 linear excitation node at the tip of the blade and (d) 40 internal modes. The first mode shapes of the blade/disk assembly are displayed for illustration on Fig. 2 and the value of the first eigenfrequencies are provided in Tab. 1. Rayleigh damping is used for this study. The damping matrix is $C = a\mathbf{K} + b\mathbf{M}$ with $a = 5.488 \cdot 10^{-7} (\text{rad/s})^{-1}$ and $b = 1.822 \text{ rad} \cdot \text{s}^{-1}$, yielding a damping ratio of $\xi_{1B} = 1 \cdot 10^{-3}$ for the first bending mode.

## 3.2 Contact scenario

In this work, the dynamics of the casing is neglected at the RSI interface. Thus, the casing is considered to be a rigid mathematical profile. The static equilibrium used to set the position of the casing is obtained after performing a static nonlinear calculation where the root friction interface is preloaded due to the presence of $f_0$ in
<figure><img src="image_2.png" /><figcaption>**Figure 2.** Modal shapes of the blade.</figcaption></figure>

<table><caption>**Table 1.** First eigenfrequencies of the blade.</caption><thead><tr><th>mode •</th><th>1B</th><th>2B</th><th>1T</th></tr></thead><tbody><tr><td>ω• (rad·s<sup>-1</sup>)</td><td>1,834.7</td><td>4,366.8</td><td>5,691.9</td></tr><tr><td>f• (Hz)</td><td>292.0</td><td>695.0</td><td>905.9</td></tr></tbody></table>

Eq. (1), modeling the centrifugal forces. In nominal conditions, the system is in rotation at an angular speed of $\Omega = 1,151.9 \text{ rad·s}^{-1} = 11,000 \text{ rpm}$.

The casing is chosen to be at a distance of $c_j = 4 \cdot 10^{-4} \text{ m}$ for all nodes at the blade-tip on its static equilibrium. The values chosen reflect the same order of magnitude that was measured experimentally on the PHARE\#2 test rig [41]. During their experimental work, no RSI was observed. This study proposes to numerically trigger the RSI dynamics by applying higher level of excitations. The main focus is to study the dynamics of the system including two interfaces of different nature. This work aims to explore the parameters that can help preventing or mitigating RSI for mechanical system under extreme conditions.

In order to progressively trigger the activation of the dynamics of the two nonlinearities, excitation forces of varying amplitudes are applied to the blade. The harmonic excitation of amplitude noted $F$ is applied to a single node at the tip of the blade in the circumferential direction, thus favoring the response of bending modes. The excitation is synchronous with the rotation of the blade: it is applied for an angular frequency $\omega = \Omega$. Firstly, the excitation forces activate the friction interface and then, when the radial vibrational amplitude exceeds the gap $c_j$ on one of the blade-tip nonlinear nodes, RSI is activated as well. The RSI contact initiation scenario is illustrated in Fig. 3. This study does not include any wear effects on the casing. Indeed, the analysis of wear patterns through
frequency domain methods remains an open problem to this day.

In this work, the behavior of the system is studied around its first bending mode frequency of $f_{1B} = 292.0$ Hz. The HBM simulations are carried out by assuming that constant angular speed solutions exist. This hypothesis is supported by the fact that numerical-experimental correlations were already made with such assumptions with a time integration numerical tool [46]. This same tool was recently confronted to frequency domain methods [17] showing good agreement. Then, it is considered that the employed frequency domain strategy is relevant even if it makes an assumption about the rotation speed of the nonlinear solutions.

<figure><img src="image_3.png" /><figcaption>**Figure 3.** RSI contact initiation scenario, casing (□) facing the jth boundary node of the blade, after preloading (○), before preloading (○).</figcaption></figure>

# 4 Static results

Before going into depth into the analysis of the coupling between the dynamics of friction and RSI nonlinearities, the evolution of the static equilibrium of the blade is analyzed with respect to the friction coefficient of the blade/disk interface $\mu^{\text{fri}}$. Firstly, the normal force at the friction interface is assessed and secondly the influence of the blade/disk friction coefficient on the blade-tip gap is evaluated.

## 4.1 Normal force at the friction interface

The static state of contact is a key information when it comes to the evaluation of friction interactions: the sliding limit is directly proportional to the normal force applied to the interface. As a way of validating the developed nonlinear tool, the normal force obtained in this work through DLFT on the reduced-order model is confronted to the reference normal force obtained by the commercial software ANSYS on the full mesh. As the value of friction coefficient $\mu^{\text{fri}}$ is rarely precisely known, a comparison is here carried out to evaluate the variation of normal force at the interface with respect to $\mu$ and displayed in Tab. 2. The quantity referred to as the normal force is the absolute value of the contact force in the normal direction of the interfaces summed on both sides of the blade attachment.

<table><caption>**Table 2.** Evolution of the normal static force on the blade root with respect to the friction coefficient.</caption><thead><tr><th>root friction coefficient &mu;<sup>fri</sup></th><th>0.15</th><th>0.3</th><th>0.45</th><th>0.6</th><th>0.75</th></tr></thead><tbody><tr><td>reference (ANSYS)</td><td>25,124 N</td><td>21,268 N</td><td>18,420 N</td><td>16,176 N</td><td>14,339 N</td></tr><tr><td>nonlinear tool</td><td>24,268 N</td><td>20,187 N</td><td>17,408 N</td><td>15,313 N</td><td>13,644 N</td></tr><tr><td>relative deviation</td><td>3.4%</td><td>5.1%</td><td>5.5%</td><td>5.3%</td><td>4.8%</td></tr></tbody></table>

A good agreement is found between the two strategies and the relative deviation is relatively small. Though, these differences are attributed to (1) the fact that some information about the centrifugal force is lost during the projection on the Craig-Bampton basis and (2) significantly less contact elements are retained for computational efficiency purposes in the reduced-order model. Hence, these results are considered to be acceptable in the chosen framework.
Note that the range of root friction coefficient chosen for the study is deliberately wide. Indeed, having such a variability of this parameter allows to observe how the overall nonlinear dynamics of the blade is affected by a large variation of $\mu^{\text{fri}}$. This shows that the numerical strategy allows to have a full control over the blade's dynamics even for unrealistic friction coefficients.

These results show a relatively high decrease of the normal force, around 40%, for the studied range of friction coefficient $\mu^{\text{fri}}$. As a way of validating this observation, a simplified analytical model is built in Fig. 4a. In order to evaluate the normal force on both sides of the interface, it is necessary to study the static equilibrium of the system in the radial direction and it reads

$$f_c = (N_1 + N_2) \sin \theta + (T_1 + T_2) \cos \theta, \qquad (23)$$

where $N_1$, $N_2$, $T_1$ and $T_2$ are the magnitude of the normal and tangent forces at the interface (subscripts ●1 and ●2 refer to both sides of the attachment). Assuming the problem is symmetric ($N_1 = N_2 = N$ and $T_1 = T_2 = T$) and that the contact is fully sliding ($T = \mu N$)-as observed numerically-, it leads to

$$f_c = 2N \sin \theta + 2\mu N \cos \theta \Rightarrow N = \frac{f_c}{2(\sin \theta + \mu \cos \theta)} \qquad (24)$$

One can note that due to the symmetry assumption, contact forces balance each other out in the circumferential direction. The numerical application of Eq. (24) in the configuration of this article, where $f_c = 15,724$ N and $\theta = 30^\circ$, is displayed in Fig. 4b alongside the values obtained by the two numerical strategies. It is observed that the small

<figure><img src="image_4.png" /><figcaption>**Figure 4.** Analytical model for static normal force computation. (a) Schematic view of the nonlinear static equilibrium at the root of the blade, (b) normal force at the root of the blade (2N), reference (ANSYS) (—), nonlinear tool (—), analytic (---), analytic sliding limit $2\mu^{\text{fri}}N$ (---).</figcaption></figure>

model developed is actually accurate and the trend that was observed numerically is validated by the analytical model. It also validates the hypothesis according to which the contact is fully sliding. However, while the normal force is reduced by 40%, the sliding limit of the interface $2\mu^{\text{fri}}N(\mu^{\text{fri}})$ (---) increases significantly with respect to $\mu^{\text{fri}}$.

The properties of the deformed contact interface can be visualized locally in Fig. 5, where the friction status of all nodes obtained by DLFT is displayed for different values of $\mu^{\text{fri}}$. On this figure, the blade root is seen from above and the left part of the graph corresponds to the pressure side while the right part is the suction side. Nodes towards the leading edge (LE) are located at the bottom of the figure and at the top for the trailing edge (TE). Each marker corresponds to a contact element and the size of each marker is related to the value of nodal normal force: the bigger the marker is, the more loaded the node is. The scales of marker size are not shared among the different subfigures of Fig. 5 and must be read independently.

Firstly, Fig. 5 shows that the centrifugal forces induces a specific kinematic to the blade/disk assembly: the blade is tilting towards the pressure side for all values of $\mu^{\text{fri}}$. It can be seen by the fact that the areas where the normal
<figure><img src="image_5.png" /><figcaption>**Figure 5.** Evolution of the static state of contact with respect to the root friction coefficient $\mu^{\text{fri}}$, sliding (●), sticking (■), the size of markers is proportional to the intensity of the nodal normal force.</figcaption></figure>

force is the highest are (1) the lower part on the suction side and (2) the higher part on the pressure side. This means that the rotation induces a static bending deformation to the blade towards the direction of rotation.

Secondly, it is observed that all nodes are sliding in static for $\mu^{\text{fri}} = 0.15$ and around 90% of the nodes are sliding for higher values of $\mu^{\text{fri}}$ which validates once again the hypothesis made on the static state of contact of the interface under centrifugal forces. Note that the notion of static sliding does not imply that the interface will be sliding in dynamics as well. The coupling between harmonics 0 (static) and 1 (dynamics) induces an overall sticking behavior at low amplitudes.

## 4.2 Blade-tip displacement

After pre-loading, i.e. when the blade is in rotation, it is necessary to monitor the tip gap in order to ensure that no blade-tip/casing contacts can occur in nominal operating conditions. As a way of understanding the nonlinear dependence between root friction and the tip gap, a parametric study is led. For all nodes of the blade-tip interface, the evolution of the radial static displacement is displayed in Fig. 6 alongside its variability on the range of friction coefficients studied. On the one hand, the curves represent the variation of the radial displacement $r_j(\mu^{\text{fri}})$ for all of the 9 contact nodes. On the other hand, the strips of same color are just highlighting the interval over which each $r_j$ spans on the interval of $\mu^{\text{fri}} \in [0.15, 0.75]$. Its width is then expressed as $\Delta r_j = \left| \max_{\mu^{\text{fri}} \in [0.15, 0.75]} r_j(\mu^{\text{fri}}) - \min_{\mu^{\text{fri}} \in [0.15, 0.75]} r_j(\mu^{\text{fri}}) \right|$.

The results show that, for almost all nodes, the blade-tip displacement decreases when the friction coefficient increases. The exact levels of variation from 0.15 to 0.75 are extracted from the color strips in Fig. 6 and summed up in Tab. 3.

<table><caption>**Table 3.** Variation of the static blade-tip displacement with respect to the friction coefficient (tip gap is $c_j = 4 \cdot 10^{-4}$ m here).</caption><thead><tr><th>node</th><th>LE</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>TE</th></tr></thead><tbody><tr><td>radial disp. variability Δr<sub>j</sub> (10<sup>-5</sup>m)</td><td>3.05</td><td>3.35</td><td>4.50</td><td>5.25</td><td>6.00</td><td>6.74</td><td>7.49</td><td>8.28</td><td>9.07</td></tr><tr><td>Δr<sub>j</sub>% = 100 · Δr<sub>j</sub>/c<sub>j</sub> (%)</td><td>7.6</td><td>8.4</td><td>11.3</td><td>13.1</td><td>15.0</td><td>16.8</td><td>18.7</td><td>20.7</td><td>22.7</td></tr></tbody></table>

The variation levels obtained are very different from node to node but it goes up to almost 0.1 mm at the trailing
<figure><img src="image_6.png" /><figcaption>**Figure 6.** Radial displacement of the blade-tip nonlinear nodes with respect to the root friction coefficient $\mu^{\text{fri}}$, nodes: LE (—), 2 (—), 3 (—), 4 (—), 5 (—), 6 (—), 7 (—), 8 (—), TE (—).</figcaption></figure>

edge. To put these displacements into perspectives, the relative deviation $\Delta r_j^% = 100 \cdot \Delta r_j/c_j$ is computed and provided in Tab. 3. The results show that if the friction coefficient between the disk and the blade is not properly estimated, there can be an uncertainty on the tip gap of up to 22.7%. In the process of designing the blade and the casing, the gap consumption induced by preloading should be taken into account to ensure a safety margin on the blade-tip gap. Indeed, several mechanisms are likely to cause evolutions of $\mu^{\text{fri}}$, such as lubricant, wear and temperature. Actually, the dependence of the tip gap with the root friction coefficient is a first coupling mechanism between the two interfaces that designers should be aware of. In this work, since the casing is numerically set at a distance of $c_j = 4 \cdot 10^{-4}$ m of the pre-loaded state of the blade-tip independently of the value of $\mu^{\text{fri}}$, this does not come into play in the calculations presented afterwards.

# 5 Independent analysis of the contact interfaces

## 5.1 Configuration 1: blade-tip/casing contacts only

Before going into the coupling of friction and blade-tip contacts, an analysis of the RSI kinematics is made around the first bending mode of the assembly. The parameters of HBM calculations are set to $N_h = 10$ harmonics, $N_t^{R/S} = 500$ instants and $\varepsilon^{R/S} = 1 \cdot 10^7 \text{ N} \cdot \text{m}^{-1}$. The choice of this set of parameters results from a parametric study on the convergence of the solutions. The blade-tip friction coefficient is set to $\mu^{R/S} = 0.15$ and remains constant for all calculations made in this article. Forced responses with various excitation amplitudes (from $F = 0.25$ N to $F = 1.5$ N) and the nonlinear mode associated with the first bending mode are computed and displayed in Fig. 7. In order to be consistent with further simulations, the disk is accounted for by linearly assembling the disk and the blade at the friction interface. This means that the system without friction has the same linear eigenfrequencies than when the blade root is in a sticking state, while not being able to enter into sliding. For all NFRC presented afterwards, the amplitude of the time signal of a circumferential DoF, i.e. associated with bending, of a linear node on the blade-tip is represented.

Firstly, the results displayed in Fig. 7 show that the numerical strategy used to handle RSI is able to compute the dynamics response of the problem when no friction occurs at the root. The nonlinear mode of the system describes correctly the evolution of the resonance frequency with respect to the response amplitude: a strong stiffening is observed (shift towards higher frequencies). In order to simplify the analysis and since the nonlinear modal analysis corroborates the results obtained in forced response, only the solutions belonging to the nonlinear mode are analyzed hereafter. Note that since the DoF used for the representation of the NFRC is in the circumferential direction and
<figure><img src="image_7.png" /><figcaption>**Figure 7.** Forced response and nonlinear mode of the system without root friction, NFRC (—), linear response (……), nonlinear mode (■), five points along the nonlinear mode (○) associated with frames that represent the state of contact of the blade-tip nonlinear nodes. Each pie chart stands for a contact node, displaying the state of contact: separation (■), contact (○), pie charts represent a whole period and *t* = 0 corresponds to a trigonometric angle of zero.</figcaption></figure>

linear, no links can be made with the initial gap in Fig. 7 and further NFRC graphs.

The blade-tip/casing contact kinematics of the presented configuration can be read in the different frames displayed in Fig. 7. Each frame is related to a solution of the nonlinear mode, *i.e.* a modal shape at a given amplitude, where the 9 blade-tip nonlinear nodes motions are displayed in pie charts, from the leading edge on the left to the trailing edge on the right. The size of each pie chart is proportional to the maximum radial contact force that the node undergoes. The different phases of contact are displayed in different colors: separation in gray (■) and contact in red (○). The circle represents the whole temporal period and the time *t* = 0 corresponds to a trigonometric angle of 0° on the circle, meaning that the phase of the different contact states can be interpreted from the pie charts.

From the five frames displayed in Fig. 7, it can be seen that (1) at low amplitudes no contact occurs, (2) the trailing edge is the first node to enter in contact with the casing, (3) two nodes around the trailing edge are contacting the casing while the leading edge enters in contact in phase opposition, (4) the contact time increases for all nodes and the maximum load shifts from the trailing edge to the leading edge and finally (5) all nodes enter in intermittent contact for high amplitudes. The contact kinematics of the blade under blade-tip contact then appears to be an alternative tilting motion going from trailing edge to the leading edge, *i.e.* from back to front, in addition to the overall bending motion. The observation of the presented dynamics will be used to put in perspective the results obtained when both nonlinearities are triggered in Sec. 6. As a way of illustrating the tilting motion, the radial displacements *r*<sub>j</sub>(*t*) and corresponding radial contact forces *f*<sub>nl,j</sub><sup>R/S,N</sup> are displayed in Fig. 8 for the frame number 4 on the nodes 2, 8 and 9 (LE). The contact kinematics interpreted from frame 4 of Fig. 7 can be seen again in the signals of Fig. 8. The contact occurs in phase opposition between the front and the back of the blade.
<figure><img src="image_8.png" /><figcaption>**Figure 8.** Time signals of radial displacements and radial contact forces for the solution of frame 4 in Fig. 7 for nodes 2, 8 and 9 (LE), radial displacements (—), casing (■), radial contact forces (—).</figcaption></figure>

**5.2 Configuration 2: blade/disk friction only**

After analyzing the dynamics of the system under blade-tip/casing contacts, the dynamics of the friction interface
is analyzed in this section in order to provide a point of comparison for the studies where both contact interfaces
are accounted for. Similarly as in Sec. 5.1, several forced responses are computed around the first bending mode
for a single friction coefficient $\mu^{\text{fri}} = 0.45$. The underlying nonlinear mode is also analyzed in order to understand
the system's dynamics response. All the results are displayed in Fig. 9. The HBM parameters are set to $N_h = 10$
harmonics, $N_t^{\text{fri}} = 100$ and $\varepsilon^{\text{fri}} = 3.3 \cdot 10^7 \text{ N}\cdot\text{m}^{-1}$. Forced responses are computed with an excitation level starting
from $F = 2.5 \text{ N}$ to $F = 25 \text{ N}$ by steps of $2.5 \text{ N}$ with an additional response at $F = 0.25 \text{ N}$ (close to linear, sticking
case). The results underline a good agreement between the resonance frequency and the nonlinear mode, acting as a
backbone curve. As the excitation amplitude increases, the friction interface dissipates more energy and the apparent
damping increases, thus limiting the response amplitude. A softening (shift towards lower frequencies) behavior is
observed on the nonlinear response. This phenomenon is commonly observed on friction damping applications [47].
In order to analyze the dynamics of the system in depth, the resonances of three forced response curves are studied.
Similarly as in Fig. 5, the dynamic state of contact of the blade root friction interface is analyzed in Fig. 10c
by displaying all the nodes that compose it. For each node, three concentric pie charts are displayed that each
correspond to one of the three points labeled in Fig. 9: the inner ring corresponds to $F = 0.25 \text{ N}$, the middle ring
<figure><img src="image_9.png" /><figcaption>Figure 9. Forced response and nonlinear mode of the system with blade/disk friction only, NFRC (—), nonlinear mode (—), resonance points: F = 0.25 N (▲), F = 5 N (□), F = 25 N (●).</figcaption></figure>

to F = 5 N and the outer ring to F = 25 N. The size of the pie charts are proportional to the nodal static normal force that is applied to the contact element due to centrifugal forces (evaluated during the static computation). The pie charts display the sticking rate in light blue (●), the sliding rate in green (■), the separation rate in dark blue (●). In order to highlight the proportion of the different contact status, no information about the phase is provided (unlike RSI pie charts). The summed up data is also available in Tab. 4 where the sliding rate is computed such as

$$S_r^{\text{sli}} = 100 \cdot \frac{\sum_{i=1}^{N_t^{\text{fri}}} \sum_{j=1}^{n^{\text{fri}}} \delta_{B^{i,j}}^{A^{i,j}}}{N_t^{\text{fri}} n^{\text{fri}}} \quad (25)$$

with δ the Kronecker symbol, $A^{i,j} = \mu^{\text{fri}} |f_{nl,j}^{\text{fri,N}}(t_i)|$ and $B^{i,j} = \|f_{nl,j}^{\text{fri,T}}(t_i)\|$. The separation rate is computed following

$$S_r^{\text{sep}} = 100 \cdot \frac{\sum_{i=1}^{N_t^{\text{fri}}} \sum_{j=1}^{n^{\text{fri}}} \delta_0^{A^{i,j}}}{N_t^{\text{fri}} n^{\text{fri}}}. \quad (26)$$

The sticking rate $S_r^{\text{sti}}$ is then computed such that $S_r^{\text{sti}} = 100 - S_r^{\text{sli}} - S_r^{\text{sep}}$.

<table><caption>Table 4. Overall dynamic state of contact at the three resonances for blade/disk friction only.</caption><thead><tr><th></th><th>sticking</th><th>sliding</th><th>separation</th></tr></thead><tbody><tr><td>F = 0.25 N</td><td>96.3 %</td><td>3.7 %</td><td>0 %</td></tr><tr><td>F = 5 N</td><td>77.7 %</td><td>22.1 %</td><td>0.2 %</td></tr><tr><td>F = 25 N</td><td>59.2 %</td><td>40.6 %</td><td>0.2 %</td></tr></tbody></table>
<figure><img src="image_10.png" /><figcaption>**Figure 10.** Dynamic state of contact along the nonlinear mode with $\mu^{\text{fri}} = 0.45$ without RSI. Dynamic state of contact (c): sticking (sticking symbol), sliding (sliding symbol), separation (separation symbol), inner ring: solution 1, middle ring: solution 2, outer ring: solution 3; contact force (a) (b) (d) (e): $\mu^{\text{fri}}|f_{\text{nl},j}^{\text{fri,N}}(t)|$ (left), $|f_{\text{nl},j}^{\text{fri,T}}(t)|$ (right).</figcaption></figure>

The results presented in Fig. 10 show the evolution of the state of friction on the three resonances labeled in Fig. 9. In Fig. 10c, it is shown that going from low amplitudes to high amplitudes, i.e. from the inner to the outer ring, sliding increases overall. This result is expected and constitutes a reference result for the rest of the study. It is worth noting the high degree of asymmetry of the sliding pattern between pressure and suction sides. This can be attributed to the complex geometry of the blade that favors asymmetric responses. Figures 10a, 10b, 10d and 10e illustrate the time-period evolution of the friction forces $\|f_{\text{nl},j}^{\text{fri,T}}(t)\|$ inside the Coulomb cone $\mu^{\text{fri}}|f_{\text{nl},j}^{\text{fri,N}}(t)|$ for four nodes identified in Fig. 10c. The contact force signals correspond to the outer ring, thus the resonance for $F = 25 \text{ N}$ (●) of Fig. 9. When $\|f_{\text{nl},j}^{\text{fri,T}}(t)\| = \mu^{\text{fri}}|f_{\text{nl},j}^{\text{fri,N}}(t)|$, the node is sliding and the background is highlighted in green (■). When $\|f_{\text{nl},j}^{\text{fri,T}}(t)\| < \mu^{\text{fri}}|f_{\text{nl},j}^{\text{fri,N}}(t)|$, the node is sticking and the background is highlighted in light blue (■). Finally when $\mu^{\text{fri}}|f_{\text{nl},j}^{\text{fri,N}}(t)| = 0$, the node is separated and the background is highlighted in dark blue (■). The rate at which the three states occur is what is reflected on the pie charts graphs of Fig. 10c. In Fig. 10a, only a sticking behavior is reached. The bottom left node associated with Fig. 10b, however, undergoes all three different states during one period: sliding, sticking and separation. The size of its associated pie chart shows that this node has a very low initial normal force which favors such changes in contact status. For Fig. 10d and Fig. 10e, an intermittent sticking and sliding motion is observed on both nodes.
# 6 Configuration 3: coupled analysis of blade/disk friction and blade-tip/casing contact

Now that the contact kinematics both at the blade root and the blade tip have been analyzed independently in Sec. 5, the coupled response—with both nonlinearities activated at once—is investigated in depth. As a first step, the dynamics of the system is studied for a single value of $\mu^{\text{fri}}$. A local analysis at the nodes of the blade/disk interface is then proposed.

## 6.1 Analysis of forced responses

First of all, for illustration, one forced response is computed for an excitation amplitude of $F = 1.5$ N for each of the three configurations. The corresponding NFRC are displayed in Fig. 11 and include: a response with RSI only (configuration 1), with blade/disk friction only (configuration 2) and a response with both interfaces activated (configuration 3). The results highlight the fact that both forced responses that include blade/disk friction (...) (—) configuration 2 (with blade/disk friction only) (...) configuration 3 (with RSI and blade/disk friction) (—).

<figure><img src="image_11.png" /><figcaption>**Figure 11.** NFRC for $F = 1.5$ N: configuration 1 (with RSI only) (—) configuration 2 (with blade/disk friction only) (...) configuration 3 (with RSI and blade/disk friction) (—).</figcaption></figure>

can be compared since they come with the same level of damping: linear and nonlinear due to blade/disk friction. Indeed, when RSI is not yet triggered (low amplitudes), the responses are indissociable. However, a difficulty arises when looking at the RSI only response: since blade/disk friction is not accounted for, the apparent level of damping is significantly lower and thus contact interactions at the tip of the blade are notably more intense (amplitudes three times higher) for the same level of forcing $F$. This observation highlights the damping capabilities of the blade/disk interface, but prevents the direct comparison of the three forced responses at once. This is why, in the rest of this article, the comparisons between configurations 1 and 3 are only conducted on the nonlinear modes. This way, the mode shapes can be compared at a given response amplitude, independently of $F$.

In order to be consistent with the calculations made in Sec. 5.2, forced responses are computed for $\mu^{\text{fri}} = 0.45$ with forcing amplitudes going from $F = 2.5$ N to $F = 25$ N by steps of 2.5 N and with a low amplitude of $F = 0.25$ N. The associated nonlinear mode is also computed. The responses of configurations 2 and 3 are depicted in Fig. 12. Moreover, the nonlinear mode of the configuration 1 is displayed for comparison. First of all, it is observed, as expected, that the blade's response in configuration 3 is identical to configuration 2 as long as the response amplitude is not high enough to trigger the activation of RSI. This can be observed both on the forced responses and on the nonlinear modes that coincide up to an amplitude of around $3 \cdot 10^{-3}$ m. However, when the amplitudes reach sufficiently high values ($\ge 3 \cdot 10^{-3}$ m) the RSI is triggered and a stiffening effect appears at first. Unlike what is commonly observed on RSI forced responses, in Fig. 7 or in the literature [13, 15, 48], the system undergoes successively a stiffening then a softening effect. It seems that the nonlinear mode of configuration 3 (—) is dominated
<figure><img src="image_12.png" /><figcaption>**Figure 12.** Forced responses and nonlinear modes of the system with blade/disk friction and RSI with $\mu^{\text{fri}} = 0.45$, configuration 3: NFRC (—), nonlinear mode (—), resonances: $F = 0.25$ N (▲), $F = 5$ N (□), $F = 25$ N (⊙), points along the nonlinear mode (⊙); configuration 2: NFRC (····), nonlinear mode (---), resonances: $F = 0.25$ N (▲), $F = 5$ N (□), $F = 25$ N (⊙), points along the nonlinear mode (⊙); configuration 1: nonlinear mode (—), points along the nonlinear mode (⊙); frames represent the state of contact of the blade-tip nonlinear nodes: separation (), contact (■), pie charts represent a whole period and $t = 0$ corresponds to a trigonometric angle of zero.</figcaption></figure>

by RSI effects at first since it is tangent to the nonlinear mode of configuration 1 (—), explaining the stiffening effect. Then, around a response amplitude of $6.5 \cdot 10^{-3}$ m, a softening behavior appears and prevails over the stiffening effect. This fact is attributed to the strong activation of blade/disk friction that makes the coupled nonlinear mode (—) (and its nonlinear resonance frequency) tend towards the blade/disk friction only (■■) mode. This observation on the nonlinear modes is also confirmed by looking at the shapes of the forced response curves that effectively reflects this fact: firstly the resonance frequency is shifted towards high frequencies and then goes back to lower frequencies at high forcing amplitudes. Moreover, on the highest forcing amplitudes displayed in Fig. 12, the shapes of the curves do not reflect the responses usually observed for this kind of interactions. Indeed, the blade/disk friction interface seems to introduce strong damping in the system, to the point where friction appears to dominate the overall nonlinear response and the NFRC look flattened.

## 6.2 Analysis of the dynamics at the contact interfaces

In order to analyze the response of the system in greater detail, this section describes how the contact kinematics are affected at each interface when accounting for both nonlinearities at once. The study is firstly done on the RSI interface and then secondly on the blade/disk interface. An extensive analysis is led on the influence of RSI contacts on the contact-related quantities at the root of the blade.

### 6.2.1 Blade-tip/casing contact interface

When looking at the blade-tip/casing contact kinematics in frames 3 and 4 in Fig. 12, the dynamics of the blade for solutions with RSI only (green background) and with both interfaces accounted for (white background) seem to be almost identical. The 2D trajectories of each contact node in the circumferential-axial plane (seen from above) are displayed in Fig. 13 for both configurations. The displayed solutions in Fig. 13a and Fig. 13b correspond to the nonlinear modes' solutions of frame 3 and frame 4 respectively. The trajectories of the nodes when RSI is not activated are also displayed for reference.
<figure><img src="image_13.png" /><figcaption>**Figure 13.** 2D trajectories of the RSI contact nodes, configuration 1 (RSI only) (—), configuration 2 (friction only) (—), configuration 3 (RSI + friction) (—), contact nodes at rest (○).</figcaption></figure>

At low amplitudes (for frame 3 in Fig. 13a), trajectories are quite similar and the blade experiences the same kinematics in both configurations 1 (—) and 3 (—). However, for higher amplitudes (frame 4 in Fig. 13b), the trajectories are notably different and the tilting behavior associated with the RSI only configuration (—) near the trailing edge is attenuated. Indeed, while the blade-tip experiences a complex 2D trajectory in the configuration with RSI only, the motion becomes almost unidirectional when blade/disk friction is introduced. Friction interactions seem to simplify the kinematics of the blade to a more classical bending motion going back and forth from the pressure side to the suction side, as for configuration 2 (—). Thus, while being relatively remote from the blade-tip, the introduction of the blade/disk friction interface is responsible for a significant change in kinematics of the blade-tip during contact interactions.

### 6.2.2 Blade/disk friction interface

As a way of investigating the behavior observed on forced responses in Sec. 6.1, the dynamics of the blade/disk friction interface are analyzed in depth in this section. Firstly, the evolution of the dynamic state of contact of the blade/disk friction interface for the coupled configuration 3 is given for the two levels of excitations at the resonance: 5 N and 25 N. As a way of visualizing the comparison between the two dynamic states, the two configurations are displayed on the graph in Fig. 14a for $F = 5$ N and in Fig. 14b for $F = 25$ N. In both Fig. 14a and Fig. 14b, the inner ring corresponds to configuration 2 (without RSI) while the outer ring stands for configuration 3. This way, the comparison between configurations 2 and 3 can be easily read by comparing the two concentric pie charts on each node.

For both resonances displayed in Fig. 14, even if small discrepancies can be observed on several nodes, the overall contact kinematics are similar. This means that the RSI interface does not fundamentally modify the blade/disk friction kinematics even for high level of forcing. This property can be more easily seen in Tab. 5 that summarizes
<figure><img src="image_14.png" /><figcaption>**Figure 14.** Comparison of both dynamic states, inner ring: configuration 2, outer ring: configuration 3, sticking (●), sliding (■), separation (◆).</figcaption></figure>

the different states of friction for the two configurations. The data presented for the configuration without RSI is actually the data from Tab. 4, recalled here for the sake of comparison.

<table><caption>**Table 5.** Overall dynamic state of contact at the three resonances with both blade/disk friction and RSI.</caption><thead><tr><th rowspan="2"></th><th colspan="2">sticking</th><th colspan="2">sliding</th><th colspan="2">separation</th></tr><tr><th>RSI</th><th>no RSI</th><th>RSI</th><th>no RSI</th><th>RSI</th><th>no RSI</th></tr></thead><tbody><tr><td>F = 0.25 N</td><td>96.3%</td><td>96.3%</td><td>3.7%</td><td>3.7%</td><td>0%</td><td>0%</td></tr><tr><td>F = 5 N</td><td>79.7%</td><td>77.7%</td><td>20.1%</td><td>22.1%</td><td>0.2%</td><td>0.2%</td></tr><tr><td>F = 25 N</td><td>58.9%</td><td>59.2%</td><td>40.7%</td><td>40.6%</td><td>0.4%</td><td>0.2%</td></tr></tbody></table>

For both configurations, the kinematics are the same at $F = 0.25$ N since RSI is not yet activated. Then, sliding is slightly decreased (–2%) for $F = 5$ N when RSI is activated. Finally, the qualitative observation made on Fig. 14b is confirmed by the overall sliding rate that is almost identical for both cases when $F = 25$ N (0.1% difference).

Sliding rate is correlated to nonlinear damping, however it is not directly proportional and other parameters should be studied. Thus, further investigations are led on the kinematics of the blade/disk friction interface. A study on the intensity and direction of friction forces is carried out. Firstly, the normal contact force is analyzed. The total normal forces at the blade/disk interface are plotted with respect to $\omega$ for responses with $F = 5$ N and $F = 25$ N respectively in Fig. 15a and Fig. 15b.

Two distinct cases are observed for the forced responses in Fig. 15: (1) when the RSI is weakly activated (Fig. 15a) the normal force of configuration 3 is higher than configuration 2 and (2) when the RSI is strongly activated (Fig. 15b)
<figure><img src="image_15.png" /><figcaption>Figure 15. Normal force at the blade/disk friction interface on forced responses, mean normal force: with RSI (—), without RSI (—), alternating normal force: with RSI (■), without RSI (■).</figcaption></figure>

the normal force starts to decrease and becomes similar to the case with blade/disk friction only (configuration 2). Consequently, at low forcing amplitudes, the sliding limit is higher when RSI is activated and therefore the magnitude of friction forces is higher when sliding occurs.

Finally, the direction of friction is examined in Fig. 16. As a way of evaluating how the direction of friction is influenced by the RSI kinematics, the evolution of the friction forces inside the friction cone is displayed. For each node, the green circle (○) represents the sliding limit. The inside of the circle, displayed in light blue (■), thus stands for the sticking area. The tangent forces are normalized by the sliding limit $\mu^{\text{fri}}|f_{\text{nl},j}^{\text{fri,N}}|$ in order to allow for schematic representation. These two-directional forces relative to the sliding limit are then plotted inside these circles in blue lines (—). By doing so, the direction of friction and the state of contact can be read along one period of vibration in Fig. 16 for each node. The plotted quantities are

$$ \bar{f}_{\text{nl},j}^{\text{fri,T1}} = \frac{f_{\text{nl},j}^{\text{fri,T1}}}{\mu^{\text{fri}}|f_{\text{nl},j}^{\text{fri,N}}|} \quad (27a) $$

$$ \bar{f}_{\text{nl},j}^{\text{fri,T2}} = \frac{f_{\text{nl},j}^{\text{fri,T2}}}{\mu^{\text{fri}}|f_{\text{nl},j}^{\text{fri,N}}|} \quad (27b) $$

$$ \text{with } \mathbf{f}_{\text{nl},j}^{\text{fri,T}} = \begin{bmatrix} f_{\text{nl},j}^{\text{fri,T1}} \\ f_{\text{nl},j}^{\text{fri,T2}} \end{bmatrix}. $$

Firstly, it is observed on Fig. 16a for configuration 2 (without RSI) that the tangent friction forces oscillate from sliding back and forth between the pressure side and the suction side (T1 direction) on the majority of nonlinear nodes. Moreover, the axial component (from LE to TE, T2 direction on the figures) of friction forces do not seem to contribute as much as the circumferential direction (except on nodes with low static normal force, located near edges, where 2D orbits are witnessed, refer to Fig. 5b). This motion is characteristic of the first bending mode and is thus expected in this configuration. When RSI is activated (Fig. 16b), the motion is very similar and the influence of the LE-TE tilting motion, induced by RSI, does not seem to have a significant influence on the direction of friction. Indeed, the friction forces still oscillate alternatively from the pressure side to the suction side with few contribution of the T2 direction. It comes that the dynamics of the blade seems to be dominated by the bending motion even when the RSI interface is accounted for.
<figure><img src="image_16.png" /><figcaption>Figure 16. Friction directions for F = 25 N, sliding limit $\mu^{fri} |f_{nl,j}^{fri,N}|$ (—), sticking (■), 2D-tangent force relative to sliding limit (—).</figcaption></figure>

## 6.3 Energy analysis

As a way of going further in the analysis of the blade/disk contact, the dissipated energy at the interface is investigated. The dissipated energy of each jth node over one period is computed by means of

$$E_{\text{dis},j}^{\text{fri}} = \int_{0}^{T} \langle \dot{\mathbf{x}}_{\text{rel},j}^{\text{fri}}, \mathbf{f}_{\text{nl},j}^{\text{fri}} \rangle \mathrm{d}t \quad (28)$$

and displayed in Fig. 17 for both configurations 2 and 3 at the resonances for F = 5 N and F = 25 N. For the case F = 5 N in Fig. 17a and Fig. 17b, while Tab. 5 highlighted a difference in terms of sliding rate, the total energy dissipated at the friction interface is equivalent for both configurations. Though, the dissipated energy distribution reveals that the energy is not dissipated on the same nodes from one configuration to the other. On the one hand, for configuration 2, the energy dissipation is concentrated on two areas: at the center of the suction side and near the trailing edge at the pressure side. On the other hand, for configuration 3, a new area of dissipation appears near the leading edge on the pressure side and the two areas already identified for configuration 2 dissipate less energy (darker colors). This demonstrates that even if the levels of dissipated energy over a period are equivalent, the distribution of dissipated energy at the interface is still affected when considering blade-tip/casing contacts. In accordance with Fig. 14, the loci of energy dissipation are concentrated on the nodes experience significant sliding phase while being highly statically loaded.

On the contrary, the observations made for F = 5 N are not applicable to the case F = 25 N in Fig. 17c and Fig. 17d, the total energy dissipated over one period decreases of around 14% when accounting for the RSI. However, in order to compare the behavior of both configurations, it is necessary to look at relative quantities. This is done
<figure><img src="image_17.png" /><figcaption>Figure 17. Energy dissipation at the blade/disk interface at different resonances ((a), (b) F = 5 N, (c), (d) F = 25 N) for configurations 2 and 3 (with and without RSI). "s" stands for suction side and "p" for pressure side.</figcaption></figure>

through the analysis of the percentage of the vibrational mechanical energy that is dissipated over one period by the blade/disk interface $E_{dis}^{fri,\%}$. It is computed following

$$E_{dis}^{fri,\%} = 100 \cdot \frac{\sum_{j=1}^{n_{dis}^{fri}} E_{dis,j}^{fri}}{E_{m}}. \qquad (29)$$

with the vibrational mechanical energy being

$$E_{m} = \frac{1}{2} \mathbf{x}^{\top} \mathbf{K} \mathbf{x} + \frac{1}{2} \dot{\mathbf{x}}^{\top} \mathbf{M} \dot{\mathbf{x}} - \frac{1}{2} \mathbf{x}_{0}^{\top} \mathbf{K} \mathbf{x}_{0} \qquad (30)$$

where the strain energy stemming from the statically preloaded displacement field $\mathbf{x}_0$ is subtracted.

In configuration 3, the vibrational mechanical energy of the blade is lower because of the introduction of RSI that causes the blade to experience lower vibrational amplitudes than configuration 2 for a fixed level of excitation F (refer to Fig. 12). Then, even though the absolute dissipation in joules is lower for configuration 3 (0.782 J for configuration 2 and 0.670 J for configuration 3), the values of $E_{dis}^{fri,\%}$ are almost equivalent for both configurations. Since both percentages of total energy dissipation are similar (39.1% for configuration 2 and 38.5% for configuration 3), the efficiency of the energy dissipation at the blade root is still equivalent between both configurations. Moreover, both energy dissipation distributions are similar and only two main areas are subject to energy dissipation. This shows that there is a complex interaction between the kinematics at both interfaces that induces distinct behaviors according to the level of energy in the system, directly related to the intensity of excitation forces.

As a way of going further, $E_{dis}^{fri,\%}$ is computed along the nonlinear modes corresponding to configurations 2 and 3 in Fig. 18. In addition, the same quantity is computed at the resonance of each response displayed in Fig. 12 and superimposed with the data associated with the nonlinear modes.

The results displayed in Fig. 18 highlight the fact that for a given mechanical energy (vertical line on the graph), $E_{dis}^{fri,\%}$ is higher for configuration 3 (with RSI) than for configuration 2 along the nonlinear modes. Superimposing the points associated with the different resonances in Fig. 18 demonstrates that this property is also true when considering forced responses. Thus, it shows that Fig. 18 can not be read at constant mechanical energy: for a fixed value of the forcing amplitude F both configurations do not respond with the same level of mechanical energy.
<figure><img src="image_18.png" /><figcaption>**Figure 18.** Percentage of the vibrational mechanical energy dissipated by the blade/disk interface, configuration 2: nonlinear mode (—) $F = 0.25$ N (o), $F = 5$ N (□), $F = 25$ N (●), other values of $F$ (▲); configuration 3: nonlinear mode (—) $F = 0.25$ N (o), $F = 5$ N (□), $F = 25$ N (●), other values of $F$ (▲), same value of $F$ (···).</figcaption></figure>

Hence, it is more appropriate to compare the responses for similar values of $F$. In the end, two effects going in opposite directions are identified: configuration 2 dissipates more energy but vibrational orbits are more energetic. These two effects seem to almost perfectly cancel out for this model and set of data: for every value of $F$, the ratio of the dissipated energy along one period over the vibrational mechanical energy is similar for both configurations (points with the same value of $F$ are linked by a black dotted line (···)). This underlines that the blade/disk friction interface dissipates the same relative amount of energy in both configurations, and the introduction of RSI neither enhances nor degrades the efficiency of the dissipation taking place of the blade/disk interface.

Moreover, the nonlinear damping ratios obtained along the complex nonlinear modes of the system with and without RSI are plotted in Fig. 19. The trend of evolution of both damping ratios is exactly following the trend

<figure><img src="image_19.png" /><figcaption>**Figure 19.** Nonlinear damping ratio of the nonlinear modes in configuration 3 (—), configuration 2 (—), configuration 1 (—), linear damping $\xi_{1B} = 1 \cdot 10^{-3}$ (---).</figcaption></figure>

identified in Fig. 18 on the percentage of vibrational mechanical energy dissipated at the blade/disk interface. This
shows that the damping of the system is almost entirely driven by the dynamics of the friction interface. This also
shows that the intermittent friction occurring at the blade-tip does not contribute significantly to the damping of the system when considering a blade rotating at a constant rotation speed. This is also reflected by the fact that the nonlinear damping ratio of the mode with RSI only stays close to the linear damping value of $\xi_{1B} = 1 \cdot 10^{-3}$.

# 7 Evolution of contact characteristics with respect to the friction coefficient and excitation level

After analyzing in depth the dynamics of the blade for a given set of parameters, a parametric study is conducted as a way of providing a broader view of the system's behavior. In this section, both the value of the blade attachment friction coefficient $\mu^{\text{fri}}$ and the level of excitation $F$ on the system are varied as a way of determining how the global response of the system evolves with respect to these first order parameters. Since RSI may have a very detrimental impact on the system's dynamics response [27], the determination of design areas where the RSI is weakly activated (if at all activated) is capital. This study aims to provide guidelines regarding the design of such systems when accounting for both blade/disk friction and RSI. At first, a subset of the responses computed for this study are presented in Sec. 7.1, then the contact characteristics of both interfaces are assessed in Sec. 7.2 and Sec. 7.3 respectively on the whole range of parameters.

The range of forcing amplitude [0.125 N, 5 N] is swept with a step of 0.125 N alongside the blade/disk friction coefficient $\mu^{\text{fri}}$ on the range [0.15, 0.75] by steps of 0.025. Both configurations 2 and 3 (blade/disk friction with and without RSI) are treated for comparison purposes. It yields a matrix of 40 amplitudes, 25 friction coefficients and 2 configurations for a total of 2000 forced responses.

## 7.1 Forced responses

Before analyzing in details the contact characteristics of the interfaces, some of the forced responses used for this study are presented in Fig. 20. It displays all 40 forced responses (all excitation levels) for four friction coefficients $\mu^{\text{fri}} \in \{0.15, 0.30, 0.60, 0.75\}$ at the blade attachment.

All the forced responses displayed in Fig. 20 reveal a first trend: the lower is the blade attachment friction coefficient for $\mu^{\text{fri}}$, the higher the forcing amplitude $F$ should be to activate RSI contacts. For instance, for $\mu^{\text{fri}} = 0.15$ in Fig. 20a, RSI is hardly activated for the selected range of $F$. Even when it is activated at high forcing amplitude $F \approx 4$ N, the response amplitudes for configuration 3 remain of the same order of magnitude than for configuration 2. Then, as $\mu^{\text{fri}}$ increases, the resonance frequencies of responses for configuration 3 are drifting further and further away. This effect is closely related to the intensity of RSI contacts: it is then expected that RSI contacts are stronger if the $\mu^{\text{fri}}$ is high. This coupling mechanism is investigated in the following sections.

## 7.2 Analysis of RSI-related quantities

Firstly, the RSI characteristics are evaluated on the configuration that includes both RSI and blade/disk friction. The number of contacting nodes at the resonance of the forced responses is displayed in Fig. 21 as a function of $\mu^{\text{fri}}$ and $F$.

The figure shows that for low values of the friction coefficient and for low levels of forcing amplitude, RSI is not even triggered which is consistent with what is expected. Indeed, low friction coefficients promote the occurrence of sliding phases during the vibration and thus dissipate more energy, bringing more damping to the system. This additional damping prevents any node of the blade-tip from entering into contact. Though, the value of $F$ that is necessary for the blade to enter blade-tip/casing contacts decreases fast with respect to $\mu^{\text{fri}}$. For instance, after reaching $\mu^{\text{fri}} = 0.3$, the value of $F$ that causes the triggering of RSI is below $F = 1$ N. Then, the number of contacting nodes at the resonance increases when $\mu^{\text{fri}}$ and $F$ increase, which is also consistent. Moreover, Fig. 21 points out the fact that for a known and constant level of forcing $F$ of up to around 3.5 N, if an uncertainty on the friction coefficient of the blade/disk interface exists on the range [0.15, 0.75], then the blade could experience no contact for low values of $\mu^{\text{fri}}$, or contacts from 1 up to 4 nodes. It means that all intensities of RSI can be encountered if the $\mu^{\text{fri}}$ is not well known and accounted for. As already pointed out in Sec. 4.2, a variability of the friction coefficient caused by thermal or wear effects could radically change the nonlinear response of the blade and may bring up RSI in extreme cases.
<figure><img src="image_20.png" /><figcaption>Figure 20. Forced responses of $\mu^{\text{fri}}$, configuration 3: NFRC (—) configuration 2: NFRC (—).</figcaption></figure>

In parallel with Fig. 21, the mapping of the maximum resultant radial contact force withstood by the blade is displayed in Fig. 22. It allows to evaluate the intensity of the blade-tip/casing contacts. The mapped quantity $f_{\text{nl,max}}^{\text{R/S,N}}$ corresponds to the maximum over time of the sum of the radial contact forces on all nodes at the resonance, computed following

$$f_{\text{nl,max}}^{\text{R/S,N}} = \max_i \left( \sum_{j=1}^{n^{\text{R/S}}} f_{\text{nl},j}^{\text{R/S,N}}(t_i) \right). \qquad (31)$$

The evolution of $f_{\text{nl,max}}^{\text{R/S,N}}$ with respect to $F$ and $\mu^{\text{fri}}$ is fully in accordance with what was observed in Fig. 21: the intensity of blade-tip/casing contacts increases with both $F$ and $\mu^{\text{fri}}$. Furthermore, at a fixed value of $F$ the maximum resultant force goes from 0 (or close to 0) up to values of several hundreds of newtons for $f_{\text{nl,max}}^{\text{R/S,N}}$, depending
<figure><img src="image_21.png" /><figcaption>Figure 21. Number of contacting nodes at the blade-tip with respect to $\mu^{\text{fri}}$ and the excitation amplitude $F$ at the resonance of the forced responses.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure 22. Maximum resulting radial contact force applied on the blade with respect to $\mu^{\text{fri}}$ and the excitation amplitude $F$.</figcaption></figure>

on the value of $\mu^{\text{fri}}$. This observation underlines the importance of determining with precision the friction coefficient at blade attachments in design phases.

## 7.3 Blade/disk friction related characteristics

The blade/disk friction characteristics are then evaluated on the whole mapping and compared for both configurations. As evidenced in Sec. 6.3, the relevant quantity allowing to evaluate the damping characteristics of the blade/disk interface in forced responses is the percentage of the vibrational mechanical energy dissipated over one period at the resonance. The energy dissipation at the blade/disk interface is, among others, directly related to the sliding rate
and the normal force at the interface. In order to understand how the dynamics of the blade evolve across the range of parameters chosen, the difference between both configurations 2 and 3 in terms of sliding rate of the interface and the normal force variation are studied at first. As a way of displaying the discrepancies between the two more clearly, the difference between the sliding rates of configuration 2 (reference) $S_r^{\text{sli,ref}}$ and configuration 3 $S_r^{\text{sli,R/S}}$, noted $\Delta S_r^{\text{sli}} = S_r^{\text{sli,R/S}} - S_r^{\text{sli,ref}}$, is mapped in Fig. 23.

<figure><img src="image_23.png" /><figcaption>**Figure 23.** Difference in sliding rate $\Delta S_r^{\text{sli}}$ of the nodes belonging to the blade/disk friction interface at the resonances for both configurations. Positive values (■) refer to the case where configuration 3 (with RSI) experiences more sliding than configuration 2 (without RSI), and negative values (■) is the other way around.</figcaption></figure>

The mapping of Fig. 23 reveals a complex evolution of the sliding rate between both configurations. Four areas are uncovered. (1) An area where RSI is not triggered and thus the responses are identical and there is no difference between the two (□). (2) For low forcing amplitudes (RSI is weakly triggered), configuration 2 experiences a slight higher sliding rate than configuration 3. (3) Then, configuration 3 experiences more sliding (■) than configuration 3 in an area that corresponds partly to the area where only one contacting point is detected along the blade tip (□) in Fig. 21. (4) Finally, in the right-top corner, corresponding to high forcing amplitudes and friction coefficients, configuration 2 enters into sliding more than configuration 3 of up to 4%: this particular area seems also related with the specific state of contact of the blade-tip/casing interface where 4 nodes are contacting the casing along the period (□) in Fig. 21. The transition from area (3) to (4) seems related to the shift from 1 (□) to 4 (■) contacting nodes in Fig. 21. It suggests that the dynamics of the blade/disk interface is influenced by the kinematics of the blade-tip/casing contacts. Indeed, strong correlations between Fig. 21 and Fig. 23 are evidenced on this application.

In order to investigate further on the dynamics of the system, the relative variation of the mean normal force with respect to the reference configuration 2, noted $\Delta N^\%$, is mapped in Fig. 24. It is computed following

$$ \Delta N^\% = 100 \cdot \frac{N^{\text{R/S}} - N^{\text{ref}}}{N^{\text{ref}}} \quad \text{with} \quad \left\{ \begin{array}{ll} N^{\text{R/S}} = \sum_{j=1}^{n_{\text{fri}}} \left| \mathbf{f}_{\text{nl},j}^{\text{fri,N}}(t) \right| & \text{for configuration 3,} \\ N^{\text{ref}} = \sum_{j=1}^{n_{\text{fri}}} \left| \mathbf{f}_{\text{nl},j}^{\text{fri,N}}(t) \right| & \text{for configuration 2.} \end{array} \right. \qquad (32) $$

Like for the sliding rate, the normal force mapping seems related to the blade-tip/casing kinematics highlighted in Fig. 21. The results highlight that for the area where there are four contacting node in Fig. 21, the normal force is higher for configuration 3, thus yielding less sliding. Actually, the evolution of $\Delta N^\%$ and $\Delta S_r$, in Fig. 24 and Fig. 23 respectively, is opposite: if the sliding rate is higher for a given configuration than the other then the normal
<figure><img src="image_24.png" /><figcaption>**Figure 24.** Normal force difference at the blade/disk interface at the resonances for both configurations. Positive values (■) refer to the case where configuration 3 (with RSI) is more loaded at the blade/disk interface than configuration 2 (without RSI), and negative values (■) is the other way around</figcaption></figure>

force is lower and vice versa. It shows that both quantities are affected by the dynamics of the blade-tip/casing contacts, as evidenced in the discussion about Fig. 23.

However, the quantity that is actually of interest for designers is the nonlinear damping ratio of the whole system. It was shown in Sec. 6.3 that the damping ratio of the system in forced response is directly driven by the proportion of the vibrational mechanical energy that is dissipated over one period by the blade/disk interface $E_{dis}^{fri,%}$. Thus, the difference between this quantity on configuration 2 $E_{dis,fri,%}$ and on configuration 3 $E_{dis,R/S,%}$, computed following $\Delta E_{dis}^{fri,\%} = E_{dis,R/S,\%} - E_{dis,fri,\%}$, is mapped in Fig. 25.

<figure><img src="image_25.png" /><figcaption>**Figure 25.** Difference in energy dissipation at the blade/disk interface at the resonances between both configurations. Positive values (■) refer to the case where configuration 3 (with RSI) dissipates more than configuration 2 (without RSI), and negative values (■) is the other way around.</figcaption></figure>
The mapping of $\Delta E_{\text{dis}}^{\text{fri,\%}$ displayed in Fig. 25 reveals that the four areas identified in Fig. 23 for the sliding rate are practically opposite in terms of energy dissipation. Indeed, while the four areas identified in Fig. 23 are also present in Fig. 25, when configuration 3 experiences more sliding than configuration 2, it also dissipated a lower proportion of energy. Hence, the direct link drawn between the blade-tip/casing kinematics and the sliding rate can also be extended to the energy dissipation percentage $E_{\text{dis}}^{\text{fri,\%}$. In order to explain this dependence between sliding rate and energy dissipation, it is necessary to analyze it in the light of the mean normal force variation of Fig. 24.

Following what was observed in Fig. 24, the mean normal force evolves in the opposite direction as the sliding rate. Then, when sliding occurs, friction forces are also higher since the sliding limit (proportional to the normal force) has increased. This shows that the variation in normal force is actually driving the energy dissipation more than the sliding rate for this case study. Indeed, the mapping related to the energy dissipation in Fig. 25 is more closely related to the normal force mapping. Though, this analysis explains only part of how the energy dissipation is evolving in Fig. 25: it is coherent for high forcing and high friction coefficients. Nevertheless, it seems that more complex behaviors occur at low forcing amplitudes that cannot be directly correlated using Fig. 24 and Fig. 23 only.

In the end, all these observations underline the close link that exists between the blade-tip/casing contacts and the local indicators of the blade/disk friction interface that are the energy dissipated, sliding rate and normal force. Overall, the sweep of both $F$ and $\mu^{\text{fri}}$ allowed to uncover specific coupling mechanisms between the interfaces. Moreover, the importance of evaluating precisely the friction coefficient of the blade/disk interface and the excitation level was highlighted. Indeed, on relatively restrained ranges for both parameters, a wide number of different behaviors can occur and threaten the blade integrity. Performing the type of calculations presented in this article in the design phase of turbomachines appears to be ever more crucial as a way to improve the safety and efficiency of aircraft engines. This section underlines that the dynamics nonlinear response of such systems can be managed by using the parameters associated with the nonlinear interfaces as new dimensions to the design space. Indeed, in the present case it is shown that controlling carefully the blade/disk friction coefficient can lead to the mitigation of blade-tip/casing contacts.

## 8 Conclusion

This paper explores the dynamics response of a state-of-the-art fan stage featuring two nonlinear interfaces: blade-tip/casing contacts and blade/disk friction. Firstly, the static state of the blade is analyzed as a function of the friction coefficient of the blade/disk interface. It reveals that a first static coupling mechanism lies in the dependence of the blade-tip displacement to the friction coefficient. The dynamics responses of the blade featuring a single of the two interfaces are analyzed in order to provide points of reference for the coupled case. These reference points are of particular importance because the dynamics reflected by the configuration that accounts for both interfaces is significantly different from each one of the individual configurations. The coupled configuration, accounting for both interfaces simultaneously, is characterized through multiple forced responses and nonlinear modes calculations. A complex evolution of the resonance frequency with respect to the response amplitude is highlighted. Analyzing the kinematics of the blade-tip contact reveals that accounting for blade/disk friction strongly modify the trajectories of the contact nodes. Then, the load of the blade/disk interface for the configuration with blade-tip/casing contact is compared to the reference case in terms of state of contact, normal force and friction direction. Opposing effects are detected on the sliding rate and the normal force of the interface, resulting in the fact that the proportion of energy dissipated by friction remains similar for both configurations for a given excitation level. Though, the location of dissipation is modified by the specific kinematics induced through the introduction of rotor/stator interaction. The nonlinear damping ratio is also shown to be directly related to the amount of energy dissipated at the blade attachments, making it the main source of damping of the system. Finally, the robustness of the strategy is demonstrated by performing a large number of forced responses on ranges of values for the friction coefficient and the excitation level. A strong correlation link is found between the kinematics at the blade-tip (number of contacting nodes over a period) and local indicators of the blade/disk interface (energy dissipated, normal force, sliding rate). Specific motions of the blade-tip induce more or less energy dissipation at the root depending on the number of contact nodes that enter into contact.

Among the perspectives of this work, the characterization of a bladed disk in cyclic symmetry would allow to make more representative simulations, but a theoretical lock still remains concerning the treatment of highly
nonsmooth nonlinear interfaces such as the blade-tip/casing interface in cyclic symmetry. Moreover, using a real aerodynamic dynamic loading on the blade could enhance the fidelity of the results. Furthermore, the model of the casing could be enriched to account for its flexibility and the characteristic shapes encountered in real life engines. The modeling of an abradable coating would also allow to take a step further in the fidelity of the solutions, but it remains an open problem for frequency domain methods.

## 9 Acknowledgments

The authors are grateful for the financial support of the ANR (project ANR-22-CPJ2-0061-01) and the European Union's Horizon HE-ART project (Hybrid Electric propulsion system for regional AiRcrafT, grant agreement ID: 101102013).

## References

[1] International Energy Agency. *Net Zero by 2050, a roadmap for the global energy sector*. 2021.

[2] M. Krack and J. Gross. Harmonic Balance for Nonlinear Vibration Problems. 2019. DOI: 10.1007/978-3-030-14023-6.

[3] S. Nacivet, C. Pierre, F. Thouverez, and L. Jezequel. "A dynamic Lagrangian frequency-time method for the vibration of dry-friction-damped systems". J. Sound Vib. Vol. 265, No. 1 (2003), pp. 201-219. DOI: 10.1016/S0022-460X(02)01447-5. hal:01635272.

[4] E. P. Petrov and D. J. Ewins. "Analytical Formulation of Friction Interface Elements for Analysis of Nonlinear Multi-Harmonic Vibrations of Bladed Disks". J. Turbomach. Vol. 125, No. 2 (2003), pp. 364-371. DOI: 10.1115/1.1539868. hal:04433199.

[5] M. Krack, L. Panning-von Scheidt, and J. Wallaschek. "On the Interaction of Multiple Traveling Wave Modes in the Flutter Vibrations of Friction-Damped Tuned Bladed Disks". J. Eng. Gas Turbines Power Vol. 139, No. 4 (2017), p. 042501. DOI: 10.1115/1.4034650.

[6] C. Pierre, A. A. Ferri, and E. H. Dowell. "Multi-Harmonic Analysis of Dry Friction Damped Systems Using an Incremental Harmonic Balance Method". J. Appl. Mech. Vol. 52, No. 4 (1985), pp. 958-964. DOI: 10.1115/1.3169175.

[7] D. Charleux, C. Gibert, F. Thouverez, and J. Dupeux. "Numerical and Experimental Study of Friction Damping Blade Attachments of Rotating Bladed Disks". Int. J. Rotating Mach. Vol. 2006 (2006), pp. 1-13. DOI: 10.1155/IJRM/2006/71302.

[8] L. Salles, L. Blanc, F. Thouverez, A. M. Gouskov, and P. Jean. "Dual Time Stepping Algorithms With the High Order Harmonic Balance Method for Contact Interfaces With Fretting-Wear". J. Eng. Gas Turbines Power Vol. 134, No. 3 (2012), p. 032503. DOI: 10.1115/1.4004236. hal:00975699.

[9] T. Heinze, L. Panning-von Scheidt, and J. Wallaschek. "Global detection of detached periodic solution branches of friction-damped mechanical systems". Nonlinear Dyn. Vol. 99, No. 3 (2020), pp. 1841-1870. DOI: 10.1007/s11071-019-05425-4.

[10] C. M. Firrone, S. Zucca, and M. M. Gola. "The effect of underplatform dampers on the forced response of bladed disks by a coupled static/dynamic harmonic balance method". Int. J. Non Linear Mech. Vol. 46, No. 2 (2011), pp. 363-375. DOI: 10.1016/j.ijnonlinmec.2010.10.001.

[11] B. Yang, M. Chu, and C. Menq. "Stick-slip-separation analysis and non-linear stiffness and damping characterization of friction contacts having variable normal load". J. Sound Vib. Vol. 210, No. 4 (1998), pp. 461-481. DOI: 10.1006/jsvi.1997.1305.

[12] T. M. Cameron and J. H. Griffin. "An Alternating Frequency/Time Domain Method for Calculating the Steady-State Response of Nonlinear Dynamic Systems". J. Appl. Mech. Vol. 56, No. 1 (1989), pp. 149-154. DOI: 10.1115/1.3176036. hal:01333697.
[13] Y. Colaïtis and A. Batailly. "The harmonic balance method with arc-length continuation in blade-tip/casing contact problems". J. Sound Vib. Vol. 502 (2021), p. 116070. DOI: 10.1016/j.jsv.2021.116070. hal-03163560.

[14] Y. Colaïtis and A. Batailly. "Stability analysis of periodic solutions computed for blade-tip/casing contact problems". J. Sound Vib. Vol. 538 (2022), p. 117219. DOI: 10.1016/j.jsv.2022.117219. hal-03764770.

[15] T. Vadcard, F. Thouverez, and 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. Eng. Vol. 419 (2024), p. 116641. DOI: 10.1016/j.cma.2023.116641. hal-04318743v1.

[16] T. Vadcard, F. Thouverez, and A. Batailly. "Computation of isolated periodic solutions for forced response blade-tip/casing contact problems". J. Eng. Gas Turbines Power Vol. 146, No. 4 (2024), p. 041011. DOI: 10.1115/1.4063704. hal-04298741v1.

[17] T. Vadcard, Y. Colaïtis, A. Batailly, and F. Thouverez. "Assessment of Two Harmonic Balance Method-Based Numerical Strategies for Blade-Tip/Casing Interactions: Application to Nasa Rotor67". J. Eng. Gas Turbines Power Vol. 144, No. 12 (2022). DOI: 10.1115/1.4055416. hal-03775621.

[18] 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 Vol. 139, No. 2 (2016). DOI: 10.1115/1.4034353.

[19] N. M. Newmark. "A Method of Computation for Structural Dynamics". J. Eng. Mech. Div. Vol. 85, No. 3 (1959), pp. 67-94. DOI: 10.1061/JMCEA3.0000098.

[20] N. J. Carpenter, R. L. Taylor, and M. G. Katona. "Lagrange constraints for transient finite element surface contact". Int. J. Numer. Meth. Eng. Vol. 32, No. 1 (1991), pp. 103-128. DOI: 10.1002/nme.1620320107. hal-01389918.

[21] L. Salles, B. Staples, N. Hoffmann, and C. Schwingshackl. "Continuation techniques for analysis of whole aeroengine dynamics with imperfect bifurcations and isolated solutions". Nonlinear Dyn. Vol. 86, No. 3 (2016), pp. 1897-1911. DOI: 10.1007/s11071-016-3003-y.

[22] D. Laxalde and F. Thouverez. "Complex non-linear modal analysis for mechanical systems: Application to turbomachinery bladings with friction interfaces". J. Sound Vib. Vol. 322, No. 4-5 (2009), pp. 1009-1025. DOI: 10.1016/j.jsv.2008.11.044. hal-00343494v3.

[23] M. Krack. "Nonlinear modal analysis of nonconservative systems: Extension of the periodic motion concept". Comput. Struct. Vol. 154 (2015), pp. 59-71. DOI: 10.1016/j.compstruc.2015.03.008. arXiv:2101.00949.

[24] T. Vadcard, F. Thouverez, and A. Batailly. "Nonlinear Normal Modes-Related Isolated Branches of Subharmonic Solutions for Forced Response Blade-Tip/Casing Contact Problems". Vol. 147. 5. 2024, p. 051022. DOI: 10.1115/1.4066713. hal-04812090v1.

[25] C. Touzé, A. Vizzaccaro, and O. Thomas. "Model order reduction methods for geometrically nonlinear structures: a review of nonlinear techniques". Nonlinear Dyn. Vol. 105, No. 2 (2021), pp. 1141-1190. DOI: 10.1007/s11071-021-06693-9. hal-03283647v2.

[26] M. Krack, L. Salles, and F. Thouverez. "Vibration Prediction of Bladed Disks Coupled by Friction Joints". Arch. Comput. Methods Eng. Vol. 24, No. 3 (2016), pp. 589-636. DOI: 10.1007/s11831-016-9183-2. hal-01825517.

[27] A. Millecamps, A. Batailly, M. Legrand, and F. Garcin. "Snecma's Viewpoint on the Numerical and Experimental Simulation of Blade-Tip/Casing Unilateral Contacts". Proceedings of the ASME Turbo Expo 2015. 2015. DOI: 10.1115/GT2015-42682. hal-01223582.

[28] A. Benaarbia, Y. Rae, and W. Sun. "Unified viscoplasticity modelling and its application to fatigue-creep behaviour of gas turbine rotor". Int. J. Mech. Sci. Vol. 136 (2018), pp. 36-49. DOI: 10.1016/j.ijmecsci.2017.12.008.
[29] E. Delhez, F. Nyssen, J.-C. Golinval, and A. Batailly. "Numerical study of bladed structures with geometric and contact nonlinearities". J. Sound Vib. Vol. 544 (2023), p. 117382. DOI: https://doi.org/10.1016/j.jsv.2022.117382. hal-03838701v1.

[30] F. Mashayekhi and S. Zucca. "Modal Derivatives for Efficient Vibration Prediction of Geometrically Nonlinear Structures with Friction Contact". Appl. Sci. Vol. 14, No. 9 (2024). DOI: 10.3390/app14093936.

[31] E. Petrov. "Method for direct parametric analysis of nonlinear forced response of bladed disks with friction contact interfaces". J. Turbomach. Vol. 126, No. 4 (2004), pp. 654–662. DOI: 10.1115/1.1776588.

[32] S. Quaegebeur, B. Chouvion, and F. Thouverez. "Nonlinear dynamic analysis of three-dimensional bladed-disks with frictional contact interfaces based on cyclic reduction strategies". Int. J. Solids Struct. Vol. 236-237 (2022), p. 111277. DOI: 10.1016/j.ijsolstr.2021.111277. hal-03429160.

[33] P. Wriggers and T. A. Laursen. Computational contact mechanics. 2006. DOI: 10.1007/978-3-540-32609-0.

[34] L. Woiwode, N. N. Balaji, J. Kappauf, F. Tubita, L. Guillot, C. Vergez, B. Cochelin, A. Grolet, and M. Krack. "Comparison of two algorithms for Harmonic Balance and path continuation". Mech. Syst. Sig. Process. Vol. 136 (2020), p. 106503. DOI: 10.1016/j.ymssp.2019.106503. hal-02424746.

[35] D. Laxalde, F. Thouverez, and J.-P. Lombard. "Forced response analysis of integrally bladed disks with friction ring dampers". J. Vib. Acoust. Vol. 132, No. 1 (2010). DOI: 10.1115/1.4000763. hal-00453467.

[36] T. Vadcard, A. Batailly, and F. Thouverez. "On Harmonic Balance Method-based Lagrangian contact formulations for vibro-impact problems". J. Sound Vib. Vol. 531 (2022), p. 116950. DOI: 10.1016/j.jsv.2022.116950. hal-03665624.

[37] L. Salles, L. Blanc, F. Thouverez, A. M. Gouskov, and P. Jean. "Dynamic Analysis of a Bladed Disk With Friction and Fretting-Wear in Blade Attachments". Proceedings of the ASME Turbo Expo 2009. 2009, pp. 465–476. DOI: 10.1115/GT2009-60151. hal-02537384.

[38] L. Pesaresi, J. Armand, C. Schwingshackl, L. Salles, and C. Wong. "An advanced underplatform damper modelling approach based on a microslip contact model". J. Sound Vib. Vol. 436 (2018), pp. 327–340. DOI: 10.1016/j.jsv.2018.08.014. hal-02369256.

[39] G. Von Groll and D. Ewins. "The harmonic balance method with arc-length continuation in rotor/stator contact problems". J. Sound Vib. Vol. 241, No. 2 (2001), pp. 223–233. DOI: 10.1006/jsvi.2000.3298. hal-01333704.

[40] V. Pagès, P. Duquesne, S. Aubert, L. Blanc, P. Ferrand, X. Ottavy, and C. Brandstetter. "UHBR Open-Test-Case Fan ECL5/CATANA". International Journal of Turbomachinery, Propulsion and Power Vol. 7, No. 2 (2022). DOI: 10.3390/ijtpp7020017. hal-03765815v2.

[41] A. P. Schneider, A.-L. Fiquet, B. Paoletti, X. Ottavy, and C. Brandstetter. "Experiments on Tuned UHBR Open-Test-Case Fan ECL5/CATANA: Performance and Aerodynamics". J. Turbomach. Vol. 146, No. 8 (2024). DOI: 10.1115/1.4064231. hal-04066082.

[42] A.-L. Fiquet, X. Ottavy, and C. Brandstetter. "UHBR Open-Test Case Fan ECL5/CATANA: Non-Linear Analysis of Non-Synchronous Blade Vibration at Part-Speed Conditions". J. Turbomach. Vol. 146, No. 7 (2024), p. 071003. DOI: 10.1115/1.4064841. hal-03800486.

[43] J. Al-Am, V. Clair, A. Giauque, J. Boudet, and F. Gea-Aguilera. "Aeroacoustic analysis of the tip-leakage flow of an ultrahigh bypass ratio fan stage". Phys. Fluids Vol. 35, No. 4 (2023), p. 047104. DOI: 10.1063/5.0146143.

[44] K. Billon, L. Sanchez, G. Bouvard, C. Gibert, L. Blanc, and F. Thouverez. "Experiment on tuned open-test-case fan ECL5/Catana : structural characterization under vacuum conditions". Proceedings of the ASME Turbo Expo 2024. 2024. DOI: 10.1115/GT2024-127539. hal-04679347.

[45] R. R. Craig and M. C. C. Bampton. "Coupling of substructures for dynamic analyses". AIAA J. Vol. 6, No. 7 (1968), pp. 1313–1319. DOI: 10.2514/3.4741. hal-01537654.

[46] A. Batailly, M. Legrand, A. Millecamps, and F. Garcin. "Numerical-Experimental Comparison in the Simulation of Rotor/Stator Interaction Through Blade-Tip/Abradable Coating Contact". J. Eng. Gas Turbines Power Vol. 134, No. 8 (2012), p. 082504. DOI: 10.1115/1.4006446. hal-00746632.
[47] C. M. Firrone and S. Zucca. "Modelling Friction Contacts in Structural Dynamics and its Application to Turbine Bladed Disks". *Numerical Analysis*. 2011. Chap. 14. DOI: 10.5772/25128.

[48] E. P. Petrov. "Multiharmonic Analysis of Nonlinear Whole Engine Dynamics With Bladed Disc-Casing Rubbing Contacts". *Proceedings of the ASME Turbo Expo 2012*. 2012, pp. 1181–1191. DOI: 10.1115/GT2012-68474.
## A Linear condensation procedure

In order to apply the linear condensation procedure, one should partition Eq. (3) according to the type of DoF: linear DoF, subscripted ●, are not subject to nonlinear forces, and nonlinear DoF, subscripted ●, that undergo nonlinear forces, such as

$$ \begin{bmatrix} \mathbf{Z}_{l,1} & \mathbf{Z}_{l,nl} \\ \mathbf{Z}_{nl,1} & \mathbf{Z}_{nl,nl} \end{bmatrix} \begin{bmatrix} \tilde{\mathbf{x}}_l \\ \tilde{\mathbf{x}}_{nl} \end{bmatrix} = \begin{bmatrix} \tilde{\mathbf{f}}_{\text{ex},1} \\ \tilde{\mathbf{f}}_{\text{ex,nl}} \end{bmatrix} + \begin{bmatrix} \tilde{\mathbf{f}}_{0,1} \\ \tilde{\mathbf{f}}_{0,nl} \end{bmatrix} + \begin{bmatrix} \mathbf{0} \\ \tilde{\mathbf{f}}_{nl,nl} \end{bmatrix}. \quad (33) $$

Expressing $\tilde{\mathbf{x}}_l$ as a function of $\tilde{\mathbf{x}}_{nl}$ then reads

$$ \tilde{\mathbf{x}}_l = \mathbf{Z}_{l,1}^{-1} (\tilde{\mathbf{f}}_{\text{ex},1} + \tilde{\mathbf{f}}_{0,1} - \mathbf{Z}_{l,nl} \tilde{\mathbf{x}}_{nl}). \quad (34) $$

Substituting Eq. (34) into Eq. (33) then leads to

$$ \mathbf{Z}_{nl} \tilde{\mathbf{x}}_{nl} = \tilde{\mathbf{f}}_{\text{ex},nl} + \tilde{\mathbf{f}}_{0,nl} + \tilde{\mathbf{f}}_{nl,nl} \quad (35) $$

where

$$ \mathbf{Z}_{nl} = \mathbf{Z}_{nl,nl} - \mathbf{Z}_{nl,1} \mathbf{Z}_{1,l}^{-1} \mathbf{Z}_{l,nl} \quad (36a) $$

$$ \tilde{\mathbf{f}}_{\text{ex,nl}} = \tilde{\mathbf{f}}_{\text{ex,nl}} - \mathbf{Z}_{nl,1} \mathbf{Z}_{1,l}^{-1} \tilde{\mathbf{f}}_{\text{ex},1} \quad (36b) $$

$$ \tilde{\mathbf{f}}_{0,nl} = \tilde{\mathbf{f}}_{0,nl} - \mathbf{Z}_{nl,1} \mathbf{Z}_{1,l}^{-1} \tilde{\mathbf{f}}_{0,1} \quad (36c) $$

## B Rigid body modes handling for harmonic 0

This appendix explains how the relative condensation of Sec. 2.2 is performed while the blade is subject to rigid body modes. Indeed, it makes its stiffness matrix singular and a workaround should be used. All the presented calculations here are only made for harmonic 0, corresponding to the static problem where the rigid body modes problem intervenes. Thus, the vector of excitation forces is not accounted for in the calculations since it is null on harmonic 0. This procedure inherently conducts to some numerical errors but these turn out to be negligible such that the results obtained through this condensation are not affected. The stiffness matrix of the blade is supposed to be invertible during all steps of the calculation and the use of its estimated value (evaluated by a procedure described in appendix B.2) is only used when absolutely necessary.

### B.1 Condensation procedure

The condensed system resulting from the linear condensation is firstly partitioned in terms of interface as follows

$$ \begin{bmatrix} \mathbf{Z}_b^{\text{fri}} & \mathbf{0} & \mathbf{Z}_b^{\text{fri,R/S}} \\ \mathbf{0} & \mathbf{Z}_d^{\text{fri}} & \mathbf{Z}_d^{\text{fri,R/S}} \\ \mathbf{Z}_b^{\text{R/S,fri}} & \mathbf{Z}_d^{\text{R/S,fri}} & \mathbf{Z}_d^{\text{R/S}} \end{bmatrix} \begin{bmatrix} \tilde{\mathbf{x}}_b^{\text{fri}} \\ \tilde{\mathbf{x}}_d^{\text{fri}} \\ \tilde{\mathbf{x}}^{\text{R/S}} \end{bmatrix} = \begin{bmatrix} \tilde{\mathbf{f}}_{0,b}^{\text{fri}} \\ \tilde{\mathbf{f}}_{0,d}^{\text{fri}} \\ \tilde{\mathbf{f}}_0^{\text{R/S}} \end{bmatrix} + \begin{bmatrix} \tilde{\mathbf{f}}_{nl,b}^{\text{fri}} \\ \tilde{\mathbf{f}}_{nl,d}^{\text{fri}} \\ \tilde{\mathbf{f}}_{nl}^{\text{R/S}} \end{bmatrix}. \quad (37) $$

After multiplying them by $(\mathbf{Z}_b^{\text{fri}})^{-1}$ and $(\mathbf{Z}_d^{\text{fri}})^{-1}$ respectively, the two first lines of Eq. (37) as a system read

$$ \tilde{\mathbf{x}}_b^{\text{fri}} + (\mathbf{Z}_b^{\text{fri}})^{-1} \mathbf{Z}_b^{\text{fri,R/S}} \tilde{\mathbf{x}}^{\text{R/S}} = (\mathbf{Z}_b^{\text{fri}})^{-1} \tilde{\mathbf{f}}_{0,b}^{\text{fri}} + (\mathbf{Z}_b^{\text{fri}})^{-1} \tilde{\mathbf{f}}_{nl,b}^{\text{fri}} \quad (38) $$

$$ \tilde{\mathbf{x}}_d^{\text{fri}} + (\mathbf{Z}_d^{\text{fri}})^{-1} \mathbf{Z}_d^{\text{fri,R/S}} \tilde{\mathbf{x}}^{\text{R/S}} = (\mathbf{Z}_d^{\text{fri}})^{-1} \tilde{\mathbf{f}}_{0,d}^{\text{fri}} + (\mathbf{Z}_d^{\text{fri}})^{-1} \tilde{\mathbf{f}}_{nl,d}^{\text{fri}}. \quad (39) $$

Subtracting Eq. (39) to Eq. (38) and accounting for the reciprocity of contact forces at the friction interface $(\tilde{\mathbf{f}}_{nl,b}^{\text{fri}} = -\tilde{\mathbf{f}}_{nl,d}^{\text{fri}} = \tilde{\mathbf{f}}_{nl}^{\text{fri}})$ leads to

$$ \tilde{\mathbf{x}}_b^{\text{fri}} - \tilde{\mathbf{x}}_d^{\text{fri}} + ((\mathbf{Z}_b^{\text{fri}})^{-1} \mathbf{Z}_b^{\text{fri,R/S}} - (\mathbf{Z}_d^{\text{fri}})^{-1} \mathbf{Z}_d^{\text{fri,R/S}}) \tilde{\mathbf{x}}^{\text{R/S}} = (\mathbf{Z}_b^{\text{fri}})^{-1} \tilde{\mathbf{f}}_{0,b}^{\text{fri}} - (\mathbf{Z}_d^{\text{fri}})^{-1} \tilde{\mathbf{f}}_{0,d}^{\text{fri}} + ((\mathbf{Z}_b^{\text{fri}})^{-1} + (\mathbf{Z}_d^{\text{fri}})^{-1}) \tilde{\mathbf{f}}_{nl}^{\text{fri}}. \quad (40) $$
By multiplying Eq. (40) by $Z_{rel} = ((Z_b^{fri})^{-1} + (Z_d^{fri})^{-1})^{-1}$, the following form is obtained

$$Z_{rel} \tilde{x}_{rel}^{fri} + Z_{rel} Z_{rel,a} \tilde{x}^{R/S} = \tilde{f}_{rel} + \tilde{f}_{nl}^{fri} \quad (41)$$

with

$$\tilde{f}_{rel} = Z_{rel} ((Z_b^{fri})^{-1} \tilde{f}_{0,b}^{fri} - (Z_d^{fri})^{-1} \tilde{f}_{0,d}^{fri}) \quad (42)$$

and

$$Z_{rel,a} = ((Z_b^{fri})^{-1} Z_b^{fri,R/S} - (Z_d^{fri})^{-1} Z_d^{fri,R/S}). \quad (43)$$

The calculation of $Z_{rel}$ involves the inverse of $Z_b^{fri}$, that is singular due to the presence of rigid body modes on the blade. The procedure allowing to evaluate $Z_{rel}$ and an estimate of $(Z_b^{fri})^{-1}$ is explained in appendix B.2.

Once the relative equation is condensed in the form Eq. (41), it is necessary to express the third line of Eq. (37) as a function of the new unknown $[\tilde{x}_{rel}^{fri}, \tilde{x}_{R/S}]^T$. This equation reads

$$Z_b^{R/S,fri} \tilde{x}_b^{fri} + Z_d^{R/S,fri} \tilde{x}_d^{fri} + Z^{R/S} \tilde{x}_{R/S} = \tilde{f}_0^{R/S} + \tilde{f}_{nl}^{R/S}. \quad (44)$$

The variables $\tilde{x}_b^{fri}$ and $\tilde{x}_d^{fri}$ are reexpressed thanks to Eq. (38) and Eq. (39), such that

$$\tilde{x}_b^{fri} = (Z_b^{fri})^{-1} (\tilde{f}_{0,b}^{fri} + \tilde{f}_{nl}^{fri} - Z_b^{fri,R/S} \tilde{x}^{R/S}) \quad (45)$$

$$\tilde{x}_d^{fri} = (Z_d^{fri})^{-1} (\tilde{f}_{0,d}^{fri} - \tilde{f}_{nl}^{fri} - Z_d^{fri,R/S} \tilde{x}^{R/S}) \quad (46)$$

Moreover, the nonlinear forces at the relative interface read

$$\tilde{f}_{nl}^{fri} = Z_{rel} \tilde{x}_{rel}^{fri} + Z_{rel} Z_{rel,a} \tilde{x}^{R/S} - \tilde{f}_{rel}, \quad (47)$$

it leads to

$$\tilde{\mathbf{x}}_{\mathrm{b}}^{\mathrm{fri}} = (\mathbf{Z}_{\mathrm{b}}^{\mathrm{fri}})^{-1} (\tilde{\mathbf{f}}_{0,\mathrm{b}}^{\mathrm{fri}} + \mathbf{Z}_{\mathrm{rel}}\tilde{\mathbf{x}}_{\mathrm{rel}}^{\mathrm{fri}} + \mathbf{Z}_{\mathrm{rel}}\mathbf{Z}_{\mathrm{rel,a}}\tilde{\mathbf{x}}^{\mathrm{R/S}} - \tilde{\mathbf{f}}_{\mathrm{rel}} - \mathbf{Z}_{\mathrm{b}}^{\mathrm{fri,R/S}}\tilde{\mathbf{x}}^{\mathrm{R/S}}) \quad (48)$$

$$\tilde{\mathbf{x}}_{\mathrm{d}}^{\mathrm{fri}} = (\mathbf{Z}_{\mathrm{d}}^{\mathrm{fri}})^{-1} (\tilde{\mathbf{f}}_{0,\mathrm{d}}^{\mathrm{fri}} - \mathbf{Z}_{\mathrm{rel}}\tilde{\mathbf{x}}_{\mathrm{rel}}^{\mathrm{fri}} - \mathbf{Z}_{\mathrm{rel}}\mathbf{Z}_{\mathrm{rel,a}}\tilde{\mathbf{x}}^{\mathrm{R/S}} + \tilde{\mathbf{f}}_{\mathrm{rel}} - \mathbf{Z}_{\mathrm{d}}^{\mathrm{fri,R/S}}\tilde{\mathbf{x}}^{\mathrm{R/S}}). \quad (49)$$

By injecting these expressions into Eq. (44), it reads

$$\begin{aligned} & Z_b^{R/S,fri} (Z_b^{fri})^{-1} \left( + Z_{rel} \tilde{x}_{rel}^{fri} + Z_{rel} Z_{rel,a} \tilde{x}^{R/S} - \tilde{f}_{rel} - Z_b^{fri,R/S} \tilde{x}^{R/S} \right) + \\ & Z_d^{R/S,fri} (Z_d^{fri})^{-1} \left( - Z_{rel} \tilde{x}_{rel}^{fri} - Z_{rel} Z_{rel,a} \tilde{x}^{R/S} + \tilde{f}_{rel} - Z_d^{fri,R/S} \tilde{x}^{R/S} \right) + Z^{R/S} \tilde{x}_{R/S} \\ & = \tilde{f}_{0,rel}^{R/S} + \tilde{f}_{nl}^{R/S}. \quad (50) \end{aligned}$$

Grouping the terms depending on $\tilde{x}_{rel}^{fri}$ and $\tilde{x}_{R/S}$ respectively, the following equation is obtained

$$\overbrace{\left[ Z_b^{R/S,fri} (Z_b^{fri})^{-1} - Z_d^{R/S,fri} (Z_d^{fri})^{-1} \right]}^{Z_{a,rel}} Z_{rel} \tilde{x}_{rel}^{fri} + \overbrace{\left[ Z_b^{R/S,fri} (Z_b^{fri})^{-1} (Z_{rel} Z_{rel,a} - Z_b^{fri,R/S}) + Z_d^{R/S,fri} (Z_d^{fri})^{-1} (-Z_{rel} Z_{rel,a} - Z_d^{fri,R/S}) + Z^{R/S} \right]}^{Z_{a,a}} \tilde{x}_{R/S} = \overbrace{Z_b^{R/S,fri} (Z_b^{fri})^{-1} (\tilde{f}_{rel} - \tilde{f}_{0,b}^{fri}) - Z_d^{R/S,fri} (Z_d^{fri})^{-1} (\tilde{f}_{rel} + \tilde{f}_{0,d}^{fri}) + \tilde{f}_0^{R/S} + \tilde{f}_{nl}^{R/S}}^{\tilde{f}_a} \quad (51)$$
Finally, the condensed equation is

$$ \underbrace{\begin{bmatrix} \mathbf{Z}_{\text{rel}} & \mathbf{Z}_{\text{rel}}\mathbf{Z}_{\text{rel},\text{a}} \\ \mathbf{Z}_{\text{a,rel}}\mathbf{Z}_{\text{rel}} & \mathbf{Z}_{\text{a,a}} \end{bmatrix}}_{\mathbf{Z}_r} \underbrace{\begin{bmatrix} \tilde{\mathbf{x}}_{\text{rel}}^{\text{fri}} \\ \tilde{\mathbf{x}}^{\text{R/S}} \end{bmatrix}}_{\tilde{\mathbf{x}}_r} = \underbrace{\begin{bmatrix} \tilde{\mathbf{f}}_{\text{rel}} \\ \tilde{\mathbf{f}}_a \end{bmatrix}}_{\tilde{\mathbf{f}}_r} + \underbrace{\begin{bmatrix} \tilde{\mathbf{f}}_{\text{nl}}^{\text{fri}} \\ \tilde{\mathbf{f}}_{\text{nl}}^{\text{R/S}} \end{bmatrix}}_{\tilde{\mathbf{f}}_{\text{nl},r}} \quad (52) $$

## B.2 Computation of $Z_{\text{rel}}$ and estimation of $(Z_b^{\text{fri}})^{-1}$

The matrix $Z_{\text{rel}}$ reads

$$ Z_{\text{rel}} = ((Z_b^{\text{fri}})^{-1} + (Z_d^{\text{fri}})^{-1})^{-1}, \quad (53) $$

however $(Z_b^{\text{fri}})^{-1}$ cannot be evaluated directly. Thus, an equivalent expression is found by assuming the properties of inverse matrices still apply, and it leads to

$$ Z_{\text{rel}} = ((Z_d^{\text{fri}})^{-1} + (Z_b^{\text{fri}})^{-1})^{-1} \quad (54a) $$

$$ = ((Z_b^{\text{fri}})^{-1} Z_b^{\text{fri}} (Z_d^{\text{fri}})^{-1} + (Z_b^{\text{fri}})^{-1} Z_d^{\text{fri}} (Z_d^{\text{fri}})^{-1})^{-1} \quad (54b) $$

$$ = ((Z_b^{\text{fri}})^{-1} (Z_b^{\text{fri}} + Z_d^{\text{fri}}) (Z_d^{\text{fri}})^{-1})^{-1} \quad (54c) $$

$$ = Z_d^{\text{fri}} (Z_b^{\text{fri}} + Z_d^{\text{fri}})^{-1} Z_b^{\text{fri}}. \quad (54d) $$

This new expression of $Z_{\text{rel}}$ does not suffer from the non-invertibility of the matrix $Z_b^{\text{fri}}$, since the only matrix that is actually inverted is $Z_b^{\text{fri}} + Z_d^{\text{fri}}$.

Finally, since a value of $(Z_b^{\text{fri}})^{-1}$ is required for the rest of the calculations, it is provided by using Eq. (54) such as

$$ (Z_b^{\text{fri}})^{-1} = Z_{\text{rel}}^{-1} Z_d^{\text{fri}} (Z_b^{\text{fri}} + Z_d^{\text{fri}})^{-1}. \quad (55) $$

Once again, this value can be estimated since its expression only requires the inversion of non singular matrices. This procedure induces numerical errors, however it was shown that it does not significantly affect the results.