# Geometric design of friction ring dampers in blisks using nonlinear modal analysis and Kriging surrogate model

Yekai Sun, Enora Denimal, Jie Yuan, Loïc Salles

► To cite this version:

Yekai Sun, Enora Denimal, Jie Yuan, Loïc Salles. Geometric design of friction ring dampers in blisks using nonlinear modal analysis and Kriging surrogate model. Structural and Multidisciplinary Optimization, 2022, 65 (3), pp.98. 10.1007/s00158-021-03093-w. hal-03591512

HAL Id: hal-03591512

https://inria.hal.science/hal-03591512

Submitted on 28 Feb 2022

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Geometric design of friction ring dampers in blisks using nonlinear modal analysis and Kriging surrogate model

Yekai Sun<fnref n="1" /> Ⓒ · Enora Denimal<fnref n="1" />,<fnref n="2" /> · Jie Yuan<fnref n="1" />,<fnref n="3" /> · Loïc Salles<fnref n="1" />,<fnref n="4" />

Received: 25 February 2021 / Revised: 27 September 2021 / Accepted: 17 October 2021

© The Author(s) 2022

## Abstract

Integrally bladed disks (blisk) have been widely used in the turbo-machinery industry due to its high aerodynamic performance and structural efficiency. A friction ring damper (FRD) is usually integrated in the system to improve its low damping. However, the design of the geometry of this FRD become complex and computationally expensive due to the strong nonlinearities from friction interfaces. In this work, we propose an efficient modelling strategy based on advanced nonlinear modal analysis and Kriging surrogate models to design and optimize the geometry of a 3D FRD attached to a high fidelity full-scale blisk. The 3D ring damper is parametrised with a few key geometrical parameters. The impact of each geometric parameter and their sensitivities to nonlinear dynamic response can be efficiently assessed using Kriging meta-modelling based on a few damped nonlinear normal modes. Results demonstrate that the damping performances of ring dampers can be substantially optimized through the proposed modelling strategy whilst key insights for the design of the rings are given. It is also demonstrated that the distribution of the contact normal load on the contact interfaces has a strong influence on the damping performances and can be effectively tuned via the upper surface geometry of the ring dampers.

**Keywords** Friction ring damper · Damped nonlinear normal mode · Geometric study · Surrogate model · Compressor blisk · Nonlinear vibration

## 1 Introduction

Industrial turbomachines are very important components in power generators and aircraft engines. The vibrational stress of turbomachines must be maintained under a certain level to avoid the structural failures of a rotating component due to the high cycle fatigue. The high modal density

and the complicated operating environment make the resonances unavoidable. Therefore, a good level of damping is necessary to maintain the operability of turbomachines. The friction damping is regarded as the major damping source compared to the material damping and the aerodynamic damping. Several types of friction dampers have been designed for the bladed disks in the literature (Pesaresi et al. 2017; Zucca et al. 2012; Brake et al. 2018). In classic bladed disks, the frictional joints between different components are the main source of friction damping that dissipates the vibration energy. However, these friction contact interfaces are absent in blisks, and so they have a low level of internal damping. To improve their damping, one solution is the use of friction ring dampers (FRDs), which are held by the centrifugal force and located under the groove of the blisks. These types of FRDs were firstly used in labyrinth seals (Niemotka and Ziegert 1995), then in gears (Zucca et al. 2010; Ye et al. 2017). The dynamic analyses of FRDs for blisks have been investigated by several researchers, with lumped parameter models (Laxalde et al. 2007; Sun et al. 2018), full-scale finite element (FE) models (Laxalde et al. 2010; Tang and Epureanu 2017; Wu et al. 2020; Sun et al.

Responsible Editor: Jianbin Du.

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

Yekai Sun
ys5113@imperial.ac.uk

Jie Yuan
jie.yuan@strath.ac.uk

<footnote n="1">¹ Department of Mechanical Engineering, Imperial College London, London SW7 2AZ, UK</footnote>

<footnote n="2">² Univ. Gustave Eiffel, Inria, COSYS/SII, I4S, Campus de Beaulieu, 35042 Rennes, France</footnote>

<footnote n="3">³ Aerospace Centre of Excellence, University of Strathclyde, Glasgow G1 1XQ, UK</footnote>

<footnote n="4">⁴ Center for Design, Manufacturing and Materials, Skolkovo Institute of Science and Technology, Moscow 121205, Russia</footnote>
2021b, c) and experimental testing (Laxalde et al. 2008). These works demonstrated the efficiency of FRDs to provide damping for blisks.

Designing friction dampers can significantly improve the damping performances and should be considered in the design process (Pollini et al. 2018). The shape and the geometry of friction dampers have a strong impact on the dynamic of blades as they modify its mass (Laxalde et al. 2010), the contact geometry (Panning et al. 2002; Sun et al. 2020) and the contact position (Niemotka and Ziegert 1995; Sun et al. 2020). A strong focus has been put so far by academics and industrials on the study of the impact of friction dampers geometry by designing and comparing different shapes, as wedge dampers (Krack et al. 2014; Pesaresi et al. 2017; Yuan et al. 2020), cylindrical dampers (Panning et al. 2002), conical dampers (Denimal et al. 2021b), asymmetric dampers (Gastaldi et al. 2018; Hüls et al. 2018), FRDs (Tang and Epureanu 2019; Sun et al. 2021b) or others. All results demonstrate without ambiguity that the shape of such component has a major impact on the dynamic of bladed discs. However, the design of such components is complex due to the presence of the frictional interfaces making the system strongly nonlinear (Salles et al. 2011; Mijar and Arora 2000; Krack et al. 2016). Thus, the first works for the design of friction dampers were based on large parametric studies, as they also allow to get insights in the influence of different parameters. In Sanliturk et al. (2001), a lumped model is used and solved with the harmonic balance method (HBM) to study the impact of the damper mass. The full nonlinear frequency response function (FRF) are computed each time for the different cases. In Panning et al. (2002), the mass and the platform angle are also studied based on the HBM for the computation of the nonlinear FRF and on a parametric study. Finally, in Sun et al. (2021b) a full 3D FE model of a blisk and of a FRD is studied and three different geometries of the FRD are compared. The nonlinear response is obtained using the damped nonlinear normal modes (dNNMs). A few works have attempted to optimise the damper properties. Thus, in Krack et al. (2014), the mass of a rigid underplatform damper is optimised by taking into consideration uncertainties in some contact parameters. The approach using dNNMs was also employed for the computation of the nonlinear properties of the system. In Yuan et al. (2020), the geometry of an underplatform damper of a two blades system is done by using surrogate modelling to speed up the simulations. The nonlinear dynamic response is computed using the HBM and continuation method, which represents a large computational time. In Denimal et al. (2021a), the topological optimisation of a 2D underplatform damper is achieved. The nonlinear response is also computed using the HBM representing a large computational time. To speed up the process, in Denimal et al. (2021c) a phase quadrature is added to the HBM formulations. To

focus on the case of FRDs, Tang and Epureanu (2019) tested several different cross sections of the FRDs using the nonlinear FRF. Whereas, results in Sun et al. (2021b) show that changing the cross section of the FRDs has limited effects on its damping performance (Sun et al. 2021b). However, the impact of the circumferential geometry of the FRDs has never been considered, which is the main objective of the present study. Indeed, previous studies have all demonstrated that a change in the contact normal load will have a strong influence on the damping properties. So it is expected that changing the circumferential geometry will have a strong effect on the damping properties as it directly influences the initial pressure. Thus, the objective of this paper is to study how a geometrical modification of a FRD influence the nonlinear dynamic response of a blisk. A full 3D FE model of a blisk will be considered, as well as for the FRD. Using 3D FE model for the FRD allows to compute the initial contact normal load and to take into consideration the inertial and elastic forces in the FRD.

In nonlinear FRF, the resonant solutions are identified by computing the dynamic responses with a certain excitation forcing level and a range of excitation frequency. To obtain the information of the resonant solutions under different excitation forcing levels, several repeated computations of the nonlinear FRF are required. This is usually done with the well established HBM with continuation method (Cameron and Griffin 1989; Allgower and Georg 2003) and used for study of the friction dampers (Pesaresi et al. 2017; Sanliturk et al. 2001; Yuan et al. 2020; Petrov 2004). It represents a large computational time, which makes unaffordable large parametric or optimisation for real industrial FE model of blisks with friction dampers. Therefore, the nonlinear modal analysis got some attention over the last two decades as it can capture the nonlinear dynamic responses at resonances, which can directly provide the critical information at resonance, i.e., the resonant solution, the resonant frequency, the damping ratio and the corresponding excitation forcing level (Krack 2015; Sun et al. 2021d). Nonlinear normal modes can be seen as the solution of an autonomous nonlinear system (Rosenberg 1960; Shaw and Pierre 1991). The concept has been extended to non-conservative nonlinear system (such as structures with friction dampers) with dNNMs (Sarrouy 2019). There are two concepts for the computation of the dNNMs for a full-scale structure with friction damping, which are complex nonlinear mode (Laxalde and Thouvez 2009) and extended periodic motion concept (Krack 2015). However, the use of dNNMs in an engineering problem is limited by that the relation between the dNNMs and the resonant solutions in nonlinear FRF cannot be directly addressed. In Sun et al. (2021a), the extended energy balance method (E-EBM) was proposed to predict the resonant solutions of the FRF from the dNNMs. The E-EBM is achieved by a balance of energy between the solutions in dNNMs
and resonant solutions in the nonlinear FRF. Similar energy analysis can be found in Hill et al. (2015), Vakakis (2001) and Grenat et al. (2018). Thanks to the dNNMs coupled with the E-EBM, the computation of the nonlinear FRF is not required to obtain the resonant solution, reducing drastically the numerical cost. From author's knowledge, this numerical methodology has never been used for the design of an industrial problem and will be used in the present work to demonstrate its ability to be used for the design of real-world systems.

Even though the dNNMs coupled with the E-EBM is efficient, it still represents an important computational time when large model are studied, especially when parametric studies on meshes are considered. To cope with this issue, a strategy based on surrogate modelling is adopted. A surrogate model, also called meta-model, has the following basic idea: from a small number of evaluations of an expensive model, a mathematical function which gives the same output for the same inputs is constructed (Cheng et al. 2020). Different methods allow to do so as Kriging (Sacks et al. 1989; Paulson and Ragkousis 2015), polynomial chaos (Wiener 1938; Xiu and Karniadakis 2002; Totis and Sortino 2020), support vector regression (Smola and Schölkopf 2004), etc. Using surrogate models has proven to be efficient in different cases for the design of friction dampers (Yuan et al. 2020; Denimal et al. 2021a). Kriging has shown good results in terms of computational cost and requires a low number of points to create an accurate surrogate model on large ranges of variation of design parameters (Denimal et al. 2021a). For this reason, it is adopted in the present work. For the first time, Kriging is used to surrogate the nonlinear modal properties at resonance of such systems through the dNNMs for nonlinear resonance computation. Once constructed and validated, these surrogate models are used to perform a large parametric study to investigate the impact of the geometric shape of the FRDs on the system nonlinear dynamics. The low computational cost of the surrogate models allows to perform a sensitivity analysis based on Sobol indices and an extensive analysis on the influence and the role of the FRD geometry can be performed.

