Computation and investigation of mode characteristics in nonlinear system with tuned/mistuned contact interface

This study derived a novel computation algorithm for a mechanical system with multiple friction contact interfaces that is well-suited to the investigation of nonlinear mode characteristic of a coupling system. The procedure uses the incremental harmonic balance method to obtain the nonlinear parameters of the contact interface. Thereafter, the computed nonlinear parameters are applied to rebuild the matrices of the coupling system, which can be easily solved to calculate the nonlinear mode characteristics by standard iterative solvers. Lastly, the derived method is applied to a cycle symmetry system, which represents a shaft-disk-blade system subjected to dry friction. Moreover, this study analyzed the effects of the tuned and mistuned contact features on the nonlinear mode characteristics. Numerical results prove that the proposed method is particularly suitable for the study of nonlinear characteristics in such nonlinear systems.


Introduction
Shaft, disk and blade are the key components of rotating machinery. In operating conditions, the deformation and vibration of subcomponents are constantly encountered and grouped. Many studies have analyzed the dynamic couplings among the shaft, disk and blades of a linear system [1][2][3][4][5][6]. Moreover, some studies [7][8][9][10] have adopted pre-twisted, thin-walled rotating blades to analyze their nonlinear vibration characteristics under different excitation condition. However, friction dampers are commonly designed and applied to attenuate the response levels and prevent high cycle fatigues. Consequently, the coupling behavior of structures that are made of components assembled by means of joints may be highly nonlinear. Over the years, many scholars have focused on the dynamic behavior of blades with friction damper.
The most common friction damper in a blade system is the blade-disk interface. Petrov and Ewins [11] developed an approach to analyze the multi-harmonic forced response of large-scale finite element modes of bladed disks by considering the nonlinear forces acting at the contact interface of the blade roots. Thereafter, the multi-harmonic vibrations for systems with friction and gaps based on analytically derived contact interface elements were analyzed using a proposed approach [12]. Ciğeroğlu and Özgüven [13] proposed a multi-degree-of-freedom (multi-DOF) model of bladed disk system subjected to dry friction dampers for the efficient vibration analysis of turbine blades with dry friction. Peeters et al. [14,15] applied a shooting method to calculate the nonlinear normal modes of a system with cyclic symmetry, thereby exposing the similar and nonsimilar normal modes and localization phenomena for some nonlinear normal modes. Zucca et al. [16] presented a method to compute the friction forces that occur at the blade root joints to evaluate their effects on blade dynamics. Another research used the harmonic balance method (HBM) to conduct a parameter study of the non-linear aero-elastic phenomena of a bladed disk for aeronautical application in the presence of friction contact by a one-way coupled method [17]. Li et al. [18] used a finite element model (FEM) to propose a dynamic model to analyze the nonlinear characteristics of a flexible blade with dry friction. Joannin et al. [19] introduced a novel reduced-order modeling technique well-suited to the study of nonlinear vibrations in large FEMs. Apart from FEM, the lump parameter model was also established to study the nonlinear vibration of bladed disks with dry friction dampers of blade roots [20]. Joannin et al. [21] computed the steady-state forced response of nonlinear and dissipative structure by presenting an extension to classic component mode synthesis methods, which was proven by previous mistuned cycle model.
Underplatform damper is another friction damper that is commonly used to reduce vibration amplitude. Petrov and Ewins [22] developed an advanced structural model for wedge and split underplatform dampers (UPDs) and proposed and realized an approach for using the new damper models in the dynamic analysis of large-scale FEMs of bladed disks. Firrone et al. [23] proposed a novel method to compute the forced response of blade disks with UPDs. Berruti et al. [24] designed a static test rig called "Octopus" to validate the numerical model and its nonlinear dynamic response of a bladed disk with UPDs. Zhang et al. [25] described an efficient method to predict the nonlinear steady-state response of a complex structure with multi-scattered friction contacts. Thereafter, blades with UPDs were used as an example to validate the approaches by calculating the steady-state response of an FEM with numerous DOFs. Pesaresi et al. [26] used an updated explicit damper model as basis to perform a nonlinear analysis and evaluated the results against a newly developed UPD test rig. Pesaresi et al. [27] used a modified Valanis model as basis to propose a 3D microslip element to describe the contact interface. This new approach allows to implicitly account for the microscale energy dissipation and the pressure-dependent contact stiffness caused by the roughness of the contact surface. In a mechanical system, the subcomponents are often connected by a single or multiple friction contact interfaces, thereby resulting in a nonlinear coupling system. To compute the nonlinear coupling mode characteristics, this study proposed an algorithm to obtain the coupling mode information of such a mechanical system. Thereafter, the representative shaft-disk-blade (SDB) coupling structure of this mechanical system was used as an example to adopt the derived method to obtain and investigate the nonlinear coupling vibration characteristics of the lumped SDB system with tuned/mistuned contact interfaces. The algorithm used in this research is well suited in analyzing the coupling vibration characteristics of such a mechanical system, thereby easily and efficiently obtaining the nonlinear modes similar to a linear system. Lastly, the effects of the tuned/mistuned contact interface on the nonlinear coupling mode characteristics of the cycle SDB structure was discussed in detail.

