Phase difference and stability of a shaft mounted a dry friction damper: Effects of viscous internal damping and gyroscopic moment

Abstract

This study investigates the stability and phase difference of a shaft mounted a dry friction damper with effects of viscous internal damping and gyroscopic moment. The equations of the system with the vibration reduction effect of the dry friction damper on the shaft are derived in the form of the rectangular coordinate and polar coordinate in the vicinity of critical speed. The phase difference characteristics in the rub-impact process and its physical mechanism are analyzed by mathematical derivation. The characteristic equation is studied to investigate the stability of the periodic solution. Effects of different parameters of the system, especially viscous internal damping of the composite shaft and gyroscopic moment on the phase difference and stability regions are presented in detail by analytical and numerical simulation based on a helicopter tailrotor driveline. The experimental investigation is conducted in a test rig to validate theoretical formulas and simulation analysis. The analysis results show that rub impact delays the change of phase difference, viscous internal damping improves the stability of synchronous full annual rub solution, and gyroscopic moment affects the increase of the phase difference.

Introduction

The driveshaft in the supercritical regime has been commonly used in many helicopter tailrotor drivelines. The lateral amplitude of the transmission shaft increases significantly when crossing the first critical speed because of the inevitable eccentricity. ¹ In an attempt to alleviate shaft vibration, dampers or shock absorbers are set in the drivelines. ² Dry friction damper, as shown in Figure 1, is used in the certain helicopter tailrotor drivelines as an effective damper with the advantages of simplicity, reliability, and so on.^3,4 The most important component of the damper is a damping ring, which is confined to the base by springs and has radial clearance between itself and the sleeve fixed on the shaft. The whirl of the shaft is restricted by the damping ring, which leads to the relationship between the shaft and the dry friction damper is similar to the rub-impact of the rotor-stator system.

Figure 1.

Schematic diagnosis of shaft and dry friction damper.

There are two patterns of rub-impact motions, one is the synchronous full annular rub, that is, the sleeve slides and rotates along the inner surface of the ring at spin speed. The second is the partial rub, the sleeve contacts and separates with the inner surface of the ring intermittently. ⁵ Comparing with the first pattern, the partial rub is more likely to cause wear and tear due to repeated impact. ⁶ These two patterns can be distinguished through the stability evaluation of the periodic solution of the system. There is extensive research work on rub-impact phenomenon and stability,^7–9 mainly including the stability of parameter regions and bifurcation boundaries analysis of rub-impact solution,^10–12 bifurcation and chaos characteristics of the rotor-stator system,^13,14 and analysis about the synchronous full annular rub and partial rub.^15,16 Xu et al. ¹⁷ studied the full annular rub motion of a Jeffcott rotor and its dynamic stability by linear perturbation method and by the Routh-Hurwitz criterion, finding that the full annular rub motions exist in a wide spinning region of the rotor. Shang et al. ¹⁸ studied the global response characteristics of a general rotor-stator system. The dry friction effect is considered in the process of rubbing, which is regarded as the main factor for the self-excited dry friction backward whirl. Zhang et al. ¹⁹ studied the nonlinear dynamic characteristics of a rotor-bearing system with rub-impact, investigating effects of parameters such as eccentricity, stiffness of the shaft, and radial stiffness of stator on the system through the numerical calculation. Wang et al. ²⁰ proposed a novel method of rubbing fault diagnosis based on variational mode decomposition (VMD) for vibration signal analysis of rotary machinery. Hu et al. ²¹ investigated the rub, crack and rub-crack coupled fault through experiments on a rotor test rig and analyzed the instantaneous frequency signatures of the faults. Wang et al. ²² developed a dual-rotor system with inter-shaft bearing able to describe the mechanical vibration caused by unbalance and fixed-point rub-impact, to understand the mechanism of fixed-point rub-impact in aero-engine.

Previous investigations are mainly based on the Jeffcott rotor, which is a typical lumped mass model But the internal damping effect due to the material of the shaft cannot be considered in this system. The application trend of composite transmission shaft in the helicopter tail rotor drive system is more extensive. This effect in the composite materials is more significant than metallic materials, so it cannot be ignored to ensure the accuracy of the analysis results.^23,24 Even though the internal damping effect of shaft materials has been studied by some researchers, such as Montagnier et al. ²⁵ investigated the stability region about the free movement of the rotating shaft. Desmidt ²⁶ found that free whirl instability occurs at shaft speeds above the first critical speed and is related to the internal and external-damping ratio. Ren et al. ²⁷ studied the dynamics of a rotating thin-walled composite shaft with internal damping analytically, showing the effect of design parameters on the instability thresholds of shafts. Monajemi and Mohammadimehr ²⁸ examined the effects of viscous and hysteretic damping on the vibrational behavior and stability of a spinning Timoshenko micro-shaft. Ghasabi et al. ²⁹ investigated the dynamic bifurcation of a viscoelastic micro rotating shaft, showing that the internal and external damping coefficients influence the critical speed, amplitude, and phase of a non-trivial solution, and radius of the limit cycle. Ben Arab et al. ³⁰ introduced an Euler–Bernoulli shaft finite element formulation including the hysteretic internal damping of composite material and transverse shear effects, then used it to evaluate the influence of various parameters, including instability thresholds. Hosseini ³¹ modeled the shaft as an in-extensional spinning beam, which includes the effects of nonlinear curvature, inertia, and the internal damping, using center manifold theory and the method of normal form, analytical expressions are obtained, which describe the behavior of the rotating shaft in the neighborhood of the bifurcations. These researches are mainly about the instability of the free whirl or flutter of supercritical shafts, however, to the best of the authors’ knowledge, no research has involved the stability analysis for a rotating system with internal damping and undergoes unbalanced force and rub-impact.

The rub-impact phenomenon with non-linear behaviors is one of the classic research fields of rotor dynamics.^32–34 Nevertheless, one can recognize that the phase difference between excitation and the response during rub-impact has not attracted enough attention in the current study. In fact, it can be regarded as another key point besides the amplitude for the dynamic characteristics of the system. The experimental results show that there is a great difference between the phase difference with and without rub impact. ³⁵ There is a certain relationship between the change of phase difference, parameter conditions for the appearance of saddle-node bifurcation, and the transition from rubbing response to non-rubbing periodic response. ³⁶ The specific phase difference characteristics in the rub-impact process and its mechanism need for in-depth and clear research, especially, the changes of phase difference after the application of composite transmission shaft.

In order to provide parameter guidelines for the design of dry friction damper used for slender flexible transmission shaft and reveal the mechanism of the phase difference, this paper focuses on the stability analysis and the change of phase difference about the shaft and dry friction damper system with viscous internal damping and gyroscopic moment under coupling of unbalanced force and non-linear rub-impact excitation. After this introduction, rub impact model between a slender shaft and a dry friction damper is established in Section 2. The relationship between the angular velocity of phase difference change and system parameters and criteria for stability and bifurcation are addressed by mathematical derivation in Section 3. Characteristics and mechanism of phase difference, effects of parameters, especially internal viscous damping and gyroscopic moment are studied by simulation in Section 4. The test rig is designed and experiments are conducted to verify simulation results in Section 5. Finally, conclusions are drawn in Section 6. A flowchart where the work flow of the present study is illustrated in Figure 2.

Figure 2.

Flowchart of the present study.

Mathematical formulation and condition of contact occurrence

The dynamic model of the system

Figure 1 shows the schematic diagram of the shaft and dry friction damper. The shaft is rotating about the X-axis, relative to the inertia-fixed coordinate frame XYZ. The X_rY_rZ_r is a rotating frame coordinate with X_r coinciding with the rotation axis. Base on the governing equations of the slender shaft in Montagnier and Hochard, ²⁵ Desmidt, ²⁶ viscous internal damping (rotating damping) of shaft generated by the variation of the dissipative energy function, and external damping (nonrotating damping) provided by the bearing and block are distinguished and introduced in the governing equations. Then we obtain equations for a rotating shaft with disc:

$$\begin{array}{c}\rho A\frac{{\partial }^{2}v(x,t)}{\partial t^2}-\rho I\frac{{\partial }^{4}v(x,t)}{\partial z^2\partial t^2}+c_s\mathrm{EI}\frac{{\partial }^{5}v(x,t)}{\partial t\partial z^4}\hfill \\ -2\rho I\mathrm{\Omega }\frac{{\partial }^{3}w(x,t)}{\partial z^2\partial t}+\mathrm{EI}\frac{{\partial }^{4}v(x,t)}{\partial z^4}\hfill \\ -c_s\mathrm{EI}\mathrm{\Omega }\frac{{\partial }^{4}w(x,t)}{\partial z^4}+\sum _{i=1}^n\gamma (x-L_i)\hfill \\ [m_{h,i}\frac{{\partial }^{2}v(L_{\mathrm{hi}},t)}{\partial t^2}-{\rho }_{h,i}{I}_{h,i}\frac{{\partial }^{4}v(L_{\mathrm{hi}},t)}{\partial z^2\partial t^2}-2{\rho }_{h,i}{I}_{h,i}\mathrm{\Omega }\frac{{\partial }^{3}w(L_{\mathrm{hi}},t)}{\partial z^2\partial t}]\hfill \\ =\rho \mathrm{Ae}(x){\mathrm{\Omega }}^2\cos (\mathrm{\Omega }t)+\sum _{i=1}^n\gamma (x-L_i)m_{e,i}{e}_{h,i}{\mathrm{\Omega }}^2\cos (\mathrm{\Omega }t)\hfill \\ -c_n\frac{\partial v(L_n,t)}{\partial t}\gamma (x-L_n)-k_{n}v(L_n,t)\gamma (x-L_n)\hfill \\ \rho A\frac{{\partial }^{2}w(x,t)}{\partial t^2}-\rho I\frac{{\partial }^{4}w(x,t)}{\partial z^2\partial t^2}+c_s\mathrm{EI}\frac{{\partial }^{5}w(x,t)}{\partial t\partial z^4}\hfill \\ +2\rho I\mathrm{\Omega }\frac{{\partial }^{3}v(x,t)}{\partial z^2\partial t}+\mathrm{EI}\frac{{\partial }^{4}w(x,t)}{\partial z^4}\hfill \\ +c_s\mathrm{EI}\mathrm{\Omega }\frac{{\partial }^{4}v(x,t)}{\partial z^4}+\sum _{i=1}^n\gamma (x-L_i)\hfill \\ [m_{h,i}\frac{{\partial }^{2}w(L_{\mathrm{hi}},t)}{\partial t^2}-{\rho }_{h,i}{I}_{h,i}\frac{{\partial }^{4}w(L_{\mathrm{hi}},t)}{\partial z^2\partial t^2}+2{\rho }_{h,i}{I}_{h,i}\mathrm{\Omega }\frac{{\partial }^{3}v(L_{\mathrm{hi}},t)}{\partial z^2\partial t}]\hfill \\ =\rho \mathrm{Ae}(x){\mathrm{\Omega }}^2\sin (\mathrm{\Omega }t)+\sum _{i=1}^n\gamma (x-L_i)m_{e,i}{e}_{h,i}{\mathrm{\Omega }}^2\sin (\mathrm{\Omega }t)\hfill \\ -c_n\frac{\partial w(L_n,t)}{\partial t}\gamma (x-L_n)-k_{n}w(L_n,t)\gamma (x-L_n)\hfill \end{array}$$

(1)

where the parameters in equation (1) are shown in the Notation Table. The Galerkin method is used to solve the above partial differential equations, the shaft transverse deflections and eccentricity are written in terms of the following modal expansion:

$$\begin{array}{c}\multicolumn{1}{c}{v(x,t)=\sum _{r=1}^{+\infty }{\mathrm{\Phi }}_r(x)Q_{v,r}(t),w(x,t)}\\ =\sum _{r=1}^{+\infty }{\mathrm{\Phi }}_r(x)Q_{w,r}(t),e(x)=\sum _{r=1}^{+\infty }{\mathrm{\Phi }}_r(x)e_r\hfill \end{array}$$

(2)

By substituting equation (2) into equation (1), multiplying the two sides of the equation by ${\mathrm{\Phi }}_r(x)$ , integrating x over X-axis, then using the orthogonality of the modal function, the differential equations are obtained,

$$\begin{array}{c}{M}_r=\rho A{\int }_0^l{\mathrm{\Phi }}_r^2(x)\mathrm{dx}-\rho I{\int }_0^l{\mathrm{\Phi }}_r(x){\mathrm{\Phi }″}_r(x)\mathrm{dx}\hfill \\ +\sum _{i=1}^n[m_{h,i}{\mathrm{\Phi }}_r(L_{\mathrm{hi}})-\rho I_{h,i}{\mathrm{\Phi }″}_r(L_{\mathrm{hi}})]\hfill \\ G_r=2\rho I{\int }_0^l{\mathrm{\Phi }}_r(x){\mathrm{\Phi }″}_r(x)\mathrm{dx}+2\rho I_{\mathrm{sz}}{\mathrm{\Phi }″}_r(L_{\mathrm{s}})K_r=\mathrm{EI}{\int }_0^l{\mathrm{\Phi }}_r(x){\mathrm{\Phi }}_r^{(4)}(x)\mathrm{dx}+k_n{\mathrm{\Phi }}_r(L_s)\hfill \\ m_{e,r}=\rho A{\int }_0^l{\mathrm{\Phi }}_r^2(x)\mathrm{dx}+\sum _{i=1}^n[m_{\mathrm{eh},i}{e}_{h,i}{\mathrm{\Phi }}_r(L_{\mathrm{hi}})]\hfill \end{array}$$

(3)

$$\begin{array}{c}{K}_{s,r}=C_{s,r}=c_s\mathrm{EI}{\int }_0^l{\mathrm{\Phi }}_r(x){\mathrm{\Phi }}_r^{(4)}(x)\mathrm{dx}\hfill \\ C_{n,r}=c_n{\mathrm{\Phi }}_r(L_{\mathrm{s}})\hfill \end{array}$$

(3)

Ref. has proved that taking the first two modes of flexible shaft system will guarantee the accuracy of the numerical simulation around the first critical speed to fourfold. While only the first critical speed is the interest region where damper works in this work, thus taking the first mode from can meet the requirement of accuracy in this work,

$$\begin{array}{c}\left[\begin{array}{cc}{M}_1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& M_1\hfill \end{array}\right]\left\{\begin{array}{c}Z″\hfill \\ \multicolumn{1}{c}{Y″}\end{array}\right\}\hfill \\ +(\left[\begin{array}{cc}{C}_{n,1}+C_{s,1}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& C_{n,1}+C_{s,1}\hfill \end{array}\right]+\mathrm{\Omega }\left[\begin{array}{cc}0\hfill & -G_1\hfill \\ \multicolumn{1}{c}{G_1}& 0\hfill \end{array}\right])\hfill \\ \left\{\begin{array}{c}Z\prime \hfill \\ \multicolumn{1}{c}{Y\prime }\end{array}\right\}\hfill \\ +(\left[\begin{array}{cc}{K}_1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& K_1\hfill \end{array}\right]+\mathrm{\Omega }\left[\begin{array}{cc}0\hfill & K_{s,1}\hfill \\ \multicolumn{1}{c}{-K_{s,1}}& 0\hfill \end{array}\right])\left\{\begin{array}{c}Z\hfill \\ \multicolumn{1}{c}{Y}\end{array}\right\}\hfill \\ =\left\{\begin{array}{c}{m}_{e,r}{e}_1{\mathrm{\Omega }}^2\cos \mathrm{\Omega }t\hfill \\ \multicolumn{1}{c}{m_{e,r}{e}_1{\mathrm{\Omega }}^2\sin \mathrm{\Omega }t}\end{array}\right\}\hfill \end{array}$$

(4)

where $C_n,C_s,K_s,G$ are external damping, internal damping, skew-symmetric stiffness, and gyroscopic moment, respectively. The damper ring restricts and reduces the whirl of the shaft to an appropriate range, as shown in Figure 1. The motion of the shaft and dry friction damper is similar to the rub-impact in the Jeffcott rotor, as shown in Figure 3, the rub-impact force is expressed as:

$$\begin{array}{c}{p}_n=k_c(\eta -\delta )\hfill \\ p_t=\mu ·p_n\hfill \end{array}$$

(5)

Figure 3.

Rub impact between the shaft and damper.

where $p_n$ is the radial impact force and $p_t$ is the tangential rub force. $\delta$ , $\mu$ , and $k_c$ are radical clearance, friction coefficient, and impact stiffness between the damper ring and the sleeve respectively. Then, they are decomposed in the Z–Y direction of the vibration,

$$\begin{array}{c}{F}_z=\mathrm{\Theta }{k}_c(Z-\mu Y)(1-\delta /\eta )\hfill \\ F_y=\mathrm{\Theta }{k}_c(Z\mu +Y)(1-\delta /\eta )\hfill \end{array}$$

