Effects of centrifugal stiffening and spin softening on nonlinear modal characteristics of cyclic blades with impact–friction coupling

This paper aims to interpret the coupling modal properties of cyclic blades under impact–friction interactions and their evolution mechanism versus operating points. Therefore, a coupling analytical model of cyclic blades is developed based on a Lagrange method and the assumed mode method (AMM), after considering centrifugal stiffening, spin softening, stagger angle, and twist angle. Then a mixed modal analysis method (MMAM) for this analytical model is extended by combining the extended periodic motion concept (EPMC) with AMM. Wherein a classic alternating frequency/time method (AFT) and the continuation method are employed to overcome the numerical divergence problem. Then damped nonlinear normal modes (dNNMs), including eigenfrequencies, modal damping ratios, and mode shapes, of the coupling system with shroud joints are finally computed and discussed under different excitation levels and contact conditions through a modal synthesis algorithm. After that, the influence laws of centrifugal stiffening and spin softening on the dNNMs are explored to reveal its evolution mechanism versus operation speeds. Finally, the Campbell diagrams of dNNMs are successfully obtained to discuss the effects of the impact–friction coupling on critical speeds (CSs) of the shrouded blades system.


Introduction
The severe load environment will induce a high vibration response of blade components [1,2]. Therefore, the application of special devices, such as shroud structure, under-platform damper, dovetail joints, strip dampers, and lacing wire [3][4][5][6][7], can dissipate vibration energy and reduce the dynamic responses to acceptable levels is usually essential. Among all mitigation measures, the shroud friction damper between adjacent blades is widely used in engineering practice, which can absorb vibration energy and greatly reduce the vibration of rotating blades. And the effectiveness of reducing forced response levels by using shroud dampers will completely depend on the choice of the damper parameters. For obtaining the optimal selection of these parameters, it is necessary to exploit the complex nonlinear vibration properties of shrouded blades from the perspective of dynamics, especially at operating speed.
To accurately predict the forced response characteristics of the shrouded blade system, many scholars have made great efforts. Siewert et al. [8] presented a reduced-order model of a shrouded blade system coupled with a state-of-the-art three-dimensional contact model, which is solved by the multi-harmonic balance method and the alternating frequency-time method. Next, the forced properties of the system and its variation versus the initial gap were shown and discussed. A good friction model is a foundation for obtaining the damping characteristics of shrouded blades. Thus, Xu et al. [9] developed a new friction model after considering the effects of the topography of contact surfaces by Hertz contact theory and fractal geometry. The proposed contact model and trajectory tracking method were applied to compute the forced response of a real shrouded blade. And analysis results indicate that the resonance amplitude especially depends on initial normal load and contact surface roughness. In the FEM model of a shrouded blade, the contact interface is usually discretized into a series of contact pairs, inducing the problem of large-scale nonlinear computation. As a consequence, Afzal et al. [10] formulated a new expression of the 3D contact model in the frequency domain, which can successfully predict the friction damping characteristics of the shroud contact interface and greatly save the computation time. After that, the method was used to compute the friction damping characteristics of a bladed disk system with shroud contact and strip damper [11]. Numerical results reveal that multiple friction contact interfaces exhibit a better damping effect in contrast to the case with a single friction interface. In addition to the efforts in the treatment of friction forces, Zucca et al. [12] proposed two reduction strategies for shrouded bladed disks by utilizing the cyclic symmetry. The comparison results of their accuracy and computational effort prove that both methods ensure high accuracy, and the ratio of the number of reserved nodes will especially influence the computational efficiency. Additionally, for overcoming the limitation of computation efficiency of the finite element model, some scholars simplified the shrouded blades as a lumped parameter [13,14] or continuum model [15,16] to carry out the vibration suppression mechanism of shroud joints, such as the influence of contact stiffnesses, rotating speed, shroud position, and excitation parameters on the vibration responses of the shrouded group blades.
The nonlinear modal analysis is also an important aspect of facilitating the understanding of the vibration attenuation mechanism of the shroud damper and improving the capability of the vibration reduction optimization design. Unlike linear systems, the superposition principle and orthogonal conditions are no longer valid in nonlinear systems. Therefore, the development of nonlinear normal mode (NNM) has to extend the concept of linear mode to nonlinear systems. At present, there are two main definitions for nonlinear modes: (i) The first one is based on the periodic motion, which was first presented by Rosenberg [17], and the other one is the definition of invariant manifold, firstly put forward by Shaw and Pierre [18,19]. In definition (i), all degrees of freedom in a system only generate proportional synchronous vibration without phase differences. This classic definition is for a nonlinear conservative autonomous system, which is extended to the periodic motion of an undamped system and called damped nonlinear normal mode (dNNM) [20]. However, this definition is not applied to the non-conservative nonlinear system due to its autonomous aperiodic solutions, e.g., the system with dry friction dampers. The advantages of this definition are that common solution algorithms for solving periodic solutions can be used to calculate NNM, such as the harmonic balance method [21,22] and shooting method [23]. In definition (ii), the NNM is defined as the invariant manifold in high-dimensional phase space. The invariant manifold here can be understood as a hypersurface in the high-dimensional phase space, which is tangent to the equilibrium point of the linear modal subspace corresponding to the degenerate linear system [24]. The definition (ii) can break through the limitation of the conservative system to the non-conservative system by both analytical [18] and numerical methods [25]. But the process of determining the entire manifold simultaneously greatly increases the computational cost. Moreover, the above analytical and numerical methods for calculating invariant manifolds cannot be applied to systems with discontinuous nonlinearity.
The comparative analysis indicates that the periodic motion definition of the NNM in a non-conservative nonlinear is easier, which only calculates a single periodic orbit at the same time. Therefore, most researchers deal with the nonlinear modal analysis of systems with dry friction based on the definition of periodic motion through three different techniques.
The first technique was called nonlinear complex mode (CNM), proposed by Laxalde and Thouverez [26]. This method extends the complex mode to the nonlinear system and describes the non-conservation of energy in the system through the modal damping ratio. Therefore, the CNM of a system with friction dampers can be expressed by a quasi-periodic solution and attenuation factor. Secondly, Krack [27] proposed an extended periodic motion concept (EPMC), in which virtual viscous damping is added to balance the energy dissipation, induced by friction dampers. Hence, the dNNM can be represented by a periodic solution and modal damping ratio (i.e., virtual viscous damping). Then the nonlinear modal characteristics of a turbine bladed disk model with shroud contact were analyzed to demonstrate the effectiveness of the solution method [28]. Recently, Sun [29] has tried to add virtual hysteretic damping to achieve a similar effect to EPMC, which may improve calculation efficiency in engineering practice. Wherein the comparative analysis of CNM and EPMC shows that CNM can only be solved by HBM, while EPMC is also applicable to the shooting method in the time domain method.
At present, the previous methods have been widely recognized to evaluate the modal characteristics of dry friction dampers through nonlinear modal analysis. But these methods are developed based on the assumption that the dynamic behavior is only dominated by an isolated resonance [24], limiting its application in systems with dense modes, e.g., the shaft-disk-blade structure. The most typical feature of such structures is that the coupling effects between sub-components will activate the eigenfrequency veering and mode exchange phenomena under various rotating speeds or other system parameters [30][31][32][33]. Therefore, to compute the dNNM of such a system with multiple friction dampers, She et al. [34] developed a new algorithm based on the harmonics equivalent technique and incremental HBM, called E-IHBM here. Then a shaft-disk-blade system subjected to dry frictions is taken as an example to compute and analyze the nonlinear modes, proving the validity of the E-IHBM method. Then based on the E-IHBM and AMM, She et al. [35] developed a representative model of a bladed disk system with dovetail joints and successfully computed the Campbell diagram of such structure. Moreover, this study clearly explained the eigenvalues of veering and merging phenomena, induced by the coupling between the bladedominated and disk-dominated families of modes.
The E-IHBM is proposed and developed based on an assumption that only the fundamental harmonic dominates the dynamics, neglecting super-harmonics. This may reduce the accuracy of modal results, especially in the system with a strong nonlinearity, e.g., shrouded blade system. Therefore, the EPMC serves as an important inspiration for the method developed in the present thesis. Firstly, a coupling model of the cyclic blades system with shroud joints is established by the analytical method, considering the centrifugal stiffening, spin softening, twist angle, and stagger angle, see Sect. 2. The dynamic model also considers the cyclic symmetry of the system and the nonlinear impact-friction coupling behavior between the adjacent blades. Then in Sect. 3, this study combines the extended periodic motion concept (EPMC) and AMM, developing a mixed modal analysis method (MMAM) for this analytical model. Finally, the dNNMs of coupling systems under different operation speeds are successfully computed and discussed utilizing the derived analysis strategy in Sect. 4, especially revealing the variation mechanism of CSs versus operation speeds and contact conditions. To exploit the nonlinear dynamics of a cyclic blade system with shroud interfaces, multi-Timoshenko beams with varied cross sections are proposed, as depicted in Fig. 1. The coupling model consists of N b blades and is assumed to be perfect cyclic symmetry. Therefore, the response of shrouded blades is also cyclic symmetric, and only a single sector with cyclic boundary conditions is considered in the present investigation, cf. Fig. 1a. Nonlinear boundary conditions of the shroud interfaces are simulated by using a nonlinear impact-friction model, introduced in the subsequent sections. The schematic of a single segment is reported in Fig. 1b, where the O-XYZ is the global coordinate system, fixed on the center of the rigid disk, and o-xyz is the local coordinate system at the blade root. Moreover, for representing the elastic deformations of a pre-twist blade, two local reference frames (o 1 -x 1 y 1 z 1 , o 2 -x 2 y 2 z 2 ) are also displayed, cf. Fig. 1c.
Note that o 1 -x 1 y 1 z 1 is fixed on an arbitrary cross section of the ith blade after rotating b(x) = b L x/L b around the x-axis, and o 2 -x 2 y 2 z 2 locates on the blade tip. u, v, and w represent the displacements of any point on the rotating blade in the radial, flexural, and swing directions, respectively. u L , v L , and w L represent the displacements of the blade tip in the radial, flexural, and swing directions, respectively. Additionally, the width and thickness of the tapered blade at the root are b 0 and h 0 , while b L and h L represent the width and thickness at the blade tip. Moreover, the stagger angle b 0 at the blade root and twist angle b L at the blade tip are also considered to depict the twisted surface of the shrouded blade. Consequently, the twisted angle of any cross section is b( the twisted angle of the blade tip relative to the blade root and b 0 is the stagger angle at the blade root.
2.1 Derivation of the ith blade's energy expression

Derivation of the kinetic energy expression
In the present analysis, a Timoshenko beam is used to simulate the bending deformation of the cyclic blades. According to Ref. [36], the cross-sectional and side views of the bending deformation are depicted in Fig. 2a and b, respectively. Note that points Q 0 and Q 1 represent the reference point before deformation and after deformation.
According to the deformation diagram in Fig. 2, the displacements of the arbitrary point Q 1 in the local coordinate o 2 -x 2 y 2 z 2 can be expressed as follows: Cyclic boundary condition β L β 0 Fig. 1 Diagram of a cyclic blades system with shroud interfaces: a cyclic blades model coupled by shroud interfaces, b a single segment of the ith shrouded blade, c the coordinate system of the ith blade where R D is the outer radius of the disk; u is the radial displacement due to the centrifugal force, g and n are the local positions of the arbitrary point Q 0 ; v and w are the bending displacements; and u v , and u w are the rotation due to bending. The rotation due to bending, u v and u w , are small so it is assumed that sinu w % u w , sinu v % u v . Using these assumptions, Eq. (1) can be rewritten as follows: Therefore, according to the reference frame in Figs. 1 and 2, the displacements of an arbitrary point Q 1 on the ith blade can be stated as follows: where b(x) = b 0 ? b L x/L b , b 0 is the stagger angle and b L is the pre-twist angle at blade's tip cross section; and # i = Xt ? (i -1)2p/N b is the position of the ith blade in the circumferential direction, in which X and N b are the rotating speed and blade number. The kinetic energy expression of an arbitrary ith blade is as follows: where q b is the blade density and L b is the blade length. After substituting Eqs. (3) into (4), the specific expression of kinetic energy can be obtained as follows: Fig. 2 Schematic diagram of the blade's deformation: a cross-sectional view, b view of bending deformation where A b is the blade cross-sectional area. While I z and I y are the blade area moment of inertia, see ''Appendix A'' for specific expressions.

Derivation of the kinetic energy expression
Based on the derivation of Eq. (3), knowing that the location vectors of the reference point before and after deformation in the local reference frame o 2 -x 2 y 2 z 2 are represented by r 0 and r 1 , which can be expressed as, (a) Before deformation (reference point Q 0 ): (b) After deformation (reference point Q 1 ): where i, j, and k are the unit vectors in the x 2 , y 2 , and z 2 directions, respectively. Thus, dr 0 and dr 1 can be given by The components of dr 0 and dr 1 can be given by dx 0 ¼dx ð8aÞ Then the classical tensor e ij can be obtained using the equilibrium equation below: e xx e xg e xn e gx e gg e gn e nx e ng e nn As for the long slender beam, the radial strain e xx is dominant over the transverse normal strains e gg and e nn . Additionally, the shear strain e gn is two orders smaller than the other shear strains e xg and e xn . Therefore, after neglecting the high-order terms, the elements of the strain tensor e ij can be simplified as: Based on previous derivations, the bending potential energy U bti of the ith blade is given by Then the final form of the shear potential energy of the ith blade, U bsi , can be stated as where j b and G b are the shear correct factor and shear elastic modulus, respectively. G b J b is the torsion rigidity. Therefore, the total potential energy of the ith rotating blade is as follows:

Discretization and equation of motion
Based on the previous sections, the total energy of the SDB unit is obtained as follows: where N b is the blade number of the coupling system. Then the assumed mode method (AMM) and Lagrange equation are applied to discretize the energy functions and obtain the motions of the coupling system. In the present analysis, the mode function of a Timoshenko beam is used to discretize the energy functions of the flexible blade, which can be expressed by where q uj (t), q vj (t), q wj (t), q wvj (t), and q wwj (t) are the generalized coordinates, respectively. In addition, U j (x) and V j (x) are the mode function of longitudinal displacement and lateral displacement, and / vj (x) and / wj (x) are the cross-sectional rotation, which are as below: where a j = (2j -1)p/2/L b . Finally, for obtaining the motion equation of the coupling system, the Lagrange equation is applied where L = T Total -U Total is system Lagrangian and q is generalized coordinates. After substituting Eqs. (15) and (16) into Eq. (14) and applying the Lagrange equation, the motion equations of the coupling system are obtained as follows: where M BB and K e are the mass and stress stiffness, respectively. K c is the centrifugal stiffening matrix, while K s originates from the spin softening effect. q BB = [q u ; q v ; q wv ; q w ; q ww ] is the generalized coordinate. A dot above q BB denotes derivative with respect to time, t. Therefore, after further considering the terms of the external force and linear damping, the nonlinear dynamic equations of the cyclic blade system are given by Þ is a nonlinear impact-friction force of the shroud interface, which will be formulated in the following sections. q nl t ð Þ is the generalized coordinate of nonlinear DOFs. A dot above q nl t ð Þ denotes derivative with respect to time, t. f ex t ð Þ is a traveling wave excitation [34], which can be expressed as where i represents the imaginary unit and m = 1, 2, …, N b . / is the phase shift between the adjacent blades, defined as / = 2p/N b . Under this periodic excitation, the periodic response of the blade can be obtained easily.

Validity of the coupling model
To validate the developed model in previous sections, a finite element model (FEM) of the blade with tapered and pre-twisted cross section, obtained by ANSYS, is established, taking account into the centrifugal stiffness and spin softening effects. Then the Campbell diagram and harmonic response results obtained by presented mode and FEM are compared, respectively. Wherein the FEM of the blade is meshed by 20 BEAM188 elements, and the node of the blade root is completely constrained. Additionally, all the parameters of the rotating blade are shown in Table 1. Figure 3 displays the effects of rotational speed on the first three modal frequencies of rotating blades, which are obtained by the present model (solid line) and FEM (dotted line), respectively. It can be figured out that the solid line and dotted line basically coincide, especially for the 1st modal frequency. The tiny error between the 3rd modal frequency of the present model and FEM can be attributed to the reason that the high-order mode may be affected by the modal truncation number. Altogether, the frequency curves of the two models agree well, validating the accuracy of the present model.
The harmonic responses of the present model and FEM are plotted in Fig. 4 to further validate the

Modal analysis method
The nonlinear modes of the cycle blades with shroud joints, i.e., a typical non-conservative nonlinear system, are sometimes referred to as damped nonlinear normal modes (dNNMs) [20], which are no longer periodic. For applying the classic methods for periodic problems to predict dNNMs, this study combines the extended periodic motion concept (EPMC) [25] and AMM, developing a mixed modal analysis method (MMAM) for the analytical model of Sect. 2, as plotted in Fig. 5.

Prediction of modal characteristics
Equation (19) is the general dynamic equation for computing the nonlinear forced response by using the classic multi-harmonic balance method (MHBM) in the frequency domain [22,34]. Nevertheless, dNNMs are the free vibration solutions of the system, i.e., not taking account into the linear damping and external force. In such a case, the relative motion of the shroud interface will dissipate energy, resulting in a negative artificial damping effect. Therefore, based on the Fig. 5 Flowchart of the nonlinear modal analysis method EPMC method, an artificial damping term C d ¼ À2fx 0 M BB is employed to balance the energy loss, Herein, f and x 0 are the nonlinear modal damping ratio and eigenfrequency of the considered nonlinear mode. Note that dNNMs will be periodic if we can find a proper value of 2fx 0 , since MHBM can be used. Therefore, the generalized coordinate vector q BB (t) and nonlinear impact-friction force f nl q nl t ð Þ; _ q nl t ð Þ; t ð Þcan be truncated after N h harmonics as and Herein, a is the modal amplitude. W k is the eigenvector of kth harmonics of the generalized coordinate vector. The vector F nl,k is the kth harmonic coefficient of the impact-friction force. It can be figured out that only W 1 = 0 and F nl,1 = 0 if the contact interface is under the full-sticking or completely separate states, similar to the linear case. Substituting Eqs. (22)-(23) into Eq. (21), the nonlinear dynamic equations in the frequency domain are rewritten as Þ is the dynamic stiffness, in which diag is a diagonal operator.
And F nl (x 0 , aW 1 , …, aW N h ) = [F T nl;0 , …, F T nl;N h ] T is the vector of impact-friction force, which will be computed by using the alternating frequency/time method in the following section. In Eq. (25), a, x 0 , and f are both unknown, leading to Eq. (25) being under-determined. To obtain a unique solution of W, additional constraints should be provided. Thus, mass and phase normalizations are used in this method,

Alternating frequency/time method
For the impact-friction behavior of the contact interface, many scholars [37][38][39][40] have developed lots of representation models. This study chooses a 2D contact model, which is widely used due to its fewer model parameters. As depicted in Fig. 6a, the full-2D impact-friction formulation of a shroud interface depends on the relative motions of adjacent blades and the motion of the damper, connected by two linear springs (K n and K t ), which represent the normal and tangential contact stiffness, respectively. According to Coulomb's law, the shear and normal relative displacements (X Ai , Z Ai ) are discretized into N points for one motion period for tracing the displacements of the damper (X Di , Z Di ). Therefore, X Ai (s j ) and Z Ai (s j ) [j = 1, 2, …, N] represent the relative displacement of the shroud interface at the jth discrete point. Therefore, the impact load F ni (s j ) and the friction force F ti (s j ) can be given by where l is the friction coefficient and N 0 denotes the initial normal load. Note that the contact forces are the function of normal (k n ) and tangential (k t ) contact stiffnesses; hence, its estimation is of vital interest. However, literature research shows that no reliable method is available for the contact stiffness estimation, particularly for the shroud joint. Consequently, the experimental results in Refs. [41][42][43] are referred to and chosen, as displayed in Table 1. Additionally, the selection of contact parameters should also consider the iteration efficiency while solving nonlinear algebraic Eq. (24) with the alternating frequency/time method.
Due to the existence of shroud angle a and pre-twist angle b(L b ) (see Fig. 6b), the local coordinate system of the shroud interface is not always identical to that of the blade tip. Consequently, the relative displacements between the (ith) and (ith ? 1) blades at the shroud position should be transformed into the local coordinate system of the shroud interface, which can be stated as where c = b(L b )-a is the angle difference between shroud interface and pre-twist direction; thus, the shroud interface and pre-twist surface are parallel for c = 0°. v i , v i?1 , and w i , w i?1 are the displacements in flexural and swing directions, respectively. Similarly, the normal impact force and tangential friction force of the arbitrary ith blade should also be transformed into the local coordinate system of the blade at shroud position. Therefore, it can be described as follows: where the subscript (i) represents blade number; thus, F nl,i is the nonlinear impact-friction force of the ith blade.  Fig. 6 Schematic of coupling behavior between adjacent shrouded blades: a the impact-friction model of a shroud joint, b a full-2D impact-friction model periodic relative displacement of shroud interfaces. Nevertheless, the nonlinear algebraic Eqs. (24) are developed in the frequency domain using MHBM; thus, a common AFT method [21] is implemented to solve the limit. The full cyclic of the AFT flow diagram is displayed in Fig. 7 and can be mathematically summarized as where FFT and iFFT denote the forward and inverse discrete fast Fourier transform.
At the first stage of the AFT procedure (time domain), the time-discrete damper motion and nonlinear contact forces are determined as formulated in Eqs. (27)-(30) using the time-discrete relative displacements. Then the discrete nonlinear forces are converted into the frequency domain by the FFT algorithm, obtaining corresponding harmonic coefficients of F nl,k (k = 0, 1,…, N h ) in Eq. (25). At the second stage of the AFT procedure (frequency domain), the harmonics of F nl,k of contact forces are employed to solve the algebraic Eqs. (25). Consequently, a new estimate Q k = aW k (k = 0, 1, …, N h ) is determined by using MHBM. Then new solved Q k (k = 0, 1, …, N h ) are converted into the time domain using the iFFT algorithm, finishing a full cyclic of the AFT flow diagram.
Additionally, the assumed mode method (AMM) is implemented during the model formulation stage, shifting the physical displacement vector, e.g., [v i , w i ] of contact DOFs to the modal displacement vector under the generalized coordinate system. Therefore, to compute the relative displacements of contact DOFs between the adjacent blades and evaluate the nonlinear force, an algorithm of displacement vector transformations between physical and generalized coordinate systems is taken account into. As depicted in Fig. 7, when a new estimate Q ¼ Q 0 ; ...; Q N h À Á (k = 0, 1, …, N h ) is determined, it is converted into the physical coordinate system. The specific expression can be stated as where Q v i and Q w i are the harmonic coefficients under the physical coordinate system. H v i and H w i are the conversion vectors of flexural and swing directions:  by FFT under the physical coordinate system, a coordinate inversion algorithm is necessary: where F nl ¼ F nl;0 ; . . .; F nl;N h À Á are the harmonics of nonlinear forces under the generalized coordinate system, which are implemented to solve the eigenproblem in Eq. (24).

Continuation method
Based on the previous description, the nonlinear modal algebraic equations can be summarized as: For solving the numerical solutions, the classical Newton-Raphson method is also employed, and the arbitrary iterative step of algebraic Eqs. (28) can be derived as where X (p?1) and X (p) represent the Fourier coefficients of modal solutions at the pth and pth ? 1 step, respectively. R(X (p) ) is the residual vector at the pth step. qR(X (p) )/qX (p) is the analytical expression of Jacobian matrix at the pth step, which can significantly improve the solution efficiency. Wherein obtaining the derivation of nonlinear impact-friction force F nl with respect to an arbitrary variable x is the main challenge, which can be addressed by the AFT method and stated as, where detail derivations of the nonlinear impactfriction force can refer to Ref. [44].
During the prediction stage of modal characteristics, the full-2D impact-friction force with an initial gap or preload of shroud interfaces can induce the turning point bifurcation, limiting the use of Newton-Raphson method. For overcoming this convergence limitation, a predictor-corrector continuation procedure of the arc-length method, applied by Ferreira [22], is used to compute the solutions of the nonlinear algebraic equation system. Due to the dependence of modal results on modal amplitude a, the continuation method is applied for an interval of [a min , a max ], as illustrated in Eq. (36).

Modal synthesis
The nonlinear modal results within the vibration level of interest including the eigenfrequencies x 0 , modal damping ratio f, and eigenvectors W can be computed based on previous methods. It should be noted that: on the one hand the impact-friction nonlinearity will result in multi-harmonics of the nonlinear modes, different from the linear system; on the other hand, the eigenvalue problem in Eq. (36) is under the generalized coordinate system due to the appliance of AMM. Therefore, these factors necessitate the synthesis procedure for mode shapes of the cycle blades with shroud joints.
where W vi and W wi are the nonlinear eigenvectors of the ith blade under the physical coordinate system.

Dynamic analysis of the shrouded blade system
In the following sections, the dNNMs of a cyclic blades system with shroud interfaces are presented and analyzed. Due to the existence of the nonlinear impact-friction model, one of the most important parameters is the retained number of harmonics N h , which involves the convergence of the nonlinear solutions. Thus, for obtaining accurate mode properties, the results on the subject of the forced response convergence as a function of the retained number of harmonics N h is reported in Fig. 8. And a traveling excitation of Eq. (20) for F ex = 2000 N and model parameters in Table 1 are performed. Note that the forced response of the blade tip rapidly converges toward a stable value for N h C 15. Since the retained harmonic terms play a key role to determine the nonlinear dynamics accurately, fifteen harmonics are retained to keep a compromise between accuracy and computational efficiency in the following investigations.

Nonlinear modal properties of cyclic blades with shroud interfaces
The steady-state response of the cyclic blade system under a harmonic force in the frequency range around the eigenfrequency of blade bending in flexural direction was investigated for varying normal load N 0 (N 0 = k n Ágap) in the shroud joint. In the following analysis, all numerical results are computed for X = 300 rad/s, and the effects of centrifugal stiffening, and spin softening due to rotating speed on the modal properties are discussed in the next section. Firstly, the frequency response curves of an initial gap case (gap = -2 mm) are illustrated in Fig. 9a under a wide range of excitation levels. Note that the excitation frequency is normalized with respect to the first eigenfrequency of blade bending under free state. The results indicate that the tip displacement of blade bending increases with the varied excitation level. When the displacement of the blade tip is higher than the initial gap of the shroud interface, the frequency response curves bend toward  the right, inducing a multi-valued response for some of the curves. This phenomenon is generally called overhanging branch, representing the local unstable region. That's why a frequency multi-harmonic balance method is applied in this thesis. Moreover, the backbone curves, which correspond to the frequency-energy characteristic plotted in the dotted line, show good agreement with the maxima of the amplitude-frequency response curves. These phenomena can be attributed to the impact-friction behavior of the shroud interface for a high response amplitude of the rotating blade. Next, the frequency response of the blade tip for preload case (gap = 2 mm) was studied under an increasing excitation level. As depicted in Fig. 9b, the shroud interface is permanently sticking for a low excitation level; thus, the resonance frequency keeps constant, resulting in a linear behavior. However, for high excitation levels, contact DOF can slip freely and the softening effect deduced from the modal properties is exhibited in the forced response. Therefore, amplitude-frequency response and backbone curves gradually bend toward the left and exhibit overhanging branches, which are locally unstable.

Initial gap case
Previous investigations indicate that the dynamics of the cycle blade system with shroud joints greatly depend on the normal load of the contact interface. Thus, understanding the impact of normal load N 0 on nonlinear mode characteristics is crucial in the optimization design of the shroud joint. Firstly, a detailed investigation focuses on the initial gap (N 0 \ 0) of shroud interface on the evolution of dNNMs of the system, as reported in Fig. 10. In the case of low vibration level, there exists a linear vibration region, where the shroud interface is permanently liftoff and no energy is dissipated. Thus, the normalized eigenfrequency of the first bending mode maintains constant at 1, and the modal damping ratio is equal to that of the free state. As the vibration level rapidly increases, the shroud interface experiences a contact state, i.e., impact-friction behavior occurs. Therefore, the modal damping ratio significantly increases. In addition to that, the blade system becomes stiffer, inducing an increase in the eigenfrequency of the blade's first bending mode. As soon as the vibration level tends toward infinity, the modal characteristics rapidly converge toward a stable value, agreeing well with the values of initial gap = 0 mm. Additionally, note that a decreased process of modal damping ratio for higher vibration level occur, which can be explained as: For a viscous damping source, the dissipated energy grows quadratically with amplitude, which leads to a constant modal damping ratio. The energy dissipated in the dry friction contact, in contrast, increases only linearly with amplitude, thus leading to a decreasing modal damping ratio for large vibration amplitudes.
For better understanding modal characteristics of the cycle blade system with shroud interfaces, the 3D mode shapes and backbone curve for initial gap = -  Fig. 11. In the case of point A, note that the eigenfrequency remains in the linear regime. Thus, the 3D mode shape of the cycle blade system is identical to that of a linear system. However, for the nonlinear regime of point B, the modal deflection shape of the rotating blade exhibits a peculiar form, which distinctly differs from the first linear bending mode. This indicates that the contact behavior of shroud joints will make the 3D mode shape of the system significantly deviate from the linear case. Next, the time series, phase portraits [45][46][47][48], and frequency contents for the blade tip deflection are illustrated in Figs. 12 and 13. Note that, for linear point A, the frequency content is only predominated by the fundamental harmonic; hence, the phase portrait of the blade tip exhibits a smooth elliptic orbit. Times series show that the shroud interface is always separated during one vibration period. However, for the case of nonlinear point B, the time series exhibit apparent nonlinear features, and one vibration cycle experiences a stick, slip, and liftoff state, respectively. Therefore, the 'wrinkles' phenomenon occurs in the phase portrait, which indicates the presence of super-harmonics. This result can be proved by the frequency content, as plotted in Fig. 13c. This also demonstrates that the changes Fig. 11 a Backbone curve, b 3D mode shape, c mode shape comparison of the shrouded blade system

Preload case
Next, the variation of modal damping ratio and eigenfrequency of the blade first bending mode versus    Fig. 14. Numerical results indicate that dNNMs of shrouded blades with fixed system parameters completely depend on vibration level, which can quantitatively describe the kinetic energy of the system. In the case of small energies, the varied normal load of the contact interface does not exceed the initial preload N 0 . Thus, the shroud interface permanently remains sticking, causing a linear regime where the modal frequency and modal damping ratio maintain constant. Beyond a critical value of energy level, the shroud interface experiences liftoff or contact states during one cycle of oscillation, exhibiting obvious nonlinear phenomena henceforth. As a consequence, the eigenfrequency decreases with energy and gradually deviates from its corresponding value in the linear regime, cf. Fig. 14a. Meanwhile, the modal damping ratio also significantly increases with the vibration level, as depicted in Fig. 14b. Finally, the modal characteristics will rapidly converge toward a stable value of preload N 0 = 0. In addition to that, a larger preload N 0 needs a higher energy level to converge to a stable value.
Moreover, the smaller the preload N 0 is, the more the modal damping ratio increases. This finely explains the results of forced response in Fig. 9.
Next, the 3D mode shape, times series, phase projection, and the frequency content of the blade tip are plotted in Figs. 15, 16 and 17 for interface gap = 2 mm, i.e., the preload case (N 0 [ 0). In the case of linear point C, the shroud interface undergoes a sticking state; thus, the 3D mode shape of the cycle blades system significantly differs from that of the free blades system. Its time series, phase portrait, and frequency content in Fig. 16 indicate that the periodic motion of the shrouded blade is only dominated by a fundamental harmonic component. For the case of high vibration level (point D), the shroud interface will occur slipping, even liftoff states. Consequently, as plotted in Fig. 17, 'wrinkles' and super-harmonics components appear in the forced response of the blade tip. Furthermore, the 3D mode shape of the system exhibits nonlinear characteristics with larger values of modal deflections. Moreover, its modal deflections are greatly higher compared to the linear regime (sticking state). These findings highlight that the constraint capability of the shroud joint will be especially affected by the motion state of contact interfaces.

Effects of centrifugal stiffening and spin softening on modal properties
One of the contents of this article is to exploit the influence mechanism of the centrifugal stiffening and spin softening factors on the nonlinear mode properties of the shrouded blade system. Thus, after considering inter-blade coupling by impact-friction behavior, centrifugal stiffening, and spin softening effects, the present study developed a mechanism model in Sect. 2. And the forced response of blade tip for varied rotational speed under F ex = 2500 N is computed to reveal the influence of centrifugal stiffening and spin softening effects, as depicted in Fig. 18. Note that increasing rotational speed will encourage the increase of resonance frequency under both initial gap cases and preload cases. However, the variation law of resonance amplitude along with rotational speed is associated with the initial load of the shroud interface. For the initial gap case (see Fig. 18a), the resonance amplitude of the shrouded blade increases firstly and then decreases gradually with rotational speed. However, in the case of preload case (see Fig. 18b), the higher the rotational speed is, the smaller the forced response decreases, suggesting that stronger centrifugal stiffening and spin softening effects will stiffen the shrouded blade. That also explains why the occurrence of nonlinear phenomena needs a larger vibration level for higher rotational speed.
The results reported in Fig. 19 highlight the impact of centrifugal stiffening and spin softening effects on the modal characteristics of the shrouded blade system. The simulation parameters are taken as gap = -2 mm or gap = 2 mm, and the shroud is at the tip of the blade, and other parameters are displayed in Table 1. As plotted in Fig. 19a and b, the modal frequencies of the shrouded blade significantly increase with rotational speed. Meanwhile, the influence laws of centrifugal stiffening and spin softening on modal frequencies are identical for different modal amplitudes. Besides, there exists an amplitude interval, where the variation of modal frequencies versus rotational speed maintains constant. This suggests that when the shroud interface keeps in a separation or sticking state, the influence laws of centrifugal stiffening and spin softening effects on modal frequencies will not change over the vibration level of the shrouded blade.
The results in Fig. 20 are computed under different rotational speeds to understand the impact of centrifugal stiffening and spin softening on the modal damping ratio of the shrouded blade system. In general, the damping ratio decreases with the rotational speed. This can be attributable to the fact that the modal damping ratio will especially depend on the stiffness of the shrouded blade. Consequently, in the case of low rotational speed, the modal damping ratio decreases slowly, whereas it decreases significantly with rotational speed under higher rotational speed.
This indicates that the influence of centrifugal stiffening and spin softening effects on modal damping ratio is more sensitive to high rotational speed, identical to the variation law of eigenfrequency. Moreover, for low rotational speed, the changing range of modal damping ratio versus modal amplitude is larger. These results mean that shroud interfaces exhibit better damping effects under the same conditions. This phenomenon can also explain the results in Fig. 20 that why the blade needs a higher blade vibration level to activate overhanging branches under larger rotational speed.

Critical speeds of cyclic blades with impactfriction coupling
In engineering practice, the rotating blade will experience the action of aerodynamic excitation, which significantly affects the dynamic characteristics of the shrouded blade. And its excitation frequency is generally directly proportional to the rotational speed, which determines the centrifugal forces. Therefore, the present analysis defines a mono-harmonic with a constant phase lag between adjacent blades to simulate the aerodynamic excitation. While the external excitation frequency enters the frequency range near the considered mode, the vibration level of the shrouded blade will conspicuously increase, inducing fatigue damage. In general, this operation condition is called critical speed (CS), which should be avoided and paid enough attention to during the design stage of the shrouded blade. Therefore, to interpret the variation mechanism of the CSs, the CSs are obtained by the intersection of the frequency loci in Campbell   shrouded blade increase, and rapidly converges toward a stable value. Besides, at a low rotation speed, the system exhibits a larger variation of the modal frequency. Then the CSs are obtained by the intersection of the frequency loci in Campbell diagrams and the dot-dash line x = 3X, as plotted in Fig. 21b. Note that the variation law of CSs versus modal amplitude is similar to that of eigenfrequency (see Figs. 10a and 14a). There exists a modal amplitude interval where the CSs increase significantly, and then rapidly converge toward a stable value. Here, the variation range of CSs is called the sensitive rotation speed region. Similar phenomena are also observed in the preload case, as illustrated in Fig. 22. Previous results indicate that the vibration energy governs the frequency loci in Campbell diagrams, which should be fully considered during the design stage.
Then the variation of CSs versus tip amplitude for different interface gaps (gap = N 0 /k n ) is computed in Fig. 23. The vibration level governs the frequency loci in the Campbell diagram. Therefore, it can be seen that the CSs results in Fig. 23 are in accordance with the variation law of eigenfrequencies. Additionally, a larger interface load will narrow the critical speed range.
Finally, the effects of interface load of shroud joints on the critical speeds (CSs) are further presented in Fig. 24 under tip amplitude q(i norm ) = 10 mm. Note that the Campbell diagrams of the shrouded blade under different interface loads are displayed in Fig. 24a. In general, the frequency loci rapidly increase with the initial load of the shroud interface. However, when the absolute value of the interface load is higher than the tip amplitude, the frequency loci maintain constant. As introduced previously, this can be explained by the fact that the shroud interface keeps in a sticking or separation state. Hence, as depicted in Fig. 24b, the variation process can be divided into three regions. During regimes I and III, the CSs keep in constant for liftoff or sticking states. In the case of regime III, the variation of CSs is sensitive and completely depends on the vibration level of the shrouded blade, which should be paid enough attention and well designed.

Conclusion
In the present study, a coupling model of cyclic blades system with shroud interfaces at the blade tip is developed to explore its modal characteristics induced by the impact-friction of adjacent shrouds. The developed model considers the participation of the centrifugal stiffening and spin softening effects, the inertia effects of cyclic blades, and nonlinear coupling effects between adjacent blades. After considering these factors, dNNMs of the system are investigated under two typical initial states (initial gap case, and initial preload case) through a mixed modal analysis method (MMAM). Then the variations of dNNMs, including eigenfrequencies and modal damping ratios, versus rotational speed are discussed to exploit the influence of centrifugal stiffening and spin softening effects on nonlinear modal characteristics of the cyclic blades under impact-friction coupling. Finally, the Campbell diagram of the system is plotted, and the critical speeds (CSs) of the system under different vibration levels and initial interface loads are obtained and compared to interpret its variation mechanism. According to the analysis results, the main conclusions are as follows: (1) Initial interface load and vibration energy govern the dNNMs of the cyclic blades system. And there always exists a linear energy regime, where the eigenfrequency and the modal damping ratio keep constant due to the permanent stick or liftoff. With the increase of vibration energy, the eigenfrequency and modal damping radio of the shrouded blade will experience a great change and then rapidly converge toward a stable value. Moreover, the friction contact can also induce higher harmonic vibration content, making the 3D mode shape of cyclic blades significantly deviating from the linear case. (2) The centrifugal stiffening and spin softening can lead to a significant influence on the modal frequencies and modal damping ratio of a shrouded blade. At low rotational speed, the variation of the Campbell diagram of dNNMs versus vibration level is more conspicuous. The modal damping ratio of a shrouded blade is more sensitive under a higher rotational speed. Moreover, for a low rotational speed, the modal damping ratio and its variation range versus vibration amplitude are larger; thus, shroud joint exhibits better vibration suppression ability. (3) The variation law of CSs versus modal amplitude is similar to that of eigenfrequency, and mainly depends on the vibration energy and initial interface load. There exists a sensitive interval of vibration amplitude where the CSs increase significantly, and then rapidly converge toward a stable value. Moreover, a larger interface load will narrow the critical speed range. 5 ðA:5Þ W T Wdx ðA:8Þ Table 2 Area integrals for kinetic and potential energy expression R R Ab dg dn ¼ A b R R Ab g 2 dg dn ¼ I z R R Ab n 2 dg dn ¼ I y R R Ab g 2 þ n 2 À Á dg dn ¼ I a   ðA:20Þ w T w w w cos 2 bdx ðA:26Þ