Local Sensitivity Analysis of Steady-State Response of Rotors with Rub-Impact to Parameters of Rubbing Interfaces

Local sensitivity analysis, which describes the relative importance of specific design parameters to the response of systems, is crucial for investigating dominant factors in optimal design. In this paper, local sensitivity analysis of the response of rotors with rub-impact to parameters of rubbing interfaces is carried out. The steady-state motion of the rotor is evaluated by a harmonic balance method and the sensitivity coefficients for every rotation speed over the speed range are derived analytically. Two classical models, including the Duffing oscillator and the gap model, are utilized to validate the accuracy and capability of the adopted methods and high accuracy is shown. Numerical investigations of sensitivities of steady-state response of rotors to parameters of rubbing interfaces are then carried out, based on a lumped Jeffcott rotor and a finite element model respectively. Conclusions are drawn that the response of rotors subjected to rubbing problems is more sensitive to initial clearance than other parameters of the applied friction model. With increase of initial gap, the response of rotors becomes more sensitive and the range of region subjected to rub-impact forces shrinks until the separation of rotor and stator.


Introduction
For the purpose of searching for optimal design parameters for rotating machinery, investigations of the dynamic characteristics of the system, such as stability analysis [1,2], reliability analysis [3] and optimal design analysis [4] based on specific criteria, have been performed by many scholars. Typically, a sensitivity analysis is also suggested when researching the dynamic characteristics of a system. Sensitivity analysis can provide quantitative information about to what extent does the responses of rotating system change when the design parameters varies from a base case, thus dominant impact factors can be found based on a restricted number of cases rather than enormous amount of data. In general, there are two different kinds of methods measuring sensitivities of the investigated system, i.e., local and global sensitivity analysis [5,6]. Local sensitivity analysis, which is in the topic of this paper, is performed by investigating the influence of slight deviation from nominal values of different design parameters on the output parameters, i.e., the maximum amplitudes of certain degrees of freedom.
Many scholars have contributed to the applications of sensitivity analysis to specific domains. Su [7] et al. proposed a method to assess the extent of rotor failure and investigated the frequency reliability sensitivity with respect to several stochastic parameters. Melchers and Ahammed [8] presented a structural sensitivity analysis approach based on moment and compared the results with those computed by a conventional Monte Carlo method. Lataillade [9] et al. combined sensitivity analysis with Monte Carlo methods by introducing a supplementary procedure. Petrov [10,11] et al. investigated the sensitivity of the steady response of bladed disk assemblies with respect to parameters of nonlinear friction contact interfaces. Leister [12] et al. proposed a new sensitivity analysis method for a rotor supported by two bearings. Yang and Sievert [13] studied the sensitivity of large turbine shaft system with respect to shaft unbalance. Pan [14] et al. applied sensitivity analysis to investigate the influence of the mass and the offset of disk on rotors. And results have revealed that the offset of disk influence the critical speed in a significant way. Papini [15] et al. carried out a sensitivity analysis on solid rotor induction machines. Ying [16] et al. investigated the load change sensitivities of a multi-support system utilizing virtual displacement principle. Based on the model they proposed, the influence of load variation on dynamic characteristics and sensitivities of a multi-support structure is researched. Anusonti-Inthra [17] et al. carried out sensitivity analysis and uncertainty qualification for a coaxial rotor system, where the sensitivities of hub loads with respect to ambient drift speed, rotor angular speed and center of gravity location on the blade are investigated.
Furthermore, rubbing phenomenon is expected when the clearance between rotor and stator becomes so small that friction and collision happens at interfaces. It makes the response of rotor strongly nonlinear and is sensitive to parameters of friction interfaces. Investigations of dynamic analysis for rotors with rub-impact have been performed by many researchers. Choy and Padovan [18] found that the rub-impact becomes more severe with increase of imbalance of rotor. Dai [19] et al. presented that the increase of excitation force would result in the transformation of contact status, i.e., from partial rubbing to full rubbing. Feng and Zhang [20] investigated the vibration characteristics of a rotor with rubimpact resulting from initial perturbations. Patel [21] et al. utilized full spectrum analysis to investigate the vibration characteristics of a cracked rotor rubbing with stators. The research done by Popprath [22] et al. showed that a reduction in structural damping will result in a periodic response of the rotor subjected to rub-impact to become non-periodic. Liu [23] et al. simulated the rubbing of the rotor experimentally and a serious rub chaos was observed. Ma [24] et al. found that different rubbing forms and transient forces may exist under different operation conditions. Shen [25] et al. researched the influence of different parameters on the rubbing condition and utilized the proposed algorithm to detect rubbing faults. Yao [26] et al. also proposed a method based on harmonic components to detect rubbing fault. Behzad [27] et al. analyzed the characteristics of backward rub both experimentally and numerically. The work done by Liu [28] et al. showed that high-order harmonic terms may appear in the response of rotors subjected to rub-impact.
However, few studies have been done investigating the sensitivity of the response of rotors subjected to rubbing forces. Zhang [29] et al. investigated the sensitivity of a Jeffcott rotor with rub-impact to different parameters, but his research was mainly concentrated on transient sensitivity in the time domain, which is not so practical since the periodic stable vibration is the conventional focus of engineering applications. In this paper, the sensitivities of steady-state response of rotors subjected to nonlinear rubbing forces with respect to parameters of the rubbing interfaces are investigated. Based on the multiharmonic response obtained in frequency domain, analytical expressions for sensitivity coefficients are derived by differentiating the maximum displacements with respect to the investigated parameters, such as normal stiffness and initial gap, etc. Considering the effect of buckling and bifurcation in nonlinear systems, the pseudo arc-length continuation method is applied to trace the equilibrium path of nonlinear equations.
In this section, the research background and status about sensitivity analysis and response of rotors with rub-impact are elaborated. In Sections 2.1 and 2.2, basic theories for the adopted methods are introduced, including the harmonic balance method (HBM) and the pseudo are-length continuation method. The analytical expressions for sensitivity coefficients are derived in Section 2.3 for the purpose of accuracy and performance. Then in Section 3.1, the accuracy of the adopted algorithm is illustrated utilizing two classical models. Sections 3.1 and 3.2 shows the applications of the proposed methods to two different rotor models, a lumped Jeffcott rotor and a finite element model. The summary of results and discussions are finally presented in Section 4.