Mathematical formulation
A lumped-parameter model is devised to study the coupling vibration characteristics of the SDBs subjected to dry friction at moderate computational cost. This model consists of several sectors (see Fig. 1), all of which are made of two DOFs to account for the blade body and blade root and one DOF to account for the disk. For each sector, the disk and blade root are connected by tenon-mortise. An ideal dry friction model is adopted to simulate the contact feature. In this model, M s , M i , m i1 and m i2 stand for the mass of shaft, ith disk, ith blade root and ith blade body, respectively. Moreover, k di and k bi correspond to the disk stiffness between the ith and (i + 1)th disk elements and blade bending stiffness of ith blade, respectively. In Fig. 1 Diagram of the lumped SDB model. addition, k sdi is the coupling stiffness of the shaft and disk. In the rotating SDB system, blades are often installed on the disk through contact interface. Consequently, the contact behavior is simplified and described using the nonlinear damping and nonlinear stiffness of the dry friction force. In particular, K ni and C ni represent the nonlinear damping and nonlinear stiffness of the dry friction damper, respectively.
The kinetic energy of the shaft-disk system (T r ) can be presented as follows: where n is the sector number of the disk and x s and x di are the displacements of M s and M i , respectively. The kinetic energy of the shaft-disk system (U r ) is provided as follows: where k s is the shaft stiffness. In Fig. 1, the displacement of the ith blade's root and tip in a global coordinate system can be expressed as follows: where x s is the displacement of the shaft (M s ), x di is the local displacement of m di with respect to M s , x i1 is the local displacement of the blade root (m i1 ) with respect to m di , and x i2 is the local displacement of the blade tip (m i2 ) with respect to m i1 . Therefore, the kinetic and potential energies of the blades can be presented as follows: where x i1 and x i2 are the displacements of the ith blade root and ith blade body, respectively. By substituting the energy expressions into the Lagrange equations, the differential equation of the vibration in the matrix notation is provided as follows: where M SDB , C SDB , K SDB , and q SDB are the mass matrix, damping matrix, stiffness matrix, and generalized coordinate vector, respectively, F e is the external excitation, and F D is the friction force vector. The meaning of the parameters in the formula is provided in the Appendix.

Harmonic excitation
Harmonic excitation is commonly used to simulate aerodynamic excitation. The excitation naturally arises from the static parts in the gas flow. The harmonic excitation acting on ith blade body is provided as follows: where F is the amplitude of the excitation, ω is the excitation frequency, N b is the blade number, and d is the Dirac delta function.

Dry friction model
In this analysis, an ideal dry friction model is established to simulate the nonlinear force of contact interface. The contact interface diagram is illustrated in Fig. 2. A timediscrete friction contact model is employed and the specific derivation of the contact force is introduced [28]. The two linear spring k t and k n in the tangential and normal directions, respectively, are used to model the local contact stiffness of the contact surface. Coulomb friction law is used as basis to derive the tangential force f(t) and normal force N(t) by considering the relative tangential (u(t)) and normal (v(t)) displacements, respectively, of the contact DOF. After using a predictor step, the value of the tangential contact force (T(t)) under the sticking condition is as follows: where Dt is the time step, w(t) is the slider displacement, and k t is the contact stiffness in the tangential direction. Thereafter, the actual value of tangential contact force is computed accordingly as follows: where N(t) is the positive pressure of the contact surface. Moreover, the slider displacement can be further obtained as follows: By decomposing the dry friction force of Eq. (9) with a first order Fourier, the dry friction force can be written as follows: where k eq (A r ) and c eq (A r ) are the nonlinear stiffness and damping, respectively, of the dry friction and A r is the response amplitude of the contact DOF.

