Stress recovery algorithm for reduced order models of mechanical systems in nonlinear dynamic operative conditions

Nonlinear forced response analyses of mechanical systems in the presence of contact interfaces are usually performed in built-in numerical codes on reduced order models (ROM). Most of the cases these derive from complex finite element (FE) models, resulting from the high accuracy the designers require in modeling and meshing the components in commercial FE software. In the technical literature several numerical methods are proposed for the identification of the nonlinear forced response in terms of a kinematic quantity (i.e. displacement, velocity and acceleration) associated either to the master degrees-of-freedom retained in the ROM, or to the slave ones after having expanded the reduced response through the reduction matrix. In fact, the displacement is the quantity usually adopted to monitor the nonlinear response, and to evaluate the effectiveness of a partially loose friction interface in damping vibrations, with respect to a linear case where no friction interfaces exist and no energy dissipation can take place. However, when a ROM is used the engineering quantities directly involved in the mechanical design, i.e. the strains and stresses, cannot be retrieved without a further data processing. Moreover, in the case of a strong nonlinear behavior of the mechanical joints, the distributions of the nonlinear strains and stresses over the structure is likely different than the one obtained as a superposition of linear mode shapes whose definition require a-priori assumptions on the boundary conditions at the contact interface. This means that the mentioned approximation cannot be used to predict the safety margins of a structure working in real (nonlinear) operative conditions. This paper addresses this topic and presents a novel stress recovery algorithm for the identification of the strains and stresses resulting from a nonlinear forced response analysis on a ROM. The algorithm is applied to a bladed disk with friction contacts at the shroud joint, which make the behavior of the blades nonlinear and non-predictable by means of standard linear analyses in commercial FE software.