Harmonic Balnace Method
Since traditional time integral method suffers from expensive cost when solving nonlinear equations, the HBM is applied in this paper, which mainly concentrates on the steady-state regimes of vibration response. In time domain, the equation of motion for a rotor rubbing with rigid stators can be expressed by: where M, C o and K o are, respectively, the linear mass, damping and stiffness matrices of the system; ω is the rotation speed of the system; x is the vector of displacements, and, . x, ..
x are the corresponding derivatives of displacements with respect to time t; ωC c is the Coriolis matrix; ω 2 K s is the matrix describing stiffness softening effect; f ub is the vector of unbalanced force; f nl is the nonlinear force existing in the rubbing interfaces between rotor and stator, which is dependent on the displacement vector x and λ, the parameter of contact interface used for sensitivity analysis.
According to the HBM, the time-variant displacement vector x in Equation (1) is transformed into frequency domain represented as a restricted number of Fourier series, i.e.,: where Q 0 is the vector of constant component in Fourier expansion; Q cj and Q sj are vectors of harmonic coefficients corresponding to cosine and sine terms respectively; m j is the harmonic index; τ is the non-dimensional time; n is the total number of harmonic terms retained, it should be selected carefully so that moderate number of harmonics are kept in harmonic expansion to reflect realistic vibration. The terms of nonlinear force and unbalanced force in Equation (1) are also expanded in frequency domain by sequentially multiplied by cosine terms cos m j τ and sine terms sinm j τ for all harmonic indexes and integrated over the whole vibration period T, where T = 2π/ω. Since all time-variant terms have been expanded in the frequency domain, the multi-harmonic representation of Equation (1) can be obtained by equating harmonic coefficients, i.e.,: is the vector of harmonic components containing all coefficients for displacements described in Equation (2); F n and F u are vectors of multi-harmonic coefficients for rubbing forces and unbalanced forces; Z(ω) is the blockdiagonal dynamic stiffness matrix and it has the following form: In Equation (5), K and C are simplified representation forms of generalized stiffness and damping matrices of system DOFs, with K = K o − ω 2 K s and C = C o + ωC s respectively. The left-hand side of Equation (3) is defined as the residual terms which are used to describe if the sought vector converges to the actual solution. Thus, by solving Equation (3) we can obtain the sought solution of Q and time-domain displacements of system DOFs can be easily obtained using Equation (2). It should be noted that since the interface force f n could only be calculated explicitly in time domain and then the nonlinear term F n in Equation (3) is obtained, the iterative process need to be shifted between time domain and frequency domain for several times until the residual vector R(Q) converge to a prescribed accuracy, i.e., R(Q) < ε.

