Analytical Scheme of Stability Analysis for 4-DoF Mechanical System Subjected to Friction-Induced Vibrations

The stability problem for non-conservative multi-parameter dynamical system is usually associated with labor-intensive calculations, and numerical methods do not always allow one to obtain the desired information. The presence of circulatory forces often leads to the so-called ”destabilization effect” of the system under the influence of small dissipative forces. In this regard, it seems important to develop analytical approaches that make it possible to use a simplified scheme for checking the stability conditions. When obtaining and analyzing stability conditions, the algebra of polynomials and elements of mathematical analysis are applied. To obtain a simplified scheme for checking the stability conditions, an asymptotic method is used. A mechanical system with four degrees of freedom which is under the action of dissipative, potential and non-conservative potential (circulatory) forces is considered. The stability problem of friction-induced vibrations is studying. In the case of weak damping an analytical approach is proposed that makes it possible to simplify the analysis of stability conditions, which, due to the presence of many uncertain parameters, are very cumbersome. With the help of numerical testing, the adequacy of the results obtained for the reduced conditions and full stability conditions was established. The results of the analysis make it possible to single out the ”advantageous” regions in the space of dimensionless parameters, which makes it possible to improve the design of the system to increase its reliability.


Introduction
The history of research on the stability of systems with nonconservative positional (circulatory) forces dates back to Euler, who analyzed the static buckling of elastic compressed rods and formulated a theory of buckling, which can be regarded as the first approximation in solving the stability problem of lightweight load-bearing structures. The problem was updated and attracted much attention of scientists in the field of structural mechanics in the second half of the last century, when flutter was theoretically discovered as a result of the action of so-called follower forces.
In 1952, Hans Ziegler investigated the flutter problem in aerodynamics and considered the model of a double pendulum fixed at one end and loaded with a tangential load at the other end [1]. He discovered a phenomenon with an unexpected property: there is a gap between the stability regions in the case of vanishingly small damping compared to the case when there is no damping. This problem attracted great attention of scientists in the field of structural mechanics in the second half of the last century, and numerous studies were devoted to the theoretical analysis of this problem [2][3][4][5][6]. Circulatory forces are present in many mechanical systems and often lead to undesirable effects such as self-excited vibrations in rotor dynamics [7,8], friction-induced vibrations [9,10], squeal vibration 1 3 in drum brakes [11,12], flutter in aeroelastic systems [13,14], the Levitron system [15] and others.
A wide variety of works has been devoted to studying the different aspects of this phenomena. The paper [16] proposed a physical explanation of the mechanism behind the destabilizing effect of small internal damping in the dynamic stability of Beck's column. Several theorems related to the stability problem of non-conservative linear systems were presented in a paper [17]. In paper [18], some paradoxical phenomena that arise in the study of the dynamic stability of finite-dimensional autonomous mechanical systems were discussed. In [6], systems with two neutrally stable damping levels were considered, at which the system firstly acquires stability and then loses it as the damping level increases. In [19], the Brouwer problem of a heavy particle in a rotating vessel was considered. This problem highlights some important effects that arise in the study of the dynamics of systems with dissipationinduced instabilities. In paper [20] some optimized algorithms for stability analysis of systems with circulatory forces were suggested.
Despite the extensive literature on the stability of mechanical systems with circulatory forces, the analytical approaches that may be applied to multi-parameter systems with several degrees of freedom are not well developed. Usually, the problem is solved using numerical analysis, which gives a very limited understanding on the dynamics of the system: in particular, research using grid search is ineffective due to the presence of a large number of unknown parameters. These unknown parameters may arise at the design stage, in the process of optimization the tuned parameters (stiffness and damping of the dynamic absorber for instance), unknown excitations etc. Besides, such methods are unreliable and too expensive in terms of computation time. The problem is complicated by the fact that there are no universal criteria for making a conclusion based on the general structure of the circulatory forces-a kind of analogs of the Kelvin-Chetaev theorems [21,22].
Many theoretical models have been presented in the literature to describe and predict friction-induced vibrations [15,[23][24][25][26]. Along with the use of very simple models with lumped parameters, more complex finite element models are considered, but experimental verification turns out to be ambiguous. One of the main reasons for this is the high sensitivity of such systems to changes in parameters which can occur under certain conditions and that predictions can be sensitive to parameters that are often neglected as irrelevant [27,28].
In the present paper, we consider the system with four degrees of freedom, which can be considered as a phenomenological model of multi-parameter multi-DoF system [25,26]. Our aim is to suggest an effective approach for analyzing stability conditions for such systems and elicit the mechanical parameters more important in order to enlarge the stability region in parameter space.
The paper is organized as follows. In the next section, some auxiliary notes from algebra are presented. In the following section, the mechanical system under study is described and a general form of stability conditions is given. The next section analyzes the case when the damping in joints is neglected and reduces the system of stability conditions to a single inequality. Finally, in the last section these results are expanded for a case of a weakly damped system, and numerical testing is provided.

Mathematical Annex
In 1914, Lienard and Shipard proposed a criterion [29] that simplifies the well-known Hurwitz criterion. It reads: The necessary and sufficient conditions for a real polynomial to have all roots with negative real parts can be presented in one of the following forms: Here Conditions (2)-(5) have a certain advantage over the Hurwitz criterion, since they contain approximately half as many determinant inequalities as the Hurwitz conditions. Conditions (2) were presented in original paper [29], and conditions (3)-(5) were proved by Gantmakher [30] and may be applied depending on the form of the polynomial f(z). For instance, if the last has the even degree, then conditions (2), (4) are more easy than conditions (3), (5) (there is no need to compute the determinants of highest degree n).
We will also use the criterion for the absence of complex roots in the fourth-degree polynomial.
Pokrovsky criterion [31,32]. For a real polynomial (1) f (z) = a n z n + a n−1 z n−1 + ⋯ + a 1 z + a 0 (a 0 > 0) had four different real roots, it is necessary and sufficient to satisfy the inequalities where dis(P) is discriminant of the polynomial P(z).
In the case of multiple roots, the last inequality becomes the equality.
The above two criteria lead to the following consequence: Statement 1. The polynomial (7) has four different negative roots if and only if the inequalities (8) supplemented by inequalities hold.

Mechanical Model and Motion Equations
The model consists of two mass blocks m 1 and m 2 which are coupled by viscoelastic spring, its damping and stiffness coefficients are denoted as d a and k a respectively. These blocks are held against three moving bands as it is shown in Fig. 1.
It is assumed that the belts move at a constant speed, and the coefficient of friction ̃ on the contacting surfaces is constant; for the sake of simplicity, this coefficient is supposed to be the same for all friction surfaces. Due to the pre-load applied to the system, it is assumed that the masses and bands are always in contact. The classical Coulomb's law is used. In addition, it is assumed that the directions of the tangential frictional forces F t do not change, since the relative velocities between the blocks and the belt are considered to be positive. Thus, the tangential friction forces are proportional to the normal forces F n , namely F t =̃F n . The linearized equations of motion have the following view [25] Here x = (x 1 , x 2 , x 3 , x 4 ) T describes the displacement vector, symmetric matrices M, D, K represent the mass, damping and structural stiffness matrices, respectively; the matrix K represents the influence of frictional forces. These matrices have the following form One can see that matrix K is non-symmetric neither skewsymmetric matrix and in terms of so-called MDGKN structuring [17] may be presented in the form K =̃(K + N) where matrix ̃K is classified as a part of positional forces influence, and ̃N represents the non-conservative positional or circulatory forces.   Fig. 1 Mechanical system (adopted from [25]) For the sake of simplicity, we consider the case of identical masses and two pairs of identical joints, namely Introducing the dimensionless parameters , p, q, , and time according to formulas we can write down the characteristic polynomial Here All coefficients of characteristic polynomial (14) are positive, and stability conditions are where (12) Expressions for Δ 50 and Δ 70 are much more cumbersome and are presented in Appendix.

Stability Conditions for the Case " = 0
Even in considered particular case (12), the expression for Δ 7 is extremely huge. For this reason, we shall use the asymptotic expansion on the small parameter . First of all let us discuss the stability conditions for "generative" system (in mathematical meaning), i.e. we take = 0. With such assumption the characteristic polynomial P( ) becomes biquartic polynomial. The system is marginally stable if all eigenvalues are purely imaginary, in other words, with respect to variable Λ = 2 this polynomial has four negative roots. 1 We write the stability conditions according to Statement 1.
where The inequalities (9) are satisfied for all admissible values of incoming parameters.
Inequalities (18)- (20) determine the region of stability in the space of dimensionless parameters (the region below the surface shown in Fig. 2). Physically, these conditions determine the range of safe values for friction coefficient . This interval can be significantly increased with a suitable (18) choice of coefficients characterizing the stiffness and damping of the hinge connecting the masses m 1 and m 2 .

Remark 1
In fact, if p ≠ p + ,p, then two real roots of polynomial U 3 ( ) are positive. This fact may be proven thoroughly using Sturm theorem or similar mathematical tools, although, technically, such proof is tedious enough and we refer here to image of 3-D surface U 3 ( , , p) = 0 presented in Fig. 2.
Consider now three resultants for each pair of polynomials U j ( ) (j = 1, 2, 3). The corresponding expressions are given by the following formulas (expressions for two others are presented in Appendix). The equality of each resultant to zero may be interpreted geometrically as part of the line belonging to the first quadrant of the parametrical plane p (the curve of fourth order and two curves of eighth order). These curves are shown in Fig. 3.
Let us take a look at Fig. 3a ( , p) is negative.) Now, bringing together these results, we can make a conclusion about hierarchy of all six roots of polynomials U j ( ) with the aim to find the solution of system (18)- (20). Results are presented in Table 1 and Fig. 4.
Let's take a look at the situation in zone 4 for example. Remind, that inequality U 1 > 0 is equivalent to restriction < 1 . The inequality U 2 > 0 is equivalent to condition ∈ (0, 21 )U( 22 , 0) if dis 2 ( , p) > 0 2 and is fulfilled otherwise. Finally, the inequality U 3 > 0 is equivalent to condition ∈ (0, 31 )U( 32 , +∞). Hence, in zone 4 in order to satisfy all three inequalities (18)-(20) we have restriction < 1 . But on subinterval [ 21 , 1 ) the second inequality is broken, and on subinterval [ 31 , 21 ) the third inequality is not valid. At the same time, all three inequalities are fulfilled if ∈ (0, 31 ). In the same manner one can see that in all other zones we also have this limitation on . In other words, system (18)-(20) is equivalent to a single inequality < 31 .