Introduction
Blades are the critical components of turbine and compressor assemblies in aircraft engines. Due to the occurrence of severe mechanical and aerodynamical loads in operative conditions, blades are known to suffer high cycles fatigue (HCF) damages, which is considered the major cost, safety and reliability issue for gas turbine engine [1]. As so much is dependent on the reliability of these components, the tendency of manufacturers would be to over-design them to largely cope the safety specifications. On the other hand, limited weights are necessary to achieve high efficiency in the latest generation gas turbine engines. This aspect unavoidably leads to design much slender blades and thinner disks, making them prone to mechanical vibrations. Such a conflict requires the best possible trade-off made by the designers to keep the static and dynamics stresses within the safety margins. Therefore, a reliable simulative prediction of blades stresses in the design phase is necessary to prevent structural failures in service and reduce the maintenance costs.
Industrial practices for the design of a bladed disk prescribe FE analyses for the stress assessment at the critical locations of the assembly [2]. A reliable numerical identification of the stresses requires a deep knowledge of actual boundary conditions on the real component, such as the applied loads, the kinematic constraints, and the presence of contact interfaces for assembly. Nowadays commercial FE software are flexible enough to handle a wide variety of structural loading and constraints conditions, giving to the designers the possibility to customize the numerical simulation according to the actual, or expected, behavior of the components. However, an unaware usage of the software for complex FE analyses could lead to results whose interpretation might be ambiguous and misleading. For this reason, on one side industrial practices aim to be robust, since these must guide the engineers to the optimum design by providing well-established procedures resulting from the legacy gained in the development of similar structures in the past. On the other side, these become strict, limiting the designers to conventional assessments where the uncertainties on the boundary conditions are somehow considered in larger safety factors.
Besides the assessment against static failure, the design of low-pressure and high-pressure bladed disks must be verified against HCF. Therefore, designers must be aware of the static and dynamic loading conditions in order to carefully predict the stress distribution in service. Examples of static loads in bladed disks are the centrifugal force, the temperature field, and steady pressure distribution on the blades' airfoils. The dynamic loads are instead represented by the unsteady component of the gas flow resulting in traveling wave engine order (EO) excitation, which are the major responsible for forced vibrations in blade arrays Fig. 1 Campbell/Waterfall diagram of a turbine bladed disk: the disk natural frequencies are plotted as a nearly horizontal lines whose shape depends on the rotor speed . This is due to the blade stiffening effect caused by the increasing centrifugal force. The straight lines starting from the axis origin denote the EO traveling wave excitation, which represent the harmonic content of the unsteady pressure distribution exciting the bladed disk [3][4][5]. The occurrence of forced vibrations can be visualized in the schematic Campbell/Waterfall diagram in Fig. 1.
Among all the crossings between the mode shapes and the EO lines, the resonance condition occurs if the following relationship is satisfied: where N s is the number of disk's sectors (i.e. the number of blades) and h is the harmonic index of the mode shape, i.e. the number of nodal diameters of the excited mode shape [5]. Although some critical crossings can be moved outside from the operative range by slightly modifying the disk design (grey resonances in Fig. 1), some others cannot be avoided and additional damping besides the hysteretic one of the material is necessary to mitigate the effects of dangerous vibrations. This practice is crucial to avoid unacceptable level of dynamic stresses that would drastically reduce the fatigue life of the blades. In the last four decades plenty of research papers have addressed the problem of the blade's vibration mitigation in turbomachinery. The most exploited and widely accepted solution is the design optimization of the mechanical joints used in the assembly, to maximize the friction damping taking place at the contact interfaces due to the occurrence of macro and micro-slip phenomena. Typical applications where such a solution is successfully applied are the blade root and shroud joints for rotating bladed disks [6][7][8], the interlocking and hook joints for stator vane segments [9]. In this regard, linear FE analyses for the prediction of the static and dynamic behavior of the structure can just be used to get a general idea on the behavior of the system, without providing insights on how the actual stress distribution is effected by the presence of nonlinearities at the contact interfaces of the joints.
The interaction between two surfaces in contact in commercial FE software is simulated by means of contact elements implementing specific contact laws, e.g. node-to-node and surface-to-surface. Contact elements can be used in static and dynamic FE analyses to estimate the effects of friction on the static and steady-state response of a structure. In the static analysis the nonlinear deflection as well as the contact status at the joints' interfaces are identified in a one-shot run, while the evaluation of a nonlinear frequency response functions (FRF) requires to solve several time domain analyses in order to collect the steady-state responses of the system for each excitation frequency. To overcome the expensive computational costs typical of a direct time integration, the scientific community have developed advanced numerical techniques, usually implemented in purposely developed codes, for the prediction of the nonlinear FRF of structures with friction contact in the frequency domain [2]. The key idea behind them is the estimation of the contact forces by means of appropriate contact models [10][11][12] in the time domain and the solution of a set of nonlinear equations of motion (EQM) by using Harmonic Balance Method (HBM) [10] or the Multi-harmonic Balance Method (MHBM) [7,8]. In order to make the FRF prediction faster or even feasible, reduced order methods are used to condense the dynamics of large FE models (FEM) to few master DOFs. This is done by means of user-friendly routines available in commercial FE software that allow to extract the reduced order matrices (mass, stiffness and nodal loads) of the components according to wellestablished Component Mode Synthesis (CMS) methods, e.g. fixed-interface and free-interface methods. If on one side the solution of the EQM in a reduced fashion allows to achieve a remarkable reduction of the computational costs, on the other side the analysis results can just be obtained in the form of kinematic quantities (i.e. displacement, velocity, and acceleration) at the master DOFs, while the response at the slave DOFs as well as the other engineering quantities (e.g. the strains and stresses) cannot be directly retrieved without a further data processing. In this regard, stress recovery algorithms (SRA) allow to expand the reduced solution back to the full FEM [13], so that stresses and strains at the critical locations of the structure can be retrieved from a nonlinear FRF predicted in the frequency domain on a ROM. This paper proposes a novel stress recovery algorithm from a reduced nonlinear FRF of a bladed disks in the presence of friction contacts at the shroud. The recovered nonlinear stresses distribution at the blade airfoil is compared to the one of the mode shape resulting from a modal analysis where the behavior of the shroud is linearized assuming the full compatibility of nodal displacements (i.e. no relative displacements) at the contact interface. This comparison highlights the relevance of the proposed methodology for the identification of the critical location at the blade airfoil, since it changes with variation of the boundary conditions occurring at the contact interface.

Forced response of shrouded bladed disks
In this section the reduced nonlinear EQM of a bladed disk in cyclic symmetry (CS) conditions are obtained from the FEM of the disk's basic sector. Furthermore, it will be shown how to compute a nonlinear forced response in the presence of friction contact assuming the HBM hypothesis.