Pseudo Arc-Length Continuation Method
For tracing the solution curve of nonlinear equations like Equation (3), the Newton-Raphson method is usually applied. However, the Newton iterative method usually fails when buckling or bifurcation happens, so the pseudo arc-length continuation method [30] is used to solve Equation (3). The detailed procedure is described in the following paragraphs.
In the pseudo arc-length continuation method, an extra parameter s is introduced to control the solution searching process. Equation (3) can be rewritten as the following expression: where ψ is represented to be related to ω because ω also varies during the iterative process. At the first iterative step, two initial solutions away from strongly nonlinear region is obtained using Newton method based on Equation (3), i.e.,: where m is the number of previous iterative step. Then the first arc length s 1 for outer iteration is calculated in the following form: where ω 1 , Q 1 and ω 2 , Q 2 are the results obtained using Newton method; ζ is the weight coefficient. Then the inner iteration could be executed but an additional vector is required firstly. It is given as: where , ∂ψ ∂ω and the superscriptˆdenotes elimination of the corresponding column. Obviously, vector v(Q, ω) is the tangential vector of the tracing curve with respect to Q and ω. The unit tangential vector of solution curve is derived as τ(Q, ω) = v v where v represents the Euclidean norm of vector v. For the mth outer iteration, the constraint equation of s m can be expressed as follows: where ∆Q(s) = Q(s) − Q m and ∆ω(s) = ω(s) − ω m ; ∆s is related to the arc length for next step of outer iteration, which can be expressed as ∆s = s m+1 − s m . Then the extra constraint equation Equation (10) can be reorganized into Equation (6) and new equations take the form: The Jacobian matrix of ψ * can be derived based on Equation (11), i.e.,: After all those derivations, the process of inner iteration takes the form: And the initial solution searching point for next outer iteration is obtained by: It should be noted that the arc length s needs to be updated after several steps of outer iteration to obtain better performance in tracing the solution curve near the resonance region.

Sensitivity Analysis
In this paper, a local sensitivity analysis is performed by evaluating the local gradients of maximum displacements with respect to parameters of rubbing interfaces. Typically, what we need is the sensitivity of the maximum displacement in time domain, x max , with respect to different parameters, i.e., the frictional coefficient, the normal stiffness k y , etc. The local sensitivities, defined as slight changes of the response caused by a shift of individual variable, are computed as follows: where s 1 = ∂x max ∂p and s 2 = ∂ 2 x max ∂p 2 are the first-order and second-order sensitivities of timedomain displacements with respect to rubbing parameter p. They need to be derived based on the results obtained in frequency domain. By differentiating Equation (3) with respect to p, the analytical expression of first-order derivatives of harmonic components with respect to the investigated parameter is obtained as: Thus the first-order sensitivity in frequency domain can be expressed as: Equation (17) is a generalized form of first-order sensitivity. If p is selected as a parameter in the friction model, considering that both Z and F u are linear parts independent of nonlinear forces, then the corresponding terms ∂Z ∂p and ∂F u ∂p are equal to zero. Thus Equation (17) can be rewritten as: where J = Z + ∂F n ∂Q is the Jacobian matrix. Since the nonlinear force is usually explicitly derived in most cases, the term ∂F n ∂p can be calculated easily with a Fourier expansion described in Section 2.1. The derivatives of harmonic coefficients of rubbing force with respect to harmonic coefficients of system DOFs can be formulated by: where the superscripts i, j = 1,· · · N refer to the i-th harmonic component of F n and the jth component of Q, respectively. To derive the detailed analytical expressions in Equation (19), the analysis in time domain is conducted at first. By differentiating nonlinear forces with respect to the displacements in local coordinate frame, one can obtain that: (20) so the expression for derivatives of rubbing force with respect to global displacements can be derived as: where T is the transformation matrix between local frame and global frame. For fixed local frame, it's the identity matrix and for rotating local frame it has the following form: This transformation matrix, T, makes the formulation of normal rubbing force and the tangential rubbing force easily to be completed when transferred to local coordinate. Additionally, define a transformation vector Γ and Ψ as follows: , cos(τ), sin(τ), · · · , cos(nτ), sin(nτ) (23) ψ = [1, cos(τ), sin(τ), · · · , cos(nτ), sin(nτ)] (24) Obviously, N is equal to 2n + 1. Then Equation (19) can be expressed as: where ⊗ represents the Kronecker product; the subscript n implies that only those terms corresponding to nonlinear forces and nonlinear DOFs are nonzero elements; the terms τ j and τ j+1 are the time instants when contact status begins and ceases respectively; n τ is the total number of contact status during a time period. Again by differentiating Equation (16) with respect to the concerned parameter we can obtain the expression of second order derivative as: Since no quadratic terms exist in the expression of rubbing force, the term ∂ 2 F n ∂p 2 is equal to zero. The term ∂ 2 F n ∂p∂Q can be obtained by differentiating Equation (25) with respect to the parameter p. Considering only the matrix D* in Equation (25) where H = 1 π τ j+1 τ j Γ T Ψdτ. The components in Equation (27) can be computed by: The lowercase letter q represents the displacements in time domain. Then the maximum displacement in time domain can be determined by the following expression: Substituting Equation (29) into Equation (16) and Equation (20) yields: The extra term ∂τ max ∂p is introduced in consideration that the non-dimensional time instant for maximum displacement is dependent on the parameter p. Thus the first-order sensitivity s 1 and the second-order sensitivity s 2 are obtained. For different parameters which possess different units, the results may be misunderstanding. So Equations (30) and (31) can be refined as: In this cases, the sum of s * 1 and s * 2 can be expressed as the variation of maximum amplitudes under slight change of specific parameters. When p + 0 = (1 + 1%)p 0 , terms in Equation (32) describe the variation of maximum amplitudes under 1 percent change of the investigated parameters, which is the main focus in the following discussion.