So the objective of this work is to study the impact of FRD geometry, and more specifically the modification in of circumferential geometry, on the nonlinear modal properties of blisk. A full 3D FE model of a blisk system is studied. The FRD is modelled with a full 3D FE model, which is parametrised to change its shape. To reduce the computational time, the nonlinear modal properties of the system are predicted with Kriging meta-models, constructed based on a number of evaluations obtained with the dNNMs and E-EBM framework.

In Sect. 2, the FE model under study as well as the model-
ling strategy for the computation of the nonlinear dynamic
analysis are presented based on the dNNMs, and E-EBM are

presented. Then, in Sect. 3, the impact of the FRD geometry
is investigated. The parametrisation of the mesh is intro-
duced and Kriging meta-models constructed and validated in
the first part of Sect. 3. The validated Kriging meta-models
are then exploited to get insights in the role of the different
parameters. Finally, specific attention is put on the optimal
geometry to highlight the gain one can have by designing
the geometry of FRD.

## 2 Nonlinear dynamic analysis of blisk with friction ring damper

This section is devoted to the presentation of the model
of the blisk and of the FRD, and to the numerical strategy
employed for the computation of the nonlinear dynamic
response. First, the FE model of the academic blisk with
a FRD is shown in Sect. 2.1. Second, a general overview
of the numerical strategy for the computation of the non-
linear dynamic response of the full system is introduced in
Sect. 2.2, and then each step is detailed between Sects. 2.3
and 2.5.

2.1 Finite element model

A 3D FE model constructed with the commercial soft-
ware Altair HyperMesh is used to represent the blisk with
its FRD and is considered as a testcase for this study. The
full annulus and a sector of the model is shown in Fig. 1.
The blisk is composed of 30 sectors and cyclic symmetric
boundary conditions are applied on an elementary sector
(see Fig. 1b) to model the full annulus. Linear bricks ele-
ments are used and homogenous titanium with a Young’s
modulus of 115 GPa is used as the material of the blisk.
The material of the FRD is assumed to be identical with
the blisk structure. In reality, a softer and cheaper mate-
rial is expected for the friction damper compared with the

<figure><img src="image_3.png" /><figcaption>Fig. 1 A three-dimensional model of the blisk with friction ring damper: a the full annulus, and b an elementary sector</figcaption></figure>
<table><caption>**Table 1** Dimensions of the finite element model</caption><thead><tr><td>Description</td><td>Value (mm)</td></tr></thead><tbody><tr><td>Thickness of the blade</td><td>3</td></tr><tr><td>Height of the blade</td><td>87.5</td></tr><tr><td>Chord of the blade</td><td>52.5</td></tr><tr><td>Thickness of the ring damper</td><td>3.9</td></tr><tr><td>Width of the ring damper</td><td>9</td></tr><tr><td>Radius of tip of the blade</td><td>443.8</td></tr><tr><td>Radius of contact interface</td><td>321.2</td></tr><tr><td>Radius of lower surface of the disk</td><td>177.1</td></tr></tbody></table>

material used in the blisk. Different materials used in the FRD will lead to a different dynamic response because of the changes in the damper's mass and stiffness. The present work focuses on the geometric effects of the FRD and the influence from the material is not considered. One sector of the jointed structure is composed of 1344 elements. There is one contact interface between the underneath of the blisk and the upper surface of the FRD. The structure is built in a global Cartesian coordinate system (X-Y-Z system), where the Z-axis is the blisk rotational axis. The mesh convergence study is not presented here for the sake of consistency, but

the current grid size is small enough to capture the main features of mode investigated, and also gives acceptable computational times. The dimensions of the FE model are listed in the Table 1. To achieve a better demonstration, the cylindrical coordinate system ($\hat{r}-\hat{t}-Z$) is used for a sector instead of the global Cartesian coordinate system (X-Y-Z, please see Fig. 1). $\hat{r}$-Axis is the radial direction, whereas $\hat{t}$-axis represents the circumferential direction.

## 2.2 Outlines of the modelling process

The general overview of the workflow for the modelling process and computation of the nonlinear dynamic response of the blisk with a FRD is introduced as shown in Fig. 2 and has already been introduced in details in Sun et al. (2021b). Each step will be detailed in the following. After the creation of the FE model presented in Sect. 2.1, the computation of the nonlinear dynamic properties for one FRD can be divided into three main steps:

- step 1 identification of the resonant case based on rotor dynamic analysis and a Campbell diagram (see Sect. 2.3 for details),

<figure><img src="image_4.png" /><figcaption>**Fig. 2** Flowchart for modelling process for Sect. 2</figcaption></figure>

<figure><img src="image_5.png" /></figure>
- step 2 computation of the nonlinear static equilibrium and the corresponding contact conditions at the friction interface between the blisk and the FRD for the chosen resonant case (see Sect. 2.4 for details),

- step 3 a nonlinear dynamic analysis based on the dNNMs and the E-EBM is performed for the chosen resonant case and the given contact conditions. The vibration amplitude $|q_{res,z}|$, the modal damping ratio $\zeta$ and the resonant frequency $\omega_0$ are obtained (see Sect. 2.5 for details). This nonlinear dynamic analysis is achieved using an in-house Python script.

## 2.3 Rotor dynamic analysis and determination of the resonance

In a real operating situation, the blisk vibrates in the form of travelling wave mode under the rotational framework.

<figure><img src="image_6.png" /><figcaption>Fig. 3 Normalized natural frequency against nodal diameter from linear modal analysis: the 1st in-plane bending, the 1st torsion, the 2nd in-plane bending and the 1st disk mode. (Color figure online)</figcaption></figure>

The vibration of the blisk can be also described in a stationary framework by considering the vibration mode as a standing wave. The centrifugal effect caused by the engine rotation is included by adding a centrifugal load on the blisk, whereas the Coriolis forces are usually ignored in the dynamic analysis of the blisk. The effect of the Coriolis force has been proven to be small enough in the current compressor blisk (Ruffini 2016).

The linear modal analysis is completed using commercial FE solver ABAQUS. The nodal diameter (ND) diagram is given in Fig. 3, which indicates the first four families of modes for this FE blisk model. The natural frequencies of the 1st disk mode (1D), the 1st in-plane bending mode (1B), the 1st torsion mode (1T) and the 2nd in-plane bending mode (2B) are provided in Fig. 3 and the modeshapes are shown in Fig. 4. From the figures, the modal families dominated by the blade motion (i.e., 1B, 2B and 1T) do not significantly change with the ND. Whereas, for the disk dominated modes (1D), the natural frequency increases with the ND and a constant value can be seen for higher NDs. As the FRD aims to reduce the vibration amplitude for those disk dominated modes (Laxalde et al. 2007), the mode 1D with ND 3 is selected and used for the following study. The normalized natural frequency of this selected mode is 2.192.

The engine order (EO) excitation is a common source of excitation in the industrial turbomachinery due to the non-uniform pressure around the full annulus. It is caused by the nonuniformity of the intake, like the presence of vanes and a rotating blisk. To identify resonances caused by the EO excitation, the Campbell diagram of the blisk model is shown in Fig. 5. This is achieved by a rotor dynamic analysis using the commercial FE solver ABAQUS. The change of

<figure><img src="image_7.png" /><figcaption>Fig. 4 The modeshapes of the first four modes with nodal diameter 3: a the 1st disk dominated mode, b the 1st in-plane bending mode, c the 1st torsion mode, and d the 2nd in-plane bending mode. (Color figure online)</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Fig. 5** Campbell diagram from rotor dynamic analysis: natural frequency against engine rotational speed (red curves), engine order excitation (grey curves), excited engine order excitation (black curves), resonance (blue circle), selected resonance (blue star). (Color figure online)</figcaption></figure>

natural frequencies of the modal family 1D with respect to the engine rotational speed is demonstrated by red curves. The EO excitations are represented by the grey (unexcited) and black curves (excited). All possible resonance cases are highlighted by the markers. When the engine rotating speed is 2081 rpm and the 27th EO excitation is applied, the selected mode with a normalised resonant frequency $\omega_0 = 2.195$ is excited. This resonance case is marked by a star in Fig. 5 and used for the following analysis.

## 2.4 Nonlinear static analysis

The initial contact status of the different contact pairs the contact interface between the blisk and the FRD is obtained from a nonlinear static analysis using the commercial FE solver ABAQUS. The contact normal load ($N_0$) and contact gap ($\Delta$) distributions within the contact interface are extracted from the nonlinear static analysis. These initial contact $N_0$ and $\Delta$ will be used as the initial contact status for the different contact pair during the nonlinear modal analysis, as it will be explained later. For the default geometry (see Fig. 6a), the $N_0$ distribution and the $\Delta$ distribution are provided in Fig. 6. From the results, the nodes in the centre

<figure><img src="image_9.png" /><figcaption>**Fig. 6** Initial contact status from nonlinear static analysis: a contact normal load $N_0$ distribution, and b contact gap $\Delta$ distribution. (Color figure online)</figcaption></figure>

have a greater normal load, whereas nodes near boundaries show a larger gap.

## 2.5 Nonlinear dynamic analysis

To determine the nonlinear vibration amplitude at resonance and for a given forcing level $\gamma$, the dNNMs $Q$ of the blisk with its FRD are computed and the E-EBM is applied to relate the computed dNNM with the resonant solution of the forced nonlinear system. This part is devoted to the presentation of this approach. First, the reduction strategy applied to the structural matrices (mass and stiffness) is introduced. Second, the nonlinear contact forces are presented and finally the dNNM framework and the E-EBM are presented.

### 2.5.1 Three-step model reduction

