A Method for the Design and Optimization of Nonlinear Tuned Damping Concepts to Mitigate Self-Excited Drill String Vibrations Using Multiple Scales Lindstedt-Poincar é

: In downhole drilling systems, self-excited torsional vibrations caused by the bit-rock interactions can affect the drilling process and lead to the premature failure of components. Especially self-excited oscillations of higher-order modes lead to critical dynamic loads. The slim drill string design and the naturally limited drilled borehole diameter limit the installation space, power supply and lead to numerous potentially critical self-excited torsional modes. Consequently, small and robust passive damping concepts are required. The variety of possible downhole boundary conditions and potential damper designs necessitates analytical solutions for effective damper design and optimization. In this paper, two nonlinear passive damper concepts are investigated regarding design and effectiveness to reduce self-excited high-frequency torsional oscillations in drill string dynamics. Based on a ﬁnite element model of a drill string, a suitable minimal model based on the identiﬁed critical mode is generated and solved analytically using the Multiple Scales Lindstedt-Poincar é (MSLP) method. The advantages of MSLP compared to conventional MS methods are shown for this example. On the basis of the analytical solution, parameter inﬂuences are determined, and design equations are derived. The analytical results are transferred to self-excited drill string vibrations and discussed using time domain simulations of the drill string model. analytical results are simulatively validated using self-excited drill string models. The inﬂuence of the nonlinear stiffness on the stability and shape of the resulting attractors is shown, and three representative dynamic responses of the self-excited structure were discussed. The equations derived from the analytical solutions with the aim of mitigating critical drill string vibrations. The procedure described above is transferable to other ﬁelds of engineering where similar conditions, such as low reactive effects, occur. The presented method and the derived analytical solutions will further be used for the efﬁcient design and optimization of nonlinear dampers. Furthermore, the complex nature of self-excitation due to bit-rock interaction with various uncertainties can be addressed in more detail in the future.


Introduction
In deep drilling applications, unwanted vibrations can reduce reliability, lead to premature failures of components, reduce the drilling efficiency and increase nonproductive time. Around one-third of all failures in drilling applications are related to vibrations [1,2]. The occurring vibrations are divided into axial, lateral and torsional vibrations according to their operating direction. Since the late 1960s, torsional vibrations have been the research focus [3,4]. Torsional vibrations are divided into low-and high-frequency vibrations. Lowfrequency vibrations like stick/slip correspond to the first torsional mode of the drilling system, with frequencies less than 1 Hz. These self-excited vibrations are distinguished according to their origins in bit and string-induced vibrations. Bit-induced vibrations originate from the bit-rock cutting process, while string-induced vibrations are caused by the contact between drilling tools and formation [5]. The bit rock interaction during drilling is a complex process with many uncertainties [6], which can be considered, in various ways, e.g., with hybrid uncertainties [7]. It has been shown that, for the torsional dynamic of the drill string, the energy input can be modeled using a nonlinear torque characteristic at the bit [8]. Due to improved and new measuring tools, high-frequency oscillations have been identified as the cause of numerous drill string failures and have been intensively investigated over the last years with studies [9,10], simulations [8,11] and experiments [12,13]. Especially drilling in hard and dense formations leads to highfrequency torsional oscillations (HFTO) with frequencies between 50 Hz and 500 Hz [14,15].
By assuming a resistance characteristic at the bit (torque characteristic) that decreases with the angular speed (revolutions per minute RPM) (Figure 1), a predictive criterion based on linearization at the operation point ( dTorque dRPM ) is derived [8], wherein D k represents the modal damping, ω 0,k the natural angular frequency and ϕ k the deflection of the mass normalized eigenvector of mode k.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 2 of 23 have been intensively investigated over the last years with studies [9,10], simulations [8,11] and experiments [12,13]. Especially drilling in hard and dense formations leads to high-frequency torsional oscillations (HFTO) with frequencies between 50 Hz and 500 Hz [14,15]. By assuming a resistance characteristic at the bit (torque characteristic) that decreases with the angular speed (revolutions per minute RPM) (Figure 1), a predictive criterion based on linearization at the operation point ( ) is derived [8], wherein represents the modal damping, , the natural angular frequency and the deflection of the mass normalized eigenvector of mode k. Additionally, a maximum amplitude at the bit (φ ) that depends on the revolutions per minute (RPM) in the operating point is derived in [16] and agrees with the observations in [14] that no backward rotation occurs at the bit. This maximum amplitude is valid when one HFTO mode dominates the dynamic of the entire drill string [17,18]. In [17], the interaction between stick/slip and HFTO is observed and [19] analyzed in detail. Using the information about stick/slip, HFTO and the torque characteristic ( Figure 1) stability maps are derived [19] and new strategies to reduce vibrations by adjusting operational parameters like the drilling velocity (RPM) and the downhole force on the bit (in drilling industry: weight on bit (WOB)) are shown [5].
As neither the parameters nor the design of the drill string can be changed arbitrarily, new downhole tools must be developed to reduce critical vibrations. An isolator tool based on the principle of mechanical low-pass filtering is analyzed in testing, simulation and operation [13] and reduces the amplitude in specific areas of the BHA. As the isolator does not mitigate the critical vibration along the entire drill string, further mitigation strategies are necessary. Increasing the damping of a system is a common approach to reduce self-excited vibration [20,21]. In [22], the effects of various friction dampers on the stability of critical HFTO modes are analyzed. The derived analytical results show the importance of damper placement regarding the multiple instable modes and multiple damper positions within the BHA. Due to the large number of variable parameters such as the critical mode (frequency and mode shape), the position of the damper and the design of the damper, analytical solutions with regard to the effective damping range of the damper are necessary for an effective design and optimization of dampers in drill string dynamics. Therefore, as in [22], analytical and semi-analytical solutions regarding the effective damping range of the dampers must be determined. For nonlinear dampers, various methods to approximate is derived in [16] and agrees with the observations in [14] that no backward rotation occurs at the bit. This maximum amplitude is valid when one HFTO mode dominates the dynamic of the entire drill string [17,18]. In [17], the interaction between stick/slip and HFTO is observed and [19] analyzed in detail. Using the information about stick/slip, HFTO and the torque characteristic ( Figure 1) stability maps are derived [19] and new strategies to reduce vibrations by adjusting operational parameters like the drilling velocity (RPM) and the downhole force on the bit (in drilling industry: weight on bit (WOB)) are shown [5]. As neither the parameters nor the design of the drill string can be changed arbitrarily, new downhole tools must be developed to reduce critical vibrations. An isolator tool based on the principle of mechanical low-pass filtering is analyzed in testing, simulation and operation [13] and reduces the amplitude in specific areas of the BHA. As the isolator does not mitigate the critical vibration along the entire drill string, further mitigation strategies are necessary. Increasing the damping of a system is a common approach to reduce selfexcited vibration [20,21]. In [22], the effects of various friction dampers on the stability of critical HFTO modes are analyzed. The derived analytical results show the importance of damper placement regarding the multiple instable modes and multiple damper positions within the BHA. Due to the large number of variable parameters such as the critical mode (frequency and mode shape), the position of the damper and the design of the damper, analytical solutions with regard to the effective damping range of the damper are necessary for an effective design and optimization of dampers in drill string dynamics. Therefore, as in [22], analytical and semi-analytical solutions regarding the effective damping range of the dampers must be determined. For nonlinear dampers, various methods to approximate nonlinearities exist. The harmonic balance method [23][24][25] is one possibility to approximate nonlinear parts of a differential equation that has been used for drill string dynamic [22]. Especially for tuned systems like nonlinear tuned mass dampers, methods like the Multiple Scales (MS) [26] or Lindstedt-Poincaré (LP) [27] are suitable approximation options.
Multiple Scales Lindstedt-Poincaré (MSLP) is the rather new combination of MS and LP to solve nonlinear systems [28].
In the following, the MSLP method is used to determine suitable analytical solutions for various nonlinear tuned damper concepts to reduce high-frequency torsional vibrations within the drill string. First, a suitable finite element (FE) model of a drill string is developed and reduced to one modal degree of freedom representing a critical torsional mode. This minimal model is extended by two tuned nonlinear damper concepts. Similar to a tuned mass damper, the basic structure of both nonlinear dampers consists of an inertia mass that is connected by a linear spring and linear damper to the structure. One nonlinear damper has an additional cubic nonlinear stiffness, and the other, a cubic nonlinear stiffness and a friction contact. Second, the resulting nonlinear models are solved using MSLP. It is shown that in this specific case MSLP is more accurate than other methods, like the MS method. Third, the parameter influences are determined, and the analytical solutions are transferred to a self-excited drill string vibration to realize a robust and optimized design. Finally, the analytical results are compared with time domain simulations of self-excited drill string vibrations using a reduced FE drill string model. Furthermore, influences and design specifications for drill string vibrations are discussed.

