Responses of a two degrees-of-freedom system with uncertain parameters in the vicinity of resonance 1:1

We analyze the dynamics of a nonlinear mechanical system under the influence of an external harmonic force. The system consists of a linear oscillator (primary mass) and attached nonlinear dynamic absorber. It is supposed that the frequency of the external force is close to the natural frequency of the main mass. Assuming that the parameters of the system are uncertain, the stability conditions of the stationary regimes of the averaged equations are obtained analytically; these regimes correspond to the quasi-periodic motions of the original input system. An analytical approach to the problem of selecting the parameters of a dynamic absorber is proposed in order to reduce the amplitude of oscillations of the main system. The results obtained are compared with the results of the numerical integration of the equations of the motion with different initial conditions and parameter values.


Introduction
In recent decades, a number of researchers have made significant efforts to solve the problem of dissipating the excessive vibration energy of mechanical systems. Vibration isolation and control is one of the most important areas of engineering research, which aims to prevent the transmission of unwanted vibrations of a basic structure to neighboring systems or to eliminate or reduce excessive vibration in the main system in order to avoid possible damage. Elimination, reduction or isolation of oscillations is the main problem of various industrial and technical practices. The use of a device to reduce the resonant vibrations firstly has been patented by Frahm [1], and later have been presented in more formalized and detailed form by Ormondroyd [2], Den Hartog [2,3] and Brock [4]. Such vibration absorber represents a dashpot consisting of a mass connected with the main system by spring with viscous damping. Through the proper tuning of the absorber (stiffness of the spring and damping coefficient), it is possible to minimize the frequency response in the vicinity of the target resonant frequency. Such a device is addressed in the literature as dynamic vibration absorber (DVA or LDVA) or tuned mass damper (TMD).
As many researchers have noted later on, linear vibration absorbers are effective only in a very narrow band of excitation frequencies, in other words, their frequency robustness is low. To overcome this obstacle and increase the frequency bandwidth, different nonlinear spring characteristics have been introduced. Roberson [5] proposes a dynamic vibration absorber with a weakly nonlinear spring characteristic, i.e., a linear spring with small cubic nonlinearity. It was shown that this nonlinear absorber reduces the vibrations over a larger frequency bandwidth compared to its linear counterpart.
More recently, absorbers with a strongly nonlinear spring characteristic (a spring with cubic nonlinearity) have been investigated. In the work of Gendelman [6], the redistribution of the energy of free oscillations of a system with two degrees of freedom consisting of a coupled linear and nonlinear oscillator was investigated. Vakakis and Gendelman showed [7] that energy transfer between weakly coupled linear and nonlinear oscillators is due to transient resonance capture on a resonant 1:1 manifold. Zhu et al. [8] have studied stability and bifurcations in 2-DOF vibration system with nonlinear damping and nonlinear spring. Gendelman and Staroswetsky [9] have demonstrated that the quasi-periodic response regimes are very typical for a periodically forced linear oscillator with the nonlinear energy sink (NES) attached. Jing and Lang theoretically studied [10] the cubic nonlinear damping in the frequency domain through a dimensionless vibration system model actuated by a harmonic input. A qualitative tuning methodology, which imposes the frequency-energy dependence of the absorber to be identical to that of the nonlinear primary system, was developed by Vigui e and Kerschen [11] to mitigate the impulsive response of an SDOF nonlinear oscillator. Petit et al. [12] have focused on the analysis of energy thresholds in 2-DOF system, below which (thresholds) no efficient vibration reduction is possible. The authors proposed a reformulation of the dynamics of the system, leading to the introduction of a new parameter threshold, and their analysis pointed out the regions in parameter space where efficient vibration reduction can be obtained. A novel NES for energy harvesting has been developed by Kremer and Liu [13]. Its spring possesses a strong nonlinear stiffness with a minimum linear component and low mechanical damping. It has been shown that the system is capable of energy localization as well as broadband vibration absorption and energy harvesting. The article of Yang et al. [14] provided an assessment of the dynamic characteristics of a nonlinear absorber attached to a nonlinear primary oscillator from the point of view of the oscillatory power flow. The power factor and kinetic energy of the nonlinear oscillator were used as performance indices. Various combinations of cubic nonlinearities of damping and stiffness in an oscillator and an absorber were investigated.
In the work of Casalotti and Lacarbonara [15] the method of multiple scales was adopted to investigate the 1:1 internal resonance arising in a two-DOF system composed of a nonlinear oscillator coupled with a nonlinear DVA exhibiting hysteresis. Taleshi et al. [16] investigated the targeted energy transfer from a harmonic-excited nonlinear plate to nonlinear and linear attachments with comparison those both. Cirillo et al. [17] have applied singularity theory for nonlinear resonance of a two-degree-of-freedom mechanical system.
In the present paper, we investigate the dynamics of a 2-DOF system with a nonlinear elastic characteristic and uncertain parameters. The emphasis is placed on the development of an analytical procedure for determining the parameters of the absorber, contributing to the maximum decrease in the amplitude of oscillations of the main system. Numerical calculations are used to verify the results obtained.
The article is structured as follows. In Sect. 2 the equations of motion of a mechanical system are given, which are then reduced to a form suitable for analysis. The method of averaging [45] is applied to simplify the mathematical model and stationary points are found. In Sect. 3, we analyze the stability of equilibria of averaged equations, which correspond to the periodic motions of the mechanical system under study. The necessary and sufficient conditions for asymptotic stability are obtained on the basis of analysis of the linearized system. The joint effect of the damping coefficient and the nonlinear component of the stiffness of DVA onto appearance of the flutter instability is investigated. Section 4 is devoted to studying the influence of the NDVA parameters on the magnitude of the oscillations of the primary mass in the vicinity of resonant frequencies. An analytical approach is proposed for localizing the values of these parameters (tuning the absorber). Section 5 conducts numerical experiments to verify the results obtained. The results of integration of averaged equations and complete equations of motion are compared. Finally, some concluding remarks are provided.