To ensure an affordable numerical cost, the three-dimensional FE model is reduced. Generally, blisks are assumed to be perfectly tuned structures. Therefore, the full annulus of the blisk and the FRD can be modelled using a sector with cyclic symmetric boundary conditions (see Fig. 1). The mass and stiffness matrices of the blisk sector and the sector of the FRD are calculated using the commercial FE solver ABAQUS. At this stage, the number of degrees of freedom (DoFs) is equal to 6057. A three-step model reduction is then applied here to reduce the size of the system. First, the internal DoFs are eliminated using the well-known Craig-Bampton method (1968). After this first reduction, the number of DoFs is equal to 1146. Second, Mead's approach is applied to project the DoFs corresponding to the left cyclic symmetric boundaries onto the DoFs related to the right cyclic symmetric boundaries (Mead 1975), reducing the number of DoFs to 783. Finally, the retained DoFs related to cyclic symmetric boundaries are further reduced using partial interface mode method (Tran 2009). A similar reduction process is applied on both the sector of blisk and the FRD respectively. At the end, the final number of DoFs is equal to 510 and one gets the reduced mass matrix $\underline{M}$ and the reduced stiffness matrix $\underline{K}$ of the blisk and FRD.

### 2.5.2 Contact model

A classic 3D contact element is used to model the frictional contact (Yang and Menq 1998). A combination of two Jenkins elements (in the contact plane) and a unilateral spring in the normal direction is used to calculate the contact forces, as shown in Fig. 7. This classic 3D contact element has been widely used for the nonlinear dynamic analyses of structures with friction joints (Zucca et al. 2010; Tang and Epureanu 2019; Pesaresi et al. 2017). In this contact model, $k_n$ is the
<figure><img src="image_10.png" /><figcaption>**Fig. 7** 3D contact model: two Jenkins elements in the contact plane and a unilateral spring in the normal direction</figcaption></figure>

normal contact stiffness; $k_{t1}$ and $k_{t2}$ are the tangential contact stiffness. A prescribed gap $\Delta$ or normal load $N_0$ corresponds to the initial contact status for each contact pair. It obtained from the nonlinear static analysis described in Sect. 2.4. Node pairs with an initial gap $\Delta$ are separated at the static equilibrium position, whereas nodes with a normal load $N_0$ are in contact. With this contact element, there are three basic contact conditions: sticking, sliding and separation. The use of such friction model for the computation of the nonlinear FRF to study the structure with friction damping has been well validated by numerous experimental and numerical works in the field of nonlinear dynamics and turbomachinery industry (Zucca et al. 2012; Panning et al. 2002; Petrov 2004; Pesaresi et al. 2017; Yuan et al. 2021; Laxalde et al. 2008; Schwarz et al. 2020).

### 2.5.3 Nonlinear modal analysis

The dynamic equation of the bladed disk under forcing is defined by Eq. 1, where $\underline{M}$ is the reduced mass matrix, $\underline{K}$ is the reduced stiffness matrix, $\underline{C}$ is the structural reduced damping matrix (from material dissipation for example), $\underline{F_c}$ is the vector of the friction force within the contact interface and $\underline{F_c}$ is the excitation force with $\gamma$ the forcing amplitude, $\varphi$ the absolute phase of the excitation force and $\Omega$ the excitation frequency. The underline symbol is used to represent the vector and the time variable $t$ does not appear explicitly in all equations in this work for the sake of clarity.

$$ \underline{M} \ddot{\underline{Q}}_f + \underline{C} \dot{\underline{Q}}_f + \underline{K} \underline{Q}_f + \underline{F}_c(\underline{Q}_f, \dot{\underline{Q}}_f) = \underline{F}_c(\gamma, \varphi, \Omega, t). \quad (1) $$

On the other side, the dNNMs of the bladed disk are defined by Eq. 2 where $\underline{Q}(t)$ is the modal displacement of the dNNM for the blisk with FRD. The dNNMs can be seen as the solutions of a free vibration problem of such systems, and the external excitation force and viscous damping terms are usually not considered (Krack 2015; Laxalde and Thouverez 2009).

$$ \underline{M} \ddot{\underline{Q}} + \underline{K} \underline{Q} + \underline{F}_c(\underline{Q}, \dot{\underline{Q}}) = 0. \quad (2) $$

With the E-EBM framework, the solution of (1) can be obtained from the solutions of (2). The first step consists in the computation of the dNNMs, and then in equating some energetic terms. They are detailed in the following.

#### 2.5.3.1 Computation of damped nonlinear normal modes
Unlike linear normal modes, the nonlinear modes and other modal properties are highly energy-dependent, i.e. the modeshapes and the resonant frequency vary with the level of system energy. To capture this energy-dependency, a new parameter, namely the modal amplitude $\alpha$, is introduced into the system to quantify the modal energy. Then, modal displacements are expressed as: $\underline{Q} = \alpha \cdot \underline{Q}_0$, where $\underline{Q}_0$ is the mass normalized modal displacement. Furthermore, the absolute phase of the mass normalized solution $\underline{Q}_0$ is arbitrary for an autonomous system. To compute the mass normalized solution, two extra constraints are applied, which are the mass normalisation and the phase normalisation. Readers can refer to Krack (2015) and Sun et al. (2021d) for further details.

To calculate the dNNM as a periodic solution for a non-conservative nonlinear system, an artificial mass proportional damping term is introduced into the system to balance the energy dissipated by the friction forces (Krack 2015; Sun et al. 2021d). The artificial damping matrix $\underline{\underline{C}}$ is expressed as $\underline{\underline{C}} = -2\omega_0\zeta\underline{\underline{M}}$ where $\zeta$ is the modal damping ratio and $\omega_0$ is the resonant frequency. The periodic solutions of the non-conservative nonlinear system are obtained by solving Eq. 3. This artificial damping compensates the friction damping generated by the FRD and so makes it possible to quantify it.

$$ \alpha \cdot \underline{\underline{M}} \ddot{\underline{Q}}_0 + \alpha \cdot \underline{\underline{C}} \dot{\underline{Q}}_0 + \alpha \cdot \underline{\underline{K}} \underline{\underline{Q}}_0 + \underline{\underline{F}}_c(\alpha \cdot \underline{\underline{Q}}_0) = 0. \quad (3) $$

To solve the nonlinear dynamic equation 3, the classic HBM (Cameron and Griffin 1989) is used. Readers are invited to refer Krack and Gross (2019) for further details. The modal displacement $\underline{\underline{Q}}_0$, the resonant frequency $\omega_0$ and the modal damping ratio $\zeta$ can be easily computed for a given modal amplitude $\alpha$ (Krack 2015; Sun et al. 2021d). The evolution of the modal properties against the vibration amplitude are obtained with the well-known continuation method (Allgower and Georg 2003) to track the evolution of the dNNM over a range of modal amplitudes.

#### 2.5.3.2 Extended energy balance method
The dNNMs can be used to represent the resonant solution in a forced system. The E-EBM is a numerical method to provide a direct relation between the dNNM and the resonant solutions under forcing, i.e. the solution $\underline{\underline{Q}}_f$ of Eq. 1. The E-EBM assumes that the resonant solution of the nonlinear FRF and that of the dNNMs are similar for a small enough level of damping and forcing. It has been demonstrated in Sun et al. (2021a) and Yuan et al. (2022) that for friction dampers, this condition is verified and so the E-EBM can be employed. Therefore, these solutions can be assumed to be identical, and are named as 'resonant shared solutions'. The solution of an
autonomous system bears resemblance with that of the non-linear FRF at resonance, so when the excitation frequency $\Omega$ equals the resonant frequency $\omega_0$, the resonant shared solution $\underline{\psi}$ is given in Eq. 4.

$$ \underline{\psi} = \underline{Q} \approx \underline{Q}_f, \quad \text{when} \quad \Omega = \omega_0. \qquad (4) $$

For each level of the modal amplitude $\alpha$, the resonant shared solution $\underline{\psi}$, the resonant frequency $\omega_0$ and the modal damping ratio $\zeta$ are known from the computation of the dNNMs, as presented in previously. From the E-EBM, if the viscous damping term $\underline{C}$ and the vector of the excitation force are given are given, then the excitation forcing level $\gamma$ and the absolute phase of the excitation force $\varphi$ can be determined. In other words, the E-EBM makes possible to find the correspondence between the dNNMs and a modal amplitude $\alpha$, and the corresponding excitation level $\gamma$ and forcing phase $\varphi$.

Recalling the Eqs. (1)–(2) and assuming that the solution of the nonlinear FRF equals the resonant shared solution $\underline{\psi}$, the subtraction of Eq. (3) to Eq. (1) must hold in its weak form (Sun et al. 2021a). So, the energy exchange caused by the forcing term $E_f(\gamma, \varphi)$ can be computed by integrating the excitation force times the velocity over one vibrational period, as shown in Eq. 5. The energy transfer caused by damping terms $E_{damp}$ can be obtained in a similar way. $E_{damp}$ is constant for a certain level of modal amplitude $\alpha$, whereas the energy exchanged by the external excitation $E_f(\gamma, \varphi)$ varies with the forcing phase $\varphi$ and the excitation forcing level $\gamma$. The values of $\varphi$ and $\gamma$ are determined by finding the single intersection between the curve of $E_f(\gamma, \varphi)$ and the constant value of $E_{damp}$. Readers are invited to refer Sun et al. (2021a) and Sun et al. (2021d) for a complete and exhaustive description of the approach.

$$ E_f(\gamma, \varphi) = \int_0^{2\pi/\omega_0} \underline{F}_e(\gamma, \varphi, \omega_0, t) \cdot \underline{\psi} \cdot \underline{\delta t}, \qquad (5) $$

$$ E_{\text{damp}} = \int_{0}^{2\pi/\omega_0} (\underline{\underline{C}} + 2\omega_0\zeta\underline{\underline{M}}) \cdot \underline{\psi} \cdot \underline{\psi} \cdot \underline{\delta t}. \qquad (6) $$

After the application of the E-EBM, one gets the values of $\varphi$ and $\gamma$ that satisfy $E_f(\gamma, \varphi) = E_{damp}$. From this, it is possible to get directly the solution $\underline{Q}_f$ from the previous dNNMs computation. This means there is no need for solving Eq. 1 to get the resonant solutions of the nonlinear FRF.

## 2.6 Illustration on the default geometry

To illustrate the general workflow to compute the nonlinear dynamic response, results obtained on the default geometry are given and commented for a better understanding of the reader. To provide a better physical interpretation, the vibration amplitude is taken from the displacement in the Z direction from the response node $|q_{res,z}|$ (see Fig. 1). The evolutions of the normalised $\omega_0$, $\zeta$ and $\gamma$ against the vibration amplitude $|q_{res,z}|$ as shown in Fig. 8a, b. When the blisk vibrates at a low vibration amplitude $|q_{res,z}|$, the contact nodes are either sticking or separated resulting in a constant resonant frequency $\omega_0$ and a zero modal damping ratio $\zeta$. The $\omega_0$ at zero amplitude is 2.195, which can be seen as a linear sticking natural frequency. This linear sticking natural frequency is obtained by considering the contact pairs in either sticking or separation condition based on the pre-loading (normal load or gap), i.e., sticking contact pairs (with normal load) connected by the linear springs in three dimensions. This can be explained by the fact that