Extended Drill String Model
To investigate the effect of the nonlinear dampers on the torsional dynamics of the drill string, a generic drill string model ( Figure 2) was constructed using beam elements. The resulting linear equation of motion M ..
consists of the mass M, damping C and stiffness K matrix, an external force vector f and the angular deviations from the operating point x for a static twist and constant rotary speed.
nonlinearities exist. The harmonic balance method [23][24][25] is one possibility to approximate nonlinear parts of a differential equation that has been used for drill string dynamic [22]. Especially for tuned systems like nonlinear tuned mass dampers, methods like the Multiple Scales (MS) [26] or Lindstedt-Poincaré (LP) [27] are suitable approximation options. Multiple Scales Lindstedt-Poincaré (MSLP) is the rather new combination of MS and LP to solve nonlinear systems [28]. In the following, the MSLP method is used to determine suitable analytical solutions for various nonlinear tuned damper concepts to reduce high-frequency torsional vibrations within the drill string. First, a suitable finite element (FE) model of a drill string is developed and reduced to one modal degree of freedom representing a critical torsional mode. This minimal model is extended by two tuned nonlinear damper concepts. Similar to a tuned mass damper, the basic structure of both nonlinear dampers consists of an inertia mass that is connected by a linear spring and linear damper to the structure. One nonlinear damper has an additional cubic nonlinear stiffness, and the other, a cubic nonlinear stiffness and a friction contact. Second, the resulting nonlinear models are solved using MSLP. It is shown that in this specific case MSLP is more accurate than other methods, like the MS method. Third, the parameter influences are determined, and the analytical solutions are transferred to a self-excited drill string vibration to realize a robust and optimized design. Finally, the analytical results are compared with time domain simulations of self-excited drill string vibrations using a reduced FE drill string model. Furthermore, influences and design specifications for drill string vibrations are discussed.

Extended Drill String Model
To investigate the effect of the nonlinear dampers on the torsional dynamics of the drill string, a generic drill string model ( Figure 2) was constructed using beam elements. The resulting linear equation of motion consists of the mass , damping and stiffness matrix, an external force vector and the angular deviations from the operating point for a static twist and constant rotary speed. The mass and stiffness matrix are determined from the design of the generic bottom hole assembly (BHA) using Young's modulus, density and geometry, while the damping matrix is unknown. Previous studies show that, in many cases, a critical mode dominates the system behavior [17,18]. Therefore, the damping matrix can be determined by the modal damping ratio of the modes. Using the predictive criterion derived in [8] and determining the , value, the BHA can be reduced to a single modal degree of freedom representing the critical high frequency mode. This modal reduction of the torsional dynamic of a drill string is despite its complexity suitable, because when HFTO occur, mostly one critical mode dominates the entire torsional system dynamic of the BHA. In addition, the large but slim design of the drill string with little available installation space leads to low reactive effects of dampers on the torsional dynamics of the drill string. In [8,16,19], a The mass and stiffness matrix are determined from the design of the generic bottom hole assembly (BHA) using Young's modulus, density and geometry, while the damping matrix is unknown. Previous studies show that, in many cases, a critical mode dominates the system behavior [17,18]. Therefore, the damping matrix can be determined by the modal damping ratio of the modes. Using the predictive criterion derived in [8] and determining the S c,k value, the BHA can be reduced to a single modal degree of freedom representing the critical high frequency mode. This modal reduction of the torsional dynamic of a drill string is despite its complexity suitable, because when HFTO occur, mostly one critical mode dominates the entire torsional system dynamic of the BHA. In addition, the large but slim design of the drill string with little available installation space leads to low reactive effects of dampers on the torsional dynamics of the drill string. In [8,16,19], a similar reduction method was used to characterize downhole vibrations and in [5,19,22] 'it is used to investigate vibration reduction strategies.