Description of the model
The schematic view of a harmonically excited linear oscillator (primary system or LO) coupled with a nonlinear absorber (secondary system or NDVA) is presented in Fig. 1.
The equations of motion of this mechanical system are where x 1 ðtÞ and x a ðtÞ are the displacements of the harmonically forced primary system and of the damped absorber with a stiffness whose linear and nonlinear coefficients of rigidity are denoted by k lin a and k nonlin a ; respectively. After adding the second equation of (2) to the first one and introducing the new variable x 2 ¼ x a À x 1 the motion equations may be rewritten in simplified form The subject of our study is the case when natural frequency of the primary mass is close to the frequency of the external excitation x: Taking into account that mass of the absorber is usually much smaller than main mass, i.e., the ratio m a =m 1 may be assumed as small parameter; therefore the ratiox 1 ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi is also close to x. We believe that use of thex 1 instead of x 1 is accountable, because the natural frequencies of 2-DOF system are changed slightly, and the magnitude x ¼ x 1 is not the ''pure'' resonant frequency.
Let us introduce the dimensionless parameters and time by formulas With the aim to reduce the number of parameters, we introduce also the relative displacements by formulas and now the motion equations are in the following form Here Fig. 1 Mechanical system: nonlinear absorber attached to linear SDOF primary system where the prime means the derivative on time s. For convenience, the symbol ''s '' over x 1 ; x 2 is subsequently discarded.
With the transformation the following shape of the motion equations (5) appears cos sðu 0 þ vÞ þ sin sðv 0 À uÞ ¼ À sin su þ cos sv: Further, taking the linear combinations À sin s Á ð8:1Þ þ cos sM Á ð8:2Þ and cos s Á ð8:1Þ þ sin sM Á ð8:2Þ; we have the following equations where the right-hand side expressions are defined by formulas þ ð2s 2 À s 4 Þv 3 2 Þ; s 2 ¼ sin 2s; c 2 ¼ cos 2s; s 4 ¼ sin 4s; c 4 ¼ cos 4s: ð10Þ Equation (9) are more cumbersome compared to the original Eq. (1), but they have the advantage that the variables u; v change much more slowly over time s, in contradistinction to ''fast'' variables x; x 0 . Accordingly, with numerical integration, the calculation fault due to the accumulation of error (which may result in stuck overflow) is more likely for system (1), and dealing with the system (9) is more secure (an illustrative example is presented below). Figure 2a refers to time history for x 1 component of system (5) (l ¼ 0:236, h ¼ 0:256; c ¼ À0:04; e ¼ À0:15; , ¼ 0:08Þ. As one can see, it may appear the avalanche-like increasing of amplitude (leads to overflow). Figure b, c are related to time histories for u 1 ; v 1 components for system (9) (with the same values of parameters and corresponding initial values) exhibiting a normal behavior. Though the discussed issue for numerical integration of system (5) does not often happen, but the system (9) is more robust in this respect.
Let us assume that u; v are the slow functions about the time s. Applying the method of averaging [45], we get the simplified version of the system (9) which do not contain the time-dependent terms, as all of them have zero average value over 2p time period. Now we shall find the stationary points of the system (11) which correspond to periodic motions with respect to displacements x 1 ; x 2 . 1 For the sake of simplicity it is appropriate to introduce the complex variables z ¼ u þ iv. Then Eq. (11) takes the following form Here the following notions are introduced 1 In fact, quasi-periodic motions due to averaging procedure.
Thus, with condition z 0 ¼ 0, we have the following system of nonlinear algebraic equations and their conjugate counterparts for z 1 ; z 2 . Depending on the value of b 11 , we have two cases.
(A) b 11 ¼ 0. According to formulas (13) this case may happen only when e ¼ 0; which corresponds to ''pure'' resonant case for system with ''frozen'' absorber, i.e., x ¼ x 0 : Then, taking into account that b 12 ¼ Àl=ð1 þ l 2 Þ 6 ¼ 0; immediately from Eq. (14) we get Thus, the system (13) has the unique constant solution (15). (B) b 11 6 ¼ 0 ðe 6 ¼ 0Þ. Let us express z 1 from the first equation of (14) and substitute it into the second one. It follows from equations Subtracting the second Eq. (16) from the first one, we write down the auxiliary equality Also we subtract the second Eq. (16) with multiplier z 20 from the first equation multiplied by z 20 . As a result, we have another auxiliary equality Coming back from complex variables to real ones we get the corresponding system which is equivalent 2 to the system (14). Expressing the sum u 2 20 þ v 2 20 from the second equation and substituting it in the first one we can express the variable u 20 in terms of v 20 : which leads to cubic equation The number of the real roots of this equation is determined by the sign of discriminant of the polynomial P 1 ðvÞ. Namely, the last has three different real roots if the expression is positive, and it has one real root if D P 1 is negative. The expression (22) can be considered as polynomial of second order with respect to parameter ,; hence, the necessary and sufficient condition for D P 1 to be positive is the following double inequality Here ; and the value of b 11 may be positive or negative, so the multiplier sgnðb 11 Þ guarantees that expression in left-hand side is smaller than the expression in right-hand side. It is easy to notice that , 1 ; , 2 are positive if and only if b 11 D B [ 0 and are negative when b 11 D B \0.
In the present case study the magnitude of e is small (say, jej 0:2), and taking into account formulas (16), (23) one can see that minfj, 1 j;j, 2 jg¼ 2 27 Thus, the hypersurfaces , ¼ minfj, 1 j; j, 2 jg are located close to the hyperplane e ¼ 0 in parameter space fl; e; c; h; ,g. The typical view of the 3Dprojections of the hypersurfaces , ¼ , 1 ; , ¼ , 2 is presented in Fig. 3.

Asymptotic stability conditions
In this section we derive the stability conditions for fixed points of system (12) in analytical form. In the literature, the authors usually underline the need of analysis of the eigenvalue problem for a characteristic matrix and the analysis itself is performed numerically [8,14,19,25,37]. However, in the case of a multiparameter problem, dealing with explicit analytical conditions seems to be important, since a purely numerical analysis may be incomplete. Consideration of the system (11) and introducing small perturbations according to formulas yields the following equations and three dots means the nonlinear terms.
The kmatrix corresponding to the linear part of the system (26) is written as Accordingly, the characteristic polynomial has the following form where r ¼ ,ðu 2 20 þ v 2 20 Þ; , ¼ 4 3 e ,: Note that, unlike the matrix (27), coefficients a j ðj ¼ 0; 4Þ do not include values u 20 ; v 20 separately, but only the combination u 2 20 þ v 2 20 . Therefore, it makes sense to get the equation directly for r instead of the determining Eq. (21) for v 2 : This can be done by excluding the variable v 2 from the system which leads to the condition The notion resðÃ 1 ; Ã 2 Þ here and below means the resultant of the polynomials Ã 1 ; Ã 2 .
Since the polynomials P 0 ; P 1 are, respectively, third and fourth degree in v 2 ; then res 1 is a seventhorder determinant. However, the required condition is quite simple Here The necessary and sufficient conditions for the fact that all roots of the polynomial P char ðkÞ lie in the left half-plane can be obtained by Lienard-Chipart criterion [46]. For a fourth-degree equation (with a 4 [ 0), these conditions are As can be seen from the formulas (28), the condition a 3 [ 0 is satisfied, and it is easy to verify that for that is, the condition a 2 [ 0 is also satisfied. Thus, the inequalities a 0 [ 0; D 3 [ 0 are necessary and sufficient conditions for the exponential stability of the zero solution of the system (26).
In case A the coefficient a 0 is obviously positive, and expression for D 30 takes the view Here the notion c 0 ¼ À 2l 2 ð1þl 2 Þ 2 is used. For given value of l the equation D 30 ¼ 0 describes the second-order surface in parameter space ðh 2 ; ,; cÞ. The Hessian matrix for (33) is indefinite; therefore this surface stands for hyperbolic paraboloid (Fig. 4a). The stability condition may be broken only if In the plane Oc, these inequalities correspond to the sets inside acute angles between straight lines presented in Fig. 4b. Note, that even in case when the inequalities (34) are valid, the region of instability is very small. As the magnitude of c is limited, the threshold level for h 2 is achieved at In other words, if the friction coefficient is not too small (h [ h min ), the motion under study is asymptotically stable.
In case B, since r is an implicit function of the five specified parameters, a ''direct'' analysis of the conditions of stability seems difficult to implement. At the same time conditions and we define the Hopf bifurcation surfaces in the region D par ðh; c; ,; e; lÞ: the passage through the hypersurface a 0 ¼ 0 ðD 3 [ 0Þ determines the divergent loss of stability, and through the hypersurface D 3 ¼ 0 ða 0 [ 0Þ determines the flutter loss of stability. The equation of the hypersurface a 0 ¼ 0 has a fairly simple form Accordingly, the condition is a sufficient condition for the instability of the zero solution of the system (26). Note that after performing the replacement we obtain the equation of a plane curve of the third order The graph of this curve is shown in Fig. 5.

Bifurcation points for averaged equations
Note that polynomial (31) can be used to determine bifurcation points, since it can be considered as the resulting function g of the Lyapunov-Schmidt reduction [47] for system (11) with the variable q ¼ jz 2 j 2 and the bifurcation parameter e: gðq; e; l; h; c; ,Þ ¼ À, The condition for the singularity is the equality og=oq ¼ 0. The following cases are possible: Fig. 5 The curve of third order (42) (A) o 2 g oq 2 ¼ 0 ðhysteresisÞ. 3 Then the following equalities hold Suppose that the mechanical parameters of the main system and the absorber (excluding the nonlinear stiffness) are specified. At the same time the exact values of the frequency and amplitude of the external influence are uncertain; however the intervals of their possible values are known, i.e., x 2 ½x ð1Þ ; x ð2Þ ; F 0 2 ½F 1 ; F 2 . Since it is usually desirable to avoid the appearance of bifurcations for the normal operation of the system, it is possible to determine the interval of ''safe values'' for k nonlin a ; in which the formulas (49) cannot be fulfilled. We illustrate this with an example. Let the system parameters be determined according to Table 1.
The proposed procedure is easy to explain using a geometric representation. Figure 7a shows the curves C 1 : f 1 ¼ h hyst ðeÞ; C 2 : f 2 ¼ 50, hyst ðeÞ.    (20). Finally, v 20 is the root of the cubic equation (21). The roots of this equation can be written down explicitly using the Cardano formula. However, this solution contains the roots of the second and third degrees and depends on the five parameters l; e; h; c; ,: Therefore, it is completely unsuitable for further analysis. For this reason, we suggest the following procedure.
With fixed values of the parameters h ¼ h Ã ; c ¼ c Ã ; , ¼ , Ã ; l ¼ l Ã ; e ¼ e Ã ; the value of v 20 is defined as a root of the polynomial P 1 and, being substituted [taking into account the formulas (20), (14)] into (50), gives the corresponding value n Ã : Since n can be considered as a polynomial P 4 with regard to v 2 and with coefficients depending on the parameters then n Ã satisfies the system P 1 ðv 2 Þ ¼ 0; P 4 ðv 2 Þ ¼ 0, that is, is the root of the following polynomial After substitution of expressions for b 11 ; b 12 ; b 22 [formulas (13)] the resulting formula is rather cumbersome to be given here. Since, by assumption, the ratio X 0 ¼ x=x 1 is close to 1, we can assume that e ( 1; and the specific magnitudes of m 1 ; k 1 ; as well as the ranges of varying of x; F 0 ; are known when setting the task for a given mechanical system (1). It is also obvious that the efficiency of the absorber is the higher, the greater its mass. At the same time, due to constructive considerations, there is usually a top limit on the value m a , that is, we can assume that the value of the parameter l is known (for a particular system). Thus, the task is to find such parameters c; h; , (tuning of the absorber), at which the maximum on e of the function (46) takes the smallest value.
In principle, this problem can be solved numerically, but such a solution can have serious errors, since the limits of variation of various parameters differ significantly. As will be shown below, there is a ''very small'' order parameter (10 À5 ), small order (10 À1 À10 À2 ) parameters and a large order (10 3 À10 4 ) parameter. Therefore, to find the optimal configuration of the absorber, we propose the following combined numerical-analytical approach. For a given value of l as a first approximation for h; c we take the same values as in the case of a linear absorber, using known formulas ( [3,48] or the like). For instance, according to [3] we have the following formulas which imply the following estimations: c DH % À0:01; h DH % 0:12 ðl ¼ 0:1Þ. Constructing the surface P 5 ðn; e; ,Þ ¼ 0, we can make a decision about the appropriate range of the parameter ,. Figure 8a displays the parameter , range is taken ½À 0:1; 0:1. It is clear that such range is too large, because in the middle of the surface there is a notch in which the values of n are much less than 4000. In Fig. 8b this range is reduced of a hundred times, and, obviously, it should be further reduced in order for the peaks at the edges of the pattern to disappear. In Fig. 8c there is no strong change in the height of the upper edge, and a rough estimate of the size can be made as , 2 ½À 0:00007; À 0:00005. Now, taking as a first approximation c ¼ À0:01; , ð1Þ ¼ À0:6 Â 10 À4 , we vary h 2 within such limits so as n not to exceed value % 210 À 220. In such a way we determine the amendment to the taken value: hh 1 ¼ ðh ð1Þ Þ 2 % 0:015 (Fig. 9a) instead of h 2 0 . Substituting hh 1 in (52), we vary c (Fig. 9b). After that, the values found are refined from the condition of equal peaks at the lowest possible value of n (Fig. 9c). Thus, the values found are c ¼ À0:0153; hh ¼ 0:015; , ¼ À0:52 Â 10 À4 .
Verifying the correctness of the result can be easily performed using the expression (52). If there is a set of values c H ; h H ; , H for which n takes a value less than that found (for example, n 1 ¼ 201:1), then the straight line n ¼ 201 intersects the curve P 5 ðn; c H ; h H ; , H ; l; eÞ ¼ 0 at four points, that is, the polynomial P 5 ðn; eÞ j n¼201 has four real roots on the interval ½e 1 ; e 2 . If this equation has two roots e ð1Þ ; e ð2Þ for any  Remark 4.2 If there are two real roots, but, e ð1Þ Á e ð2Þ \0, this means that two peaks are involved in, then a set of parameters ðc H ; h H ; , H Þ has been found for which we have the equality of the peaks. Such a situation is possible from a mathematical point of view, although it can be numerically realized only with some accuracy limit specified in advance.
Summarizing the content of this section, we believe that methodology suggested allows to reduce the range for search of the appropriate values for nonlinear component of stiffness. A possible alternative approach is to solve the optimization problem by a grid search. However, in a multidimensional model space, such approach will be computationally more expensive and scales as N d , where d stands for the dimensionality of the problem. For example, the given problem has four free parameters to optimize ðd ¼ 4Þ. Given the high sensitivity of the function to the parameters, the number of evaluations per grid direction N will be in the range of at least several tens. Thus, the total number of function evaluations required to achieve the desired optimization accuracy using the grid search can be in the range of hundreds of thousands.

Numerical analysis
In this section, the case studies are inspected to check the performance of the absorber attached to primary oscillator. The responses of the system obtained from analytical approximations, which are based on the averaging method, are compared with numerical integrations based on a fourth-order Runge-Kutta method. In order to facilitate the comparison with the results of Sects. 3 and 4, we use dimensionless parameters and variables introduced according to the formulas (3). (9) and (11) First, we compare the results of the numerical integration of the equations of the motion and the averaged ones. Namely, we check how well the information obtained from the analysis of the averaged equations corresponds to the behavior of the solutions of the original mechanical system. The values of the parameters were taken as follows: two values for l were taken in turn, namely 0.1 and 0.14, parameters c; h; , are ranging from À 0:2 to 0.2, from 0.05 to 0.5, and from À 0:001 to 0.001, respectively. 4 The Fig. 9 Refinement of the first approximation frequency ratio was varying from 0.5 to 1.5, though the peaks of the responses are situated closer to 1. For chosen values of parameters c; h; , the distance between peaks is proportional to l and at first approximation is close to 2l. Figures 10 and 11 show projections of typical phase trajectories on the planes Ou 1 v 1 (the primary system) and Ou 2 v 2 (the absorber). Figure 10a exhibits the projections of the phase trajectory of the system (9) (solid curve) and averaged Eq. (11) (dot curve). The following initial values are taken: u 1 ð0Þ ¼ 8; v 1 ð0Þ ¼ À3, u 2 ð0Þ ¼ À2:5; v 2 ð0Þ ¼ 3; and the time interval is s 2 ½0; 50. In Fig. 10b, c, the same trajectories are shown separately within the time interval s 2 ½50; 150. Fig. 11 shows the time histories of the components u 1 ðsÞ; v 1 ðsÞ (Fig. 11a, b), as well as the behavior of the absorber (Fig. 11c). Two circumstances should be noted in this regard. The first is that the absorber calms down ostensibly faster than the main mass (compare Figs. 10b and 11c). In fact, this is not so-we recall that the variable e x 2 and, consequently, u 2 ; v 2 are taken with a multiplier l [see formula (4)], so the absorber responses are 10 times more than those shown in the figure. The second circumstance concerns the time histories of coordinates u j ; v j ðj ¼ 1; 2Þ for the complete system of Eq. (9) and the averaged curve shown in Fig. 11. As one can see, the points belonging to the dashed curve accurately reflect the behavior of the function

Comparing the behavior of the trajectories of systems
(the same situation holds for u 2 ðsÞ; v 2 ðsÞ, although they are not shown in the Figure). In the case of the variable u 1 , we see some discrepancy, i.e., the curve corresponding to the averaged equations follows the lower edge of the ''true'' curve u 1 ðsÞ. This discrepancy is explained by the fact that the first equation of the system (9) contains a term sin 2s that is of order 1, but it disappears with averaging (observe that for other equations of the system, this does not occur).
Finally, Fig. 12 shows that the trajectories of the system are tightened to the attractor. In the case of averaged equations, this attractor, according to the results of Sect. 3, is the fixed point (stable focus), and the trajectories are pulled to this point regardless of the  Fig. 10 (a, b); motion of the absorber for s 2 ½0; 100 (c) choice of initial values (of course, taking into account restrictions on the magnitude of initial perturbations). Three different paths are shown in Fig. 12a, which, as can be seen in Fig. 12b, after a certain period of time enters a small neighborhood of a point M 0 ðu 10 ; v 10 ; u 20 ; v 20 Þ. As for Eq. (11), here the limit cycle is an attractor (Fig. 12c, d), that is, the solutions from a certain neighborhood of the point M 0 are asymptotically periodic. 5 As one can see, the shape of this limit cycle is close to the cardioide with the parameter a % 0:5.

Remark 5.1
Here and further, commenting on the figures, we say ''a limit cycle'' for brevity. Of course, we are talking about projections of the limit cycle in 4-dimensional space onto the planes Ou 1 v 1 and Ou 2 v 2 : It is obvious that the limits of applicability of the averaging method (in general) and for this problem (in particular) are bounded. Consider as an illustration the case when Eq. (22) has three real roots. To find the appropriate values of the parameters, we can use the inequalities (23) (or Fig.3). Choosing l ¼ 0:1; c ¼ À0:01; h ¼ 0:01; it follows from the Fig. 3d that we can take, for example, e ¼ À0:01; , ¼ À2 Â 10 À4 . Substituting the latter values into (31), we have r 1 % À0:02957; r 2 % À 0:7709; r 3 % À 1:1017. For a smallest absolute value of sigma, dynamics is predicted by the used averaging procedure. For the averaged system (11), the trajectory tends to the limit point M 01 ðÀ 23:12; À 10:44; 10:90; 11:24Þ, whereas for the full system it is attracted to the limit orbit. However, a fundamental discrepancy already exists for the values r 2 ; r 3 . For the averaged system, the corresponding stationary point is an attractor (stable focus). For instance, with r ¼ r 2 Eq. (28) has two pairs of complex conjugate roots À 0:005 AE 0:666i; À 0:273 Â 10 À4 AE 0:898 Â 10 À2 i, respectively, and the damping rate is very weak (Fig. 13a). However, the trajectory corresponding to Eq. (9) behaves differently (Fig. 13b, c). This is because the amplitude of the oscillations of the main mass is too large (M 02 ðÀ 512:56; À 63:02; 39:89; 4:11Þ), and the amplitude of each individual pulse is also large (Fig. 13b). Figure 14 reports in more detail the initial phase of the perturbed motion. As one can see, from the very beginning the trajectory corresponding to the averaged equations does not reflect the behavior of the trajectory of the system (9) (compare Figs. 11b and 14d).

Influence of the DVA nonlinear stiffness component
This subsection describes some of the features inherent in the nonlinear component k nonlin a , which is an important part of tuning of the DVA and it has been presented in our studying through the parameter ,. As it follows from Sect. 4, in the case of a ''one-sided'' frequency range (e 1 e 2 [ 0) 6 a softening spring should be taken when X 0 \1, and a hardening spring in opposite case. This conclusion is confirmed by numerical calculations. The magnitude of the parameter , is determined according to the procedure described in Sect. 4 and is influenced by mass ratio and frequency range.
For low frequencies ratio (e\0), one deals with a softening spring, and gradual increase in the parameter , values from zero to a certain limit helps to reduce the Fig. 12 Attraction of the phase trajectories to the fixed point of system (11) for different initial perturbations: a time interval s 2 ½0; 100, b time interval [100, 130] ; attractors for system (9) with initial values: ð8; À 2:5; À 3; 3Þ for case (c) and ð6; 5:5; À 5; 0Þ for case (d) amplitude of the oscillations (Fig. 15a, ,ð1Þ ¼ 0; ,ð6Þ ¼ À0:7 Á 10 À4 ). It should be noted that for curves 5, 6 as compared with 2,3, the amplitude decreased by 30À35%. At the same time, the duration of the transition process has doubled and more. Meanwhile, a further increase in the magnitude of , leads to a change in the qualitative picture. Namely, the beats appear (Fig. 15b), and the transition time to the limiting cycle increases significantly (another 2-3 times). The further growth in , leads to instability regime. Although in the left part of the Fig. 15c the amplitude of oscillations is low, but in the right part it is higher than the counterpart shown in Fig. 15a (the curves 5, 6).
The following interesting feature has been detected: when changing the parameter ,, the shape of the limit cycle for the main mass almost does not change (or changes insignificantly with a significant (several times) increase in the parameter module). At the same time, the shape of the limit cycle for the absorber has been greatly modified. Our numerical experiments have yielded three qualitatively different forms: (A) a If the frequency interval is two-sided (e 1 e 2 \0), then it is advisable to use a hardening spring, and the j,j should be taken 8-10 times smaller in comparison with the previous case. For this reason, the amplitude of oscillations, generally speaking, increases, including for positive values of e, but this is a kind of compensation to protect the primary system against much greater growth of oscillations amplitude with possible negative values e.
Also, the following important circumstance should be taking into account. As noted above, one of the key points of this paper is based the assumption that the parameters of the main system and the external excitation are uncertain. That is, when solving the problem of setting up the device to reduce oscillations of a particular mechanical system, we can assume that the quantities m 1 ; k 1 ; e 1 ; e 2 ; (and m a ) are known. At the same time, unlike c; h, the parameter , depends significantly on F 0 . If its magnitude is known accurately enough, then the value of , is easily transformed into a value for k nonlin a , as predicted. However, a situation may appear when the interval ½F 1 ; F 2 of a possible change of F 0 (within the framework of one task) is not narrow; say F 2 is of 5-10 times more than F 1 . In this case, to determine the value of k nonlin a one should take the maximum value for F, that is, For this reason a significant decrease in ''hardening ability'' of the spring may occur, which leads to some degradation of the absorber effectiveness. However, if we take some ''more neutral'' value, for instance ðF 1 þ F 2 Þ=2, and F 2 ¼ 5F 1 then, in case F 0 ¼ F 2 we have the quadruple increase in value for j,j, and this may lead to much more increase in amplitude of the responses of primary system (see Fig. 15).

Tuning the DVA and numerical validation
Now we check the recommendations regarding the choice of absorber's parameters based on the numerical integration of Eqs. (9) and (11). We consider that the mass of the absorber and the interval of possible values of the frequency of external influence are known: l ¼ 0:14; e 2 ½À 0:15; 0:2. As can be seen from Fig. 8, in this case the surface has two ridges, that is, the amplitude-frequency curve has two peaks. Performing the appropriate calculation, we obtain the following values c % À0:022; h 2 % 0:028; , % À0:54 Â 10 À4 . The corresponding value of function n in this case is % 104. Substituting the values of the parameters and the initial values u 1 ð0Þ ¼ 7; v 1 ð0Þ ¼ À6; u 2 ð0Þ ¼ À0:5; v 2 ð0Þ ¼ 1, we build the phase portrait, and the projections of limit cycles onto the plane Ou 1 v 1 are shown in Fig. 17a, b. These cases correspond to two peaks related to the system (11). The coordinates of the point M 0 obtained by solving the system (14) are ð6:72; 7:46; À 2:56; 6:17Þ for e ¼ Fig. 16 Evolution the shape of the limit cycle for DVA under increasing value of , À0:1 (the left peak) and ðÀ 4:67; 9:10; À 3:12; À 6:68Þ for e ¼ 0:11 (the right peak).
Note that a slight change in the selected parameters does not affect the magnitude of the maximum responses too much, although it may worsen some of the secondary characteristics. For example, in Fig. 17b, the time of the transition period (the time during which the trajectory reaches a small neighborhood of a stationary point) is % 90 s (on s-scale), and in Fig. 17c-% 105 s. Recall also that c is a detuning parameter of the linear component of the stiffness of the absorber relative to the stiffness of the main mass. Therefore, a change in k lin a of 1À2% corresponds to a change of 50À100% for c, and this component is the most vulnerable link in terms of tuning accuracy. Changing the parameter h characterizing the damping, say 10%, does not matter much (Fig. 17c). The most ''vigorous'' role plays the nonlinear component of stiffness. In Fig. 17d-f there are presented limit cycles for three different values of , : weaker the sample one (À 0:3 Â 10 À4 against À 0:54 Â 10 À4 ) and two other ones for values À 10 À4 and À 2 Â 10 À4 . Even the increasing the absolute value of the magnitude of this component at four times leads to rather little increase of the amplitude of oscillations of primary system. However, such a flexible nature of the coefficient , does not apply to a change in its sign. The change of its value from À 0:3 Â 10 À4 to 0:3 Â 10 À4 makes the difference. At the frequency of the former peak e ¼ 0:11), the amplitude of oscillations does not increase (Fig. 17g). However, when frequency ratio changes a little (e ¼ 0:12), the motion becomes unstable (Fig. 17h).

On the choice of the frequency ratio
The term ''frequency ratio'', apparently, was used by Dan Hartog [3] firstly in order to denote the ratio of the frequency of an external excitation to the frequency of the main mass. This designation is natural, understandable and generally accepted. At the same time, although we are talking about the dynamics of the ''main mass ? absorber'' system, this value does not reflect the influence of the latter. For this reason, it seems to us that more natural for ''frequency ratio'' is the expression X 0 ¼ X=ð1 þ l 2 Þ, which represents the ratio of the frequency of an external force to the frequency of the main system with additional mass m a (''frozen'' absorber). Of course, usually the mass of the absorber is small compared to the main mass, and the difference between X and X 0 is quite small. However, we see the following arguments in favor of preference X 0 before X: (1) Referring to the surface view P 5 ðnÞ ¼ 0 which at each point in the region D par determines the value of the square of the amplitude of the responses of the main system. Recall that according to the formulas (3), the value corresponding to x 0 is e ¼ 0, and the value for x 1 is e ¼ l 2 =ð1 þ l 2 Þ. As we see from Fig. 18, the plane e ¼ 0 passes through the ''vale'' of the surface (Fig. 18a, b), in the plane it is very close to the minimum of the function nðeÞ (Fig. 18c). (2) Another point arises from an analytical reason.
According to Sect. 2, with e ¼ 0 the system (11) has a single stationary point for any set of parameters of the absorber. For nonzero value of e, as can be seen from the Fig. 3, there are configurations of the DVA which bring three stationary points into consideration. Thus, the zero value for e logically may be considered as starting point for various analytical constructions (the perturbation theory for instance). (3) Frequently, the authors divide the frequency range of x by using the x 1 as a border and the low frequencies (X\1) and high frequencies (X [ 1) as dividing values. The subsequent analysis is carried out taking into account this partition. For instance, in [14] the authors analyze the effect of a nonlinear absorber on a linear basic system and give the conclusion: ''In the low-frequency range of X\1, a softening stiffness absorber performs well by increasing power absorption ratio and reducing peak. ...In comparison, the hardening stiffness absorber enhances power absorption and attenuates the peak of kinetic energy in the high-frequency range of X [ 1.'' We believe, that such conclusion will be more accurate by replacing X with X 0 . We consider the following examples to validate our statements.
Example 5.1 Let the mass ratio is m a =m 1 ¼ 0:0196 ðl ¼ 0:14Þ, and the frequency range is e 2 ðÀ 0:2; 0Þ ðX 0 \1Þ. In this case a softening spring (, [ 0) is preferable. Since the amplitude curve has only one peak, it is necessary to minimize the largest of the values nðc; h; ,Þj e¼À0:2 , n max ðc; h; ,Þ; nðc; h; ,Þj e¼0 : Taking into account that on oe j e¼À0:2 [ 0, two values need to be compensated, i.e., the peak value and its counterpart at the right border of the interval.
The ranges for parameters c; h; , we determine according to approach described in Sect. 4. Thereafter the optimal values are determined by simple procedure which is illustrated in Fig. 19. Now taking c ¼ À0:094; h ¼ 0:095; , ¼ 0:9 Â 10 À4 ; we plot the frequency-amplitude curve (Fig. 20a). The maximum value is n max % 52:3. From the other side, counting the right border as e ¼ 1 À 1=ð1 þ l 2 Þ % 0:0192 and optimizing the parameters, we get another curve which gives the bigger value for n :% 65:2.
Example 5.2 Another possibility to compare two candidates x and x 0 is the testing of the debatable interval ½x 0 ; x or I e ¼ ½0; 1=ð1 þ l 2 Þ. If one choose x as a borderline, this is a low-frequency range, and a  high-frequency range with the choice of x 0 . At the first case we have to take the softening stiffness (, [ 0) and the hardening stiffness in opposite case. Figure 20 presents the frequency-amplitude curves (Fig. 20b) and the results of integrating the motion equations for both cases (limit cycles in Fig. 20c). Thus, the hardening spring gives about 25% lower amplitude of the responses.

Concluding remarks
We have studied the problem of determining the parameters of NDVA in order to reduce the maximum possible amplitude of oscillations of the main system in the conditions of uncertainties (frequency ratio and amplitude of external excitation). Particular attention is paid to the development of an analytical procedure that allows to easily algorithmize the procedure for setting up the absorber and testing the results obtained. The influence of the nonlinear elastic characteristic of DVA on the efficiency of work is studied, the recommendations are given on the choice of the spring (hardening or softening) depending on the expected frequency ranges and the amplitude of the external influence. The necessary and sufficient conditions for the asymptotic stability of the quasiperiodic motions of the mechanical system under study are found based on the use of a simplified mathematical model (the averaged equations). Technically, the proposed approach may be described by the following scheme shown in Fig. 21.
The validity of using such a model is confirmed by the results of numerical experiments. The approach used in the article can be extended to cases of a nonlinear main system and a system with many degrees of freedom.
Finally, we emphasize importance of our study. Namely, under conditions of uncertain parameters, the direct numerical optimization can be inefficient or at least costly. This is because for this problem there is no explicit expression for the objective function, moreover, it is a function of many ð ! 4Þ variables, some of which are tougher (a small change weakly affects the behavior of the system), while other parameters are very sensitive (linear and nonlinear absorber stiffnesses). In this work, to set up the absorber, an approach based on visualization is used an image in the form of a surface of the dependence of the oscillation amplitude of the main system on the frequency and nonlinear component (nl.c.) of the absorber stiffness. This allows with easy localize the interval of suitable values for nl.c. After that, the remaining parameters are refined. Localization of the interval for nl.c. seems to be an important step, since initially this parameter is completely undefined (in dimensionless parameters, the appropriate values appear of the order of 10 À5 , so iterating over possible values with a step 0.01, for instance, will be ineffective. The mathematical expression for the objective function is rather cumbersome (and the function is specified implicitly), but computer algebra helps to carry out the necessary analysis.