Numerical Testing and Discussion
Now let us consider the case when the weak damping in the joints between the masses is present. Let us turn to the condition Δ 7 > 0. For small values of , its sign is determined by the sign of the expression Δ 70 , which for q = 0 coincides with the sign of U 3 ( ). As one can see from Fig. 5, the fulfillment of condition Δ 7 > 0 entails the fulfillment of two other conditions (16) In addition, numerical experiments show that varying the parameter q within a fairly wide range has a rather weak effect on the boundary of permissible values of (Fig. 6). Thereby, we shall give our conclusion based on the analysis of inequality (20) and thereafter shall verify its correspondence to the stability conditions (16). Let us set a task of finding the optimal values for stiffness ratios , p and damping ratio q which allow to maximize the range of dimensionless parameter . 3 We shall use the asymptotical expansion for Δ 7 on , and as a first step we consider the generating case, i.e. stability conditions for = 0 which are governed by value 31 ( , p). Remind that this is the least root of the quartic polynomial U 3 ( ). Actually, it can be written in explicit form using, for instance, Ferrari's method, but this way is completely impractical, because the resulting expressions are extremely bulky (and strongly irrational), and there is a negligible chance to extract some information from their analysis. For this reason, we consider 31 as implicit function U 3 ( , , p) = 0. At the moment, we cannot separate the least root of this equation,  The possibility of a direct finding the solutions of system (24) is questionable because of high degrees of incoming polynomials. From the other side, the necessary conditions for system equalities (24) be consistent are the equalities to zero two pairs of resultants of these three polynomials. Both resultants taken on variable = are huge enough (polynomial of 28th order on variables , p ), but allow the factorization with the following multipliers and only two last of them may be equal to zero (both are present in the decomposition of resultants). According to Lemma 2, if (3 − 2)p + 2 2 − 3 + 1 = 0, then the motion is unstable (flutter instability), hereby this case may be excluded from further consideration. Finally, we have two candidates for extremum p = 1 ±  √ 2) > 0, one can see that these values have opposite signs. If p > 0, then negative root of 4 corresponds to 31 , and positive root corresponds to 32 (as 31 ≤ 32 according to our nomenclature). In case p < 0 the correspondence is on the contrary.
The key moment here is that point p = 1 − √ 2 is not the point of extremum, because both branches B 1 , B 2 of the curve U 3 ( , p) = 0 are monotonic in the vicinity of this point. Moreover, there are no other points of extremum, thus these branches are monotonic for any value of p (Fig. 7). Thereby, the curve = 31 is composed of two semi-branches-semi-branch B 1 for p ≤ 1 − √ 2 and semi-branch B 2 for p > 1 − √ 2 , and the maximal possible value for 31 is defined according to formula (27).