Algorithm of the nonlinear analysis
Nonlinear modes differ from the linear analog because the former is energy-dependent. Consequently, the computation of nonlinear modes results in more computational effort than the linear system. This study uses the incremental HBM (IHBM) as basis to derive an algorithm to efficiently yield the mode characteristics. This section provides an overview of the main steps of the calculation. The specific procedure is as follows. First, using a new time scale τ = ωt, where is the angular frequency, Eq. (6) is transformed as follows: and No is the number of degree freedom. Second, let Q 0 and ω 0 denote the state of the vibration by adding the increment DQ, Dω is as follows: By substituting Eq. (13) into Eq. (12) and disregarding all small terms of the higher order, the incremental equation is as follows: where C n and K n are the Jacobin matrices, and C n ¼ ∂FðdQ=dτ,QÞ=∂ðdQ=dτÞ and K n ¼ ∂FðdQ=dτ,QÞ=∂Q. Thereafter, let ½Δa in cosðnτÞ þ Δb in sinðnτÞ: (16) By substituting Eqs. (15) and (16) into Eq. (14), the Galerkin method is used to obtain the algebraic equation, which can be written as follows: where DA =[Da 10 , Da 11 , Db 11 , Da 12 , Db 12 , …, Da m2 ]. When the residual of DA is zero, A, C n , and K n approach the exact value. The matrices and vectors in Eq. (17) depend on the initial value of A 0 and ω 0 . Accordingly, ω is set as the active incrementing parameter. On the basis of Eq. (17), the response amplitude of the coupling system can be obtained. Lastly, by substituting the equivalent stiffness K n and damping C n into the coupling matrices, the damping and stiffness matrices in Eq. (6) can be rewritten as follows: where K n and C n are the nonlinear stiffness and damping matrices, respectively, assembled by the K ni and C ni of the different contact DOF of the blades. According to the vibration theory, it can be assumed that the solution of Eq. (6) is of the form q SDB ¼ ηe lt , then Eq. (6) can yield the following equation: where η is the undetermined coefficient vector and l is the The mode characteristics can be obtained by solving the corresponding equations. Thereafter, the nonlinear coupling mode characteristics and coupled vibration of the SDB structure can be further studied similar to a linear system.
In Fig. 3, the response amplitude A can be obtained by using IHBM to compute the nonlinear stiffness K n and nonlinear damping C n . This section uses a single sector to compute the amplitude-frequency response of the blade root, nonlinear stiffness and damping of the dry friction contact feature. In Fig. 4, the nonlinear stiffness and nonlinear damping under different excitation frequency and rotational speed are plotted with the parameters listed below: Contact parameters: The numerical results indicate that nonlinear stiffness and nonlinear damping are constant when DOF of the blade root is in a stationary sticking state and gradually shift toward a significant change when slipping occurs. The effect of rotational speed on nonlinear damping and stiffness reveals the difficulty of changing for high rotational speed. This phenomenon can be explained through the hard activation of nonlinearity under high rotational speed case. Therefore, the nonlinear stiffness and damping of the contact surface remain constant.

Numerical results and discussion
The contact force can be calculated approximately using the predictor-corrector strategy [28]. This study aims to investigate the coupling vibration characteristics under the effect of contact interface. Consequently, the variation of the normal force is assumed to be constant. Moreover, the normal force N is generated by the centrifugal force of the blade system and the specific expression is This section summarizes the coupling modal properties by using the algorithm provided in the previous section based on a lumped SDB model with the 12 previously derived identical sectors. The results consist of the natural frequencies and mode shapes in the stick-slip areas, respectively. Moreover, the variation of the coupling modal characteristics versus the excitation frequency is highlighted. The model parameters in this study are illustrated in Table 1.