<figure><img src="image_11.png" /><figcaption>Fig. 8 Damped nonlinear normal modes using nonlinear modal analysis: a normalised resonant frequency $\omega_0$ against amplitude $|q_{res,z}|$, b modal damping ratio $\zeta$ against amplitude $|q_{res,z}|$, and c excitation forcing level $\gamma$ against amplitude $|q_{res,z}|$</figcaption></figure>
the nonlinear dynamic response lies on the linear sticking response before the sliding occurs. By continuing to increase the $|q_{res,z}|$, more and more contact nodes experience sliding. A softening effect is observed and so the $\omega_0$ decreases (see Fig. 8a). The total energy of the blisk is dissipated due to the rubbing motion of those sliding nodes. Therefore, a positive $\zeta$ is obtained. The maximum $\zeta$ occurs with a value of 0.15% with a $|q_{res,z}|$ of 0.164 mm and a $\omega_0$ of 2.187.

The results obtained from the E-EBM are illustrated in Fig. 8c, where for each vibration amplitude $|q_{res,z}|$ the corresponding excitation forcing level $\gamma$ is provided. For a given point on this curve, the excitation forcing level is predicted when the resonant solution of the forced system has the similar solution with the dNNM. For a better interpretation, an extensive example is provided in Fig. 9 where the full FRF is computed and compared to the E-EBM results (see Fig. 9a) where the left part shows the $\gamma-|q_{res,z}|$ plot obtained from the E-EBM and the right part presents the normalised resonant frequency $\omega_0-|q_{res,z}|$ plot obtained from the computation of the FRF. The nonlinear FRF with different values of $\gamma$ are computed ($\gamma = 20$ N, 60 N, 80 N, and 101 N), and for each the vibration amplitude at resonance is reported on the $\gamma-|q_{res,z}|$ plot of the left. An excellent agreement is found with the resonant solutions in the nonlinear FRF according to a similar value in $\omega_0$ and $|q_{res,z}|$.

Finally, the energy dissipation ($E_d$) for each contact pairs is also calculated for the predicted resonance from the dNNMs and resonance from the nonlinear FRFs, which show a great agreement. Only the $E_d$ distributions within the contact interface from the dNNMs are shown in Fig. 9b. One can notice that the nodes in the middle can dissipate more energy than the other nodes.

Previous works for the design of friction dampers are based on the computation of the full nonlinear FRF to extract the resonant information (Yuan et al. 2020; Panning et al. 2002; Denimal et al. 2021a). This is extremely numerically very expensive as it requires many simulations. On

the opposite, using the dNNMs and the E-EBM, as here, the values of $\omega_0$, $\zeta$ and $|q_{res,z}|$ at the resonance can be efficiently predicted as the computation of the nonlinear FRFs is not required. To get a rough idea of the simulation time, using dNNMs and the E-EBM requires about 36 min for one simulation of the blisk and the computation of full nonlinear FRF with HBM and continuation takes about 48–60 min. If one wants to use the nonlinear FRF to construct the similar curves in Fig. 8, at least 80 nonlinear FRFs are required and more than 60 computational hours are expected. Therefore, the use of dNNM coupled with E-EBM is a very efficient way to predict the resonant solutions in the present case study, and makes the geometric design of the blisk with FRD achievable. The high computational time still remains critical due to the model sizes. To cope with this issue a strategy based on meta-modelling, and more precisely on Kriging, is employed to investigate the geometric effects of the FRD on the dynamic behaviour of the blisk.

### 3 Geometric exploration of ring damper

In Sect. 2, the numerical strategy of the nonlinear dynamic analysis for a blisk with a FRD is introduced and illustrated on the default FRD geometry. In this section, the influence of the damper geometry on the nonlinear dynamic behaviour is studied. As reminder, with such a nonlinear analysis based on the dNNMs and E-EBM, the modal damping ratio $\zeta$, the resonant frequency $\omega_0$ and the vibration amplitude at resonance $|q_{res,z}|$ are obtained efficiently without the computation of the full FRF. Even if this approach gives competitive computational time compared to others numerical approaches, it remains important and critical to consider parametric studies. To cope with this issue a strategy based on meta-modelling, and more precisely on Kriging, is employed to investigate the geometric effects of the FRD on the dynamic

<figure><img src="image_12.png" /><figcaption>**Fig. 9** a Numerical validation of the $\omega_0-|q_{res,z}|$ plot and $\gamma-|q_{res,z}|$ plot using computed nonlinear frequency response function, and b energy dissipation $E_d$ distribution for the resonant solutions. (Color figure online)</figcaption></figure>
behaviour of the blisk. First, the FRD is parametrized with a few key design parameters and their effects are illustrated for some limits cases. Second, Kriging surrogate models are constructed and validated. They are then employed to predict efficiently the nonlinear behaviour of the blisk for different FRD geometries, and to perform a global sensitivity analysis to get deep insights in the impact of the FRD geometry on the dynamic behaviour of the blisk. Finally, the focus is put on the best geometry. The general framework is summarized in Fig. 10.

## 3.1 Damper geometry parametrisation

### 3.1.1 Design parameters

The coordinate in radial direction of the FRD can vary (in the $\hat{r}$-axis), but the coordinates in the rotational axis (in the Z-axis) and the circumferential direction (in the $\hat{t}$-axis) are unchanged (see Figs. 4, 11 for axis names). The geometry of the FRD is parametrized using three design parameters,

$\kappa_1, \kappa_2$ and $\kappa_3$. The first one describes the upper surface, when the two others describe the lower surface. The sector of the FRD is supposed to be symmetric in its centre about the Y-axis. Each parameter corresponds to the height of a bump (or a hole for a negative height) in different locations of the damper. The geometric modification of the FRD comes from several considerations:

- Modifying the initial contact status from Sun et al. (2021c), one can notice that the initial contact status (either the normal load or initial gap) has a significant influence on the damping ability of the damper. Thus, the aim of the geometric design of the FRD is to improve the damping performances through a modification of the initial contact status.

- Modifying the geometry in the circumferential direction from Fig. 9b, one can see that contact nodes under the blade (in the middle) dissipate more energy (higher $E_d$) than the other nodes. To improve the dissipative ability of the FRD, one would expect to increase the dissipative ability of the other contact nodes (the

<figure><img src="image_13.png" /><figcaption>Fig. 10 Flowchart for geometric exploration for Sect. 3</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Fig. 11** Different geometries of a sector of ring damper: **a** default shape, and **b-g** other six geometries when varying one the parameter, left: $\hat{r}$ – $\hat{i}$ view of geometries, middle: contact normal load $N_0$ distributions, right: contact gap $\Delta$ distributions. (Color figure online)</figcaption></figure>

contact nodes that are not under the blade), so changing the shape of the FRD in the circumferential direction ($\hat{r}$ axis) is considered instead of in the Z-axis (the cross section).

- *Modifying the upper surface* from the results in Sun et al. (2020) and Panning et al. (2002), different contact angles can lead to a distinctive damping performance. The modification of the upper surface of the FRD is considered, because it can lead to different local contact angles.

- *Modifying the lower surface* in Sun et al. (2018) and Laxalde et al. (2007, 2010), different mass of the FRD can result in a different normal load leading to a different damping ability. To change the mass of the FRD, the geometric modification of the lower surface of the FRD is included.

- *Adding bumps or concaves* a continuous and smooth upper surface of the FRD is desired based on the consideration of easy fabrication. Therefore, only one bump

or a concave shape is considered for the upper surface. Whereas for the lower surface, the idea is to have a wavelet shape of the FRD. A wavelet shape of the FRD has a great potential to further improve the damping ability by treating the additional mass (bump) as a tune vibration damper (Lupini et al. 2019). In this work, only three bumps/concaves are included to illustrate the potential of this idea.

For the upper surface, only one bump over the full surface is considered, whereas three different bumps can be considered for the lower surface: one in the middle and one on each side, these outer two are similar to ensure the symmetry of the damper. Each parameter varies between -1 and 1, and the default geometry corresponds to the set [0, 0, 0] which is represented in Fig. 11a. This corresponds to a parametrization with three parameters: one for the upper bump, one for the lower central bump and one for the lower side bumps.
Seven extreme geometries are obtained by only changing one of the design parameters to its maximal or minimal and keeping the rest unchanged. These seven geometries are used to illustrate the geometric parametrisation as shown in by Fig. 11. To go into more details:

- $\kappa_1$ Controls the upper surface (i.e. the contact interface and contact geometry) of the FRD using a spline function. A bump in the centre can be obtained when $\kappa_1 > 0$ (see Case A Fig. 11b), whereas a negative $\kappa_1$ leads to a concave surface (Case B Fig. 11c). When $\kappa_1$ varies, the mass of the damper as well as the local coordinate at the contact interface are impacted. It is worth emphasizing that the shape of the groove underneath the blisk always changes with the FRD geometry to ensure that the contact points are always overlapped.

- $\kappa_2$ Controls the two outer bumps of the lower surface, similarly to $\kappa_1$ it controls the geometry with a spline function and positive values lead to two concave shapes, when negative values lead to two bumps (one on each side). See Case C and D in Fig. 11d, e for graphical illustrations. Deforming the lower surface can change the flexibility of the ring and the mass distribution along the circumferential direction but has no impact on the local coordinate within the contact interface.

- $\kappa_3$ Controls the middle bump of the lower surface, it controls the geometry with a spline function and positive values lead to a concave shape, when negative values lead to a bump. See Case E and F in Fig. 11f, g for graphical illustrations. Effects on the structure are similar to those for $\kappa_2$.

When one of the parameters $\kappa_i$ with $i \in \{1, 2, 3\}$ is modified, the location of the geometrical points of the mesh are modified as well. Each parameter is normalized to an interval $[-1, 1]$. The physical values $\delta r$ for each normalized parameters are given in Table 2. The physical values $\delta r$ represent the change of the radial coordinate (r-axis), leading to a bump (positive $\delta r$) or a concave (negative $\delta r$). The geometry of the upper surface and the lower surface is defined using a quadratic spline function by fixing the radial coordinate of the central node with its given physical value $\delta r$. To keep a consistent mesh quality, the number of nodes and elements of the FRD are unchanged when the FRD is deformed. During this process, the physical coordinates of

<table><caption>**Table 2** Range of variation of the geometrical parameters</caption><thead><tr><th>Parameter</th><th>Normalized values</th><th>Physical values δr (mm)</th></tr></thead><tbody><tr><td>κ<sub>1</sub></td><td>[-1, 0, 1]</td><td>[-1.5, 0, 2.5]</td></tr><tr><td>κ<sub>2</sub></td><td>[-1, 0, 1]</td><td>[-2.0, 0, 1.5]</td></tr><tr><td>κ<sub>3</sub></td><td>[-1, 0, 1]</td><td>[-1.5, 0, 1.0]</td></tr></tbody></table>