Remark 2
There is no need to do such verification for point p = 1 + √ 2 , because, according to Lemma 3, this case is related to the roots 33 , 34 which are complex everywhere except this point (the imaginary part turns to zero here). The branches B 1 , B 2 are still monotonic (Fig. 8d-f).
Expression for (0) from formula (27) has a minimum at = 2∕3, , respectively, taking into account the constraint < √ 2∕2, we see that the limit value for increases with decreasing value of . Due to the design features of the system, the value for has some lower constraint ≥ 0 ( cannot be "too" small), therefore, the optimal choice will be = 0 , p = p − = 1 − √ 2 0 . On the other hand, if the parameter p is "moderately tunable", that is, due to the features of design may take values from some limited range [p 1 , p 2 ] which does not contain point p − , then, the optimal choice depends on whether the value of exceed 2∕2 and p − < p 1 , then values which are close to p 1 are preferable-the line = 31 (p) decreases slowly (Fig. 8). If p − > p 2 , p 1 > p Q , then the value p = p 2 is our choice. If p 1 < p Q < p 2 < p − , then quantities 31 (p 1 ) and 31 (p 2 ) should be compared.
In case ≥ √ 2∕2 the choice is more evident, because the function = 31 (p) is monotonic, thus values from the right side of interval [p 1 , p 2 ] are likely to apply (Fig. 8f). Of course, it is necessary to take into account, that this is true for values of which are not close to the value 2/3, because in the vicinity of this point the threshold for is drastically lower comparatively to values ≥ 2 (Fig. 8d, e). Thus, to increase the range of permissible values for friction coefficient ̃, we can offer two options for the structure of the matrix K ∶ (a) the stiffness of the hinges 12, 22 (Fig. 1) is 3-5 times less than the stiffness of the hinges 11, 21, and the stiffness in the joint between the masses is approximately equal to the last one ( p ≈ 1 ); (b) "reverse" distribution-the stiffness of the hinges 12, 22 (Fig. 1) is 3-5 times higher than the stiffness of the hinges 11, 21, and p ≥ 1 ). At the same time, choosing the values for from interval ≈ [0.7, 1.3] is a bad choice, because the threshold for is about ten times lower.