Models for Verification
Firstly, two classical one-dimensional nonlinear models, including the Duffing oscillator and the gap model, are utilized to verify the accuracy and efficiency of the methods adopted.

The Duffing Oscillator
A Duffing oscillator is utilized in this section, as shown in Figure 1, with stiffness coefficient k = 40 and damping coefficient c = 0.1 √ 40 respectively. The nonlinear force is expressed in the form of cubic function, and the friction coefficient k τ is the investigated parameter in sensitivity analysis.

The Duffing Oscillator
A Duffing oscillator is utilized in this section, as shown in Figure 1, with coefficient k = 40 and damping coefficient c = 0.1 40 respectively. The nonlinear expressed in the form of cubic function, and the friction coefficient k τ is the inve parameter in sensitivity analysis.  ..
where the excitation force f e has the form f e = −100 sin(ωt). Although Equation (33)  The sensitivity coefficients are specifically defined as the variation of maximum displacements when the parameters of interfaces change by one percent, as shown in Figure 5. Results reveal that with increase of friction coefficient, the sensitivity of the response will decrease. An abrupt change of first-order sensitivity happens at around ω = 6.5 rad/s. In this case, the first-order sensitivity is extremely sensitive to the variation of friction coefficient thus the values of second-order sensitivity rise, as shown in Figure 5b.   The sensitivity coefficients are specifically defined as the variation of maximum displacements when the parameters of interfaces change by one percent, as shown in Figure  5. Results reveal that with increase of friction coefficient, the sensitivity of the response    The sensitivity coefficients are specifically defined as the variation of maximum displacements when the parameters of interfaces change by one percent, as shown in Figure  5. Results reveal that with increase of friction coefficient, the sensitivity of the response will decrease. An abrupt change of first-order sensitivity happens at around ω = 6.5 rad/s. In this case, the first-order sensitivity is extremely sensitive to the variation of friction coefficient thus the values of second-order sensitivity rise, as shown in Figure 5b.

The Gap Model
As plotted in Figure 6, a similar oscillator subjected to pure-collision nonlinear force is applied in this section. Since no friction force is introduced in this system, it is referred to as the gap model

The Gap Model
As plotted in Figure 6, a similar oscillator subjected to pure-collision nonlinear force is applied in this section. Since no friction force is introduced in this system, it is referred to as the gap model.