each node for the FRD are automatically updated according to the defined upper and lower surface.

### 3.1.2 Effects on the contact status

For each geometry, before computing the systems dynamic response, the $N_0$ and $\Delta$ distributions must be determined to initialize the different contact properties. This is done through a nonlinear static analysis where centrifugal effects are considered using the commercial FE solver ABAQUS as described in Sect. 2.4. In order to illustrate the impact of the damper geometry on such contact properties, the $N_0$ and $\Delta$ for the seven geometries are given as an example in Fig. 11 in the middle and right columns, respectively.

When $\kappa_1$ is positive, an increasing number of inner nodes are in contact and less outer nodes are separated (see Case A in Fig. 11b). Whereas when $\kappa_1$ is negative, it leads to more separated nodes with larger gaps (see Case B in Fig. 11c). The values of the $N_0$ stay nevertheless similar. When $\kappa_2$ is equal to 1, the $N_0$ is concentrated in the middle under the blade with higher values (up to 2e3 kPa, see Case C in Fig. 11d). It also brings the largest gap near boundaries amongst all seven geometries. A more dispersive $N_0$ distribution can be seen when $\kappa_2$ is -1. The highest $N_0$ is observed when $\kappa_3$ is -1 and, the $N_0$ located under the blade, it also brings large gaps for outer nodes. From these different cases, one can see that by changing of three parameters, a large variety of initial pressure and gaps distributions can be reached.

### 3.1.3 Effects on damped nonlinear normal modes

To get a quick overview on the impact of the damper geometry on the system dynamics, the evolution of the dNNMs, the resonant frequency $\omega_0$, the modal damping ratio $\zeta$ and the predicted excitation forcing level $\gamma$ with respect to the vibration amplitude $|q_{res,z}|$ are computed for the seven extreme geometries using the modelling process in Sect. 2 and Fig. 2. The results are demonstrated in the form of $\omega_0|q_{res,z}|$ plots, $\zeta|q_{res,z}|$ plots and $\gamma|q_{res,z}|$ plots (see Fig. 12). The dNNMs for these seven geometries are used to initially demonstrate the effects of each design parameter.

The results of the default geometry are represented by the black curves. Comparing the results for the different geometries (see Fig. 12), it appears that the $\omega_0$, the $\zeta$ and the $\gamma$ are very sensitive to $\kappa_1$ since the largest variation in the dynamic response is observed by comparing Case A and Case B with the default case. As described above, a positive $\kappa_1$ (Case A) leads to a more dispersive $N_0$ distribution and less separated contact nodes. In nonlinear modal analysis, a lower $\omega_0$ and a better damping performance are observed. Indeed, the $\omega_0$ equals 2.185 at a low level of $|q_{res,z}|$ and the $\zeta$ reaches 0.28%. Moreover, for the same level of $\gamma$, the $|q_{res,z}|$ are always lower
<figure><img src="image_15.png" /><figcaption>Fig. 12 Damped nonlinear normal modes for all seven geometries (Fig. 11): a normalised resonant frequency $\omega_0$ against amplitude $|q_{res,z}|$, b modal damping ratio $\zeta$ against amplitude $|q_{res,z}|$, and c excitation forcing level $\gamma$ against amplitude $|q_{res,z}|$</figcaption></figure>

for Case A. Less variations are observed for $\kappa_2$ and $\kappa_3$, especially on the $\omega_0$ and the $|q_{res,z}|$. However, large variation can be observed on the $\zeta$. Conclusively, changing $\kappa_1$ can improve the damping performance greatly with a significant shift of $\omega_0$ and allows for large variations of the dynamic response. Whereas, different $\kappa_2$ and $\kappa_3$ can ensure a better damping performance with a negligible shift in $\omega_0$.

## 3.2 Kriging meta-models

### 3.2.1 Kriging: some fundamentals

Considering the high computational cost associated to the computation of the dNNMs, the use of a surrogate model is considered here. From a set of inputs and outputs of the nonlinear analysis, a surrogate model is built and then used to predict the dNNM. This approach has the main advantage to be non-intrusive. Kriging is used here as it is the best linear unbiased estimator and has proved to reduce drastically the numerical simulation time. From the knowledge of the authors, it has never been used yet for the prediction of the dNNMs.

In the present study, the input parameters are the three design parameters, namely $\kappa_1$, $\kappa_2$ and $\kappa_3$. For each modal amplitude $\alpha$, four Kriging meta-models are built: one for the vibration amplitude $|q_{res,z}|$, one for the resonant frequency $\omega_0$, one for the modal damping ratio $\zeta$ and one for the excitation forcing level $\gamma$. In the following, $\underline{x} = [\kappa_1, \kappa_2, \kappa_3]$ will denote the vector of inputs, and $\gamma$ will

denote one of the considered output. A set $\{(\underline{x}^{(i)}, y^{(i)})\}$ of N points are used to train the surrogate models. A Kriging meta-model $\hat{y}$ approximates a linear function $y$ as (Sacks et al. 1989; Jones et al. 1998):

$$y(\underline{x}) \approx \hat{y}(\underline{x}) = \mu + \mathcal{R}(\underline{x}), \quad (7)$$