Conclusions
The article considers the stability problem for a multiparameter mechanical system with four degrees of freedom, considered earlier in works [25,26]. On the basis of the algebraic Lienard-Shipard criterion, stability conditions are obtained and analyzed. In the case of weak damping, these conditions can be significantly simplified, which makes it possible to give recommendations about convenient and inconvenient system configurations. In particular, the ratios between the components of the stiffness matrix are indicated, which make it possible to significantly expand the range of permissible values for friction coefficient ̃.
As numerical calculations show, the moderate variation of the damping ratio q has little effect on the maximum allowable value of the parameter . At the same time, the weakening of the friction coefficient in the joint between the masses increases the rate of damping of system oscillations. For instance, at q = 1 the maximum Lyapunov exponent in dimensionless time is −0.0123 ( = 0.2, p = 1, = 0.4), and at q = 0.1 it is equal to −0.0147, i.e. the damping rate is 20 percent higher. In this regard, the future work will be related to the assessment of the Lyapunov exponents of the system (10) and the study of the influence the structure of the matrix D on the rate of damping of system oscillations.
set S 1 = {( , p) ∶ u 21 < 0} has common point (0, 0) with the closure of the set S 2 = {( , p) ∶ dis 2 < 0}, we conclude that one of these sets is the subset of another one.
To determine which one is smaller, let us take the point belonging to the boundary of one of these sets. We take S 1 , because it determines the part of ellipse (more simple then curve of fourth order S 2 ). For instance, let it be point M 1 (0.1, 0.1). Calculating the corresponding value for dis 2 , we have −0.2275 < 0. In addition, along the line segment = p, p ∈ (0, 0.1) the sign of u 21 is opposite to the sign of expression p(3p − 1 + p) + 6p 2 = 10p 2 − 1 < 0. As a consequence, the boundary of set S 1 as well as the set itself belong to the closure of the set S 2 . ◻ Eventually, if the polynomial U 2 ( ) has real roots 21 < 22 , then the coefficient u 21 is negative, hence both these roots are positive. It is obvious, that in the domain dis 2 < 0 the condition (19) is fulfilled for any values of the parameter .

Proof of Lemma 2
To prove the statement of the Lemma 2 it is sufficient to show that fulfillment of the condition U 1 > 0 leads to inequality U 3 < 0. Since parameter p is positive, possible values for are limited by intervals (0, 1/2) and (2/3, 1). Substituting relation (a1) into expressions for U 1 , U 3 and introducing an auxiliary parameter 1 = (2 − 3 ) 2 we have   where Function (a2) does not have any extremum, because the system 3 ∕ = 0, 3 ∕ = 0 is inconsistent. Indeed, the resultant of these partial derivatives may be presented in the form and is positive (thereby, not equal to zero), because the polynomial of ninth order has four pairs of complex conjugate roots and one negative root. Thus, if the expressions in brackets in formula (a4) is not equal to zero, then dis(U 3 ) < 0, and the statement of the Lemma takes place. Now, suppose that p = 1 ± √ 2 . Substituting into expression for U 3 , we can directly find the roots. There is a double positive root and a pair of complex conjugate roots In the particular case p = p + ( ) ≜ 1 + √ 2 there are four positive roots (the last root is of multiplicity two).
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature. The authors did not receive support from any organization for the submitted work.

Conflict of interest
On behalf of all authors, the corresponding author states that there is no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will

3
need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.