Tuned system
The forced response and mode characteristics are computed in the case of the tuned system. The tuned SDB model with contact interface in this study is an ideally periodic structure, thereby enabling all sectors to show the identity characteristics. Therefore, only the results of Sector 1 are shown in this section. In this study, the coupling dynamic behavior of contact interface depends on a few key parameters, namely, rotational speed, contact parameters, excitation level and load frequency. Therefore, the following study selects some representative parameters as examples for the mechanism analysis of the coupling vibration. Figure 5 illustrates the amplitude-frequency response curves of the blade in Sector 1 under different excitation levels for W = 200 rad/s. Further analysis of the nonlinear characteristics points out that the resonance frequency is constant at low excitation level. When the excitation level  is sufficiently high, the nonlinearity is activated and the resonance frequency gradually shifts toward a low frequency. This result indicates that the contact interface in this study is softening nonlinearity. Figure 6 shows the variation of the nonlinear stiffness and nonlinear damping versus excitation frequency. In Fig. 6, nonlinear stiffness and nonlinear damping remain constant at low excitation levels. The values of the nonlinear stiffness and nonlinear damping are equal to that of a linear system with frictionless interface. However, complex mode characteristics may occur during the stickslip transitions when nonlinearity is activated at high excitation amplitudes. In the slipping area, the nonlinear stiffness of the contact interface initially decreases and eventually increase to the value of a linear system with frictionless interface. However, nonlinear damping exhibits an opposite change rule. Although the friction DOFs are in a sipping state, nonlinear damping significantly increases. Moreover, a broad frequency range of slipping state exists at high excitation levels.
The variations of the natural frequencies of the blade body DOFs and blade root DOFs against the excitation frequency at F = 500 N are plotted in Fig. 7. Note that the coupled natural frequencies are constant and equal to those of the underlying linear system with bonded DOFs at the contact surface. However, the natural frequencies of the coupling system decrease significantly with excitation frequency when frictional DOFs are in a slipping state. This result indicates that nonlinearity has a significant impact on the natural frequencies. Moreover, multiplicity of the coupling frequencies is present. Further analysis is conducted on the nonlinear mode shapes of the cyclic structures. The mode shapes of Points A (ω = 200 Hz) and B (ω = 250 Hz) (see Fig. 7), which are in the stick-and slipareas, respectively, are plotted in the following section. The natural frequency of the blade in Sector 12 exhibits a conspicuous difference with other blades. The single frequency reveals that it is associated with the shaft mode. Therefore, the coupled effect can ideally account for the distinct natural frequency.
The natural frequencies of the coupling system at ω = 200 Hz are shown in Table 1 and the mode shapes, which are predominated by DOFs of the blade body, are exhibited in Fig. 8. The natural frequencies of the coupling SDB  system can be divided into three frequency ranges, namely, blade body frequency, blade root frequency and disk frequency ranges. The natural frequencies of the blade body show inconspicuous difference, except for 262.5346 Hz. This result can be explained by the coupling vibration among the shaft and blade. However, the wide frequency range of the disk can be explained by the strong coupling effects among the subcomponents. Furthermore, repeated frequencies occur and the single frequencies are induced by the coupling vibration of the shaft.
For a clear interpretation of the nonlinear characteristics of the coupling system, the mode shapes for the sticking and slipping states are given in Figs. 8 and 9, respectively. The corresponding mode characteristics of the tuned system can provide a reference for the mistuned system in the following analysis. Therefore, the mode property for  the sticking state (ω = 200 Hz) is first computed and illustrated (see Fig. 8). In such a case, the nonlinear damping and stiffness of different sectors are constant and equal to one another. Given that no slipping occurs between the contact DOFs at the interface, the mode shapes are identical to those of a linear system with bonded DOFs at the interface. The majority of the mode shapes of the cyclic symmetric structure in the sticking state are in harmonic form and the subcomponents exhibit identical mode properties. This result indicates that the cyclic symmetric structure exhibits a synchronous vibration phenomenon. The results indicate that disk flexibility does not induce the split of the contact interface in different sectors. The majority of the modes are repeated except for the last two modes. This result can be explained by the flexibility of the shaft.
Further analysis of the nonlinear coupling mode is shown in Fig. 9 for ω = 250 Hz (Point B in Fig. 7). Although the increase of the excitation frequency activates the sticking state of the contact DOFs, the distribution pattern of the coupled frequencies remains constant. However, the natural frequencies of the blade DOFs significantly decrease in the slip area, which is why friction is commonly referred to as a softening nonlinearity. However, the natural frequencies of the disk's DOFs exhibits a slight decrease owing to the strong coupling between the disk's DOFs. The mode shapes of the coupling system with respect to the slip-area (Point B in Fig. 7) are depicted in Fig. 9. Although the frictional DOFs are in a slipping state, the mode shapes remain similar to the stationary sticking state.
The computation of the nonlinear mode shows that the contact interface exhibits conspicuous influence on the coupling system. In a linear system, the natural characteristics are known to be constant. However, the natural characteristics of a nonlinear system completely depend on the motion of the contact DOFs. When the friction DOFs are in a stationary sticking state, the natural characteristics of the coupling system are identical to those of the underlying linear system with linear connection. Although the slipping state occurs, the natural frequencies of a nonlinear system significantly decrease and is lower than that of the sticking state. Therefore, a series of parameters, such as contact parameters, can alter the modal characteristics of a nonlinear system with the contact interface. Accordingly, the representative contact parameters are considered as examples for the mechanism analysis of the coupling vibration.