(6)

$$\mathrm{\Theta }=\{\begin{array}{cc}1\hfill & \eta \ge \delta \hfill \\ \multicolumn{1}{c}{0}& \eta <\delta \hfill \end{array}\begin{array}{c}\eta =\sqrt{Z^2+Y^2}\hfill \end{array}$$

The dimensionless equations are shown as

$$\begin{array}{c}z″+2({\overline{\zeta }}_s+{\overline{\zeta }}_n)z\prime -\omega \overline{g}y\prime +z+\omega {\overline{k}}_{s}y\hfill \\ +\mathrm{\Theta }{\overline{k}}_c(1-1/\rho )(z-\mu y)=\overline{e}{\omega }^2\cos (\omega \tau )\hfill \\ y″+2({\overline{\zeta }}_s+{\overline{\zeta }}_n)y\prime +\omega \overline{g}z\prime +y-\omega {\overline{k}}_{s}z\hfill \\ +\mathrm{\Theta }{\overline{k}}_c(1-1/\rho )(\mu z+y)=\overline{e}{\omega }^2\sin (\omega \tau )\hfill \end{array}$$

(7)

where non-dimensional variables are defined as

$$\begin{array}{c}{\mathrm{\Omega }}_s^2=\frac{K_1}{M_1}\omega =\frac{\mathrm{\Omega }}{{\mathrm{\Omega }}_s}{\overline{\zeta }}_s=\frac{C_{s,1}}{2M_1{\mathrm{\Omega }}_s}{\overline{\zeta }}_n=\frac{C_{n,1}}{2M_1{\mathrm{\Omega }}_s}\hfill \\ \overline{g}=\frac{G_1}{M_1}{\overline{k}}_s=\frac{K_s}{{\mathrm{\Omega }}_s{M}_1}{\overline{k}}_c=\frac{K_c}{K}\tau ={\mathrm{\Omega }}_{s}tz=\frac{Z}{\delta }\hfill \\ y=\frac{Y}{\delta }\overline{e}=\frac{e}{\delta }\rho =\frac{\eta }{\delta }=\sqrt{z^2+y^2}\hfill \end{array}$$

(8)

Solution existence region

We assume that the steady-state periodic solutions of equation (7) are the form as follow:

$$\begin{array}{c}z=A\cos (\omega \tau -\varphi )\hfill \\ y=A\sin (\omega \tau -\varphi )\hfill \end{array}$$

(9)

The unbalanced response increases significantly with the spin speed of the shaft after tailrotor drivelines start running. Noting that the rubbing force is zero until the shaft contacts the damper ring, and $2{\zeta }_s={\overline{k}}_s$ , the amplitude of unbalanced response can be obtained as follow by substituting equation (9) into equation (7) and set $\mathrm{\Theta }$ as zero.

$$A=\frac{\overline{e}{\omega }^2}{\sqrt{{(1-{\omega }^2-\overline{g}{\omega }^2)}^2+4{{\overline{\zeta }}_n}^2{\omega }^2}}$$

(10)

Differentiating $A^2$ with respect to ${\omega }^2$ , it obtains as:

$$\frac{d{A}^2}{d{\omega }^2}=\frac{2{\overline{e}}^2{\omega }^2[1+{\omega }^2(2{{\overline{\zeta }}_n}^2-1-\overline{g})]}{{[4{\omega }^2{{\overline{\zeta }}_n}^2+{({\omega }^2+\overline{g}{\omega }^2-1)}^2]}^2}$$

(11)

If ${\overline{\zeta }}_n>\sqrt{(1+\overline{g})/2}$ , equation (11) is always greater than zero at any angular velocity, $A^2$ is monotonically increasing. However, in practical systems, this situation is not considered because of the small value of non-rotating damping, so ${\overline{\zeta }}_n<\sqrt{(1+\overline{g})/2}$ is assumed. Under this assumption, $A^2$ increases first and then decreases, when $\omega =1/\sqrt{1+\overline{g}-2{{\overline{\zeta }}_n}^2}$ , $A^2$ reaches its peak value $A_{max}=\overline{e}/2{\zeta }_n\sqrt{\overline{g}+1-{\overline{\zeta }}_n^2}$ . If $A_{max}<1$ , the shaft will never touch the damper ring. So, if the damper can work near the critical speed, it must satisfy the parametric constraints:

$${\overline{e}}^2>4{\overline{\zeta }}_n^2(\overline{g}+1-{\overline{\zeta }}_n^2)$$

(12)

When $A=1$ , the sleeve is just in contact with the damper ring. The expression of boundary rotation speeds at which rub occurs can be yielded by

$$\begin{array}{c}[{(\overline{g}+1)}^2-{\overline{e}}^2]{\omega }^4+(4{{\overline{\zeta }}_n}^2-2-2\overline{g}){\omega }^2+1\hfill \\ =[{(\overline{g}+1)}^2-{\overline{e}}^2]({\omega }^2-{\omega }_1^2)({\omega }^2-{\omega }_2^2)=0\hfill \end{array}$$

(13)

In this paper, $\overline{e}$ , $\overline{g}$ is much smaller than 1 in magnitude, so $(\overline{g}+1{)}^2-{\overline{e}}^2>0$ , two positive boundary speeds ${\omega }_1$ , ${\omega }_2$ ( $0<{\omega }_1<{\omega }_2$ ) can be derived by equation (13). Rub occurs during $\omega \ge {\omega }_1$ and dry friction damper works. Setting $\mathrm{\Theta }=1$ and substituting equation (9) into equation (7), then its synchronous full annual rub solution about amplitude is yielded by:

$$a_1{\kappa }^2+a_2\kappa +a_3=0$$

(14)

where

$$\begin{array}{c}\kappa =A-1\hfill \\ a_1=({\overline{k}}_c\mu +2{\overline{\zeta }}_n\omega {)}^2+({\overline{k}}_c-\overline{g}{\omega }^2-{\omega }^2+1{)}^2\hfill \\ a_2=2{\overline{\zeta }}_n\omega ({\overline{k}}_c\mu +2{\zeta }_n\omega )\hfill \\ -(\overline{g}{\omega }^2+{\omega }^2-1)({\overline{k}}_c-\overline{g}{\omega }^2-{\omega }^2+1)\hfill \\ a_3=[{(\overline{g}+1)}^2-{\overline{e}}^2]{\omega }^4+(4{{\overline{\zeta }}_n}^2-2-2\overline{g}){\omega }^2\hfill \\ +1=[{(\overline{g}+1)}^2-{\overline{e}}^2]({\omega }^2-{\omega }_1^2)({\omega }^2-{\omega }_2^2)\hfill \end{array}$$

Viscous internal damping does not enter the equation about synchronous full annual rub solution. $\kappa$ must be a real number and only $\kappa >0$ to assure that actual steady-state periodic solutions exist. a₁ is always positive, and $(\overline{g}+1{)}^2>{\overline{e}}^2$ since $\overline{e}$ is negligible contrasting with 1 in magnitude. The expression of a₃ is the same as equation (13), so the following two cases of roots are given:

1. If $a_3<0$ , then $a_2^2-a_1{a}_3>0$ , ${\kappa }_1>0$ , ${\kappa }_2<0$ , so there have only one solution ${\kappa }_1$ , and the condition for $a_3<0$ can be obtained from equation (13), that is, ${\omega }_1<\omega <{\omega }_2$ .

2. If $a_3>0$ , $-a_2-\sqrt{a_2^2-a_1{a}_3}>0$ and $a_2^2-a_1{a}_3\ge 0$ , then ${\kappa }_1>{\kappa }_2>0$ , there have two solutions ${\kappa }_1,{\kappa }_2$ . It can be simplified into two conditions $a_3>0$ , $a_2^2-a_1{a}_3\ge 0$ . The condition to satisfy $a_3>0$ is that $\omega <{\omega }_1$ or $\omega >{\omega }_2$ , but $\omega <{\omega }_1$ should be abandoned. Aim to analyze the second condition $a_2^2-a_1{a}_3\ge 0$ , express it as the function of ${\omega }^2$ .

$$\begin{array}{c}{f}_1({\omega }^2)=a_2^2-a_1{a}_3={\overline{e}}^2{\omega }^4\hfill \\ [{({\overline{k}}_c\mu +2{\overline{\zeta }}_n\omega )}^2+{({\overline{k}}_c-\overline{g}{\omega }^2-{\omega }^2+1)}^2]\hfill \\ -{[2{\overline{k}}_c\omega {\zeta }_n+{\overline{k}}_c\mu (\overline{g}{\omega }^2+{\omega }^2-1)]}^2\ge 0\hfill \end{array}$$

(15)

In this work, the eccentricity of the rotor system is strictly controlled, and the gyroscopic moment is very small because of the structure of the shaft itself, while the impact stiffness is relatively large by orders of magnitude, moreover, this study mainly concentrates on the first critical speed, so $\overline{e}$ , $\overline{g}$ is much smaller than ${\overline{k}}_c$ , and ${\omega }_2\approx 1$ , the following inequality holds.

$${{\omega }_2}^2<\frac{1+{\overline{k}}_c}{1+\overline{g}}$$

Substituting ${\omega }_2$ into $f_1({\omega }^2)$ , then $a_3=0$ , it makes $f_1({\omega }_2^2)>0$ . Substituting $(1+{\overline{k}}_c)/(1+\overline{g})$ into $f_1({\omega }^2)$ , it makes

$$f_1(\frac{1+{\overline{k}}_c}{1+\overline{g}})=({\overline{k}}_c-\overline{g}{\omega }^2-{\omega }^2+1{)}^2({\overline{e}}^2{\omega }^4-{\overline{k}}_c^2)<0$$

(16)

Because of the continuity of $f_1({\omega }^2)$ , there is a point ${\omega }_3<(1+{\overline{k}}_c)/(1+\overline{g})$ to makes $f_1({\omega }_3^2)=0$ . Therefore, during ${\omega }_2<\omega <{\omega }_3$ , $f_1({\omega }^2)>0$ . To sum up, ${\kappa }_1$ , that is, $A_1$ exists in the speed range ${\omega }_1<\omega <{\omega }_3$ , while ${\kappa }_2$ , that is, $A_2$ exists in the speed range ${\omega }_2<\omega <{\omega }_3$ .

In this section, the range of parameters in which dry friction dampers can work has been derived, the expression for rub-impact steady-state periodic solutions have been given, and the region of the existence of the solutions have been analyzed and discussed.

Phase difference and stability analysis

Mathematical derivation of phase difference

The phase difference of shaft, which is defined as the angle between the eccentric excitation and the response, is quite different with or without a dry friction damper. Once the rub between the sleeve and damper happens, the change of the phase difference is suppressed and maintained for a wide speed range. The equation governing motion in the Z–Y direction of the system cannot directly reflect the phase relationship and phase migration process, so the response of the axis trajectory is expressed in polar coordinates as follows

$$\begin{array}{c}z=\rho \cos \psi \hfill \\ y=\rho \sin \psi \hfill \end{array}$$

(17)

Substituting it into equation (7) and converting it into polar coordinates, the dimensionless governing equation in polar coordinate get

$$\begin{array}{c}\rho ″-\rho {\psi \prime }^2+2({\overline{\zeta }}_s+{\overline{\zeta }}_n)\rho \prime -\omega \overline{g}\rho \psi \prime +\rho \hfill \\ =\overline{e}{\omega }^2\cos (\omega \tau -\psi )-P_n\hfill \\ 2\rho \prime \psi \prime +\rho \psi ″+2({\overline{\zeta }}_s+{\overline{\zeta }}_n)\rho \psi \prime +\omega \overline{g}\rho \prime -\omega {\overline{k}}_s\rho \hfill \\ =\overline{e}{\omega }^2\sin (\omega \tau -\psi )-P_t\hfill \\ P_n=\mathrm{\Theta }{\overline{k}}_c(\overline{\rho }-1)\begin{array}{c}\mathrm{\Theta }=\hfill \end{array}\{\begin{array}{c}1\hfill \\ \multicolumn{1}{c}{0}\end{array}\begin{array}{c}\overline{\rho }\ge 1\hfill \\ \multicolumn{1}{c}{\overline{\rho }<1}\end{array}\hfill \end{array}$$

(18)

The non-dimensional variables are defined in equation (8). $\varphi =\omega \tau -\psi$ is the phase difference between eccentric excitation and the whirl response, as shown in Figure 3. During the rub occurring, $\overline{\rho }″=\overline{\rho }\prime =\psi ″=0$ , $\psi \prime =\omega ,\overline{\rho }=A,\mathrm{\Theta }=1$ , equation (18) becomes:

$$\begin{array}{c}(1-{\omega }^2-{\omega }^2\overline{g})A+{\overline{k}}_c(A-1)=\overline{e}{\omega }^2\cos \varphi \hfill \\ \omega [2({\overline{\zeta }}_s+{\overline{\zeta }}_n)-{\overline{k}}_s]A+\mu {\overline{k}}_c(A-1)=\overline{e}{\omega }^2\sin \varphi \hfill \end{array}$$

(19)

By eliminating the amplitude A from the above formula, one gets

$$d_1\sin \varphi -d_2\cos \varphi -d_3=0$$

(20)

where

$$\begin{array}{c}{d}_1=({\overline{k}}_c-\overline{g}{\omega }^2-{\omega }^2+1)\overline{e}{\omega }^2\hfill \\ d_2=({\overline{k}}_c\mu +2{\overline{\zeta }}_n\omega )\overline{e}{\omega }^2\hfill \\ d_3=2\omega {\overline{k}}_c{\zeta }_n+{\overline{k}}_c\mu (\overline{g}{\omega }^2+{\omega }^2-1)\hfill \end{array}$$

Substituting $\vartheta =\varphi /2$ into equation (20) and doing some calculation, it is yielded by

$$(d_2-d_3){\tan }^2\vartheta +2d_1\tan \vartheta -(d_2+d_3)=0$$

(21)

Then, the corresponding expression for the tangent of the phase difference is obtained as

$$\tan {\vartheta }_1,\tan {\vartheta }_2=\frac{-d_1\pm \sqrt{d_1^2+d_2^2-d_3^2}}{d_2-d_3}$$

(22)

If $\tan \vartheta$ has real number solutions, it must satisfy

$$\begin{array}{c}{d}_1^2+d_2^2-d_3^2={\overline{e}}^2{\omega }^4\hfill \\ [{({\overline{k}}_c\mu +{\overline{\zeta }}_n\omega )}^2+{({\overline{k}}_c-\overline{g}{\omega }^2-{\omega }^2+1)}^2]\hfill \\ -{[2{\overline{k}}_c\omega {\overline{\zeta }}_n+{\overline{k}}_c\mu (\overline{g}{\omega }^2+{\omega }^2-1)]}^2\ge 0\hfill \end{array}$$

(23)

This expression is the same as $f_1({\omega }^2)$ , so $\tan \vartheta$ has real solutions in the speed range ${\omega }_1<\omega <{\omega }_3$ .

The critical position where the phase difference changes, that is, $\varphi =\pi /2$ , is crucial. Therefore, we first analyze the case of $\varphi <\pi /2$ , that is, $\vartheta <\pi /4$ . Therefore, substituting boundary conditions $0<\tan \vartheta <1$ to equation (22) and doing some calculation, then it can be divided into two cases:

1. If $d_2-d_3<0$ , then it must satisfy $0>\sqrt{d_1^2+d_2^2-d_3^2}-d_1>d_2-d_3$ to ensure $0<\tan \vartheta <1$ . After simplification, the condition becomes $(d_1-d_3)(d_2-d_3)<0$ and $(d_2+d_3)(d_2-d_3)<0$ . It is evident that $d_2>0,d_3>0$ , so $d_1-d_3>0$ is the final simplified condition.

2. If $d_2-d_3>0$ , similarly, it can get $d_1-d_3>0$ .

Therefore, no matter under any conditions, $d_1-d_3>0$ is necessary and sufficient condition to satisfy $0<\tan \vartheta <1$ , that is, $\varphi <\pi /2$ . The condition $d_1-d_3>0$ is studied and $d_1-d_3$ is expressed as the function of $\omega$ as follow:

$$\begin{array}{c}\multicolumn{1}{c}{f_2(\omega )=d_1-d_3=-\overline{e}(\overline{g}+1){\omega }^4}\\ +[\overline{e}({\overline{k}}_c+1)-{\overline{k}}_{c}f(\overline{g}+1)){\omega }^2]-2{\overline{k}}_c{\overline{\zeta }}_n\omega +{\overline{k}}_c\mu \hfill \end{array}$$

(24)

$$\begin{array}{c}\multicolumn{1}{c}{f{\prime }_2(\omega )=-4\overline{e}(\overline{g}+1){\omega }^3+2\omega [\overline{e}({\overline{k}}_c+1)}\\ -{\overline{k}}_c\mu (\overline{g}+1)]-2{\overline{k}}_c{\overline{\zeta }}_n\hfill \end{array}$$

(25)

$${f″}_2(\omega )=-12\overline{e}(\overline{g}+1){\omega }^2+2\overline{e}({\overline{k}}_c+1)-2{\overline{k}}_c\mu (\overline{g}+1)$$

(26)

Through equations (24)–(26), the function curve of $f_2(\omega )$ , $f{\prime }_2(\omega )$ , ${f″}_2(\omega )$ can be roughly obtained. Three real roots of $f{\prime }_2(\omega )$ are ${\omega }_a<0<{\omega }_b<{\omega }_c$ , and $f{\prime }_2(0)<0$ . Study on all references of rotor dynamics, the phase difference of the rotor without rub-impact satisfies $\varphi <\pi /2$ within the speed range $\omega <1$ , namely, it satisfies $d_1-d_3>0$ . ${\omega }_1<1$ and ${\omega }_1$ is the boundary speed of the non-rub and the rub-impact, so $\varphi ({\omega }_1)<\pi /2$ , that is, $f{\prime }_2({\omega }_1)>0$ . The range of interest here is $\omega >{\omega }_1$ , which has three cases:

1. If ${\omega }_1>{\omega }_c$ , in the range $\omega >{\omega }_1$ , $f_2(\omega )$ decreases progressively with the increase of $\omega$ until there is $\omega ={\omega }_4$ to make $f_2({\omega }_4)=0$ .

2. If ${\omega }_b<{\omega }_1<{\omega }_c$ , in the range $\omega >{\omega }_1$ , $f_2(\omega )$ increases first and then decreases, there still has $\omega ={\omega }_4$ to make $f_2({\omega }_4)=0$ .

3. If ${\omega }_1<{\omega }_b$ , in the range $\omega >{\omega }_1$ , $f_2(\omega )$ decreases in ${\omega }_1<\omega <{\omega }_a$ , increases in ${\omega }_a<\omega <{\omega }_b$ , and then decreases again in $\omega >{\omega }_c$ , deriving ${\omega }_b$ from equation (24) and then substituting it into equation (25), if $f_2({\omega }_b)>0$ , there is $\omega ={\omega }_4$ as well to make $f_2({\omega }_4)=0$ .

To sum up, in the range of ${\omega }_1<\omega <{\omega }_4$ , if $f_2({\omega }_b)>0$ , then $0<\tan \vartheta <1$ , the phase difference less than $\pi /2$ . When the speed exceeds ${\omega }_4$ , then $\tan \vartheta >1$ , the phase difference greater than $\pi /2$ , but rub-impact can still be maintained until $\omega ={\omega }_3$ , then it converts to no-rub motion.

Parameter criteria for stability and bifurcation

There are two patterns of rub-impact motions. The synchronous full annular rub will turn into partial rub after the occurrence of instability. Therefore, the stability of the solution is studied in this section. For the stability of the nonlinear dynamical system mentioned above, the equation can be approximated by the linear-approximation around its periodic solution, that is, the approximate linear system can be used to analyze the stability of the original system by its eigenvalues. Introducing state vector $\mathbf{X}=[x_1,x_2,x_3,x_4{]}^T$ , where $x_1=z,x_2=y,x_3=z\prime ,x_4=y\prime$ , the state equations are used to describe the governing equations as:

$$\mathbf{X}\prime =\mathbf{f}(\mathbf{X},\tau )=\mathbf{LX}+\mathbf{N}$$

(27)

Where L and N are the coefficient matrix and the nonlinear vector,

$$\begin{array}{c}\mathbf{L}=\left[\begin{array}{cccc}0\hfill & 0\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill \\ \multicolumn{1}{c}{-1}& -\omega {\overline{k}}_s\hfill & -2({\overline{\zeta }}_s+{\overline{\zeta }}_n)\hfill & \omega \overline{g}\hfill \\ \multicolumn{1}{c}{\omega {\overline{k}}_s}& -1\hfill & -\omega \overline{g}\hfill & -2({\overline{\zeta }}_s+{\overline{\zeta }}_n)\hfill \end{array}\right]\hfill \\ \mathbf{N}=\left[\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{-{\overline{k}}_c(1-\epsilon /\eta )(x_1-\mu x_2)+\overline{e}{\omega }^2\cos (\omega \tau )}\\ \multicolumn{1}{c}{-{\overline{k}}_c(1-\epsilon /\eta )(\mu x_1+y_1)+\overline{e}{\omega }^2\sin (\omega \tau )}\end{array}\right]\hfill \end{array}$$

The steady-state periodic solution ${\mathbf{X}}_0=[x_{10},x_{20},x_{30},x_{40}{]}^T$ is assumed, applying a perturbation $\delta \mathbf{X}$ to steady-state periodic solution, the stability of the solution ${\mathbf{X}}_0$ can be transformed into the stability analysis of $\delta \mathbf{X}$ . Linear approximation is used to expand the equation (27) around the periodic solution ${\mathbf{X}}_0$ , yields by:

$$\mathbf{\delta }\mathbf{X}\prime ={\mathbf{D}}_{\mathbf{x}}\mathbf{f}(\mathbf{X},\mathbf{\tau }){|}_{\mathbf{X}={\mathbf{X}}_0}\mathbf{\delta }\mathbf{X}$$

(28)

where $\mathbf{J}={\mathbf{D}}_{\mathbf{x}}\mathbf{f}(\mathbf{X},\tau ){|}_{\mathbf{X}={\mathbf{X}}_0}$ . $\mathbf{J}$ is the Jacobian matrix, which is composed of first-order partial derivatives, $x_{10}=A\cos \theta$ $x_{20}=A\sin \theta$ is assumed for simplification of equation (9), after calculating, $\mathbf{J}$ can be obtained as

$$\mathbf{J}=\left[\begin{array}{cccc}0\hfill & 0\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill \\ \multicolumn{1}{c}{b_1}& b_2\hfill & -2({\overline{\zeta }}_s+{\overline{\zeta }}_n)\hfill & \overline{g}\omega \hfill \\ \multicolumn{1}{c}{b_3}& b_4\hfill & -\overline{g}\omega \hfill & -2({\overline{\zeta }}_s+{\overline{\zeta }}_n)\hfill \end{array}\right]$$

(29)

where

$$\begin{array}{c}{b}_1=-1-{\overline{k}}_c+{\overline{k}}_c/A(1-{\cos }^2\theta +\mu \sin \theta \cos \theta )\hfill \\ b_2=-\omega {\overline{k}}_s+{\overline{k}}_c\mu -{\overline{k}}_c/A(\mu -\mu {\sin }^2\theta +\sin \theta \cos \theta )\hfill \\ b_3=\omega {\overline{k}}_s-{\overline{k}}_c\mu +{\overline{k}}_c/A(\mu -\mu {\cos }^2\theta -\sin \theta \cos \theta )\hfill \\ b_4=-1-{\overline{k}}_c+{\overline{k}}_c/A(1-{\sin }^2\theta -\mu \sin \theta \cos \theta )\hfill \end{array}$$

According to the Floquet theory, the stability of the solution $\mathbf{\delta }\mathbf{X}$ , which also reflects the stability of ${\mathbf{X}}_0$ , can be determined through the examination of the sign of the real parts of the Floquet exponents. ¹⁸ The synchronous full annular rub solution is stable if all the real parts of the eigenvalues of $\mathbf{J}({\mathbf{X}}_0)$ are less than zero. However, the equation (29) is a periodic coefficient matrix, which cannot be directly used to judge the stability. It’s observed that the shaft is relatively static in the rotating frame coordinate. The rotation transformation matrix should be introduced as:

$$\mathbf{R}=\left[\begin{array}{cccc}\cos \theta \hfill & -\sin \theta \hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{\sin \theta }& \cos \theta \hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \cos \theta \hfill & -\sin \theta \hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \sin \theta \hfill & \cos \theta \hfill \end{array}\right]$$

(30)

Then, make the following transformation:

$$\mathbf{\delta }\mathbf{X}=\mathbf{R}\mathbf{\delta }\mathbf{T}$$

(31)

Equation (31) is substituted into equation (28), and yields

$$\mathbf{\delta }\mathbf{T}\prime ={\mathbf{J}}_{\mathbf{R}}\mathbf{\delta }\mathbf{T}$$

where ${\mathbf{J}}_{\mathbf{R}}={\mathbf{R}}^{-1}(\mathbf{JR}-\mathbf{R}\prime )$ , $\mathbf{R}\prime =\partial \mathbf{R}/\partial \tau$

$$\begin{array}{c}{\mathbf{J}}_{\mathbf{R}}=\hfill \\ \left[\begin{array}{cccc}0\hfill & \omega \hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{-\omega }& 0\hfill & 0\hfill & 1\hfill \\ \multicolumn{1}{c}{-{\overline{k}}_c-1}& -\omega {\overline{k}}_s-{\overline{k}}_c\mu (1/A-1)\hfill & -2({\overline{\zeta }}_s+{\overline{\zeta }}_n)\hfill & (\overline{g}+1)\omega \hfill \\ \multicolumn{1}{c}{\omega {\overline{k}}_s-{\overline{k}}_c\mu }& -1+{\overline{k}}_c(1/A-1)\hfill & -(\overline{g}+1)\omega \hfill & -2({\overline{\zeta }}_s+{\overline{\zeta }}_n)\hfill \end{array}\right]\hfill \end{array}$$

(32)

${\mathbf{J}}_{\mathbf{R}}$ is independent of time, so the eigenvalue of ${\mathbf{J}}_{\mathbf{R}}$ can be used to judge the stability of $\mathbf{\delta }\mathbf{T}$ , in addition, the stability of $\mathbf{\delta }\mathbf{X}$ . The characteristic equation can be expressed as

$$c_4{\lambda }^4+c_3{\lambda }^3+c_2{\lambda }^2+c_1\lambda +c_0=0$$

(33)

with

$$\begin{array}{c}{c}_4=1\hfill \\ c_3=4({\overline{\zeta }}_s+{\overline{\zeta }}_n)\hfill \\ c_2=2+4({\overline{\zeta }}_s+{\overline{\zeta }}_n{)}^2+2{\overline{k}}_c+({\overline{g}}^2+1{)}^2{\omega }^2-\frac{{\overline{k}}_c}{A}\hfill \\ c_1=(4({\overline{\zeta }}_s+{\overline{\zeta }}_n)-4{\overline{k}}_s-2\overline{g}{\overline{k}}_s){\omega }^2+2{\overline{k}}_c\mu (2+\overline{g})\omega \hfill \\ +4(1+{\overline{k}}_c)({\overline{\zeta }}_s+{\overline{\zeta }}_n)\hfill \\ -\frac{{\overline{k}}_c\mu (2+\overline{g})\omega +2({\overline{\zeta }}_s+{\overline{\zeta }}_n){\overline{k}}_c}{A}\hfill \\ c_0=(\overline{g}+1{)}^2{\omega }^4+[4{{\overline{\zeta }}_n}^2-2({\overline{k}}_c+1)(\overline{g}+1)]{\omega }^2\hfill \\ +4{\overline{k}}_c\mu {\overline{\zeta }}_n\omega +(1+{\overline{k}}_c{)}^2+{{\overline{k}}^2}_c{\mu }^2\hfill \\ +\frac{{\overline{k}}_c(\overline{g}+1){\omega }^2-2{\overline{k}}_c\mu {\overline{\zeta }}_n\omega -({{\overline{k}}^2}_c+{{\overline{k}}^2}_c{\mu }^2+{\overline{k}}_c)}{A}\hfill \\ \begin{array}{c}=d_1\hfill \end{array}+\frac{d_2}{2A}\hfill \end{array}$$

If ${\overline{k}}_s,{\overline{\zeta }}_s,\overline{g}$ are removed, the characteristic equation without gyroscopic moment, rotating damping, and skew-symmetric stiffness can be obtained, and it is the same as that in Xu et al. ¹⁷ Consequently, it is a general equation composed of the parameters of the system and can also be applied to the classical Jeffcott rotor. According to the stability criterion of Routh–Hurwitz, the necessary and sufficient conditions for the stability of rub-impact solutions are shown as

$$\begin{array}{c}{c}_n>0,c_3{c}_2-c_4{c}_1>0,c_3{c}_2{c}_1-c_4{c}_1^2-c_0{c}_3^2\hfill \\ =c_1(c_3{c}_2-c_1{c}_4)-c_0{c}_3^2>0\hfill \end{array}$$

(34)

By substituting two periodic solutions $A_1,A_2$ into the first condition of equation (34) then obtain $c_0(A_2)=-\sqrt{{(d_2/2)}^2-d_3{d}_1}/A_2<0$ , and $c_0(A_1)=\sqrt{{(d_2/2)}^2-d_3{d}_1}/A_1>0$ , so $A_2$ is not a steady-state periodic solution under any parametric conditions. Therefore, only the stability of the solution $A_1$ is analyzed in the next section of the work. It is known from equation (33) that $A>1$ leads to $c_1>0$ , $c_2>0$ , and $c_3>0$ is obvious in the case of non-underdamped, combining the latter two conditions of equation (34), the stability of periodic solutions is determined by the following unique conditions:

$$c_1(c_3{c}_2-c_1{c}_4)-c_0{c}_3^2>0$$

(35)

According to Zhang et al., ⁵ the boundary of Hopf bifurcation reveals the transition from synchronous full annular rub motion to partial rub motion. Therefore, the stability of the solution and Hopf bifurcation condition can be investigated by combining equation (35), equation (33), and equation (14).

Saddle-node bifurcation occurs if one of the roots of the characteristic equation equals zero, that is

$$c_0=0$$

(36)

Combining equation (36) and equation (14), we get the parameter conditions for the change from the synchronous full annular rub solution to the saddle-node bifurcation, where the rub-impact motion disappears and the response converts to the non-rub periodic response. The rotational speed derived from this parameter condition corresponds to the ${\omega }_3$ described in the previous section.

In this section, the characteristic equation and bifurcation including viscous damping and the gyroscopic moment are got. But what’s more concerned and interesting is the stability regions and boundaries with the change of bifurcation parameters. The variation and setting of these parameters are the key investigative factors to the design of the dry friction damper.

Simulation and discussion

The simulation of this work is based on the practical parameters of the dry friction damper-shaft system on the tailrotor drivelines of a helicopter, as shown in Table 1 by dimensionless treatment. While gyroscopic moment and viscous internal damping are used as control parameters to investigate their effect on the phase characteristics and periodic solutions in detail employing the stability analysis method proposed above.

Table 1.

Dimensionless parameters of the tailrotor drivelines of a helicopter.

Properties	Value
${\overline{k}}_c$	0.433
$\overline{e}$	0.188
${\overline{\zeta }}_n$	2.64e-2
${\overline{\zeta }}_s$	5.28e-3
${\overline{k}}_s$	$2{\overline{\zeta }}_s$
$\overline{g}$	−1.56e-3
$\mu$	0.15

Characteristics and mechanism of phase difference

Figures 4 and 5 illustrate amplitude-frequency and corresponding the phase difference characteristics of the synchronous full annular rub solution derived from equation (14) and equation (20) based on parameters in Table 1. The dotted curves represent there is no damper, and solid curves represent the shaft equipped with the damper under excitation of the eccentricity and rub-impact. The amplitude and the phase difference of the shaft increase slowly when the rotational speed is far from the critical value. Once approaching the critical speed, both increases abruptly until the sleeve contacts the damper ring, where $\omega ={\omega }_1$ . The synchronous full annular rub solution $A_1$ and the phase difference ${\varphi }_1$ appear since parameters satisfy equation (35). Thereafter the growth of amplitude and the phase difference decreases obviously due to the restriction of the damper ring, as red solid curves. Unlike the non-rub motion, the amplitude does not reverse and the phase difference does not reach $\pi /2$ even though crossing the critical speed. Another rub solution $A_2$ appears when $\omega ={\omega }_2$ and varies along the green solid curves, but it is not a steady-state periodic solution as analyzed in Section 3.2. The phase difference ${\varphi }_1$ continues to increase with the speed until $\omega ={\omega }_4$ , where ${\varphi }_1=\pi /2$ . And it is greater than $\pi /2$ after $\omega >{\omega }_4$ .

Figure 4.

Response characteristics for $(\kappa +1)-\omega$ .

Figure 5.

Response characteristics for $\varphi -\omega$ .

To validate the theoretical formula, we plot the change of the value of $f_1(\omega ),f_2(\omega ),f{\prime }_2(\omega )$ with ω in Figure 6. $f_1({\omega }^2)>0$ is the necessary condition for $f_2({\omega }^2)>0$ , so ${\omega }_4<{\omega }_3$ , as mention in Section 3.1. $f{\prime }_2(\omega )$ has three real roots. However, only ${\omega }_c$ is included in Figure 5 due to the length limitation of it. ${\omega }_1>{\omega }_c$ is conformed with case 1, $f_2({\omega }^2)$ decreases until $\omega ={\omega }_4$ , where $f_2({\omega }_4)=0$ .

Figure 6.

Value of $f_1(\omega ),f_2(\omega ),f{\prime }_2(\omega )$ with ω.

As mention in Section 2.2, $f_1({\omega }^2)>0$ until there appears a point ${\omega }_3$ to makes $f_1({\omega }_3^2)=0$ . The sleeve and damper ring are keep rubbing until $\omega ={\omega }_3$ , then rub-impact motion convert to non-rub suddenly, the amplitude and the phase difference jumps to dotted curves along the dashed line. The changing trajectory of the amplitude and the phase difference follow curve 1–2–3. Therefore, the theoretical analysis of section 4 is rigorously validated by numerical simulation.

The physical mechanism of the above phenomenon can be explained as: if there is no damper, the phase difference of shaft will migrate from the same direction (less than $\pi /2$ ) to the opposite direction (more than $\pi /2$ ) during passing through the first critical speed. However, if the shaft is equipped with a dry friction damper, the radial whirling amplitude of the shaft which increases with the rotational speed is hindered since the existence of the damper ring. The sleeve in the shaft slides and rotates along the inner surface of the ring, whirling speed can increase synchronously with the rotational speed in a wider speed range. The increase of phase difference is limited by the damper ring as well, so the rotational speed where phase difference migrates from the same direction to the opposite direction is delayed, rub-impact will be maintained for a long time after the speed is greater than the critical speed.

Effects of iViscous internal damping

The viscous internal damping of the composite material shaft is larger than the metal shaft, ¹³ so it is considered as the primary parameter when a dry friction damper is designed for the composite shaft.

Above all, the regions and boundaries in the following figures should be introduced. The amplitude of the shaft isn’t significant enough to contact the damper ring with low speed. The rub-impact occurs as the speed approaches the first critical speed. Therefore region I is the non-rub region, II is the partial rub region and III is the steady-state periodic solution region, that is, the synchronous full annual rub region. HP is the Hopf boundary between region II and III. ω₁, ω₃ are boundary between the non-rub region and the rub region. ω₁, ω₂, ω₃ correspond to the angular velocity in the previous theoretical analysis.

Figure 7 illustrates the effect of internal viscous damping and skew-symmetric stiffness on response characteristics and stability domains. Viscous internal damping of metal materials is trivial compared with composite materials, even ignored in Ozaydin and Cigeroglu. ⁴ Therefore, ${\overline{\zeta }}_s=0,$ ${\overline{\zeta }}_s=0.00538,$ ${\overline{\zeta }}_s=0.0106$ are set, corresponding to Hopf boundary curves HP₀, HP₁, HP₂, respectively. Other parameters are shown in Table 1. It’s found that region III is broadened by the increase of ${\overline{\zeta }}_s$ , namely, Even if friction coefficient increases or external damping decreases, the system can maintain a steady-state periodic solution after adding viscous damping, as shown in Figure 7(a) and (b), respectively. The increase of viscous internal damping plays a positive role in the stability of the periodic solution under coupling of unbalanced force and rub-impact excitation. Viscous internal damping does not affect the boundaries ω₁, ω₂, and ω₃, so all the boundary values coincide.

Figure 7.

Response characteristics of viscous internal damping ${\overline{\zeta }}_s=0,{\overline{\zeta }}_s=0.00538,{\overline{\zeta }}_s=0.0106$ on parameter plane: (a) $\mu -\omega$ and (b) ${\overline{\zeta }}_n-\omega$ .

Figure 7(a) shows that the stability boundary value of friction coefficient decreases at first, and then increases with the increase of frequency ratio $\omega$ , which is consistent with the Zhang et al. ⁵ Reducing the friction coefficient can make the transition from region II to region III, which is considered to be conducive to reducing impact wear. Because repeated impact in region II is easier to cause the material to fall off. ω₁, ω₂ are vertical straight lines, which also be inferred in equation (19) that they don’t change with the friction coefficient.

Conversely, the stability boundary value of external damping firstly increases first, and then decreases with the increase of frequency ratio $\omega$ , as shown in Figure 7(b). If external damping is reduced, the periodic solutions will change to Hopf bifurcation. ω₁, ω₂ are concave inward, which shows that the rub-impact area decreases with the increase of external damping.

In the upper part, the effect of viscous internal damping on periodic solutions has been analyzed, further study focuses on the influence of it on quasi-period solutions. It is difficult to deduce the expression of the solution with $\mathrm{\Theta }=0$ and $\mathrm{\Theta }=1$ together. However, the numerical integration method with variable step solver can be used to draw the dimensionless axis orbit in Figure 8, based on points $p_1,p_2$ in Figure 7(a). Comparing Figure 8(a) and (b), viscous damping does not affect the synchronous full annual rub solution, which makes the same axis orbit. Unbalance and rub-impact produce a composite synchronous excitation which leads to the disappearance of viscous damping and skew-symmetric stiffness terms. While viscous damping changes the partial rubbing response and leads to less chaos, which is denoted by the different axis orbit in Figure 8(c) and (d). The sleeve of the shaft rubs the dry friction damper intermittently, the shaft rotates in the deflected configuration and subject to deformations that change in time, which is exposed to viscous damping force like equation (19) This effect on partial rub response also affects the boundary between it and synchronous full annual rub, which explains the change of stability of the steady-state periodic solution after considering viscous internal damping.

Figure 8.

Axis orbit of shaft in Case 1 and Case 2: (a) P₁, $\mu =0.15$ , ${\overline{\zeta }}_s=0.00538$ , (b) P₁, $\mu =0.15$ , ${\overline{\zeta }}_s=0$ , (c) P₂, $\mu =0.3$ , ${\overline{\zeta }}_s=0.00538$ , and (d) P₂, $\mu =0.32$ , ${\overline{\zeta }}_s=0$ .

Viscous internal damping shows destabilizing effects on free whirl or flutter of shafts in the supercritical range according to Montagnier and Hochard. ²⁵ But the above study can be summarized that under unbalanced force and rub-impact excitation, internal damping improves the stability of periodic solutions in subcritical and supercritical range.

Effects of gyroscopic moment

Due to the structural characteristics of the slender shaft, the influence of the gyroscopic moment seems trivial in this work, but in the multi-disk rotor system, the influence of the gyroscopic torque can not be ignored. Therefore, the effects of the gyroscopic moment are amplified to observe the stability and the phase difference.

Figure 9 is used to demonstrate the effect of the gyroscopic moment on the stability of periodic solution and boundaries. Subscript 0, −, + correspond to $\overline{g}=0,\overline{g}=-0.05,\overline{g}=0.05$ , respectively. Even if $\overline{g}$ is increased or decreased at least 400% to amplify the response that it brings, and the ordinate is set in a very small range, HP₊, HP₀, HP₋ are very close to each other, which demonstrates that the gyroscopic moment has little effect on Hopf bifurcation. It’s found that reducing the value of the gyroscopic moment is beneficial to the stability of the periodic solution, but the effect is negligible. Materials with a low friction coefficient is beneficial to the stability of the periodic solution, as shown in Figure 9(a). Increasing the external damping in bearing and block is also beneficial to the stability, as shown in Figure 9(b). The gyroscopic moment affects ω₁, ω₂, and ω₃. The rub region on the plane, that is, region II and region III, moves to the high-speed area obviously by reducing it. On the contrary, the rub region moves to low-speed area if it is enlarged.

Figure 9.

Response characteristics of the gyroscopic moment $\overline{g}=0,\overline{g}=-0.05,\overline{g}=0.05$ (Corresponding subscript 0, −, +) on parameter plane: (a) $\mu -\omega$ and (b) ${\overline{\zeta }}_n-\omega$ .

The response characteristics of the gyroscopic moment on phase difference is illustrated in Figure 10. The gyroscopic moment affects both the synchronous full annual rub solution and the non-rub solution as previously studied. For the case of non-rub, the increase of the phase difference delays after reducing the gyroscopic moment. The curves move to high-speed region, and the first-order critical speed increases, as displayed by ${\omega }_{c,+},{\omega }_{c,0},{\omega }_{c,-}$ . For the case of the synchronous full annual rub, reducing gyroscopic moment have the same effects, which change the response interval of rub-impact. The increase of ω₁ and ω₃ mean that the rational speeds for getting into rub-impact and getting off rub-impact have been deferred.

Figure 10.

Response characteristics of the gyroscopic moment $\overline{g}=0,\overline{g}=-0.05,\overline{g}=0.05$ (Corresponding subscript 0, −, +) on the phase difference $\varphi$ . Sf: simulation solution without the damper; s1d: simulation solution A₁ with the damper; s2d: simulation solution A₂ with the damper.

Experimental investigation

Test rig design and experimental setup

Considering the structure of the tail drive system of the dry friction damper-shaft system on the tailrotor drivelines of a helicopter, the simulated experiment test rig has been designed by the similarity principle, as shown in Figure 11. The parameters of its components are shown in the Table 2. The dimensionless parameters of the test rig are designed to be as close to the original system as possible, so as to preserve the characteristics of the original system. Horizontal and vertical eddy current transducers are set in two locations. One is the driveshaft near the damper mounting to measure the whirl of the drive shaft, another is the motor output shaft to measure the excitation of the spin of it. In this way, the phase difference is obtained from the difference between them. Discs are designed to be small so as not to interfere with the mode of the shaft. The damper adopts similar materials and structural features as the actual dry friction damper. In order to validate the analytical formula and simulation results and ensure the repeatability of the experiment, experimental studies are conducted on the test rig adopting loading protocol as follows:.

Figure 11.

Diagram and configuration of a simulated experiment test rig for the helicopter tailrotor driveline.

Table 2.

Parameters of components.

Components	Material composition	Value
shaft	steel	Φ10 × L3 1000 mm, Elastic modulus 211 Gpa, Poisson’s ratio 0.31, internal damping 0.001 N s/m
disc	steel	Φ78 × 34 mm, L1 = L4 = 150 mm
sleeve	steel	Φ16 × 20 mm, L2 = 400 mm
damping ring	graphite, POB and PTFE	R 16.6 × T 3.6 × W 5.2 mm
bearing		Left and right: stiffness 80 KN/m, damping 55 N s/m
unbalance in disc		Left and right: 2.5g, eccentricity distance 35 mm
motor and control		constant acceleration from 10 to 2000 rpm in 20 s

Three cases are set:

1. The metal drive shaft is mounted with the disc in the initial position.

2. The metal drive shaft is mounted and two discs are moved 200 mm toward the center, respectively.

3. The composite drive shaft is mounted with the disc in the initial position.

Each case is divided into the following two groups:

1. The response data from star-up through the first critical speed of the shaft with the damper.

2. The response data from star-up through the first critical speed of the shaft without the damper.

Each group needs to collect amplitude and phase difference data separately: The amplitude is based on sampled data at uniform time interval. The phase difference is based on speed tracking and angle domain resampling to record the whirl of the drive shaft at a fixed phase point for each revolution of the motor output shaft.

Simulation validation

Firstly, the amplitude-frequency characteristics from star-up through the first critical speed are depicted in Figure 12 based on data derived from the test and simulation, respectively. The simulation is still adopted the above method based on the dimensionless parameters in Table 3. Experimental curves are plotted with all the collected data of vibration displacement in Y-axis. The amplitude variations in the experiment are similar to the simulation results, but the experimental response under the rub impact of the damper is slightly delayed comparing with the simulation. According to the parameters of the test rig, $a_3$ in equation (14) satisfies $a_3>0$ which meets the case 2 in Section 2.2. Therefore the existence region of $A_1$ is ${\omega }_1<\omega <{\omega }_3$ , which is also consistent with the experimental results. There is no amplitude curve like A₂ in the experiment because A₂ is not a steady-state periodic solution as stated in Section 3.2.

Figure 12.

Experimental response characteristics for $(\kappa +1)-\omega$ . s1d: simulation solution 1 with damper; s2d: simulation solution 2 with damper, sf- simulation solution without the damper, ef: experimental result without a damper; ed: experimental result with damper.

Table 3.

Dimensionless parameters of the test rig.

Properties	Value
${\overline{k}}_c$	2.233
$\overline{e}$	0.185
${\overline{\zeta }}_n$	5.09e-2
${\overline{\zeta }}_s$	1.02e-3
${\overline{k}}_s$	$2{\overline{\zeta }}_s$
$\overline{g}$	−3.24e-3
$\mu$	0.15

The phase-frequency characteristics is depicted in Figure 13. Smooth curves represent simulation results, the experimental results with fluctuations match with them to some extent. Key speed points ${\omega }_1,{\omega }_4,{\omega }_3$ are in the position of theoretical analysis of Section 3.1. In the same way, the phase change is delayed, especially obvious in the region of high speed. On the whole, theoretical and formula derivation, simulation results, and experimental results can be mutually verified in this work, and the above explanation of the physical mechanism gets the test and verification.

Figure 13.

Experimental response characteristics for $\varphi -\omega$ .

Effect of the gyroscopic moment and viscous internal damping validation

Moving the axial position of discs can change the gyroscopic moment of the system, which is used to examine the effect of the gyroscopic moment on response characteristics. The change of it is in a limited range since the disc mass is small. Two discs are moved 200 mm toward the center of the shaft, respectively. Then $\overline{g}$ is changing from −3.24e-3 (subscript 0) to −4.69e-3 (subscript −). Although the variation value is small, the delay of the phase difference can be distinguished as shown in Figure 14. The change of phase difference moves to high-speed region. These results are consistent with the trend in Figure 10, which revalidate response characteristics of the gyroscopic moment on the phase difference.

Figure 14.

Experimental response characteristics of gyroscopic moment $\overline{g}$  = 3.24e-3 (ef0,ed0), $\overline{g}$  = 3.24e-3 (ef-,ed-) on phase difference $\varphi$ .

Viscous internal damping comes from the material of the shaft in the rotating frame. In the above Section 4.2 has simulated and illustrated that viscous damping changes the partial rubbing response and leads to less chaos. To validate the effect of it, a flexible matrix composite shaft (bending elastic modulus 117 Gpa and Poisson’s ratio 0.307) with the same structural dimensions is installed on the test rig instead of a metal shaft, as shown in Figure 15. Viscous internal damping of carbon fiber reinforced laminate composite varies with fiber orientation and the number of layers, but the value of it is much larger than metal anyway. ²⁵ So ${\overline{\zeta }}_s$ get an increase in magnitude. Then the experimental data is used to reveal the response of different ${\overline{\zeta }}_s$ in the plane of the axis orbit.

Figure 15.

Flexible matrix composite shaft.

The axis orbit of metal and composite shafts under the restriction of dry friction damper at dimensionless critical speed is shown in Figure 16. When the composite driveshaft is installed in the test rig, there are slight partial rub between the sleeve on the shaft and the damping ring, showing a circular ring. However, when the metal shaft is installed, serious rubbing occurs, showing an irregular petal shape. The results show clearly that the vibration response changes from chaotic motion to quasi-periodic motion only by the increase of the damping. If the processing, installation error, data noise, and other factors are ignored, these characteristics correspond to Figure 8(c) and (d). The results also demonstrate that the internal damping affects the partial rub. The increasing of the internal damping can alleviate the chaotic effect of the response of the shaft under rub impact.

Figure 16.

Axis orbit of shaft of: (a) composite shaft and (b) metal shaft.

Conclusion

The stability and the phase difference of a shaft mounted a dry friction damper system with effects of viscous internal damping and gyroscopic moment are investigated in this paper. The phase difference characteristics and its mechanism are analyzed in the view of formula derivation. The stability and bifurcation of rub-impact solutions with considering viscous internal damping of the shaft, especially the composite shaft and gyroscopic moment are derived. According to the helicopter tailrotor drivelines, boundaries of synchronous full annual rub motion and partial rub and phase difference are analyzed and discussed by analytic and numerical simulation method.

Rub-impact will be maintained for a long time after the speed is greater than the critical speed because of the limit from the damping ring. The stability of rub-impact solution can be improved by choosing materials of damping ring with low friction coefficient, or increasing external damping in bearing and block. Viscous internal damping affects partial rub motion and improves the stability of synchronous full annual rub motion in the subcritical and supercritical range. While gyroscopic moment has little effect on the stability of the solution, but it affects the rub-impact boundary and the increase of the phase difference. The experimental investigation is conducted in a test rig whose characteristics are similar to the tailrotor driveline of a helicopter. Experimental results have verified the simulation results and theoretical derivation. The trend of the amplitude and the phase difference is consistent with the simulation. The increase of the gyroscopic moment causes movement of the change of phase difference to the high-speed region. The increasing of the internal damping can alleviate the chaotic effect of the response of the shaft under rub impact and proves the stability of synchronous full annual rub solution.