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

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.


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. 98 Page 2 of 25 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 , 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 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 opitmisation 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 Thouverez 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 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 nonlinar 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 modelling 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 introduced 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.

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 nonlinear 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.

Finite element model
A 3D FE model constructed with the commercial software 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 elements 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 material is expected for the friction damper compared with the 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 ( r-t-Z) is used for a sector instead of the global Cartesian coordinate system (X-Y-Z, please see Fig. 1). r-Axis is the radial direction, whereas t -axis represents the circumferential direction.

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), 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 and the resonant frequency 0 are obtained (see Sect. 2.5 for details). This nonlinear dynamic analysis is achieved using an in-house Python script.

Rotor dynamic analysis and determination of the resonance
In an real operating situation, the blisk vibrates in the form of travelling wave mode under the rotational framework.
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 nonuniform 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 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 0 = 2.195 is excited. This resonance case is marked by a star in Fig. 5 and used for the following analysis.

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 ( ) distributions within the contact interface are extracted from the nonlinear static analysis. These initial contact N 0 and 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 distribution are provided in Fig. 6. From the results, the nodes in the centre have a greater normal load, whereas nodes near boundaries show a larger gap.

Nonlinear dynamic analysis
To determine the nonlinear vibration amplitude at resonance and for a given forcing level , 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.

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 and the reduced stiffness matrix of the blisk and FRD.

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 normal contact stiffness; k t1 and k t2 are the tangential contact stiffness. A prescribed gap 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 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).