The Gap Model
As plotted in Figure 6, a similar oscillator subjected to pure-collision nonlinear force is applied in this section. Since no friction force is introduced in this system, it is referred to as the gap model The equation of motion for this model is exactly the same with Equation (33) while the nonlinear force has the following form: where g is the prescribed clearance and ky is the normal stiffness. Figure 7 shows the response of the system with different normal stiffness and initial gaps. It can be seen that an offset of the critical speed exists in both cases. With the increase of normal stiffness or decrease of initial gap, the critical speed tends to increase. And multi-solution region emerges as a result. Results have been shown to be consistent with the work done by Petrov et al. [31]. Time-domain displacement results at ω = 7.029 rad/s obtained by the HBM and TMM are compared in Figure 8 and the accuracy of the HBM can be proved. The equation of motion for this model is exactly the same with Equation (33) while the nonlinear force has the following form: where g is the prescribed clearance and k y is the normal stiffness. Figure 7 shows the response of the system with different normal stiffness and initial gaps. It can be seen that an offset of the critical speed exists in both cases. With the increase of normal stiffness or decrease of initial gap, the critical speed tends to increase. And multi-solution region emerges as a result. Results have been shown to be consistent with the work done by Petrov et al. [31]. Time-domain displacement results at ω = 7.029 rad/s obtained by the HBM and TMM are compared in Figure 8 and the accuracy of the HBM can be proved. Moreover, the sensitivities of the response with respect to normal stiffness and initial gaps are obtained, as illustrated in Figures 9 and 10. It can be seen that multi-solution regions also appear in the sensitivity curves. With an increase of normal stiffness, the multi-solution region tends to grow larger, as shown in Figure 9a. Moreover, with an increase of the initial gap, the sensitivities of the system increase and the multi-solution region shrinks, as illustrated in Figure 10a. Moreover, the sensitivities of the response with respect to normal stiffness and initial gaps are obtained, as illustrated in Figures 9 and 10. It can be seen that multi-solution regions also appear in the sensitivity curves. With an increase of normal stiffness, the multi-solution region tends to grow larger, as shown in Figure 9a. Moreover, with an increase of the initial gap, the sensitivities of the system increase and the multi-solution region shrinks, as illustrated in Figure 10a.

Sensitivity Analysisof a Jeffcott Rotor
Since the accuracy and capability of the methods adopted are demonstrated in Section 3.1, in this section the sensitivities of a rotor subjected to rubbing forces are investigated by adopting a lumped model, which is also named a Jeffcott rotor. A schematic diagram of the Jeffcott rotor is illustrated in Figure 11. It's assumed that the periodic exci-

Sensitivity Analysisof a Jeffcott Rotor
Since the accuracy and capability of the methods adopted are demonstrated in Section 3.1, in this section the sensitivities of a rotor subjected to rubbing forces are investigated by adopting a lumped model, which is also named a Jeffcott rotor. A schematic diagram of the Jeffcott rotor is illustrated in Figure 11. It's assumed that the periodic excitation force is induced by the unbalanced mass mub and the eccentric distance is e. The