Equations of motion
The design of real engineering structures often requires the generation of accurate FEMs with a high density of nodes. Nonlinear forced response analyses on such structures are usually carried out with in-house developed codes [2,[14][15][16], which operate on a limited number of DOFs in order to avoid the expensive computational costs typical of iterative solution methods (e.g. the Newton-Raphon method (NRM) [17]). Substructuring methods in the class of CMS (e.g. fixed-interface and free-interface methods [18]) can therefore be used to export for each FE component the mass and stiffness matrices, as well as the loads vector, in a reduced form [19]. In this research the CMS method used to reduce the matrices of a bladed disk FEM is the well-known Craig-Bampton (CB-CMS) [20]. The choice is due to the fact that the CB-CMS method is usually included in the commercial FE software. However, the study is not strictly dependent on the chosen reduction method, and other reduction techniques can be used (e.g. Rubin [21]). Figure 2 shows the basic sector of a bladed disk. It consists of two different FE components, i.e. the blade and the disk sector, that are assembled by enforcing the compatibility of the nodal displacements at the blade root joint. Let us denote with x the vector of DOFs, which consists on all the nodal displacements in the three directions of a cylindrical reference frame Oρθ z, having the z-axis coincident to the disk axis and pointing as shown in Fig. 2. The DOFs vector x can be partitioned as follows: where x m represents the DOFs to retain as a master in the CB-CMS reduction, while x s is the set of remaining, i.e. slave, DOFs. The subscripts m and s are also used to denote the size of the master and slave partitions, which contain n m and n s DOFs, respectively. For the purposes of the nonlinear analyses discussed in the following, the sets of master DOFs can be further partitioned as follows: where: x c r and x c l are the vectors of n c contact DOFs, at the right and left shroud interfaces respectively (i.e. the DOFs corresponding to the green nodes in Fig.  2); -x a is the vector of n a active DOFs, i.e. the DOFs where the response is monitored in terms of kinematic quantities (i.e. the DOFs corresponding to the red nodes in Fig. 2). Being this set non-uniquely defined, their number and locations must be carefully chosen, because of their influence on the definition of the fixed-interface normal modes. Thus, it is suggested to adopt smart selection strategies for the master nodes to enhance the agreement between the ROM and the full FEM [22]. -x d r and x d l are the vectors of n d DOFs, at the right and left disk frontiers for the application of cyclic symmetry constraints (i.e. the DOFs corresponding to the orange dashed lines in Fig. 2).
Note that for a reliable application of the CS constraints all the DOFs at the disk frontiers have to be retained as a master, while the definition of the sets x c r , x c l and x a might be arbitrary and depends on the application. For instance, a good reduction practice would require to define x a by choosing few active nodes on the model so that the overall kinematics could be identified with sufficient accuracy. Similarly, assuming the hypothesis of conforming meshes at the contact interface, if the number of contact DOFs is particularly large, x c r and x c l might be defined by selecting homologous nodes at the right and left shroud interface in regions where the contact is expected. Otherwise, all the contact DOFs at the shroud could be retained as a master and further reduction methods might be used to approximate the contact DOFs by means of few interface modes [23,24].
According to the CB-CMS method, the full vector of FE DOFs in Eqn. 2 can be approximated as follows: where R C B = is the CB-CMS reduction basis, = I mm T sm T is the matrix of constraint modes, = 0 mm sk ] T T is the reduced matrix of n k n s fixed interface normal modes, and η k is the reduced vector of modal coordinates corresponding to the reduced set of fixed-interface normal modes sk . If the CB-CMS reduction is carried out with a commercial FE software, the complete transformation matrix R C B is computed and stored for a further stress recovery [19], and the CB-CMS reduced matrices and load vector are finally obtained applying the following transformations: Forced responses in bladed disks occur in the form of traveling waveforms excited by traveling engine orders (EO). The number of waves of the response at resonance denotes the harmonic index h of the excited mode shape. As a consequence, the DOFs at the disk frontiers vibrate with a phase shift ϕ h that is known as inter-blade phase angle and is defined as: This involves to express the motion of one frontier in terms of the displacements occurring on the other. Assuming as independent the set x d r , x d l can be expressed as: Eqn. 7 represents the so-called CS constraints and allows the exact representation of all the mode shapes with h nodal diameters for the full bladed disk, even if these are just applied to the sector frontiers. Moreover, it can be noted that CS constraints perform a further model order reduction that halves the number of DOFs at the disk frontiers (usually very large) from 2 × n d to n d . By applying the CS constraints of Eqn. 7 the CB-CMS reduced vector can be written as: where the bar symbol denotes the CS quantities. When the number of DOFs at the independent disk's frontier is particularly large, it could be convenient to apply an interface reduction method to operate a further reduction of the set x d r . If ii is a basis of n i interface modes, which are computed for instance by using the GSI [23] or Tran [25] method, the CB-CMS vector in CS conditions becomes: where the hat symbol denotes the array that describes the kinematics of the frontier in terms of modal displacements. The mass and stiffness matrices as well as the load vector in the final reduced form can be finally obtained as follows: Note that the order reduction in Eqn. 10 is effective if n i n d . In the following, the reduced matrices of Eqn. 10 will be used to build the nonlinear EQM of the bladed disk with shroud contact, and the hat symbol will be omitted for brevity.
The nonlinear reduced EQM of a bladed disk with shroud contacts can be written as: where M, K and f are defined in Eqn. 10, f c (t, x(t),ẋ(t)) is the vector of nonlinear contact forces acting of the physical partitions x c r and x c l , while C is the structural damping matrix obtained by applying the inverse modal transformation to the modal dampingC: (12) whereˆ is the modal matrix resulting from the solution of the eigenproblem associated to the homogeneous part of Eqn. 11, while ω n j and ζ j are the eigenvalue and modal damping corresponding to the j-th mode shape.
The steady state response of the system can be obtained solving the EQM in the frequency domain using the Harmonic Balance Method (HBM) [10]. Due to the periodicity of the EO excitation, the displacements and the nonlinear contact forces at the shroud can be approximated as the sum of n h harmonic terms: where nl represent the f -th order complex amplitudes of the displacements and contact forces respectively, ω the circular frequency and n h the number of retained harmonics. The substitution of Eqn. 13 into the differential EQM in Eqn. 11 leads to the following set of nonlinear, complex, algebraic equations: where D is the f -th order dynamic stiffness matrix, and the coefficients F nl , which are typically obtained from the time domain contact forces, depend on the complex amplitudes X ( f ) for f = 1, . . . , n h . For the purposes of this research the system of EQM 14 can be truncated to the 1-st order (HBM) without loosing the general validity of the method for a MHBM case. Therefore, assuming the HBM approximation the residual vector becomes: with the superscript (1) omitted for sake of clarity. Note that the solution for the unknown amplitudes X cannot be found analytically due to the nonlinear nature of Eqn. 14. The norm of the residual vector r has therefore to be minimized numerically by using iterative solution schemes such as the NRM.