Minimal Model of a Drill String and Damper
Assuming that the reactive effect of the damper on the mode is negligible due to the limited installation space in deep drilling and the assumption that no other modes near the considered mode or higher harmonics are excited by the nonlinear forces [22], a modal single degree of freedom (sdof) was derived. The modal reduction of the extended drill string model yields to ..
consisting of the angular natural frequency ω 0,i , the modal damping D i , the modal deviation q i of the considered mode i and an external torque M acting at the j-th node (position) with a mass normalized modal amplitude ϕ i,j of the i-th mode at the j-th node. By adding the degree of freedom representing the considered damper, the minimal model yields ..
consisting of the physical degree of freedom of the damper x, the inertia of the damper J, the linear stiffness c and linear damper d connecting the damper and the structure at the j-th node. The nonlinear torque M nl is with a cubic stiffness α for the cubic nonlinearity and for the cubic nonlinearity and the additional frictional contact consisting of a normal force F N , friction radius r and coefficient of friction µ.

Adaption and Simplification of the Minimal Model
With the assumption that the additional damper does not influence the movement of the drill string itself, but only the modal amplitudeq via the energy balance, the vibration of the structure is assumed to be harmonic, q =q cos(ω 0,i t), represents a sdof-model. By introducing the difference coordinates u = x − ϕ i,jq cos(ω 0,i t) and their derivations, the amplitude and phase of the system were directly taken into account in the results due to the considered relative motion. The equation of motion for the damper with cubic nonlinearity is and with cubic nonlinearity and a friction contact is

Multiple Scales Lindstedt-Poincarè
A quite new technique for analyzing nonlinear systems is the Multiple Scales Lindstedt-Poincaré method was first introduced by Pakdemirli in 2009 in [28]. As the method's name already shows, it is a combination of both the Multiple Scales (MS) [26] and the Lindstedt-Poincaré (LP) [26,27] method. A small parameter ε is introduced to perform the MSLP. The modal natural angular frequency of the structure is now expressed as Ω, and, to simplify following equations, all parameters concerning the friction damper are expressed through the summarizing variable µ.
In analogy to the LP method, a new time variable τ = ωt is introduced and further divided according to the MS method in while the dividing of the replacement variable stays the same with and is now only dependent on τ. The introduced frequency ω now has the structure and slightly differs from those introduced in the LP method. The next step is inserting these expressions in the nonlinear differential equation of motion and separating for the orders of ε similar to LP, although, in contrast to LP ω 0 is replaced. The secular terms in the higher order equations have two unspecified expressions, one being the time derivative of the complex amplitude A c out of the MS, the other being the frequency expression out of LP. Setting these terms follows a simple rule: D n A c is set to zero, giving an expression for ω n , unless it is a complex expression; then, ω n is set to zero, giving an expression for D n A c . The main goal is to make sure no complex expressions occur for the frequency ω. The combination of the derivatives of the complex amplitude follows the same steps as those for the MS. The same applies to the combining of the frequencies regarding LP.

Applying MSLP on the Damper with Cubic Nonlinearity
In the following, the MSLP method is applied to the drill string model (Equation (10)) to find the analytical results to determine the steady-state amplitude. A small parameter ε is introduced, and Equation (10) is normalized regarding the inertia J ..
with d * = d ε 2 J , α * = α εJ and q * =q ε 2 . Introducing MSLP coordinates leads to Using the detuning parameter σ, Ω can be written as Ω = ω 1 + ε 2 σ . Now, Equation (5) becomes and simplifies the agreement in analyzing only up to O ε 2 to A separation regarding the order ε leads to Appl. Sci. 2021, 11, 1559 with D n = ∂ ∂τ n representing the n th time derivative. Equation (8) solves with A c being the complex amplitude with its complex conjugates A c and +cc denote the complex conjugate of the previous expression. The secular terms of the equation and the same for the complex conjugate expression. For secular terms in MSLP, D 1 A c is set to zero, giving Now, Equation (20) results in the particular solution recalling again that the homogenous solution is already compensated by the result Equation (10) of Equation (8). The right-hand side of Equation (7) can be written as using the detuning parameter to express Ω. By that, the secular terms of O ε 2 are Setting D 2 A equal to zero leads to a complex expression for ω n , as described above in this case, ω 2 is set to zero, leading to Using the Nayfehs method of recombination [27]: and introducing β = στ 2 − γ can rearrange the equation after a little manipulation where a split into real and imaginary parts conducted for the steady state By using mathematical identities: is obtained and yields the frequency-response curve. with Without ε the frequency-response curve is The solution of the equation of motion is obtained by using the solutions of the differential equations of O(1) and O(ε).