Sensitivity Analysisof a Jeffcott Rotor
Since the accuracy and capability of the methods adopted are demonstrated in Section 3.1, in this section the sensitivities of a rotor subjected to rubbing forces are investigated by adopting a lumped model, which is also named a Jeffcott rotor. A schematic diagram of the Jeffcott rotor is illustrated in Figure 11. It's assumed that the periodic excitation force is induced by the unbalanced mass m ub and the eccentric distance is e. The structural damping coefficient, c, is assumed to be proportional to the stiffness k and the mass m. And the relevant parameters are listed in Table 1. 12 can be derived in time domain easily. As for fn, since the values of force along normal direction will always be positive, it can be expressed by the following formula: Then the tangential friction force can also be determined in the following form: where μ is the friction coefficient at the interfaces; fτ is always negative considering that the direction of friction force is opposite to the direction of relative velocity.
Based on this friction model, firstly, the response of the Jeffcott rotor with different initial parameters of rubbing interfaces is illustrated. These parameters include the normal stiffness of stator, the friction coefficient and the initial gap. The first seven harmonic components are kept to approximate the periodic response.  Firstly, the sensitivity of the response to different parameters is investigated using the once-at-a-time (OAT) method [32], which implies changing an individual variable slightly to investigate its influence on the response. According to Figure 13a, the variation of initial gap g has a significant influence on the response curve of the Jeffcott rotor. It can not only influence the maximum displacements, but also determines the range of the  The influences of other parameters of rubbing interfaces, such as normal stiffness k y and initial gap g, are investigated. The normal values of normal stiffness k y , initial gap g and friction coefficient µ are 8 × 10 4 N/m, 1 mm, and 0.7, respectively.
The friction model applied is a classical Coulomb friction model, as shown in Figure 12. It should be noted that the model is established based on local rotating frame. The analytical expressions for tangential rubbing force f τ and normal rubbing force f n in Figure 12 can be derived in time domain easily. As for f n , since the values of force along normal direction will always be positive, it can be expressed by the following formula: Then the tangential friction force can also be determined in the following form: where µ is the friction coefficient at the interfaces; f τ is always negative considering that the direction of friction force is opposite to the direction of relative velocity.
Based on this friction model, firstly, the response of the Jeffcott rotor with different initial parameters of rubbing interfaces is illustrated. These parameters include the normal stiffness of stator, the friction coefficient and the initial gap. The first seven harmonic components are kept to approximate the periodic response.  Firstly, the sensitivity of the response to different parameters is investigated using the once-at-a-time (OAT) method [32], which implies changing an individual variable slightly to investigate its influence on the response. According to Figure 13a, the variation of initial gap g has a significant influence on the response curve of the Jeffcott rotor. It can not only influence the maximum displacements, but also determines the range of the Firstly, the sensitivity of the response to different parameters is investigated using the once-at-a-time (OAT) method [32], which implies changing an individual variable slightly to investigate its influence on the response. According to Figure 13a, the variation of initial gap g has a significant influence on the response curve of the Jeffcott rotor. It can not only influence the maximum displacements, but also determines the range of the speed region subjected to rubbing forces. Furthermore, with the increase of initial gap g, the influence of initial gap gradually becomes less significant since the variation of response curve becomes smaller. Although normal stiffness k y and friction coefficient µ don't influence the response as the way initial gap does, they can change the characteristics of multi-solution region. As shown in Figure 13b,c, with the decrease of k y or increase of µ, the multi-solution region narrows. The influence of ky and µ is much smaller than that of the initial gap g, according to the variation of the response curve. speed region subjected to rubbing forces. Furthermore, with the increase of initial gap g, the influence of initial gap gradually becomes less significant since the variation of response curve becomes smaller. Although normal stiffness ky and friction coefficient μ don't influence the response as the way initial gap does, they can change the characteristics of multi-solution region. As shown in Figure 13b,c, with the decrease of ky or increase of μ, the multi-solution region narrows. The influence of ky and μ is much smaller than that of the initial gap g, according to the variation of the response curve.  Figure 14 shows the sensitivity coefficients of the response, which is defined as the variation of maximum displacements when the parameters of rub interfaces change by one percent. In Figure 14a, first-order sensitivities with respect to all parameters of rub interfaces are equal to zero when the rotation speed range is outside of 138-195 rad/s. A similar phenomenon can be found in Figure 13a, where the response characteristics of the  Figure 14 shows the sensitivity coefficients of the response, which is defined as the variation of maximum displacements when the parameters of rub interfaces change by one percent. In Figure 14a, first-order sensitivities with respect to all parameters of rub interfaces are equal to zero when the rotation speed range is outside of 138-195 rad/s. A similar phenomenon can be found in Figure 13a, where the response characteristics of the rotor can be regarded as a linear system outside the speed range 138-195 rad/s. The abrupt jump of the sensitivity curve with respect to initial gap g represents the change of contact/separation status. Among them the initial gap g plays a dominant role since it determines whether rub-impact exists and the results is consistent with that obtained by the OAT method. Similar phenomena also exist in the second-order sensitivity curves with the initial gap being the most significant impact factor, as illustrated in Figure 14b. It should be noted that first-order sensitivities are much larger than second-order sensitivities, so our discussions are mainly based on the former case. should be noted that first-order sensitivities are much larger than second-order sensitivities, so our discussions are mainly based on the former case. In Figures 15-17 results are shown for cases when normal values of parameters are varied. As illustrated in Figure 15a,b, the first-order sensitivities decrease while friction coefficients are rising. Moreover, an obvious offset of the sensitivities is also observed, which reflects the emergence of a multi-solution region. Figure 16 shows the results of sensitivities with different normal stiffness coefficients. The variation of sensitivities is relatively much smaller in this case, which implies that the response is relatively less sensitive to ky. As is shown in Figure 17, the response becomes more sensitive with the increase of the initial gap and after it reaches a threshold the rotor becomes a linear system, as seen in the separated condition illustrated in Figure 13a.  In Figures 15-17 results are shown for cases when normal values of parameters are varied. As illustrated in Figure 15a,b, the first-order sensitivities decrease while friction coefficients are rising. Moreover, an obvious offset of the sensitivities is also observed, which reflects the emergence of a multi-solution region. Figure 16 shows the results of sensitivities with different normal stiffness coefficients. The variation of sensitivities is relatively much smaller in this case, which implies that the response is relatively less sensitive to k y . As is shown in Figure 17, the response becomes more sensitive with the increase of the initial gap and after it reaches a threshold the rotor becomes a linear system, as seen in the separated condition illustrated in Figure 13a.