Mistuned system
This section presents and discusses the effects of random mistuned contact interface on the nonlinear coupling mode. In the previous section, the mode characteristics of the cycle structure under different excitation levels have already been analyzed. Therefore, only the case of F = 500 N is used as an example to analyze the mode characteristics of the mistuned system. Accordingly, the mistuned contact interface is added to the model through the mistuned contact parameters. For a given contact parameter, the random mistuning parameter f m is obtained by the addition of a random mistuning error ε taken in the normal distribution ε~(0, σ 2 ). To provide a clear overview of the mistuned mode characteristics, the standard deviation of mistuning error is σ= 5%, which may be exaggerated in contrast to the practical case. Therefore, the mistuned parameter is as follows: where f is the tuned value of the contact parameter.

Mistuned friction coefficient
The influence of the friction coefficient mistuning is presented and investigated. The specific values of the mistuned friction coefficient are given in Table 2. Thereafter, the amplitude-frequency response curves of the different sectors are plotted in Fig. 10 using the algorithm presented in Section 4. Compared with the tuned system, the forced response of the blades no longer exhibits a single resonance peak. The multiplicity of the blade response can be explained by the mistuned contact interface. However, the split of the response only occurs in the slipping area. This result indicates that the random mistuned friction can only take effect when the friction DOFs are in a slipping state.
Further analysis of the nonlinear parameters of the contact interface is plotted in Fig. 11. In the stationary sticking state, the nonlinear stiffness and nonlinear damping of the different contact interfaces are constant and equal to one another, similar to a tuned system. However, the nonlinear stiffness and nonlinear damping of the contact interfaces in the different sectors lose their identity characteristics in the slipping state. The multiplicity of the nonlinear stiffness and nonlinear damping indicates a mistuning of the contact interfaces. Figure 12 shows the variation of the coupling frequency with the excitation frequency. Evidently, the changed rule of the natural frequency is identical to that of the tuned contact interface case.
A focus on the natural characteristics enables an improved understanding of the mode characteristics when ω = 200 Hz (Point C in Fig. 12). In the sticking state, the natural frequencies in Table 3 show no difference with those of the tuned system. The mode shapes, which are predominated by the blade body, are plotted in Fig. 13. The corresponding mode shapes show the same distribution pattern as the tuned system. Moreover, the majority of the modes occur repeatedly except for the last two modes. The equal deformation of the blade and disk DOFs reveals  that the single mode arises from the coupling vibration of the shaft. Thus, mistuning may have no impact on the natural characteristics of the cycle symmetric structure in the stationary sticking case. An analysis of the natural characteristics is performed for the slipping state (Point D in Fig. 12). The natural frequencies of the coupling system are given in Table 3. Note that repeated frequencies disappear and each DOF exhibits a peculiar frequency. The frequency split is caused by the difference in the nonlinear stiffness of the mistuning contact interface. Figure 14 illustrates the mode shapes predominated by the blade. In the mistuning system, the   mode shapes of the slipping state show a significant difference and the harmonic form disappears. The mode information exhibits an evident localization phenomenon. Therefore, the mistuned friction coefficient will result in the localization of the cycle structure when the contact DOFs are in a slipping state.