Including Friction in MSLP
Including friction damping in the analysis creates the needs for approaching the sign function that changes the sign regarding the relative angular velocity, which characterizes the friction damping, by, e.g., polynomial expressions. Many ways of this are shown in the literature-for example, by Nayfeh in [26]. One way is to use a Fourier series to approximate the friction term sgn . u . Regarding the analysis, only multiples of the angular natural frequency are considered so that the expression can be simplified using only exp(inω 0 τ 0 ). The term with n = 1 is of interest for the secular terms. Introducing a new variable φ = ω 0 τ 0 shortens Equation (40) to This expression can now be placed in the term of the corresponding order of ε. This way of approximating the sign function also shows that considering analyzing higher orders of ε makes the analysis much more complicated, as then higher terms of the Fourier series are needed. The MSLP analysis of the combination of cubic stiffness and the friction contact (Equation (1) with µ * = µ ε 2 J . By using Ω = ω 1 + ε 2 σ and MSLP coordinates, result in the same for O(1) and O(ε) as the MSLP analysis for the damper with cubic stiffness (Equations (8) and (20)) and the friction damping terms affect solutions of O ε 2 and higher. Solving O ε 2 yields the secular terms where ω 2 is set to zero; otherwise, it would be complex, giving an expression for D 2 A c The term resulting from the Fourier series must now be further analyzed. Therefore, it needs to be considered that can be separated into its real and imaginary parts, leading towards Im : Re : so that the whole expression can be written as i4µ * . This yields, similar to obtaining Equation (30) in the earlier MSLP, the combined derivative of the complex amplitude by where splitting into real and imaginary parts results for the steady state in In general, the frequency-response curve now turns out to be with Equation (22) for ω or without ε and being the same as for the MSLP of the viscous damper. When using the solutions, including a friction damper, special attention should be paid to the stick and stick-slip areas of the friction damper. The friction damper can also be included in the analysis by using its equivalent viscous damping. However, the equivalent damping ratio must be used after the MSLP coordinates are introduced with respect to the different characteristics of friction damping.

Comparing MSLP and MS
Pakdemirli [28,29] has shown that there are minor differences between MS and MSLP solutions for normally excited damped systems with cubic nonlinearities. In this specific example with difference coordinates and an excitation that is dependent on the frequency, the MSLP method is more precise. However, not only the method can affect the accuracy of the solutions, the exponent of the small parameter ε also affects the feasibility of the analysis and, thus, the accuracy. Comparing the results obtained by MS using ε in the first power for the damping and excitation term, in the second power and MSLP with ε in the second power leads in the typical Multiple Scales analysis to (57) and now analyzing only up to O(ε) due to the very complex expressions for O ε 2 ; the function for the frequency-response curve is It is important to recognize that, according to the power of ε in the equation describing the system, the expression for calculating the excitation frequency out of the detuning parameter and the angular natural frequency also changes according to Ω = ω 0 + ε n σ for MS methods or Ω = ω(1 + ε n σ) for MSLP.
The differences in the obtained solutions are shown in Figure 3a,b in comparison to the results simulated in the time domain. MSLP results are accurate even for such large nonlinearities and excitations that subharmonic resonance occurs, MS solutions meanwhile show spurious peaks that are bend to much or in the wrong direction.
It is important to recognize that, according to the power of in the equation describing the system, the expression for calculating the excitation frequency out of the detuning parameter and the angular natural frequency also changes according to Ω = + for MS methods or Ω = (1 + ) for MSLP.
The differences in the obtained solutions are shown in Figure 3a,b in comparison to the results simulated in the time domain. MSLP results are accurate even for such large nonlinearities and excitations that subharmonic resonance occurs, MS solutions meanwhile show spurious peaks that are bend to much or in the wrong direction. In Figure 3a, the solutions of MS do not cover the whole excitation branch in contrast to the MSLP solution, which matches the point representing the numerically obtained solution very accurate. For the strong excitation in Figure 3b, MS solutions even show very differently bended curves than the numerical solution or MSLP. Furthermore, the MS solution for ( ) is bent backwards and, therefore, not just quantitatively but, also, qualitatively inaccurate. Due to the dependency of the excitation force on the excitation frequency in this case, the MSLP method is more suitable to approximate the equation of motion.

Limits of the Operating Range
The deviation of analytical and numerical solutions beyond the effective frequency range of the nonlinear damper, frequencies much higher or lower than the natural frequency of the damper (Figure 3b), are rather high. This is caused in the step of expressing the excitation frequency through the detuning parameter and then sorting the terms for only up to ( ). Thereby, no attention is payed to the Ω -characteristic of the excitation and it becomes the same as a normally excited system. This happens when the right-hand side of Equation (16) is rewritten by with the use of the detuning parameter. When analyzing only up to ( ), this becomes In Figure 3a, the solutions of MS do not cover the whole excitation branch in contrast to the MSLP solution, which matches the point representing the numerically obtained solution very accurate. For the strong excitation in Figure 3b, MS solutions even show very differently bended curves than the numerical solution or MSLP. Furthermore, the MS solution for O ε 2 is bent backwards and, therefore, not just quantitatively but, also, qualitatively inaccurate. Due to the dependency of the excitation force on the excitation frequency in this case, the MSLP method is more suitable to approximate the equation of motion.

Limits of the Operating Range
The deviation of analytical and numerical solutions beyond the effective frequency range of the nonlinear damper, frequencies much higher or lower than the natural frequency of the damper (Figure 3b), are rather high. This is caused in the step of expressing the excitation frequency through the detuning parameter and then sorting the terms for ε only up to O ε 2 . Thereby, no attention is payed to the Ω 2 -characteristic of the excitation and it becomes the same as a normally excited system. This happens when the right-hand side of Equation (5) is rewritten by with the use of the detuning parameter. When analyzing only up to O ε 2 , this becomes where ω 2 , as being independent from the excitation, can simply be treated like another parameter among the others, and a "normal" excitation develops. However, as the damper has no significant effect outside the resonance area and, therefore, only the results for the resonance are of special interest; there is no such limiting factor, in this case. Furthermore, negative values for the nonlinearity parameter α can cause problematic results in the case of weakening nonlinearities by bending the peaks of the frequency-response curve unphysically. Again, this is irrelevant for the considered cases.

Parameter Dependency of the Operating Range
For engineering purposes, it might be interesting to know the maximum amplitude of the system x max =û max + q max to set the parameters respective to the available space. The maximum difference amplitude can be calculated by looking at the square root of the frequency-response Equation (36) in question. Since the maximum point is the point where the expression containing the negative square root goes over into the expression with the positive square root, the maximum point is where both expressions are equal. This means that the square root is zero. Thereby wherein only the positive square roots are being used. An additional friction element yields to a much larger expression for the maximum amplitude: which is again obtained by analyzing the square root expression of the frequency-response equation in question. Using these techniques, a three-dimensional plot can be created, showing how the damping parameter d * and the coefficient of the nonlinearity α * affected the maximum amplitudeû max for a given excitation. As the point of the maximum amplitude is characterized by the square root in Equation (36) being zero, the frequency of the maximum amplitude can be described by which, by the use of Equation (45), can be further simplified to or as an expression for the frequency shift ∂ f of the viscous damper (see Figure 4a).
which, by the use of Equation (63), can be further simplified to  Figure 4b shows again that a higher nonlinearity coefficient and a lower damping coefficient result in a higher maximum amplitude. According to the imagination, a higher nonlinearity coefficient leads to a higher frequency shift. However, an increase of the damping coefficient * decreases this frequency shift as it decreases the maximum amplitude. Thereby, the combination of a very small maximum amplitude and a very large frequency shift can be a mutual contradiction.

Influence and Normalization of the Angular Natural Frequency
Changing the angular natural frequency has a huge effect on the frequency-response curve, because its center changes to another excitation frequency , depending on . This has a large impact due to the characteristic of the right-hand side and its appearance in the argument of the cosine. Wanting to apply an obtained frequency-response curve to another excitation frequency requires not only a simple introduction of a frequency ratio but, also, an adjustment of the other parameters. Therefore, it is necessary to have a system that leads to the same resulting equation. In this angular frequency normalization, the units are used similar to Section 4. The units were considered during normalization; however, note that we omitted the unit label in the following for readability. Figure 4b shows again that a higher nonlinearity coefficient and a lower damping coefficient result in a higher maximum amplitude. According to the imagination, a higher nonlinearity coefficient α leads to a higher frequency shift. However, an increase of the damping coefficient d * decreases this frequency shift as it decreases the maximum amplitude. Thereby, the combination of a very small maximum amplitude and a very large frequency shift can be a mutual contradiction.

Influence and Normalization of the Angular Natural Frequency
Changing the angular natural frequency ω 0 has a huge effect on the frequencyresponse curve, because its center changes to another excitation frequency Ω, depending on ω 0 . This has a large impact due to the Ω 2 characteristic of the right-hand side and its appearance in the argument of the cosine. Wanting to apply an obtained frequency-response curve to another excitation frequency requires not only a simple introduction of a frequency ratio but, also, an adjustment of the other parameters. Therefore, it is necessary to have a system that leads to the same resulting equation. In this angular frequency normalization, the units are used similar to Section 4. The units were considered during normalization; however, note that we omitted the unit label in the following for readability. Using Equation (4) with a new natural frequency according to ω 0,new = nω 0,old , the excitation frequency to achieve the same plot is Ω new = nΩ old . Inserting this and equalizing Equation (4) yields ..
with u being some kind of response related to cos(Ω old nt) as discussed above. Therefore, ..  un and u new = u. So α * new = α * ω 2 0,new and d * new = d * ω 0,new can give the same structure as the original differential equation and, therefore, the same frequencyresponse curve as for ω 0 = 1s −1 . This way of adapting the parameters can also be used when varying one of the others and wanting to achieve the same frequency-response curve.
Again, a consideration of a friction damper requires a deeper understanding of the method. The way as shown above would lead to µ * new = µ * ω 2 0,new , whereas the correct solution is µ * new = µ * ω 3 0,new , because the approach via Fourier series or equivalent damping ration includes an additional division by ω 0 . This must be observed to get the same frequency-response curve.

Influence of the Normalized Parameters on the Dynamic Response of the Damper
Changing the system parameters affects the frequency-response curve and, thereby, the dynamic of the system in different ways. Changing the nonlinear stiffness of the cubic nonlinearity, for example, not only modifies the bending angle of the curve but, also, the length of the branches (Figure 5a).
can give the same structure as the original differential equation and, therefore, the same frequency-response curve as for = 1 . This way of adapting the parameters can also be used when varying one of the others and wanting to achieve the same frequency-response curve. Again, a consideration of a friction damper requires a deeper understanding of the method. The way as shown above would lead to * = * , , whereas the correct solution is * = * , , because the approach via Fourier series or equivalent damping ration includes an additional division by . This must be observed to get the same frequency-response curve.

Influence of the Normalized Parameters on the Dynamic Response of the Damper
Changing the system parameters affects the frequency-response curve and, thereby, the dynamic of the system in different ways. Changing the nonlinear stiffness of the cubic nonlinearity, for example, not only modifies the bending angle of the curve but, also, the length of the branches (Figure 5a). Changing the damping ratio * affects the length of the branch (Figure 5b). A higher damping ratio shortens the frequency-response curve and results in a smaller maximum amplitude.
The same applies when increasing the friction damping coefficient, here represented by * (Figure 6a). It might be worth mentioning that a stronger friction force not only shortens the frequency-response curve but, also, slims down the foot of the branches around the natural frequency stronger than the viscous damping. Changing the damping ratio d * affects the length of the branch (Figure 5b). A higher damping ratio shortens the frequency-response curve and results in a smaller maximum amplitude.
The same applies when increasing the friction damping coefficient, here represented by µ * (Figure 6a). It might be worth mentioning that a stronger friction force not only shortens the frequency-response curve but, also, slims down the foot of the branches around the natural frequency stronger than the viscous damping. Looking at the frequency-response curve in relation to the excitation, there are two parameters that can change. The one is the excitation frequency, and the other is the excitation amplitude. For the excitation amplitude, an increase results in an increase of the maximum amplitude; hence, the length of the curve (Figure 6b).
Looking at the frequency-response curve in relation to the excitation, there are two parameters that can change. The one is the excitation frequency, and the other is the excitation amplitude. For the excitation amplitude, an increase results in an increase of the maximum amplitude; hence, the length of the curve (Figure 6b).
The variation of the inertia J is more complex, because not only the parameters but the natural frequency of the system itself are changed (Figure 7). Looking at the frequency-response curve in relation to the excitation, there are two parameters that can change. The one is the excitation frequency, and the other is the excitation amplitude. For the excitation amplitude, an increase results in an increase of the maximum amplitude; hence, the length of the curve (Figure 6b).
The variation of the inertia is more complex, because not only the parameters but the natural frequency of the system itself are changed (Figure 7).

Transfer to Self-Excited Drill String Vibrations
Section 4.1 shows that the analytical solutions derived using MSLP correspond well with the numerical simulations of the single degree of freedom system in Equation (9) (respectively, Equations (10) and (11)). In the following, the analytical solutions derived using MSLP are transferred to self-excited drill string vibrations using the equivalent damping ratio and discuss using time domain simulations based on Equation (5). Based on this, instructions for using the results in self-excited drill string systems are derived.

Transfer to Self-Excited Drill String Vibrations
Section 4.1 shows that the analytical solutions derived using MSLP correspond well with the numerical simulations of the single degree of freedom system in Equation (9) (respectively, Equations (10) and (1)). In the following, the analytical solutions derived using MSLP are transferred to self-excited drill string vibrations using the equivalent damping ratio and discuss using time domain simulations based on Equation (5). Based on this, instructions for using the results in self-excited drill string systems are derived.

Energy Dissipation and Equivalent Damping
One possibility to transfer analytical solutions, based on the assumption of a constant harmonic amplitude of the underlying structure to self-excited structures, is the energy balance and the related equivalent damping ratio (cf. [22]). The dissipated energy of the damper can be calculated using the frequency-response curve (Equations (23) and (37)) and the equation of the damper motion (Equations (24) and (38)). In the case of a viscous damper with a cubic nonlinearity is the dissipated energy. When analyzing a system containing viscous and friction damping with a cubic nonlinearity, the total dissipated energy consists of the dissipated energy of the viscous damper E d,dis and the dissipated energy of the friction damper E f ric,dis u dt (71) over one period. As described in [22] and shown in Equation (1), the modal damping ratio is decisive for the stability and the resulting amplitude of the self-excited systems. Therefore, the total dissipated energy of the damper E ges,dis is converted by in the equivalent modal damping ratio D eq,mod .

Influence of the Critical Torsional Mode
First, the effect of the modal parameters on the damping provided by the damper is determined using the equivalent damping ratio (Equation (52)). Figure 8 shows the relative angular displacement and equivalent damping ratio over the natural frequency of the underlying structure (excitation frequency of the analytical system) and the angular displacement of the structure (excitation amplitude of the analytical system). For nonlinear systems, several possible amplitudes can occur at a single frequency (stable and unstable solutions), only the maximum stable amplitudes and the corresponding damping ratios are shown in Figure 8.
In contrast to a linear tuned mass damper (TMD), where a constant damping ratio occurs over the entire amplitude range but is only significantly high for a small frequency range, the nonlinear stiffness causes a shift of the characteristic frequency of the damper with the amplitude. At small amplitudes and, thus, small relative displacements between the structure and the damper, the dynamic motion of the nonlinear damper resembles that of the linear TMD. For higher amplitudes, the motion is affected significantly by the nonlinear stiffness and results in a frequency change. By harmonic linearization of the nonlinear stiffness [23][24][25], the linearized stiffness c hl is calculated as a function of the amplitude: resulting in the characteristic relative displacement: that corresponds to the effective damping range of the nonlinear damper (maximum in Figure 8b). When the self-excitation frequency (excitation frequency of the analytical solution) is higher than the characteristic frequency of the linear part of the damper, an increase of the relative amplitude between the damper and structure results in a change of the effective frequency range of the nonlinear damper. This means that the effective frequency range of the nonlinear damper depends on the amplitude but is larger than the frequency range of a linear damper.

Comparison with Time Domain Simulations
The damping diagram in Figure 9a is derived by determining the equivalent damping ratio (Equation (52)) from the analytical solution (Equation (23)) for one natural frequency of the underlying self-excited structure (excitation frequency of the analytical solution). The natural frequency of the underlying structure is higher than the linear natural frequency of the damper to show the results regarding the nonlinear stiffness damper. In addition to the damping ratio provided by the damper, Figure 9a shows the absolute values of three different self-excited damping ratios. For two damping ratios D i = −1% and D i = −1.6% intersections occur between the self-excited absolute damping ratios and the damping provided by the damper. The highest negative damping ratio (D i = −2.1%) exceeds the maximum additional damping provided by the damper. Figure 9b shows the relative displacement corresponding to the damper with respect to the amplitude of the underlying structure and the characteristic relative displacement (Equation (54)).

Influence of the Critical Torsional Mode
First, the effect of the modal parameters on the damping provided by the damper is determined using the equivalent damping ratio (Equation (72)). Figure 8 shows the relative angular displacement and equivalent damping ratio over the natural frequency of the underlying structure (excitation frequency of the analytical system) and the angular displacement of the structure (excitation amplitude of the analytical system). For nonlinear systems, several possible amplitudes can occur at a single frequency (stable and unstable solutions), only the maximum stable amplitudes and the corresponding damping ratios are shown in Figure 8.  For the smallest instability-thus, a linearized negative damping ratio of D i = −1%-the self-excited oscillation is limited and, thus, stabilized by the damping effect of the nonlinear damper. Due to the nonlinear stiffness and the frequency difference between the linear natural frequency of the damper and the vibration frequency of the structure, the damping effect is sufficient to stabilize the system above a certain amplitude. Therefore, at small amplitudes, the amplitude increases exponentially, as usual in self-excited systems (Figure 10a). At a certain amplitude, the damping ratio provided by the damper is sufficient to stabilize the self-excited system. The vibration amplitude of the structure is similar to the amplitude of the intersection in Figure 9a. As shown in Figure 8b the characteristic frequency and, thus, the damping is increased near the characteristic relative displacement. Due to the small, required damping ratio, the relative displacement does not exceed the characteristic relative displacement (Figure 10b). For the smallest instability-thus, a linearized negative damping ratio of = −1%-the self-excited oscillation is limited and, thus, stabilized by the damping effect of the nonlinear damper. Due to the nonlinear stiffness and the frequency difference between the linear natural frequency of the damper and the vibration frequency of the structure, the damping effect is sufficient to stabilize the system above a certain amplitude. There- fore, at small amplitudes, the amplitude increases exponentially, as usual in self-excited systems (Figure 10a). At a certain amplitude, the damping ratio provided by the damper is sufficient to stabilize the self-excited system. The vibration amplitude of the structure is similar to the amplitude of the intersection in Figure 9a. As shown in Figure 8b the characteristic frequency and, thus, the damping is increased near the characteristic relative displacement. Due to the small, required damping ratio, the relative displacement does not exceed the characteristic relative displacement (Figure 10b). For higher self-excitations and, thus, higher negative damping ratios ( = −1.6%), the intersection occurs at higher damping ratios. Figure 11a shows similar to Figure 10a at first an exponential increase of the amplitude at small relative displacements between the damper and the structure. The characteristic frequency of the nonlinear damper for For higher self-excitations and, thus, higher negative damping ratios (D i = −1.6%), the intersection occurs at higher damping ratios. Figure 11a shows similar to Figure 10a at first an exponential increase of the amplitude at small relative displacements between the damper and the structure. The characteristic frequency of the nonlinear damper for this amplitude is below the vibration frequency of the structure. Due to the self-excitation, the amplitude, as well as the relative displacement between the damper and structure, increases. This increase of the relative displacement changes the characteristic frequency of the damper and, thus, increases its damping effect. This increased damping results in a reduction of the vibration amplitude of the structure. The smaller amplitude of the structure results in a smaller relative displacement and, thus, in a reduced damping effect. This repetitive increase and decrease of the amplitude and damping ratio results in a nonperiodic attractor (Figure 11a).
Appl. Sci. 2021, 11, x FOR PEER REVIEW 18 of 23 this amplitude is below the vibration frequency of the structure. Due to the self-excitation, the amplitude, as well as the relative displacement between the damper and structure, increases. This increase of the relative displacement changes the characteristic frequency of the damper and, thus, increases its damping effect. This increased damping results in a reduction of the vibration amplitude of the structure. The smaller amplitude of the structure results in a smaller relative displacement and, thus, in a reduced damping effect. This repetitive increase and decrease of the amplitude and damping ratio results in a nonperiodic attractor (Figure 11a). As Figure 11b shows, the relative displacement amplitude occurs around the characteristic relative amplitude. Similar to the self-excitation with = −1%, the amplitude of the structure is similar to the amplitude of the intersection in Figure 9a. For even higher self-excitations, e.g., = −2.1%, the additional damping is not sufficient to stabilize the unstable mode. Thus, an exponential increase of the amplitude occurs (Figure 12a,b). In Figure 12, a change in the dynamic motion of the structure regarding the characteristic frequency is visible at ≈ 0.2 . As Figure 11b shows, the relative displacement amplitude occurs around the characteristic relative amplitude. Similar to the self-excitation with D i = −1%, the amplitude of the structure is similar to the amplitude of the intersection in Figure 9a. For even higher self-excitations, e.g., D i = −2.1%, the additional damping is not sufficient to stabilize the unstable mode. Thus, an exponential increase of the amplitude occurs (Figure 12a,b). In Figure 12, a change in the dynamic motion of the structure regarding the characteristic frequency is visible at t ≈ 0.2s. As Figure 11b shows, the relative displacement amplitude occurs around the characteristic relative amplitude. Similar to the self-excitation with = −1%, the amplitude of the structure is similar to the amplitude of the intersection in Figure 9a. For even higher self-excitations, e.g., = −2.1%, the additional damping is not sufficient to stabilize the unstable mode. Thus, an exponential increase of the amplitude occurs (Figure 12a,b). In Figure 12, a change in the dynamic motion of the structure regarding the characteristic frequency is visible at ≈ 0.2 .
(a) (b)  Figure 9a shows that the analytical results can be easily applied to self-excited structures. For self-excitations near the maximum of the damping ratio provided by the damper (e.g., −1.9% in Figure 9a), the system is unstable, although according to the analytical solution, the system should be stabilized. Reasons for this small inaccuracy are the assumptions mentioned in Sections 2 and 3, like the constant harmonic oscillation of the structure, the negligible reactive effect of the damper on the structure and the order of MSLP. One additional and very important reason is the combination of multiple possible occurring solitons in combination with the energy input and, thus, change of the amplitude in self-excited structures. This effect is shown in Figure 11. The damping effect is increased when the relative angular displacement is close to the characteristic relative angular displacement. The necessary relative displacement is achieved when the amplitude of the structural vibration increases. This is seen well at higher natural frequencies of the self-excited structure, where the influence of the nonlinearity is high. Figure 13a shows the damping diagram when the natural frequency of the structure is increased to 1500 1 s . Figure 13b shows that the instabilities of D i = −0.8% still resulted in a stable attractor due to the damping provided by the damper (similar to Figure 11). Figure 13a shows additionally to the damping diagram a principal depiction of the "closed loop" process, resulting in the attractor (Figure 13b). For a self-excited structure without significant shocks, the amplitude increases exponentially from a small initial amplitude. Due to the cubic stiffness nonlinearity of the system, the relative displacement and the corresponding damping ratio is low at small amplitudes and, thus, leads to an increase of the amplitude (Figure 13a arrow 1). At some point, an additional stable solution occurs, but the provided damping is low due to the initial small relative displacement (Figure 13a arrow 1). At higher amplitudes, the solution with small relative displacements disappears, and a rapid increase of the relative displacement between damper and structure leads to an increase of the provided damping (Figure 13a arrow 2). In Figure 13, the increased damping is sufficient to stabilize the instable mode and results in a decrease of the amplitude of the structure (Figure 13a arrow 3). This decrease of the structural amplitude leads to an increase of the provided damping and, thus, to an even faster decrease of the amplitude. Finally, the amplitude decreases below the amplitude of the maximum possible damping ratio, which causes the solution with high relative displacements and damping ratios to disappear (Figure 13a arrow 4). Due to the small relative displacement and damping ratio a repeated increase of the amplitude occurs (Figure 13a arrow 1). A slight increase of the instability to D i = −1.2% (Figure 14), although significantly below the theoretical maximum damping ratio, is unstable due to this change in amplitudes and relative displacements.
Again, beginning at the small amplitudes of the underlying structure, the relative displacement and the resulting additional damping is small (Figure 14a arrow 1). Due to the self-excitation, the amplitude increases similar to Figure 13a arrow 1 until the stable solution with small relative displacements disappears. This leads to a rapid increase of the relative displacement and provided damping (Figure 14a arrow 2). This effect is also visible in the time domain simulations in Figure 14b (t ≈ 0.3s). However, at these amplitudes, the maximum of the theoretically provided damping ratio is already exceeded, and the current damping is not sufficient to stabilize the unstable mode. Similar to Figure 12, the amplitude increases further (Figure 14a arrow 3 and Figure 14b (t > 0.3s)), while the provided damping ratio decreases. This means that, in contrast to Figure 8, the maximum amplitude does not correlate with the occurring additional damping but, rather, the smallest stable amplitude.
natural frequency of the structure is increased to 1500 . Figure 13b shows that the instabilities of = −0.8% still resulted in a stable attractor due to the damping provided by the damper (similar to Figure 11). Figure 13a shows additionally to the damping diagram a principal depiction of the "closed loop" process, resulting in the attractor (Figure 13b). For a self-excited structure without significant shocks, the amplitude increases exponentially from a small initial amplitude. Due to the cubic stiffness nonlinearity of the system, the relative displacement and the corresponding damping ratio is low at small amplitudes and, thus, leads to an increase of the amplitude (Figure 13a arrow 1). At some point, an additional stable solution occurs, but the provided damping is low due to the initial small relative displacement (Figure 13a arrow 1). At higher amplitudes, the solution with small relative displacements disappears, and a rapid increase of the relative displacement between damper and structure leads to an increase of the provided damping (Figure 13a arrow 2). In Figure 13, the increased damping is sufficient to stabilize the instable mode and results in a decrease of the amplitude of the structure (Figure 13a arrow 3). This decrease of the structural amplitude leads to an increase of the provided damping and, thus, to an even faster decrease of the amplitude. Finally, the amplitude decreases below the amplitude of the maximum possible damping ratio, which causes the solution with high relative displacements and damping ratios to disappear (Figure 13a arrow 4). Due to the small relative displacement and damping ratio a repeated increase of the amplitude occurs (Figure 13a arrow 1). A slight increase of the instability to = −1.2% (Figure 14), although significantly below the theoretical maximum damping ratio, is unstable due to this change in amplitudes and relative displacements.  Again, beginning at the small amplitudes of the underlying structure, the relative displacement and the resulting additional damping is small (Figure 14a arrow 1). Due to the self-excitation, the amplitude increases similar to Figure 13a arrow 1 until the stable solution with small relative displacements disappears. This leads to a rapid increase of the relative displacement and provided damping (Figure 14a arrow 2). This effect is also vis-

Damper Design for Self-Excited Drill String Vibrations
The analytical results in Figure 8b with the theoretical maximum damping ratio are adjusted using the information about the stable and instable solutions and amplitude change (Figures 13 and 14) to determine a realistic analytical solution ( Figure 15). While Figure 8b shows the maximum theoretic damping ratio, Figure 15 shows the damping corresponding to the smallest stable amplitude. For self-excited systems, Figure  15 is very accurate, because self-excited vibrations that are not influenced by significant shocks or other special effects show an exponential increase of the amplitude, beginning at small amplitudes. As shown in Figure 13, this means that there is an increase in the amplitude and then a rapid increase from low relative displacements to high relative displacements and, thus, from low-to-high damping ratios. Using the analytical results (Section 3), the information on self-excited drill string vibrations (Equation (1)), the maximum amplitudes of high frequency vibrations (Equation (2)) and the transfer to self-excited vibrations (Figures 13-15), it is now possible to efficiently design and optimize nonlinear dampers with respect to the critical drill string modes; the position within the drilling system and the design of the nonlinear damper in terms of inertia, damping and stiffness. The results obtained for the damper with cubic nonlinear stiffness are equally applicable to the other analytical solutions of the damper with cubic stiffness and friction. In addition, the results are transferable to other self-excited systems, where, for example, the reactive effect of the dampers on the dynamics of the structure is negligible due to, e.g., limited installation space.

Conclusions
In this paper, the Multiple Scales Lindstedt-Poincaré method was used to derive analytical solutions for two nonlinear dampers by using a specifically adapted drill string model. The advantages of the MSLP method compared to conventional methods like the Multiple Scales were discussed for this specific self-excited system. In this case, MSLP leads to exceptionally good results due to the frequency dependence of the excitation force in self-excited structure. The derived analytical solutions were transferred to self-excited drilling systems by analyzing the energy output and, thus, the stability of the considered self-excited drill string modes. In contrast to time-consuming time-domain simulations, the resulting equivalent damping ratio allows direct conclusions about the dynamic motion and stability of modes affected by dampers. Due to the variety of potentially critical modes (e.g., frequency and mode shape), of drill string and damper models (e.g., design and position) and the boundary conditions in deep drilling, such analytical solutions are essential for the effective design and optimization of dampers in deep drilling systems. In addition to the investigation of influencing parameters on the dynamic motion of the dampers, the analytical results are simulatively validated using self-excited drill string models. The influence of the nonlinear stiffness on the stability and shape of the resulting attractors is shown, and three representative dynamic responses of the self-excited structure were discussed. The equations derived from the analytical solutions with the aim of While Figure 8b shows the maximum theoretic damping ratio, Figure 15 shows the damping corresponding to the smallest stable amplitude. For self-excited systems, Figure 15 is very accurate, because self-excited vibrations that are not influenced by significant shocks or other special effects show an exponential increase of the amplitude, beginning at small amplitudes. As shown in Figure 13, this means that there is an increase in the amplitude and then a rapid increase from low relative displacements to high relative displacements and, thus, from low-to-high damping ratios. Using the analytical results (Section 3), the information on self-excited drill string vibrations (Equation (1)), the maximum amplitudes of high frequency vibrations (Equation (2)) and the transfer to self-excited vibrations (Figures 13-15), it is now possible to efficiently design and optimize nonlinear dampers with respect to the critical drill string modes; the position within the drilling system and the design of the nonlinear damper in terms of inertia, damping and stiffness. The results obtained for the damper with cubic nonlinear stiffness are equally applicable to the other analytical solutions of the damper with cubic stiffness and friction. In addition, the results are transferable to other self-excited systems, where, for example, the reactive effect of the dampers on the dynamics of the structure is negligible due to, e.g., limited installation space.

Conclusions
In this paper, the Multiple Scales Lindstedt-Poincaré method was used to derive analytical solutions for two nonlinear dampers by using a specifically adapted drill string model. The advantages of the MSLP method compared to conventional methods like the Multiple Scales were discussed for this specific self-excited system. In this case, MSLP leads to exceptionally good results due to the frequency dependence of the excitation force in self-excited structure. The derived analytical solutions were transferred to self-excited drilling systems by analyzing the energy output and, thus, the stability of the considered self-excited drill string modes. In contrast to time-consuming time-domain simulations, the resulting equivalent damping ratio allows direct conclusions about the dynamic motion and stability of modes affected by dampers. Due to the variety of potentially critical modes (e.g., frequency and mode shape), of drill string and damper models (e.g., design and position) and the boundary conditions in deep drilling, such analytical solutions are essential for the effective design and optimization of dampers in deep drilling systems. In addition to the investigation of influencing parameters on the dynamic motion of the dampers, the analytical results are simulatively validated using self-excited drill string models. The influence of the nonlinear stiffness on the stability and shape of the resulting attractors is shown, and three representative dynamic responses of the self-excited structure were discussed. The equations derived from the analytical solutions with the aim of mitigating critical drill string vibrations. The procedure described above is transferable to other fields of engineering where similar conditions, such as low reactive effects, occur. The presented method and the derived analytical solutions will further be used for the efficient design and optimization of nonlinear dampers. Furthermore, the complex nature of self-excitation due to bit-rock interaction with various uncertainties can be addressed in more detail in the future.