Sensitivity Analysis Based on Finite Element Model
Another rotor with rub-impact is discretized using the finite element method (FEM). As shown in Figure 18, the model is made up of 20 equal-length beam elements with the Young's modulus, density and Poisson's ratio equaling to 2 × 10 5 MPa, 7.8 × 10 3 kg/m 3 and 0.3, respectively. The model is supported at both ends of the shaft. The motion of each element is assumed to be within the plane vertical to axial direction. The multiplication of unbalanced mass and its eccentric distance is equal to 7.0711 × 10 −4 kg·m in this paper and the excitation force is applied based on it. Moreover, the rubbing forces are applied to node 7, with the Coulomb friction model mentioned in Section 3.2. The nominal values of parameters of rubbing interfaces are also prescribed, with stator stiffness k y , the friction coefficient µ and the initial gap g equal to 1.0 × 10 7 N/m, 0.5, 0.1 mm, respectively. For the purpose of performance and accuracy, the time sampling number is 1024 and the first seven harmonic components are kept to approximate the response. The friction model applied on elements subjected to rubbing forces is exactly the same as that of Section 3.2.
In Figures 15-17 results are shown for cases when normal values of parameters are varied. As illustrated in Figure 15a,b, the first-order sensitivities decrease while friction coefficients are rising. Moreover, an obvious offset of the sensitivities is also observed, which reflects the emergence of a multi-solution region. Figure 16 shows the results of sensitivities with different normal stiffness coefficients. The variation of sensitivities is relatively much smaller in this case, which implies that the response is relatively less sensitive to ky. As is shown in Figure 17, the response becomes more sensitive with the increase of the initial gap and after it reaches a threshold the rotor becomes a linear system, as seen in the separated condition illustrated in Figure 13a.

Sensitivity Analysis Based on Finite Element Model
Another rotor with rub-impact is discretized using the finite element method (FEM). As shown in Figure 18, the model is made up of 20 equal-length beam elements with the