where $\mu$ is the unknown process average and $\mathcal{R}$ is a stationary Gaussian process of zero mean of variance $\sigma^2$ which is to be determined. Its covariance function is given by $C(\underline{x}, \underline{x}') = \sigma^2 \mathcal{R}(\underline{x}, \underline{x}')$ where $\mathcal{R}(\underline{x}, \underline{x}')$ is the correlation function, given by the user. It is often chosen as a kernel function parametrized by a hyperparameter $\theta$ to be found. The Kriging prediction at a point $\underline{x}_0$ is given by:

$$\hat{y}(\underline{x}_0) = \mu + \underline{r}_0^T \underline{R}^{-1} (\underline{y} - \underline{1}\mu), \quad (8)$$

where $\underline{R}$ is the correlation matrix whose coefficients are $\underline{R} = C(\underline{x}^{(i)}, \underline{x}^{(i)}), \underline{r}_0$ as the vector of $C(\underline{x}_0, \underline{x}^{(i)})$ and $\underline{1}$ a vector of ones. Three parameters have to be determined, namely the process average $\mu$, the process variance $\sigma^2$ and the hyperparameter $\theta$. It is done by solving a likelihood optimisation problem. For more details about the theoretical background, the readers are invited to refer to Sacks et al. (1989), Paulson and Ragkousis (2015), Jones et al. (1998), Kleijnen (2017) and Cheng et al. (2020). In the present study, the Kriging meta-models are built based on the ordinary Kriging using pyKriging which is an open source Python Kriging toolkit (Paulson and Ragkousis 2015).
### 3.2.2 Construction and validation of the Kriging meta-models

A set of 100 points are generated for the construction of the Kriging meta-models. It means that 100 dNNMs analyses have been realized for 100 different geometries. Then, the modal amplitude is discretized and for each value, four different Kriging surrogate models are created: for the vibration amplitude $|q_{res,z}|$, for the resonant frequency $\omega_0$, for the modal damping ratio $\zeta$ and for the excitation forcing level $\gamma$. The number 100 has been chosen based on a convergence study, not shown here for the sake of consistency. This sample set includes 27 points located on the boundary of the three-dimensional design space and 77 points generated using a Latin hypercube sampling through the pyKriging package (Paulson and Ragkousis 2015).

An additional set of 20 geometries has been used to validate the results predicted from the Kriging meta-models. This comparison is done in Fig. 13a–d, where the predicted values are plotted versus the reference values from the numerical simulations (i.e. the points are expected to be on the line $y = x$). It is worth reminding that for each point, Kriging meta-models are constructed for each modal amplitude values, so it leads to a much larger number of points than 20. For $|q_{res,z}|$ and the $\gamma$, the predictions are very accurate as shown in Fig. 13a, d, where all the predictions are very close to the reference values (i.e. points are on the line with gradient of 1). However, for $\zeta$, the accuracy of the Kriging model is slightly lower at low $\zeta$ values (see Fig. 13c). By looking at Fig. 13b, the Kriging meta-models show less accurate predictions at larger $\omega_0$. For a deeper insight in the error committed, the predicted dNNMs and the reference dNNMs from the nonlinear modal analysis are compared in Fig. 13e–g. Three geometries with design parameters [−0.85, 0.65, 0.39], [−0.03, −0.41, 0.59] and [0.25, 0.45, 0.89] are selected from the validation set as examples to show the accuracy of the prediction. The predicted results are shown by the orange dashed curves whereas the reference results are represented by the black solid curves. It can be seen, the predictions are very close to the reference results, which gives confidence in the Kriging meta-models and their validation. In conclusion, these surrogate models can provide an accurate prediction of dNNMs with given design parameters.

## 3.3 Design space exploration and discussion

In this section, the surrogate models are exploited and analysed to assess the effect of each design parameter. This analysis is divided into two steps: first the influence of each parameter is studied when the two others are fixed to 0. Second, the whole design space is explored.

### 3.3.1 Influence of each design parameter

In this part, only one parameter can vary between −1 and 1 when the two others are fixed to 0. It allows to investigate the effect of each design parameter independently. Results are given in Fig. 14 for $\kappa_1$ (upper surface modification), in Fig. 15 for $\kappa_2$ and in Fig. 16 for $\kappa_3$. In each case, the evolutions of the $\omega_0-|q_{res,z}|$ plot, of the $\zeta-|q_{res,z}|$ plot and of the $\gamma-|q_{res,z}|$ plot and the considered parameter are displayed in a 3D view (the left column of the figures). For a better readability, two other views are also given: a side view that gives the envelopes of the dNNM properties versus the $|q_{res,z}|$ (the middle column in the figures) and the dNNM calculated from the default geometries [0, 0, 0] is highlighted by the black dashed curves, and a top view (right column).

- *Impact of* $\kappa_1$ from Fig. 14, one can see that large variations of the dNNM properties are observed when $\kappa_1$ varies. Indeed, the $\omega_0$ varies between 2.17 and 2.20 depending on the $|q_{res,z}|$. The $\omega_0$ tends to decrease with an increasing value of $\kappa_1$. As a modification of the upper surface leads to a variation of the number of points in contact, a strong variation for $\omega_0$ is expected as it modifies the stiffness of the global system. The $\zeta$ is also largely impacted by the variation of $\kappa_1$, and especially when $\kappa_1$ is positive. Indeed, the $\zeta$ tends to have similar evolution for negative $\kappa_1$ with one peak at the $|q_{res,z}|$ of 0.05 mm and a greater second one at 0.15 mm. But when $\kappa_1$ is positive, a large increase of the $\zeta$ for $|q_{res,z}|$ between 0.05 and 0.25 mm is observed (see Fig. 14d). The maximum $\zeta$ achieved is around 0.29% when $\kappa_1$ equals 1, with the $|q_{res,z}|$ of 0.13 mm. By looking at the $N_0$ distribution (see Case A in Fig. 11b), a more dispersive $N_0$ distribution is obtained when $\kappa_1 = 1$ and more nodes are in-contact. The evolution of the $\gamma$ follows the same trend as the $\zeta$ when $\kappa_1$ is varying. Since the total energy injected into the system by the excitation is expected to balance the total energy dissipated by the damping terms, including both the friction damping and viscous damping. Higher $\zeta$ leads to more dissipated energy due to the friction forces. Therefore, an excitation with a large forcing level is required to balance the energy dissipated due to the friction damping and viscous damping.

- *Impact of* $\kappa_2$ from Fig. 15, the $\omega_0$ at a low level of $|q_{res,z}|$ (also called as the linearized sticking frequency) decreases with an increasing value of $\kappa_2$. However, the $\omega_0$ at a large $|q_{res,z}|$ does not change with $\kappa_2$. The two outer bumps on the lower surface do not have a significant effect on the $\omega_0$. For the damping performances (see Fig. 15d–f), the maximum $\zeta$ of a value of 0.19% is achieved when $\kappa_2 = -1$, at the $|q_{res,z}|$ of 0.03 mm. By looking at the $N_0$ distribution (see Case D in Fig. 11e), a more dispersive $N_0$ distribution is obtained when $\kappa_2 = -1$
<figure><img src="image_16.png" /><figcaption>**Fig. 13** Meta-models validation: **a–d** comparison between the predicted values using Kriging against the reference values from 20 validation cases: **a** vibration amplitude $|q_{res,z}|$, **b** normalized resonant frequency $\omega_0$, **c** modal damping ratio $\zeta$, **d** excitation forcing</figcaption></figure>

level $\gamma$, and **e–g** comparison of the predicted results and the reference results using three validation cases as examples: **e** validation case 1 [−0.85, 0.65, 0.39], **f** validation case 2 [−0.03, −0.41, 0.59], **g** validation Case 3 [0.25, 0.45, 0.89]

and more nodes are in-contact with a relatively smaller $N_0$. The $\gamma$ follows a similar trend with the $\zeta$ as discussed in the previous point.

- *Impact of* $\kappa_3$ by looking at Fig. 16, the variation of $\kappa_3$ does not have a remarkable impact on the modal properties.

However, it is worth emphasizing that a lower peak can be noticed when $\kappa_3 = 0.3$ for both $\omega_0$ (at lower amplitude) and $\zeta$ (see Fig. 16c, f). In more details, the $\zeta - |q_{res,z}|$ curve shows two peaks around the $|q_{res,z}|$ of 0.03 mm (the first peak) and 0.15 mm (the second peak). An increasing
<figure><img src="image_17.png" /><figcaption>**Fig. 14** Effect of $κ_1 (-1 ≤ κ_1 ≤ 1)$ on the damped nonlinear normal mode [$κ_1$, 0, 0]: **a–c** normalized resonant frequency $ω_0$, **d–f** modal damping ratio $ζ$, and **g–i** excitation forcing level $γ$, black dashed curve: default geometry [0, 0, 0]. (Color figure online)</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Fig. 15** Effect of $κ_2 (-1 ≤ κ_2 ≤ 1)$ on the damped nonlinear normal mode [0, $κ_2$, 0]: **a–c** normalized resonant frequency $ω_0$, **d–f** modal damping ratio $ζ$, and **g–i** excitation forcing level $γ$, black dashed curve: default geometry [0, 0, 0]. (Color figure online)</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Fig. 16** Effect of $κ_3 (-1 ≤ κ_3 ≤ 1)$ on the damped nonlinear normal mode [0, 0, $κ_3$]: **a–c** normalized resonant frequency $ω_0$, **d–f** modal damping ratio $ζ$, and **g–i** excitation forcing level $γ$, black dashed curve: default geometry [0, 0, 0]. (Color figure online)</figcaption></figure>
value of $κ_3$ leads to a higher first peak and a lower second peak. The maximum $ζ$ of a value 0.19% is obtained when $κ_3 = -1$, with the $|q_{res,z}|$ of 0.17 mm.

In a nutshell, all nonlinear modal properties are very sensitive to the design parameter $κ_1$. For $κ_2$ and $κ_3$, the impact on $ω_0$ and $γ$ are negligible. In other words, damping performances can be significantly improved by varying $κ_1$, resulting in a remarkable change in $ω_0$. Whereas, by varying $κ_2$ and $κ_3$, the damping ability can be improved in a certain amount with a negligible variation of the $ω_0$. Furthermore, a dispersive $N_0$ distribution and an increasing number of contact nodes can lead to a larger $ζ$.

### 3.3.2 Whole design space exploration

Finally, the surrogate models are used to explore the whole design space. It is done by evaluating the different meta-models on a grid of 50 × 50 × 50 points. The dNNMs for all possible geometries are given in Fig. 17. The black dashed line represents the dNNM of the default geometry [0, 0, 0].

Compared to the previous cases (Figs. 14, 15, 16), larger variations of all nonlinear modal properties are obtained as shown in Fig. 17. Within the whole design space, the $ω_0$ falls in a range of 2.168–2.203 and in a percentage of –1.3% to 0.3% of the natural frequency. For $ζ$, the maximum is always reached at the $|q_{res,z}|$ that are smaller than 0.2 mm. The highest $ζ$ is around 0.35% and this geometry is defined as the optimal geometry in the following. The dNNM that corresponds to the optimal geometry is represented by the

dash line. The design parameters for this optimal geometry are [1, –1, –1]. More details about the optimal geometry are discussed in the next section.

### 3.3.3 Global sensitivity analysis

A global sensitivity analysis is achieved using a variance based approach using Sobol indices. The fundamental knowledge of Sobol indices is given in the “Appendix”. The Sobol indices are calculated using an open source software, OpenTURNS (Baudin et al. 2015). In the present study, the input parameters are the three design parameters, namely $κ_1$, $κ_2$ and $κ_3$. There are three output parameters, which are the resonant frequency $ω_0$, the modal damping ratio $ζ$ and the excitation forcing level $γ$. Four sets of Sobol indices are computed at different levels of the vibration amplitude, 0.05 mm, 0.1 mm, 0.2 mm and 0.4 mm (refer to Fig. 17).

The Sobol indices are computed within the whole design space and given in Fig. 18, where circles represent the first order Sobol indices and squares represent the total order indices. Orange markers show the Sobol indices of $ω_0$. Generally, the first order index of $κ_1$ for $ω_0$ almost reaches 1. Therefore, $κ_1$ has the most important influence on $ω_0$ regardless the level of amplitude. The interaction of $κ_1$ with the other input parameters can be ignored, since the total indices are very close to the first order index. This is in agreement with the observations done in Figs. 14, 15 and 16.

The blue markers describe the Sobol indices for the modal damping ratio $ζ$ and green for the excitation forcing level $γ$. As discussed in Sect. 3, the variation of $γ$ always shows a similar trend with $ζ$ (please refer the online

<figure><img src="image_20.png" /><figcaption>**Fig. 17** Design space exploration of the damped nonlinear normal mode: **a** normalized resonant frequency $ω_0$, **b** modal damping ratio $ζ$, and **c** excitation forcing level $γ$, black dashed curve: default geometry</figcaption></figure>

[0, 0, 0], red dashed curve: geometry with maximum modal damping ratio [1, –1, –1]. (Color figure online)
<figure><img src="image_21.png" /><figcaption>**Fig. 18** Global Sobol indices for different levels of vibration amplitude within the whole design space ($\kappa_{1,2,3} \in [-1, 1]$): **a** 0.05 mm, **b** 0.1 mm, **c** 0.2 mm, and **d** 0.4 mm. (Color figure online)</figcaption></figure>

version for the description of the color in Fig. 18). By looking at the indices for $\zeta$ (blue), it is easy to notice that the contributions for each input parameter are distinguishable for different levels of the vibrational amplitude. At a low amplitude (see Fig. 18a), the total order indices for the three inputs are 0.46, 0.39 and 0.33 and the first order indices are 0.41, 0.27 and 0.21 respectively. It shows that the three parameters have a similar influence on $\zeta$ and $\gamma$ as the Sobol indices are similar. Moreover, for $\kappa_2$ and $\kappa_3$, the difference between the first order and total order Sobol indices demonstrates that small coupling effects between the different parameters are present. At higher vibration amplitudes, the contribution of $\kappa_1$ becomes predominant and the influence of $\kappa_2$ and $\kappa_3$ becomes negligible on the modal damping ratio $\zeta$ and the excitation forcing level $\gamma$.

Therefore, from the global sensitivity analysis, one can come up with a few points. First, the direct effect of $\kappa_1$ is the most important factor for the variance three targeting outputs. The interaction of $\kappa_1$ with the two other inputs is almost negligible. $\kappa_2$ and $\kappa_3$ have a very limited effect on $\omega_0$ regardless the vibration amplitude. For a low amplitude, the influence of the three input parameters on $\zeta$ and $\gamma$ are almost similar. For a higher amplitude, the direct contribution of $\kappa_2$ and $\kappa_3$ decreases and interaction effect vanished.

## 3.4 The optimal geometry

From the exploration of the whole design space using the surrogate models, the optimal geometry of the FRD with a set of design parameters $[1, -1, -1]$ is identified. In this section, the results of the optimal geometry, including its geometry, its $N_0$ and $\Delta$ distributions, its dNNMs and the nonlinear FRFs are explained in detail. Finally, the damping performances of the optimal geometry is compared to the default geometry.

The optimal geometry shows an upper bump and three lower bumps, as shown in Fig. 19a. Compared to the other cases, it has the largest mass and the lowest flexibility. This finding seems to be contradictory with the results in Laxalde et al. (2010, 2007) and Sun et al. (2018). However, a major difference exist in terms of modelling compared to Laxalde et al. (2007, 2010) as the ring damper was modelled only by punctual masses or by Euler beams, whereas here a full FE model of the ring damper. This allows to take into consideration the real effect of the damper geometry. Moreover, in Laxalde et al. (2010, 2007) and Sun et al. (2018) the contact normal load was assumed to be uniform, whereas here the initial contact normal load distribution is obtained from a nonlinear static analysis. Results of the present work have demonstrated how a change in the contact normal load
<figure><img src="image_22.png" /><figcaption>**Fig. 19** Results for the optimal geometry: **a** optimal geometry, **b** contact normal load $N_0$ distribution, **c** contact gap $\Delta$ distribution, **d** normalised resonant frequency $\omega_0$ against amplitude $|q_{res,z}|$, **e** modal damping ratio $\zeta$ against amplitude $|q_{res,z}|$, **f** excitation forcing level $\gamma$ against amplitude $|q_{res,z}|$, **g** black curves for the default geometry and red curves for the optimal geometry, **h** numerical validation using several nonlinear frequency response function, and **i** energy dissipation $E_d$ distribution for the resonant solutions. (Color figure online)</figcaption></figure>

distribution impacts the damping properties, so a difference in terms of results is expected.

By comparing the initial pressure $N_0$ and initial gap $\Delta$ distributions between the optimal geometry (Fig. 19b) and the default geometry (Fig. 11a), one can see that more nodes are in contact and that the normal load is also higher. Apart from the effect on $N_0$, the optimal geometry is less flexible in the Y direction leading to less separation (smaller gap and less separated contact pairs) within the contact interface. Compared to the default geometry, the resonant frequency $\omega_0$ of the optimized one is reduced from 2.195 to 2.182 (0.6% of reduction). For a given value of $\gamma$, the optimal geometry shows always the lowest $|q_{res,z}|$ compared to the whole design space, but also has the lowest $\omega_0$ when the amplitude is superior to 0.05 mm. Moreover, this geometry always has the highest $\zeta$ when $|q_{res,z}|$ is superior to 0.07 mm. Finally, several nonlinear FRFs are computed for with both FRD geometries for different excitation level and displayed in

Fig. 19g, the predicted resonant solutions show an excellent agreement with the ones obtained from the nonlinear FRFs.

As discussed in the design consideration in Sect. 3.1, the initial idea of the geometric design of the FRD is to modify the initial contact status to achieve better damping performances. From Fig. 9b, one can noticed that most energy dissipation comes from the nodes in the middle under the blade. If the dissipative ability of the other contact nodes (nodes near the boundaries) can be improved, then the overall damping performance of the FRD can be enhanced. By looking the energy dissipation distributions of the optimal geometry (see Fig. 19h), it is obvious that the nodes near the boundaries can dissipate more energy in the case of the optimized geometry, which validates the initial motivation of the geometric design for the FRD.

To further demonstrate the improvement of the optimal geometry, the values of the resonant amplitude predicted using the dNNMs with the E-EBM are compared. The results are shown in Table 3. The predicted resonant
<table><caption>**Table 3** Comparison between the optimal geometry and the default geometry</caption><thead><tr><th></th><th>Resonant amplitude for the default geometry (mm)</th><th>Resonant amplitude for the optimal geometry (mm)</th><th>Reduction (%)</th></tr></thead><tbody><tr><td>Excitation with 20 N</td><td>0.069</td><td>0.049</td><td>29.0</td></tr><tr><td>Excitation with 60 N</td><td>0.191</td><td>0.114</td><td>40.3</td></tr><tr><td>Excitation with 80 N</td><td>0.301</td><td>0.145</td><td>51.8</td></tr><tr><td>Excitation with 101 N</td><td>0.439</td><td>0.215</td><td>51.0</td></tr></tbody></table>

amplitudes for both optimal and default geometry are listed in the table, as well as the reduction percentage. The nonlinear FRF for the default geometry are shown together with the optimal geometry in Fig. 19g. From Table 3 and Fig. 19g, one can notice that the resonant amplitude of the optimal geometry is always lower than the default one under a same level of the excitation force. The maximum reduction in the resonant amplitude is around 51.8%.

In summary, thanks to such an efficient numerical method (use of the dNNMs, the E-EBM and Kriging meta-models), a geometric design of the FRD has been explored. The optimal geometry identified through the design process shows a better damping performance than the default one. The maximum damping ratio is improved around 200% (default 0.12% and optimal 0.36%). The resonant amplitude can be reduced more than 50%.

### 3.4.1 Sensitivity analysis around the local optimal geometry

The optimal geometry is identified as shown in Fig. 19. To investigate the sensitivity around the optimal geometry and to investigate its robustness, the Sobol indices are calculated in a similar way. The three input parameters for the optimal design are [$\kappa_1, \kappa_2, \kappa_3$] = [1, -1, -1]. Therefore, the range for each input parameter is chosen as: $\kappa_1 \in [0.9, 1]$, $\kappa_2 \in [-1, -0.9]$ and $\kappa_3 \in [-1, -0.9]$. The local Sobol indices are shown in Fig. 20. The local Sobol indices do not show significant difference compared to the global ones. In brief, the resonant frequency $\omega_0$ is mostly driven by $\kappa_1$. Considering the modal damping ratio $\zeta$ and the excitation forcing level $\gamma$, at low amplitude, less coupling between the parameters are observed with the optimal geometry. Indeed, the first order and total order indices are equal when they

<figure><img src="image_23.png" /><figcaption>**Fig. 20** Local Sobol indices for different levels of vibration amplitude around the optimal geometry ($\kappa_{1,2,3} \in [-1, 1]$): **a** 0.05 mm, **b** 0.1 mm, **c** 0.2 mm, and **d** 0.4 mm. (Color figure online)</figcaption></figure>

<figure><img src="image_24.png" /></figure>
were slightly different on the full space. For larger amplitudes, results are similar as the Sobol indices are very low.

As a conclusion, the upper surface of the FRD has a high influence on the damper characteristics, whereas the lower surface geometry has a very limited impact. But if the damper is aimed to be used at low amplitude level, $\zeta$ can be tuned with $\kappa_2$ and $\kappa_3$. The coupling effect between $\kappa_2$ and $\kappa_3$ cannot be ignored at this case. Moreover, it appears that the three parameters do not have significant coupling effect at large amplitude, which means that they can be treated independently in a design process.

## 4 Conclusion

A comprehensive investigation on the impact of the geometry of a FRD on its damping performances was conducted in this work. A three-dimensional FE model of a blisk was used as a test case. A linear modal analysis, a rotor dynamic analysis and a nonlinear static analysis are done in the first place to identify the resonant case, the contact normal load and gap distributions. Then, the dNNM coupled with E-EBM is applied to calculate efficiently the nonlinear modal responses of the blisk with FRD.

The ring damper geometry is parametrized in order to modify the lower and upper surfaces of the latter. To study in depth the impact of the modification of the geometry, Kriging meta-models are constructed and used to predict with a low numerical cost the nonlinear modal response of the blisk with different geometries of the ring dampers. The parametric study is completed with a sensitivity analysis based on Sobol indices. The following main results are obtained:

- The modification of the ring damper geometry has a major effect on the initial contact condition at the friction interface and on the structural properties (mass and flexibility) of the ring damper. The modification of these properties impact substantially the resonant frequency and the modal damping ratio.

- The shape of the upper surface has the largest impact on the nonlinear modal properties, whereas the impact of the lower surface on the resonant frequency and on the modal damping ratio are limited or negligible. By changing the damper geometry, the contact normal load distribution over the contact surface is modified. It is shown that a better energy dissipation is reached (an so better damping performances) when the contact normal load is less localised.

- An optimum ring damper geometry is identified. The maximum modal damping ratio can be improved by about 200% and the vibration amplitude at resonance can be reduced more than 50%.

This work demonstrates how the design of the geometry of friction dampers can increase substantially the damping performances. The numerical framework presented here based on the dNNMs, extend energy balance method and Kriging meta-models makes possible to consider the design of realistic dampers modelled with industrial FE models with an acceptable computational cost. For future work, shape and/or topological optimisation of such components should be considered to consider more possible geometries. Finally, an experimental validation is expected in a near future to validate the optimal FRD.

## Appendix: Sobol indices

Sobol indices can be considered as a type of the variance based sensitivity analysis (Sobol 1993; Sobol et al. 2007). This method is a form of global sensitivity analysis and can be used for nonlinear dynamic problems. The sensitivity indices are obtained by the decomposition of the output variance into fractions due to a single input parameter or an interaction effect of more than one input parameters. If a vector of $\mathbf{\kappa} = [\kappa_1, \kappa_2, \kappa_3]$ denote the three input parameters, and $y$ will denote one of the considered outputs. The output can be seen as a function of inputs $y = f(\mathbf{\kappa})$. Sobol indices are used to investigate the direct influence of the different input parameters $\kappa_i$, as well as their coupling effects on the output variance. Therefore, the output variance $V(y)$ can be decomposed in:

$$V(y) = \sum_{i=1}^{3} V_i + \sum_{i<j}^{3} V_{i,j} + V_{1,2,3}, \qquad (A.1)$$

where $V$ denotes the variance and $E$ represents the expectation (Saltelli et al. 1999, 2002), and where the other terms are defined as:

$$ \begin{gathered} V_i = V(E(y|\kappa_i)), \\ V_{i,j} = V(E(y|\kappa_i, \kappa_j)) - V_i - V_j, \\ V_{1,2,3} = V(E(y|\kappa_1, \kappa_2, \kappa_3)) - \sum_{i=1}^{3} V_i - \sum_{i<j}^{3} V_{i,j}. \end{gathered} \qquad (A.2) $$

Then, the first three orders Sobol indices are defined as:

$$ \begin{gathered} S_i = V_i/V(y), \\ S_{i,j} = V_{i,j}/V(y), \\ S_{1,2,3} = V_{1,2,3}/V(y). \end{gathered} \qquad (A.3) $$

The total order of the Sobol indices $ST_i$ can be calculated as below where $V(E(y| - \kappa_i))$ is the total contribution to the variance of $y$ due to non-$\kappa_i$.

$$ST_i = 1 - V(E(y| - \kappa_i))/V(y). \qquad (A.4)$$
**Acknowledgements** We thank Professor David Nowell (Imperial College London) and Professor Bogdan Epureanu (University of Michigan, Ann Arbor) for many constructive suggestions for the present work.

**Funding** Dr. Yekai Sun is grateful to China Scholarship Council (File No. 201708060239) for providing the financial support for this project.

## Declarations

**Conflict of interest** The authors declare that they have no conflict of interest.

**Replication of results** The results presented in this work can be replicated by implementing the equations and the finite element model presented in this paper. All relevant equations have been included to enable readers to replicate the results. The mesh of the finite element model will be made available by request.

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

## References

* Allgower EL, Georg K (2003) Introduction to numerical continuation methods. In: Classics in applied mathematics, vol 45. SIAM, Philadelphia, p xxvi-388

* Baudin M, Dufoy A, Iooss B, Popelin AL (2015) OpenTURNS: an industrial software for uncertainty quantification in simulation

* Brake MRW, Groß J, Lacayo RM, Salles L, Schwingshackl CW, Reuß P, Armand J (2018) Reduced order modeling of nonlinear structures with frictional interfaces. Springer, Cham, pp 427-450

* Cameron TM, Griffin JH (1989) An alternating frequency/time domain method for calculating the steady-state response of nonlinear dynamic systems. J Appl Mech 56(1):149-154

* Cheng K, Lu Z, Ling C, Zhou S (2020) Surrogate-assisted global sensitivity analysis: an overview. Struct Multidisc Optim 61(3):1187-1213

* Craig R, Bampton M (1968) Coupling of substructures for dynamic analyses. AIAA J 6(7):1313-1319

* Denimal E, El Haddad F, Wong C, Salles L (2021a) Topological optimization of under-platform dampers with moving morphable components and global optimization algorithm for nonlinear frequency response. J Eng Gas Turbines Power 143(2):021021

* Denimal E, Wong C, Salles L, Pesaresi L (2021b) On the efficiency of a conical underplatform damper for turbines. J Eng Gas Turbines Power 143(2):021020

* Denimal E, Renson L, Salles L (2021c) Topological optimisation of friction dampers for nonlinear resonances mitigation. In: Nodycon 2021

* Gastaldi C, Berruti TM, Gola MM (2018) Best practices for underplatform damper designers. Proc Inst Mech Eng C 232(7):1221-1235

* Grenat C, Baguet S, Dufour R, Lamarque CH (2018) Generalization of conservative nonlinear normal modes to the non-conservative equation of motion. Paper DETC2018-85799. ASME

* Hill T, Cammarano A, Neild S, Wagg D (2015) Interpreting the forced responses of a two-degree-of-freedom nonlinear oscillator using backbone curves. J Sound Vib 349:276-288

* Hüls M, Panning-von Scheidt L, Wallaschek J (2018) Influence of geometric design parameters onto vibratory response and high-cycle fatigue safety for turbine blades with friction damper. J Eng Gas Turbines Power 141(4):041022

* Jones DR, Schonlau M, Welch WJ (1998) Efficient global optimization of expensive black-box functions. J Glob Optim 13(4):455-492

* Kleijnen JP (2017) Regression and Kriging metamodels with their experimental designs in simulation: a review. Eur J Oper Res 256(1):1-16

* Krack M (2015) Nonlinear modal analysis of nonconservative systems: extension of the periodic motion concept. Comput Struct 154:59-71

* Krack M, Gross J (2019) Theory of harmonic balance. Springer, Cham, pp 11-46

* Krack M, Tatzko S, Panning-von Scheidt L, Wallaschek J (2014) Reliability optimization of friction-damped systems using nonlinear modes. J Sound Vib 333(13):2699-2712

* Krack M, Salles L, Thouverez F (2016) Vibration prediction of bladed disks coupled by friction joints. Arch Comput Methods Eng 24:589-636

* Laxalde D, Thouverez F (2009) Complex non-linear modal analysis for mechanical systems: application to turbomachinery bladings with friction interfaces. J Sound Vib 322(4):1009-1025

* Laxalde D, Thouverez F, Sinou JJ, Lombard JP (2007) Qualitative analysis of forced response of blisks with friction ring dampers. Eur J Mech A 26(4):676-687

* Laxalde D, Gibert C, Thouverez F (2008) Experimental and numerical investigations of friction rings damping of blisks. Paper No. GT2008-50862. ASME

* Laxalde D, Thouverez F, Lombard JP (2010) Forced response analysis of integrally bladed disks with friction ring dampers. J Vib Acoust 132(1):011013

* Lupini A, Mitra M, Epureanu BI (2019) Application of tuned vibration absorber concept to blisk ring dampers: a nonlinear study. J Eng Gas Turbines Power 141(10):101016

* Mead DJ (1975) Wave propagation and natural modes in periodic systems: I. Mono-coupled systems. J Sound Vib 40(1):1-18

* Mijar A, Arora J (2000) Review of formulations for elastostatic frictional contact problems. Struct Multidisc Optim 20(3):167-189

* Niemotka MA, Ziegert JC (1995) Optimal design of split ring dampers for gas turbine engines. J Eng Gas Turbine Power 117(3):569-575

* Panning L, Sextro W, Popp K (2002) Optimization of the contact geometry between turbine blades and underplatform dampers with respect to friction damping. Paper No. GT-2002-30429. ASME

* Paulson C, Ragkousis G (2015) pyKriging: a Python Kriging toolkit (July). https://doi.org/10.5281/zenodo.21389

* Pesaresi L, Salles L, Jones A, Green J, Schwingshackl C (2017) Modeling the nonlinear behaviour of an underplatform damper test rig for turbine applications. Mech Syst Signal Process 85:662-679

* Petrov EP (2004) A method for use of cyclic symmetry properties in analysis of nonlinear multiharmonic vibrations of bladed disks. J Turbomach 126(1):175-183

* Pollini N, Lavan Oren L, Amir O (2018) Adjoint sensitivity analysis and optimization of hysteretic dynamic systems with nonlinear viscous dampers. Struct Multidisc Optim 57(6):2273-2289
Rosenberg RM (1960) Normal modes of nonlinear dual-mode systems. J Appl Mech 27(2):263-268

Ruffini V (2016) Coriolis effects in bladed discs. PhD Thesis, Imperial College, London

Sacks J, Welch WJ, Mitchell TJ, Wynn HP (1989) Design and analysis of computer experiments. Stat Sci 4(4):409-423

Salles L, Blanc L, Thouvez F, Gouskov A (2011) Dynamic analysis of fretting-wear in friction contact interfaces. Int J Solids Struct 48(10):1513-1524

Saltelli A, Tarantola S, Chan K (1999) A quantitative, model independent method for global sensitivity analysis of model output. Technometrics 41(1):39-56

Saltelli A, Tarantola S, Chan K (2002) Making best use of model evaluations to compute sensitivity indices. Comput Phys Commun 145:280-297

Sanliturk KY, Ewins DJ, Elliott R, Green JS (2001) Friction damper optimization: simulation of rainbow tests. J Eng Gas Turbines Power 123(4):930-939

Sarrouy E (2019) Phase driven modal synthesis for forced response evaluation. In: 7th International conference on nonlinear vibrations, localization and energy transfer, Marseille, France, July 2019

Schwarz S, Kohlmann L, Hartung A, Gross J, Scheel M, Krack M (2020) Validation of a turbine blade component test with frictional contacts by phase-locked-loop and force-controlled measurements. J Eng Gas Turbines Power 142(5):051006

Shaw SW, Pierre C (1991) Non-linear normal modes and invariant manifolds. J Sound Vib 150(1):170-173

Smola AJ, Schölkopf B (2004) A tutorial on support vector regression. Stat Comput 14(3):199-222

Sobol IM (1993) Sensitivity analysis for nonlinear mathematical model. Math Models Comput Exp 1(4):407-414

Sobol IM, Tarantola S, Gatelli D, Kucherenko SS, Mauntz W (2007) Estimating the approximation errors when fixing unessential factors in global sensitivity analysis. Reliab Eng Syst Saf 92:957-960

Sun Y, Yuan J, Pesaresi L, Salles L (2018) Nonlinear vibrational analysis for integrally bladed disk using frictional ring damper. J Phys Conf Ser 1106:012026

Sun Y, Yuan J, Pesaresi L, Denimal E, Salles L (2020) Parametric study and uncertainty quantification of the nonlinear modal properties of frictional dampers. J Vib Acoust 142(5):051102

Sun Y, Vizzaccaro A, Yuan J, Salles L (2021a) An extended energy balance method for resonance prediction in forced response of systems with non-conservative nonlinearities using damped nonlinear normal mode. Nonlinear Dyn 103(4):3315-3333

Sun Y, Yuan J, Denimal E, Salles L (2021b) Nonlinear modal analysis of frictional ring damper for compressor blisk. J Eng Gas Turbines Power 143(3):031008

Sun Y, Yuan J, Denimal E, Salles L (2021c) A study of the contact interface for compressor blisks with ring dampers using nonlinear modal analysis. IOP Conf Ser Mater Sci Eng 1081:012041

Sun Y, Yuan J, Vizzaccaro A, Salles L (2021d) Comparison of different methodologies for the computation of damped nonlinear normal modes and resonance prediction of systems with non-conservative nonlinearities. Nonlinear Dyn 104:3077-3107

Tang W, Epureanu BI (2017) Nonlinear dynamics of mistuned bladed disks with ring dampers. Int J Nonlinear Mech 97:30-40

Tang W, Epureanu BI (2019) Geometric optimization of dry friction ring dampers. Int J Nonlinear Mech 109:40-49

Totis G, Sortino M (2020) Polynomial chaos-Kriging approaches for an efficient probabilistic chatter prediction in milling. Int J Mach Tools Manuf 157:103610

Tran DM (2009) Component mode synthesis methods using partial interface modes: application to tuned and mistuned structures with cyclic symmetry. Comput Struct 87(17):1141-1153

Vakakis AF (2001) Inducing passive nonlinear energy sinks in vibrating systems. J Vib Acoust 123(3):324-332

Wiener N (1938) The homogeneous chaos. Am J Math 60(4):897-936

Wu Y, Li L, Fan Y, Zucca S, Gastaldi C, Ma H (2020) Design of dry friction and piezoelectric hybrid ring dampers for integrally bladed disks based on complex nonlinear modes. Comput Struct 233:106237

Xiu D, Karniadakis GE (2002) The Wiener-Askey polynomial chaos for stochastic differential equations. SIAM J Sci Comput 24(2):619-644

Yang BD, Menq CH (1998) Characterization of 3D contact kinematics and prediction of resonant response of structures having 3D frictional constraint. J Sound Vib 217(5):909-925

Ye H, Wang Y, Jiang X (2017) A method for the optimal design of split ring dampers for aviation gears. Paper DETC2017-67231. ASME

Yuan Y, Jones A, Setchfield R, Schwingshackl C (2020) Robust design optimisation of underplatform dampers for turbine applications using a surrogate model. J Sound Vib 494(6):115528

Yuan J, Fantetti A, Denimal E, Bhatnagar S, Pesaresi L, Schwingshackl C, Salles L (2021) Propagation of friction parameter uncertainties in the nonlinear dynamic response of turbine blades with underplatform dampers. Mech Syst Signal Process 156:107673

Yuan J, Sun Y, Schwingshackl C, Salles L (2022) Computation of damped nonlinear normal modes for large scale nonlinear systems in a self-adaptive modal subspace. Mech Syst Signal Process 162:108082

Zucca S, Firrone CM, Facchini M (2010) A method for the design of ring dampers for gears in aeronautical applications. Paper IMECE2010-38788. ASME

Zucca S, Firrone CM, Gola MM (2012) Numerical assessment of friction damping at turbine blade root joints by simultaneous calculation of the static and dynamic contact loads. Nonlinear Dyn 67(3):1943-1955

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.