Contact forces prediction
The problem of modeling periodical contact forces and the implementation in numerical solvers has been addressed by several authors. In the past years several node-to-node contact models for the evaluation of the contact forces between two nodes in contact have been proposed [26]. The choice of the contact element depends on the kinematics of the pair of nodes in contact. Being the shroud interface characterized by a 3-D periodic relative displacement, a contact element able to capture the normal load variation has to be preferred [27]. In order to well capture the 2-D trajectory of the relative displacement on the contact plane, two different modeling approaches are available. The first considers to combine two 1-D contact elements with normal load variation [12,28], whose tangential directions (i.e. the directions parallel to the contact surface) are orthogonal to each other. In this way the 2-D in-plane trajectory can be projected onto the two directions. The second accounts for using a 2-D contact element with normal load variation [29], for which the two orthogonal components of the tangential relative displacement on the contact plane are mutually coupled. In general, the projection of a 2-D in-plane trajectory onto two linear trajectories involves an underestimation of the friction damping [26]. Therefore, although less precise, the first modeling technique is more conservative from the point of view of the dynamic design, since it would lead to predict larger vibration amplitudes. In this paper the 1-D contact element with normal load variation is used to compute the periodic contact forces for a given periodic relative displacement. The 1-D contact element models three different contact states: stick, slip and separation. Figure 3 gives a schematic representation of the contact model, where the two dimensional relative displacement is decomposed into two perpendicular directions: one in-plane tangential components denoted by the node-to-node relative displacement u and the slider's displacement w, and one out-of-plane normal component v.
The contact model's parameters are represented by the tangential and normal contact stiffnesses k t and k n , the coefficient of friction μ and the normal preload n 0 (see Fig. 3).
At every time instant the normal contact force f n (t) is defined as: If n 0 is positive, the bodies are in contact before vibration starts, while if an initial gap g 0 is assumed, a corresponding negative n 0 can be calculated as n 0 = − g 0 k n . In the tangential direction the contact force is defined as: The stick, slip and separation states alternate each other during the vibration period according to the transition criteria reported in [11].
The contact models works only in the time domain, but the Fourier coefficients of the contact forces F nl participate to the dynamic equilibrium that is written in the frequency domain (i.e. Eqn. 11). Therefore, according to the Alternating Frequency Time (AFT) method [30], the inverse Fourier transform is used to recover the contact nodes displacements in the time domain from the Fourier coefficients X. Then, the time history of the tangential and normal relative displacements are used as inputs for the contact model in order to compute the time history of the contact forces. Finally, the Fourier coefficients of the contact forces are computed by using the Fast Fourier transform (FFT) and substituted into Eqn. 15.
The previous considerations on the contact force calculation hold for a pair of nodes sharing the same geometric location. In the case of a blade shroud, the contact is at the interface between two sectors, and the displacements used in the contact model belong to two adjacent sectors. When a single sector is used for the simulation in CS condition, the displacements of the pair of nodes of two adjacent sectors can be calculated by using the displacement of one single sector. Given a pair of homologous nodes N l and N r at the left and right shroud interface, having coordinates (ρ l , θ l , ζ l ) and (ρ r , θ r , ζ r ), it holds that: being θ s the sector angle. In this case the relative displacement to feed into the contact model must be computed using CS constraints similar to that of Eqn. 9. In particular, being x (s) the reduced vector of DOFs for the s-th sector treated in CS conditions (see Eqn. 9), the vector of DOFs for the sector s + 1 can be found as: Therefore, the relative displacement has to be computed between the node located at the right interface of the s-th sector and the node at the left interface of the (s + 1)-th sector (Fig. 4): where: The vector x rel , where the relative displacements u and v for each contact element are stored, must be used in the contact model so that the contact force f (s) cl is found. Finally, the contact forces at the left interface of the (s + 1)-th has to be shifted back in order to be applied at the left interface of the basic sector according to the following equation:

Stress recovery from the nonlinear forced response on the reduced order model
The nonlinear forced response of a bladed disks with shroud contacts is found by solving the EQM in the form of the residual function r, for each ω in a given frequency range [ω i , . . . , ω j , . . . , ω f ]. For the j-th frequency, the solution X j satisfying Eqn. 15 can be written as: an increase of the static preload n 0 makes the contact interface increasingly tight. This involves an increase of the resonance frequency from the free to the full-stick condition being X j and X j the real and imaginary part of the X j . The shroud effectiveness in terms of friction damping provided to the structure can be assessed in Fig. 5, where a typical nonlinear forced response is plotted together to the two extreme linear cases resulting from the linear behavior of the shroud.
For a given amplitude of the excitation, an increase of the preload n 0 at the shroud has the effect of shifting the resonant frequency towards higher values, moving from the free to the full-stick condition. The free condition corresponds to the absence of contact at the shroud interfaces, while the fully-stuck conditions occurs when the DOFs at the shroud interface are assembled by mean of the local tangential and normal contact stiffnesses k t and k n . Note that the relative displacement at the contact interface can be fully prevented if the CS constraints of Eqn. 7 are applied also to the shroud DOFs as follows: This condition, which is here referred to as tight condition, differs from the full-stick one for the total lack of local compliance at the contact interfaces. Nonlinear forced responses with a behavior close to the free and fully-stuck conditions can be found respectively for large negative and positive n 0 . A large negative n 0 would lead to an initial gap g 0 that never closes during the oscillation cycle. On the contrary, a large positive n 0 would prevent the occurrence of separation and sliding at the contact interface. From the behavior associated to the forced response in Fig. 5 important considerations can be inferred for two different aspects: the mechanical design of bladed disks and the identification of dynamic stresses on experimentally tested components. First, due to the absence of a source of friction damping, the prediction of dynamic stresses on the basis of the linear responses (either free or full-stick) is certainly conservative. If on the one side this implies the design of more robust structures due to the over-sizing of the components, on the other side it would involve a considerable increase in the weights, with a consequent loss of efficiency. Second, for the experimental identification of dynamic stresses on the structure, an easy and fast approach consists in the scaling of a measured quantity (e.g. the displacement) by means of modal quantities obtained from a FE modal analysis. For instance, if u P is the physical displacement measured at a point P of the blade airfoil, the distribution of the physical stresses ς on the structure can be computed as: where theς is the vector of modal stresses and K is a scaling factor defined as: This approach, hereafter denoted as a linear-mode scaling, can also be used in numerical simulations: a critical stress can be calculated even if it is not directly associated to the retained nodes in the ROM by scaling the simulated displacement in an active node. In this case K is calculated as the ratio of the simulated quantity divided by the corresponding modal quantity. However, scaling the stresses as in Eqn. 25 leads to assume for the operative deformed shape the stress distribution of the predominant linear mode shape. This assumption is acceptable if the Q-factor of the experimental forced response is consistent to the material damping, while it fails in the case of an evident nonlinear behavior of the joints. To overcome the mentioned limitation, a stress recovery algorithm (SRA) for the identification of the stress distribution on a bladed disks showing a nonlinear behavior of the partially loose shrouds is proposed. The SRA involves the re-expansion of the reduced nonlinear response at resonance to all the structure's DOFs, and allows the calculation of the corresponding stresses by means of static analyses performed in a commercial FE software. LetX res be the reduced nonlinear solution of Eqn. 15 at the resonance frequency f res = ω res /(2π) (Fig.  6). The vector of complex amplitudes of all the master DOFs retained in the CB-CMS reduction can be reconstructed by combining the transformation and reduction of Eqs. 8 and 9 respectively: BeingX res ∈ C (2×n c +n a +n i +n k ) , X C B is Eqn. 27 is a complex vector of size 2 × (n c + n d ) + n a + n k , meaning that the real and imaginary part of the expanded response to all the DOFs in the full FE model can be obtained using the transformation of Eqn. 3: Given the real (imaginary) part of the expanded response X at all the nodes of the FEM, the real (imaginary) displacement U = [U P ρ , U P θ , U P z ] T = U(P, ω res ) at a particular point P can be calculated using the shape functions η(P) for displacements [31], so that: Under the hypothesis of a linearly elastic material, the strains (P, ω res ) depends linearly on the displacements U (P, ω res ) through the strain operator S: Using the strain of Eqn. 30 and the constant material matrix E, the stresses in P can be computed as: If Eqn. 31 is applied for all the nodes coordinates in the FE model, the real and imaginary nodal stress components are found, and the global vector of complex nodal stresses can be finally obtained as follows: From a practical point of view, this methodology can be implemented in many ways. Since the FE software are graphically optimized to manage large FE models and offer functions to accurately plot and visualize detailed information related to the mesh, a good balance between custom and commercial FE software' features has to be found. In this study, the work flow for the identification of the stress distribution corresponding to a nonlinear response of the bladed disk at resonance is summarized in the following step-by-step procedure: -Mesh the bladed disk sector and extract the reduced CB-CMS mass matrix M C B , stiffness matrix K C B , and nodal force vector f C B (Eqn. 5. This step is performed in a FE software (e.g. Ansys APDL), so that the bladed disk superelement is created for further solution expansion and the reduction matrix R C B stored. -Import the CB-CMS ROM into a Matlab-based code (or equivalent), and apply the CS constraints and the interface reduction method as in Eqs. 8 and 9. Note that the CS constraint can be applied in the FE environment before performing the CB-CMS reduction. In the first case, the reduced matrices have a more general form and are suited to be treated with different CS constraint as required by the user. In the second case, the reduced matrices are specifically ready to be used with a specific harmonic index h, but other CB-CMS reductions are needed if the CS constraints must be updated for a different h. -Define the contact parameters and compute the nonlinear forced response on the ROM of Eqn. 11.
-Expand the nonlinear displacements at the resonance frequency f res to all the master DOFs, as in Eqn. 27. -Apply the complex displacements at the master nodes into the superelement model, and solve a FE harmonic analysis for the resonance frequency f res , assuming a null damping matrix. Note that this step just allows to transfer the complex amplitudes of the master DOFs to the CB-CMS superelement's master nodes. -Expand the real and imaginary parts of the solution separately to all the nodes of the full FE model using the transformation of Eqn. 28. -Given as inputs the real and imaginary displacements computed at the previous step, perform two static analyses in order to compute the real and imaginary stresses (see Eqs. 29, 30 and 31).
The SRA is here presented assuming the first-order approximation for the contact displacements and contact forces as well (Eqn. 15). When using the MHBM, it would be necessary to implement the SRA for each harmonic component retained in the approximation. The actual real stress is represented by the maximum stress value resulting from the superposition of the harmonic stress components using the inverse Fourier transform.

Application
In this section the SRA is used to recover the stresses from a set of nonlinear forced responses of a bladed disk for different values of the static preload n 0 . The Modal Assurance Criterion (MAC) for complex vectors [32] is used to compare the shape of the nonlinear stress distributions with a linear one, obtained by preventing the relative displacements at the shroud nodes in a CS modal analysis (full-stick case). Furthermore, the Von Mises stresses are computed in order to keep track of the critical location at the blade airfoil for each n 0 . The bladed disk has N s = 112 identical sectors, each of which consists on a shrouded blade and a disk sector assembled at the blade root joint (Fig. 7). Due to the limited friction occurring at blade root it was decided to neglect the damping effects caused by such joint on the nonlinear forced responses. Therefore, the total compatibility of nodal displacements was enforced in the FE model so that no relative displacement might occur at the contact interfaces between disk and blade. The sector representative of the whole wheel The CB-CMS ROM of the fundamental sector was created in the ANSYS environment performing a classic substructuring analysis. The size of the master DOFs partitions retained in the CB-CMS reduction are given in Table 1: x c l 189 x a 78 η i 50 η k 200 Fig. 8 Bending mode shape in the free (a) and tight (b) conditions. It can be noted that in the free condition the maximum vibration amplitude occurs at the mid-airfoil, while the shroud participation is less important. This leads to a frequency shift between the two modes of ∼ 10% (Fig. 9) CS constraints were applied to the {x d l } and {x d r } for the harmonic index h = 12 (Eqn. 7), and a basis of n i = 50 GSI modes [23] was used to further reduce the physical displacements at the independent disk frontier. The reduction process led to a highly compressed ROM, whose first 10 natural frequencies differ for less than 0.05% with respect to the natural frequencies of the full FE model solved in ANSYS for the same h. The size of the DOFs partitions for final ROM are listed in Table 2.
The nonlinear forced responses were performed assuming all the shroud nodes in contact. The frequency range was chosen in order to excite a bending mode in both the free and tight conditions shown in Fig. 8.
The contact stiffnesses k t and k n were set in order to be close to natural frequency of the mode with a fullystick shroud, to that of the FE model where the CS Fig. 9 Free and full-stick linear FRF (blue and red dashed plots), and nonlinear FRFs for different n 0 at the shroud contact interface (continuous plots). The vibration amplitudes and the frequency range are normalized with respect to the peak and the resonance frequency ( f res ) of the full-stick response constraints were also applied to the contact interfaces. The friction coefficient μ and the modal damping ratio ζ were assumed to be equal to 0.5 and 0.002 respectively. The magnitude of the reduced vector of nodal forces was tuned to get realistic vibration amplitudes for the linear response, i.e. ∼ 10 −4 m 0-peak at the node where the maximum displacement occurs for the mode shape with shroud in full-stick condition.
Being p 0 the reference static preload corresponding to the design interference at the shroud, the nonlinear analyses were performed for each of the following n 0 : From the point of view of the numerical simulation p 0 is the average preload obtained by solving a FE nonlinear contact analysis where the design interference is implemented as a compenetration of the shroud interfaces. Figure 9 shows the complete set of nonlinear forced responses (continuous line plots) as well as the free and full-stick linear responses (dashed line plots) for the accessory node shown in Fig. 7, in the tangential direction. Note that the linear cases were obtained by performing linear forced response analyses and are here reported to better understand the nonlinear behavior of the blade caused by the shroud.
It can be observed that for a fixed external load an increase of the static preload n 0 has the effect to increase the frequency of the resonance peak (Fig. 10), moving from a contact condition where sliding occurs (curve close to the free case), to a condition where smaller relative displacements at the shroud take place Besides the changing in the resonance frequency, which suggests a variation of the joint stiffness, varying n 0 involves a change in the amount of friction damping. In particular, a decrease of n 0 from the largest positive value (i.e. n 0 = 10 · p 0 ) up to the largest negative value (i.e. n 0 = −10· p 0 ), makes the joint increasingly loose, with a consequent increase of slippage that involves two different trends (Fig. 11).
First, decreasing n 0 from 10 · p 0 down to p 0 corresponds to decrease the vibration amplitude at resonance in a significant manner, (i.e. -89.9%). This is due to the increasing effectiveness of slippage in the generation of friction damping. Second, once the maximum damping effectiveness is achieved, a further decrease of n 0 let the relative displacements at the contact interface increase. In this region the slippage is large because the vibration amplitude is increasing too, and the joint Here it can be observed that the highest values correspond to the largest preloads, while the minimum occurs for n 0 = p 0 .
For each n 0 the nonlinear response at resonance was expanded to the full FE model by following the workflow described in Sect. 3. The real and imaginary parts of the solutions were exported in terms of equivalent Von Mises stress for all the nodes at the blade airfoil. The resulting nonlinear, complex stress distributions were then compared to the corresponding complex linear stress distribution obtained from the modal analysis in tight condition. The term "corresponding" is used to highlight that the same set of shroud nodes was used as a contact set for the nonlinear forced response, and as a coupling set for the modal analysis. An equivalent comparison could be performed for a smaller set of shroud nodes detected from a nonlinear static anal- Fig. 13 Trend of the MAC between the modal stress distribution obtained with the shroud in tight condition and all the expanded nonlinear forced responses ysis, where the blade is loaded by the rotational speed, the steady pressure, the thermal loads, and the design interference at the shroud [33].
The comparison was carried out using the MAC for complex eigenvectors [32], which is here re-arranged for the purposes of the current comparison as follows: where σ p j is the complex vector of physical Von Mises stresses expanded from the forced response corresponding to the j-th static preload (Table 3), whilẽ σ t is the complex vector of modal Von Mises stresses obtained with the shroud in tight conditions. Figure  13 shows the MAC evolution due to the n 0 variation, while Fig. 14 shows a qualitative comparison between the map of the reference Von Mises modal stressσ t (Fig. 14a), and the stress maps of the expanded nonlinear Von Mises stresses σ p 1 , σ p 5 and σ p 9 (Fig. 14b, c and d respectively). From Fig. 13 it can be noted that for 2.5 · p 0 ≤ n 0 ≤ 10 · p 0 the MAC gets close to the unity, meaning that the shape of the excited nonlinear mode is close to that of the eigenvector obtained for tight conditions at the shroud. For n 0 < 2.5 · p 0 the nonlinear mode drastically changes its shape, and the corresponding MAC suddenly drops to ∼ 0.27 when n 0 decreases down to −2.5 · p 0 . Finally, for large negative n 0 s the MAC slightly increases, approaching to the value of ∼ 0.39, which corresponds to the MAC between the FE stress eigenvectors obtained for the free and tight conditions at the shrouds. Furthermore, the change in  (Table 3). For each stress map the location where the maximum stress occurs is circled in red Fig. 15 Evolution of the critical location at the blade airfoil for each nonlinear forced response. The colored circles denote stress concentration at specific airfoil locations: lower fillet (blue, critical location N 1 ), mid-airfoil (green, critical location N 2 ) and upper fillet (red, critical location N 3 ) the shape of the excited mode leads in general to a different location where the maximum Von Mises stress occurs. In particular, let σ r N c be the absolute value of the Von Mises stress at the critical location N c for the r -th expanded forced response: The positions of σ r N c for all the expanded solutions are shown in Fig. 15, where the critical locations N 1 (blue circle) and N 3 (red circle) are the same identified for the linear modes obtained for the shroud in free and tight conditions respectively.
From Fig. 15 it can be noted that N 1 and N 3 are located at the lower and upper fillet of the blade airfoil respectively. This result is somehow expected for different reasons. First, the fillets are machined to avoid high stress concentrations due to a sudden variation of the blade geometry, as would be the case for the boundaries blade platform-airfoil (lower fillet) and airfoilshroud (upper fillet). Nevertheless, such boundaries still remain critical from the stress concentration perspective. Second, the blade with the shroud in free conditions behaves as a clamped-free beam (Fig. 8a), which leads to expect higher stresses close to the clamp, i.e. at the blade root, where the blade is clamped to the disk and the moment due to bending is higher. Third, the blade with shroud in tight conditions behaves as a clamped-clamped beam (Fig. 8b), where the higher stress due to bending is expected close to one of the clamps (either lower or upper fillet), depending on the geometry of the fillets.
It can also be noted that the critical location N 3 for the expanded nonlinear response, remains the same identified for the linear mode in tight conditions, until the corresponding MAC in Fig. 13 remains close to the unity (cases 6 to 9). For lower MACs the critical location changes as well, moving from the upper fillet (critical location N 3 ) to the lower fillet (critical location N 1 ), passing by the mid airfoil for n 0 = p 0 (critical location N 2 ). The latter case corresponds to a MAC slightly less than 0.7 (Fig. 13), which suggests the occurrence of a nonlinear stress distribution significantly different from the ones obtained for the two extreme linear cases, i.e. blade with tight shroud (MAC close to unity) and blade with free shroud (MAC close to ∼ 0.39). For this reason, it could be expected to find a critical location different than N 1 and N 3 as actually happened.
The change in the MAC as well as in the critical location makes the linear-mode scaling described in Sect. 3 (Eqn. 25) no more valid in general. In fact, although Eqn. 25 can be applied to identify the stresses for all the nonlinear shapes having MAC close to the unity, it does not hold for the cases with low MACs and different critical location with respect to that of the linear mode. For instance, if the 5-th forced response is considered, from Fig. 15 it can be noted that the critical location at the blade airfoil does not coincide to that of the linear tight mode (critical location N 3 ). If σ 5 N 3 is the stress at the location N 3 (i.e. at the critical location for the tight mode shape) for the 5-th forced response, it occurs that σ 5 N 2 /σ 5 N 3 = 1.22. Meaning that if Eqn. 25 is used to scale the modal stresses of the tight linear mode using the K factor obtained as the ratio between the maximum displacement at a certain accessory node and the corresponding (tight) modal displacement at the same node, the maximum stress would be underestimated of the ∼ 22%.

Conclusions
The application presented in the previous section highlights the need to adopt a SRA for the stress identification in a nonlinear forced response. In the technical literature the performance of a friction joint on the dynamic response of a structure is usually assessed by monitoring the maximum vibration amplitude at a specific, given location. However, a robust mechanical design against HCF requires the prediction of the static and dynamic stresses, in order to define the working point on a Haigh diagram for the identification of the safety factor.
Although the static stress assessment might be carried out by just solving a static analysis in commercial FE software, the evaluation of dynamic stresses for a nonlinear contact problem is not straightforward. In fact, due to the convenience to solve the reduced nonlinear EQM in the frequency domain, the solution, usually expressed in terms of kinematic quantities at the master DOFs, has to be re-expanded to the FE mesh for the stresses calculation. This practice allows to detect the actual stress distribution on the structure, and overcomes the limitations related to the unreliable scaling of the modal stresses in the case of a pronounced nonlinear behavior of the joint. In particular, it can be noted that scaling the modal stresses holds when a weak variation of the joint stiffness occurs. In fact, for the forced responses 6, 7, 8 and 9, which are close to the fullstick forced response, the resonance frequency does not change significantly, and the shape of the stress distribution still remains similar to that of the tight linear mode shape. Instead, for the 4-th forced response, although slight variation of the resonance frequency, the MAC drops to 0.27. This result suggests an important modification of the mode shape as well as of the critical location, which moves from the upper fillet to the lower fillet of the mid airfoil, passing by the mid trailing edge (Fig. 15). Therefore, as previously shown in Sect. 4, besides the wrong identification of the critical location, scaling the modal stresses would have led to a non-negligible underestimation of the maximum stress.
Finally, the SRA, if embedded into the design process, can provide qualitative and quantitative indications on the changing of the stress distribution due to the loss of contact interference and variation of n 0 . In fact, the frictional joint does not maintain the same assembly interference in operation due to the occurrence of wear due to fretting [34,35]. Therefore, predicting the evolution of the critical locations might improve the mechanical design in order to ensure safety operative conditions for all the loading cases.
Due to the large employment of frictional joints in a wide variety of mechanical systems, the presented methodology is suitable not only for turbomachinery applications, but also for many other engineering fields where the nonlinear dynamic behavior of the structure is usually predicted in the frequency domain on a ROM.

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

Data availability
The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.

Code availability
The code is written according to the proposed model.
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/.