Sensitivity Analysis Based on Finite Element Model
Another rotor with rub-impact is discretized using the finite element method (FEM). As shown in Figure 18, the model is made up of 20 equal-length beam elements with the 5 3 3 Similar conclusions can be drawn based on this model, i.e., the response is much more sensitive to the initial gap g and with its increase the sensitivity decreases, as in Figure 19a. Outside of the speed range 260.9 rad/s to 298.9 rad/s there is no rub-impact forces so the response is exactly the same with that of separation status. When rubbing forces exist, a multi-solution region emerges because of the nonlinearity of the system. Moreover, the increase of normal stiffness and decrease of friction coefficient results in the expansion of the multi-solution region, as seen in Figure 19b,c.  Similar conclusions can be drawn based on this model, i.e., the response is much more sensitive to the initial gap g and with its increase the sensitivity decreases, as in Figure 19a. Outside of the speed range 260.9 rad/s to 298.9 rad/s there is no rub-impact forces so the response is exactly the same with that of separation status. When rubbing forces exist, a multi-solution region emerges because of the nonlinearity of the system. Moreover, the increase of normal stiffness and decrease of friction coefficient results in the expansion of the multi-solution region, as seen in Figure 19b Similar conclusions can be drawn based on this model, i.e., the response is much more sensitive to the initial gap g and with its increase the sensitivity decreases, as in Figure 19a. Outside of the speed range 260.9 rad/s to 298.9 rad/s there is no rub-impact forces so the response is exactly the same with that of separation status. When rubbing forces exist, a multi-solution region emerges because of the nonlinearity of the system. Moreover, the increase of normal stiffness and decrease of friction coefficient results in the expansion of the multi-solution region, as seen in Figure 19b,c.  Figure 20 shows the sensitivities of the frequency-domain response to different parameters of rub-impact interfaces. Similar phenomena have been observed as obtained by the Jeffcott rotor. The responses are more sensitive to the initial gap g compared to the other investigated parameters. The sensitivity curve near the vicinity of the rotating speeds where the separation-contact status transformation happens is discontinuous, since the sensitivities are equal to zero when no rub-impact exists. An obvious offset of sensitivities is seen under this condition, as a result of the existence of a multi-solution region in the response curve.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 18 of 21 normal stiffness at g = 0.1 mm and μ = 0.5; (c) Response curves with different normal friction coefficients at ky = 1 × 10 7 N/m and g = 0.1 mm. Figure 20 shows the sensitivities of the frequency-domain response to different parameters of rub-impact interfaces. Similar phenomena have been observed as obtained by the Jeffcott rotor. The responses are more sensitive to the initial gap g compared to the other investigated parameters. The sensitivity curve near the vicinity of the rotating speeds where the separation-contact status transformation happens is discontinuous, since the sensitivities are equal to zero when no rub-impact exists. An obvious offset of sensitivities is seen under this condition, as a result of the existence of a multi-solution region in the response curve. By changing initial values of initial gap g, the sensitivities to the initial gap with different prescribed g values are investigated, as shown in Figure 21. Still, s discontinuity exists where separation-contact transformation happens. With the increase of initial gap g, the speed region subjected to rub-impact narrows and the response becomes more sensitive. Moreover, with an increase of the initial gap, the multi-solution region of sensitivities narrows as well.  By changing initial values of initial gap g, the sensitivities to the initial gap with different prescribed g values are investigated, as shown in Figure 21. Still, s discontinuity exists where separation-contact transformation happens. With the increase of initial gap g, the speed region subjected to rub-impact narrows and the response becomes more sensitive. Moreover, with an increase of the initial gap, the multi-solution region of sensitivities narrows as well.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 18 of 21 normal stiffness at g = 0.1 mm and μ = 0.5; (c) Response curves with different normal friction coefficients at ky = 1 × 10 7 N/m and g = 0.1 mm. Figure 20 shows the sensitivities of the frequency-domain response to different parameters of rub-impact interfaces. Similar phenomena have been observed as obtained by the Jeffcott rotor. The responses are more sensitive to the initial gap g compared to the other investigated parameters. The sensitivity curve near the vicinity of the rotating speeds where the separation-contact status transformation happens is discontinuous, since the sensitivities are equal to zero when no rub-impact exists. An obvious offset of sensitivities is seen under this condition, as a result of the existence of a multi-solution region in the response curve. By changing initial values of initial gap g, the sensitivities to the initial gap with different prescribed g values are investigated, as shown in Figure 21. Still, s discontinuity exists where separation-contact transformation happens. With the increase of initial gap g, the speed region subjected to rub-impact narrows and the response becomes more sensitive. Moreover, with an increase of the initial gap, the multi-solution region of sensitivities narrows as well.  Similarly, as seen in Figure 22, results show that with decrease of friction coefficient the offset of sensitivities, which is more obvious in the 300-330 rad/s speed range when µ is equal to 0.4, becomes more severe with the decrease of friction coefficients. In other words, the multi-solution region grows while µ is decreasing. This is consistent with the results obtained by the OAT method, as illustrated in Figure 19c. Similarly, as seen in Figure 22, results show that with decrease of friction coefficient the offset of sensitivities, which is more obvious in the 300-330 rad/s speed range when μ is equal to 0.4, becomes more severe with the decrease of friction coefficients. In other words, the multi-solution region grows while μ is decreasing. This is consistent with the results obtained by the OAT method, as illustrated in Figure 19c. According to Figure 23, the variation of normal stiffness ky may also change the stable sensitivity region. Moreover, the variation of ky has relatively little influence on the sensitivities, but the multi-solution region develops indeed with the increase of ky.  According to Figure 23, the variation of normal stiffness k y may also change the stable sensitivity region. Moreover, the variation of k y has relatively little influence on the sensitivities, but the multi-solution region develops indeed with the increase of k y . Similarly, as seen in Figure 22, results show that with decrease of friction coefficient the offset of sensitivities, which is more obvious in the 300-330 rad/s speed range when μ is equal to 0.4, becomes more severe with the decrease of friction coefficients. In other words, the multi-solution region grows while μ is decreasing. This is consistent with the results obtained by the OAT method, as illustrated in Figure 19c. According to Figure 23, the variation of normal stiffness ky may also change the stable sensitivity region. Moreover, the variation of ky has relatively little influence on the sensitivities, but the multi-solution region develops indeed with the increase of ky.

Conclusions
This paper investigated the sensitivity characteristics of rotors subjected to rub-impact problems to parameters of rubbing interfaces. The harmonic balance method and pseudo arc-length continuation method are utilized to obtain the steady-state response of the system. Furthermore, the analytical expressions for the sensitivity coefficients are derived for the purpose of accuracy and performance. Two classical models, namely the Duffing oscillator and the gap model, are utilized to verify the accuracies and capabilities of the adopted methods. Then numerical investigations of the sensitivities of steady-state response to parameters of rubbing interfaces are performed and conclusions are drawn as follows: (1) The response of rotors subjected to rub-impact is sensitive to the parameters of rub interfaces only when the amplitudes are large enough to reach or exceed the prescribed initial clearance. If that is not the case, the whole system is equal to a linear system under a periodical excitation force. (2) The response of the rotor is much more sensitive to the initial gap with the two other factors being relatively insignificant. The sensitivities to the initial gap g are not contiguous at rotating speeds where a separation/contact status transformation happens since it determines the existence of rub-impact. (3) With the increase of initial gap, the rotors become more sensitive. In contrast, the friction coefficient has negative impact on the sensitivities and the multi-solution region grows larger with its decrease. Moreover, the variation of normal stiffness has relatively little influence on the sensitivities but it does influence the characteristics of the multi-solution region, i.e., the region expands as it increases.