Nonlinear modal analysis
The dynamic equation of the bladed disk under forcing is defined by Eq. 1, where is the reduced mass matrix, is the reduced stiffness matrix, is the structural reduced damping matrix (from material dissipation for example), F c is the vector of the friction force within the contact interface and F e is the excitation force with the forcing amplitude, the absolute phase of the excitation force and 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.
On the other side, the dNNMs of the bladed disk are defined by Eq. 2 where 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).
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. (1) (2) Q + Q + F c (Q,Q) = 0.

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 , is introduced into the system to quantify the modal energy. Then, modal displacements are expressed as: is the mass normalized modal displacement. Furthermore, the absolute phase of the mass normalized solution 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 nonconservative 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 ̃ is expressed where is the modal damping ratio and 0 is the resonant frequency. The periodic solutions of the nonconservative 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.
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 Q 0 , the resonant frequency 0 and the modal damping ratio can be easily computed for a given modal amplitude (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.

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 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 nonlinear FRF at resonance, so when the excitation frequency equals the resonant frequency 0 , the resonant shared solution is given in Eq. 4.
For each level of the modal amplitude , the resonant shared solution , the resonant frequency 0 and the modal damping ratio are known from the computation of the dNNMs, as presented in previously. From the E-EBM, if the viscous damping term and the vector of the excitation force are given are given, then the excitation forcing level and the absolute phase of the excitation force can be determined. In other words, the E-EBM makes possible to find the correspondence between the dNNMs and a modal amplitude , and the corresponding excitation level and forcing phase . Recalling the Eqs.
(1)-(2) and assuming that the solution of the nonlinear FRF equals the resonant shared solution , 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 ( , ) 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 , whereas the energy exchanged by the external excitation E f ( , ) varies with the forcing phase and the excitation forcing level . The values of and are determined by finding the single intersection between the curve of E f ( , ) 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.
After the application of the E-EBM, one gets the values of and that satisfy E f ( , ) = E damp . From this, it is possible to get directly the solution 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.

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 0 , and 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 0 and a zero modal damping ratio . The 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

Fig. 8
Damped nonlinear normal modes using nonlinear modal analysis: a normalised resonant frequency 0 against amplitude |q res, z | , b modal damping ratio against amplitude |q res, z | , and c excitation forcing level against amplitude |q res, z | 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 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 is obtained. The maximum occurs with a value of 0.15% with a |q res, z | of 0.164 mm and a 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 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 -|q res, z | plot obtained from the E-EBM and the right part presents the normalised resonant frequency 0 -|q res, z | plot obtained from the computation of the FRF. The nonlinear FRF with different values of are computed ( = 20 N, 60 N, 80 N, and 101 N), and for each the vibration amplitude at resonance is reported on the -|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 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 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 0 , 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.

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 , the resonant frequency 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 metamodelling, and more precisely on Kriging, is employed to investigate the geometric effects of the FRD on the dynamic 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.

Design parameters
The coordinate in radial direction of the FRD can vary (in the r-axis), but the coordinates in the rotational axis (in the Z-axis) and the circumferential direction (in the t-axis) are unchanged (see Figs. 4, 11 for axis names). The geometry of the FRD is parametrized using three design parameters, 1 , 2 and 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 other contact nodes (the 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: -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 1 > 0 (see Case A Fig. 11b), whereas a negative 1 leads to a concave surface (Case B Fig. 11c). When 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. -2 Controls the two outer bumps of the lower surface, similarly to 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. -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 2 .
When one of the parameters i with i ∈ {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 r for each normalized parameters are given in Table 2. The physical values r represent the change of the radial coordinate ( r-axis), leading to a bump (positive r ) or a concave (negative 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 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 each node for the FRD are automatically updated according to the defined upper and lower surface.

Effects on the contact status
For each geometry, before computing the systems dynamic response, the N 0 and 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 for the seven geometries are given as an example in Fig. 11 in the middle and right columns, respectively. When 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 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 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 2 is −1 . The highest N 0 is observed when 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.

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 0 , the modal damping ratio and the predicted excitation forcing level 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 0 -|q res, z | plots, -|q res, z | plots and -|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 0 , the and the are very sensitive to 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 1 (Case A) leads to a more dispersive N 0 distribution and less separated contact nodes. In nonlinear modal analysis, a lower 0 and a better damping performance are observed. Indeed, the 0 equals 2.185 at a low level of |q res, z | and the reaches 0.28% . Moreover, for the same level of , the |q res, z | are always lower for Case A. Less variations are observed for 2 and 3 , especially on the 0 and the |q res, z | . However, large variation can be observed on the . Conclusively, changing 1 can improve the damping performance greatly with a significant shift of 0 and allows for large variations of the dynamic response. Whereas, different 2 and 3 can ensure a better damping performance with a negligible shift in 0 .

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 1 , 2 and 3 . For each modal amplitude , four Kriging meta-models are built: one for the vibration amplitude |q res, z | , one for the resonant frequency 0 , one for the modal damping ratio and one for the excitation forcing level . In the following, = [ 1 , 2 , 3 ] will denote the vector of inputs, and y will denote one of the considered output. A set ( (i) ), y (i) of N points are used to train the surrogate models. A Kriging meta-model ŷ approximates a linear function y as (Sacks et al. 1989;Jones et al. 1998): where is the unknown process average and Z is a stationary Gaussian process of zero mean of variance 2 which is to be determined. Its covariance function is given by C( , � ) = 2 R( , � ) where R( , � ) is the correlation function, given by the user. It is often chosen as a kernel function parametrized by a hyperparameter to be found. The Kriging prediction at a point 0 is given by: where is the correlation matrix whose coefficients are ij = C( (i) , (j) ) , 0 as the vector of C( 0 , (i) ) and a vector of ones. Three parameters have to be determined, namely the process average , the process variance 2 and the hyperparameter . 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 build based on the ordinary Kriging using pyKriging which is an open source Python Kriging toolkit (Paulson and Ragkousis 2015). Fig. 12 Damped nonlinear normal modes for all seven geometries (Fig. 11): a normalised resonant frequency 0 against amplitude |q res, z | , b modal damping ratio against amplitude |q res, z | , and c excitation forcing level against amplitude |q res, z |

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 0 , for the modal damping ratio and for the excitation forcing level . 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 , 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 , the accuracy of the Kriging model is slightly lower at low values (see Fig. 13c). By looking at Fig. 13b, the Kriging meta-models show less accurate predictions at larger 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 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.

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.

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 1 (upper surface modification), in Fig. 15 for 2 and in Fig. 16 for 3 . In each case, the evolutions of the 0 -|q res, z | plot, of the -|q res, z | plot and of the -|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 1 from Fig. 14, one can see that large variations of the dNNM properties are observed when 1 varies. Indeed, the 0 varies between 2.17 and 2.20 depending on the |q res, z | . The 0 tends to decrease with an increasing value of 1 . As a modification of the upper surface leads to a variation of the number of points in contact, a strong variation for 0 is expected as it modifies the stiffness of the global system. The is also largely impacted by the variation of 1 , and especially when 1 is positive. Indeed, the tends to have similar evolution for negative 1 with one peak at the |q res, z | of 0.05 mm and a greater second one at 0.15 mm. But when 1 is positive, a large increase of the for |q res, z | between 0.05 and 0.25 mm is observed (see Fig. 14d). The maximum achieved is around 0.29% when 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 1 = 1 and more nodes are in-contact. The evolution of the follows the same trend as the when 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 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 2 from Fig. 15, the 0 at a low level of |q res, z | (also called as the linearized sticking frequency) decreases with an increasing value of 2 . However, the 0 at a large |q res, z | does not change with 2 . The two outer bumps on the lower surface do not have a significant effect on the 0 . For the damping performances (see Fig. 15d-f), the maximum of a value of 0.19% is achieved when 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 2 = −1 and more nodes are in-contact with a relatively smaller N 0 . The follows a similar trend with the as discussed in the previous point. -Impact of 3 by looking at Fig. 16, the variation of 3 does not have a remarkable impact on the modal properties.
However, it is worth emphasizing that a lower peak can be noticed when 3 = 0.3 for both 0 (at lower amplitude) and (see Fig. 16c, f). In more details, the -|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 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 .

Whole design space exploration
Finally, the surrogate models are used to explore the whole design space. It is done by evaluating the different metamodels 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.

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 version for the description of the color in Fig. 18). By looking at the indices for (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 and as the Sobol indices are similar. Moreover, for 2 and 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 1 becomes predominant and the influence of 2 and 3 becomes negligible on the modal damping ration and the excitation forcing level .
Therefore, from the global sensitivity analysis, one can come up with a few points. First, the direct effect of 1 is the most important factor for the variance three targeting outputs. The interaction of 1 with the two other inputs is almost negligible. 2 and 3 have a very limited effect on 0 regardless the vibration amplitude. For a low amplitude, the influence of the three input parameters on and are almost similar. For a higher amplitude, the direct contribution of 2 and 3 decreases and interaction effect vanished.

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 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. (2010Laxalde et al. ( , 2007 and Sun et al. (2018). However, a major difference exist in terms of modelling compared to Laxalde et al. (2007Laxalde et al. ( , 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. (2010Laxalde et al. ( , 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 distribution impacts the damping properties, so a difference in terms of results is expected.
By comparing the initial pressure N 0 and initial gap 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 0 of the optimized one is reduced from 2.195 to 2.182 ( 0.6% of reduction). For a given value of , the optimal geometry shows always the lowest |q res, z | compared to the whole design space, but also has the lowest 0 when the amplitude is superior to 0.05 mm. Moreover, this geometry always has the highest 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 Fig. 19 Results for the optimal geometry: a optimal geometry, b contact normal load N 0 distribution, c contact gap distribution, d normalised resonant frequency 0 against amplitude |q res, z | , e modal damping ratio against amplitude |q res, z | , f excitation forcing level against amplitude |q res, z | , d-f black curves for the default geometry and red curves for the optimal geometry, g numerical validation using several nonlinear frequency response function, and h energy dissipation E d distribution for the resonant solutions. (Color figure online) 98 Page 22 of 25 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%.

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 [ 1 , 2 , 3 ] = [1, −1, −1] . Therefore, the range for each input parameter is chosen as: 1 ∈ [0.9, 1]) , 2 ∈ [−1, −0.9]) and 3 ∈ [−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 0 is mostly driven by 1 . Considering the modal damping ratio and the excitation forcing level , 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  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, can be tuned with 2 and 3 . The coupling effect between 2 and 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.

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 an single input parameter or an interaction effect of more than one input parameters. If a vector of = [ 1 , 2 , 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 ( ) . Sobol indices are used to investigate the direct influence of the different input parameters i , as well as their coupling effects on the output variance. Therefore, the output variance V(y) can be decomposed in: where V denotes the variance and E represents the expectation (Saltelli et al. 1999(Saltelli et al. , 2002, and where the other terms are defined as: Then, the first three orders Sobol indices are defined as: The total order of the Sobol indices ST i can be calculated as below where V(E(y| − i )) is the total contribution to the variance of y due to noni .