Mistuned contact stiffness coefficient
This section adds the random mistuned contact stiffness coefficient to the dry friction model in the different sectors. Table 4 illustrates the specific values of the contact stiffness in the different sectors. The forced response of the blade DOFs in the different sectors is plotted in Fig. 15 for F = 500 N. Unlike the mistuned friction coefficient case, only a slight difference in resonance amplitude occurs. However, the impact of the mistuned contact stiffness constantly exists in the overall frequency range. Moreover, note that the mistuning of the contact stiffness coefficient activates a new resonance peak in the frequency of the shaft.
The nonlinear stiffness and nonlinear damping of the contact interfaces are given in Fig. 16 for W = 200 rad/s and F = 500 N. The results show that the nonlinear stiffness splits into distinct values throughout the frequency range. However, the multiplicity of the nonlinear damping only occurs in the slipping area. Moreover, nonlinear stiffness and nonlinear damping exhibit a conspicuous shift in the frequency of the shaft. This result indicates that the random mistuning of the contact stiffness activates the transient slipping state in this frequency range. Figure 17 provides an overview of the variations of the natural frequencies.
The results show varied natural frequencies throughout the frequency range. For the frequency of the shaft, the natural frequencies of the blade DOFs exhibit a significant decrease. These results indicate that the mistuning of the contact stiffness has a conspicuous impact on the modal characteristics of the coupling system.
In this section, the natural characteristics of coupling system are investigated under the stationary sticking state (Point E in Fig. 17). Table 5 shows the natural frequencies of the coupling system when ω = 200 Hz. All DOFs show a single natural frequency, which is different from the mistuning case of the friction coefficient. A corresponding analysis of the mode shapes is plotted in Fig. 18. Although all contact interfaces are in stationary sticking, localization remains in the coupling system.
An overview of the natural characteristics of the coupling system in a slipping state (Point F in Fig. 17) is obtained for W = 200 rad/s and F = 500 N. Similar to the stationary sticking state, there exists identical mode characteristics except for the low natural frequencies. All DOFs exhibit a single and peculiar natural frequency that originates from the repeated frequency. Therefore, the repeated modes also disappear, which differ from the tuned system (see Fig. 19). Moreover, localization phenomenon also occurs in the mode shapes and is induced by the mistuning feature. The mistuned contact stiffness coefficient may constantly have an effect on the modal properties throughout the frequency range.

Conclusions
This paper uses IHBM as basis to propose a new method for calculating the coupling mode characteristics of a nonlinear system with contact interface. The proposed method is successfully applied to obtain the coupling mode characteristics of a cycle symmetry structure with tuned and mistuned contact interface. Accordingly, this method is proven to be beneficial and efficient, thereby enabling the direct determination of the nonlinear mode information  of different nonlinear structures with the contact interface. Moreover, the modal analysis of the cyclic symmetric structure with the contact feature is presented as follows. In the tuned system, all sectors show identical nonlinear characteristics in the cycle symmetrical structure. When contact DOFs are in a stationary sticking state, the nonlinear stiffness and nonlinear damping of contact interface remain constant. However, these nonlinear parameters vary significantly when the slipping state of the contact DOFs is activated. The majority of the modes occur repeatedly except for the mode coupled with shaft vibration. The mode shapes of the coupling system are in harmonic form, which is similar to a linear system with frictionless interfaces.
The mistuning of the contact interface shows a conspicuous impact on the nonlinear characteristics. The   mistuned friction coefficient can only affect the system's natural frequency in the slipping state. However, the effects of mistuned contact stiffness do not depend on the state of the friction DOFs. When the mistuned feature takes effect, the nonlinear stiffness and nonlinear damping of each contact interface exhibit a peculiar value that differs  from one another. The repeated modes in the tuned system