Nonlinear vibration of a lumped system with springs-in-series

The paper deals with the dynamics of a lumped mass mechanical system containing two nonlinear springs connected in series. The external harmonic excitation, linear and nonlinear damping are included into considerations. The mathematical model contains both differential and algebraic equations, so it belongs to the class of dynamical systems governed by the differential–algebraic system of equations (DAEs). An approximate analytical approach is used to solve the initial value problem for the DAEs. We adopt the multiple scales method (MSM) that allows one to obtain the sufficiently correct approximate solutions both far from the resonance and at the resonance conditions. The steady and non-steady resonant vibrations are analyzed by employing the modulation equations of the amplitudes and phases which are yielded by the MSM procedure.


Introduction
The mechanical systems which contain parallel or serially connected massless springs are widely investigated and discussed in the theoretical and applied mechanics. They have found applications in mechanical and civil engineering, mechatronic devices, and more recently in micromechanical systems. Various configurations of the connections between the springs, including also their spatial orientation, can lead to the complex dynamical behavior of those systems, especially when the elastic elements have the nonlinear characteristics. Such systems could exhibit a variety of interesting behaviors, sometimes even surprising which especially concern the resonance states.
Models of many real systems demand to introduce rigid body approximation where some springs and dampers are connected in various configurations. The car suspension containing systems of the parallel and serially connected springs is investigated in [1,2]. The authors showed that such connections have a great impact on the vibration transmissibility from the rough road to the car body.
Telli and Kopmaz [3] studied a one-dimensional oscillator mounted via two springs wherein one of them is linear and the second one has nonlinear features. They proposed two mathematical models for the system considered. The differential-algebraic equations on which the first approach is employed have been solved numerically. The second model based on a single differential equation, obtained using the relative displacement variables allows for achieving the approximate analytical solution. Comparing the results obtained using the two approaches they determined the values of the parameters describing elastic properties of the springs for which the high agreement is observed.
The system similar to the one analyzed by us is the subject of the paper [4], where tethered satellites are modeled by two serially connected springs. Another application of the similar system one can find in the paper [5], where the springs assembled both in the serial and parallel configurations are used by the authors as a model of the structural equivalent stiffness. In that way, they produce a powerful procedure to study structural behavior. A one-dimensional mass-spring system of two degrees of freedom is studied in [6]. Its mathematical model consists of two strongly coupled differential equations. The analytical method of nonlinear normal modes was successively applied to qualitative analysis of the behavior of that system.
Manevitch and Musienko [7] investigated, two systems with one-and two degrees-of-freedom having the serially connected springs. The equations governed the systems contain both differential and algebraic equations. The appropriately adopted multiple scales method (MSM) has been presented and successfully applied by the authors in the case of the undamped free vibration.
In this paper, a wider analysis of the one degree-offreedom system containing two serially connected nonlinear springs is presented. The damped and forced vibration are analyzed using the asymptotic approach. The paper is organized in the following way. In Sect. 2, the mechanical model is described and then the mathematical model is derived. Section 3 deals with the asymptotic solution to the set of the differential-algebraic equations obtained using MSM. Section 4 is devoted to the vibrations far from resonance, while in Sect. 5, the main resonance case is analyzed. In that part, both the steady and nonsteady states are investigated. Nonlinear damping is included into the mathematical model in Sect. 6. Concluding remarks regarding the obtained results are outlined in Sect. 7.
2 Mathematical model of the oscillator A rigid body of mass m due to the constraints imposed moves only in the horizontal direction. The body is connected to the immovable wall by two serially linked springs and the viscous damper with the damping coefficient equal to C. Moreover, the body is excited by the force of the harmonically changing value F 0 sin Xt ð Þ. The physical model of the system is shown in Fig. 1.
The position of the massless point S at which the springs are connected to each other is L 01 þ X 1 t ð Þ, and the body position equals L 01 þ L 02 þ X 2 t ð Þ, where L 01 and L 02 are the nominal lengths of the springs.
The kinetic energy of the system has the form The springs have nonlinear properties of the cubic type, thus the potential energy is as follows where k i and K i stand for the stiffness coefficients of the i-th spring. The external excitation and the damping force are considered as the generalized force Since the springs are connected serially, the following equilibrium equation for the massless point S must be satisfied The equation of motion of the system is derived using Lagrange's formalism and then transformed into a more convenient dimensionless form. Similarly, the equilibrium equation and the initial condition are also formulated in the dimensionless form. The governing equations supplemented with the initial conditions have the form Let us introduce the reference mechanical system with the analogous structure to the considered one but with the springs whose elastic properties are linear. The effective stiffness of the reference system of two spring connected in series is k e ¼ k 1 k 2 k 1 þk 2 : Taking into account various reasons of the geometric, strength, or structural nature one can determine an admissible value of the static elongation of the springs. Having denote it by D, we write D ¼ mF 0 k e ; where m ) 1: Assuming so understood parameter as the reference length, we define the dimensionless coordinates x 1 ¼ The coordinates x 1 and x 2 are functions of the dimensionless time s ¼ tx, where x ¼ ffiffiffiffiffiffiffiffiffiffi k e =m p : The other dimensionless parameters are defined as follows: The mathematical model describing the dynamics of the system under consideration contains one algebraic and one differential equation of the second order. The latter one is supplemented by the initial condition. Some special MSM approach is needed to solve the mathematical problem of such kind.
3 An analytical solution to the problem The differential-algebraic Eqs. (5)-(7) are solved analytically in the asymptotic sense using the MSM.
Let us introduce the small parameter e; where 0\e ( 1. The approximate solution, in asymptotic sense, to the problem (5)-(7) is assumed to be equal the sum of few products of the small parameter in various powers and the unknown functions n 1k and n 2k dependent on several time variables termed as the time scales s k ¼ e k s for k ¼ 0; . . .; n À 1. So, we write e k n 1k s 0 ; . . .; s nÀ1 ð ÞþO e n ð Þ; ð8Þ e k n 2k s 0 ; . . .; s nÀ1 ð ÞþO e n ð Þ: ð9Þ The differential operators take the form Limiting our considerations to the weakly nonlinear type of the system, we assume the smallness of some parameters. The smallness can be strictly expressed using the small parameter e. Introducing the notation with tilde we assume: where r i e 1; 2 f g andã 1 ;ã 2 ;c;f 0 are of the order O 1 ð Þ: To solve problem (5)-(7) two variants of MSM with two and three time scales have been applied and compared. The exponents in definitions (11) in the case of two time scales equals:

Vibration far from resonance
Relations (8)-(11) substituted into Eqs. (5)-(6) yield the equations in which the small parameter e appears in a few different powers. We order the left sides of the equations according to the powers of the small parameter. We require each of the both equations to be satisfied for any value of the small parameter e. So, omitting all terms of the order O e n ð Þ and higher ones, we equate to zero the coefficients standing at the particular powers of e. This procedure leads to the equations of the subsequent order approximation.
In the case of two time scales, the substitution of definitions (8)-(11) into Eqs. (5)- (6) and equating components at the same powers of e leads to the equations of the first and second-order approximation of the following form: kn 21 À 1 þ k ð Þn 11 ¼ kã 2 n 10 À n 20 ð Þ 3 þã 1 n 3 10 : In the case of three time scales, the first, second and third order approximation are obtained similarly as above. They are as follows: (i) equations of the approximation of the order e 0 o 2 n 20 (ii) equations of the approximation of the e 1 kn 21 À 1 þ k ð Þn 11 ¼ kã 2 n 10 À n 20 ð Þ 3 þã 1 n 3 10 ; ð19Þ (iii) equations of the approximation of the e 2 o 2 n 22 The above DAEs (12)-(15) or (16)-(21) is solved recursively, i.e. the solutions to the lower order approximation equations are introduced into the ones of the higher order.
The Eqs. (12)-(13) and (16)-(17) of the lowest order (i.e. e 0 ) have the same form when using two or three time scales. The general solution to these equations are where B and its complex conjugate B are the unknown complex-valued functions of the slower time scales. Inserting the solution of the lower order approximation to the equation of the higher order approximation implies appearing so-called secular terms, that should be eliminated due to the requirement to obtain the limited solutions. The procedure is described in detail in the paper [7]. Elimination of the secular terms leads to the solvability conditions, that are differential equations of the first order with unknown functions B and B.
It is convenient to express the complex functions B and B in the exponential form where a and w are unknown real-valued functions and stand for the amplitude and the phase of the first term of the asymptotic solutions (9). Introducing relationships (24) into the solvability conditions (not given here explicitly), then returning to the original notations according to (11) and using the definition of the differential operator (10) 1 allow us to write the modulation equations.
When two time scales are adopted, the modulation equations have the form of the following differential equations of the first order Assuming the initial conditions in the form we obtain the solution to the problem (25)-(27) as follows The motion of the body can be now described by the approximate function of the analytical form. The expressions for the amplitude and phase (28) are substituted into relations (24), and then into solutions of the first and second order. Finally, from the expansions (8)- (9), and using the original notations according to (11) we get, the dimensionless displacements of the body as follows: where a s ð Þ and w s In the case of the solution basing on the approach with three time scales, the modulation equations have the form: Equations (30)-(31), together with initial conditions (27), form the initial value problem, that cannot be solved analytically. So, the modulation problem must be solved numerically in this case. Now, the approximate analytical solution to the governing Eqs. (5)-(6) for the body displacement x 2 has the following form where a s ð Þ and w s ð Þ are solutions to the modulation problem (30)-(31).
The time histories of the body displacement obtained using the asymptotic solution (29) and (32), respectively for two and three time scales, are presented in Fig. 2. They show the body vibrations for the transient vibration stage. The solution obtained by the numerical integration of the governing Eqs. (5)-(6) are given in the same graphs to verify the accuracy of the approximate analytical solution. The relationship between initial conditions (7) and (27) has been determined using the analytical form of the solution (29) or (32) for the two and three time scales, respectively.
The results presented in Fig. 2 are obtained for the following fixed values of the parameters: a 1 ¼ 0:87; a 2 ¼ 1:21; k ¼ 0:7; f 0 ¼ 0:08; p ¼ 0:215; c ¼ 0:07; x 0 ¼ 0:1; and v 0 ¼ 0. The compatibility of the two curves is very high, which confirms the correctness of the derived analytical solutions in both cases i.e. with two and three time scales.
Moreover, we propose some objective measure of the analytical solution correctness. It is checked the fulfillment of the governing equation by the asymptotic solution. For this purpose, we rewrite the governing Eqs. (5)-(6) in the compact form where The error of the satisfying of Eqs. (5)- (6) is defined as follows where x 1a s ð Þ; x 2a s ð Þ are the approximate solutions (29) or (32), depending on the context, obtained employing the MSM, and s s , s e denote the chosen instants.
Calculating the error of the analytical solutions according to definition (34) for the results presented in Fig. 2, we assume: s s ¼ 0 and s e ¼ 1000. In the case of solution (32), obtained using three time scales, d 1 ¼ 0:0009073, d 2 ¼ 0:00004152, whereas for the more rough approximated solution (29), using only two time scales, d 1 ¼ 0:0009076, d 2 ¼ 0:00004151; which can be considered as very good results. It turns out that two time scales are enough to obtain high accuracy for non-resonant vibration. It means that the more sophisticated analysis with three time scales insignificantly improves the results.

Vibration at resonance
The main resonance appears in the system when the frequency of the external force is close to the eigenfrequency of the analogous linearized system, i.e. p % 1. To investigate this case let us introduce the detuning parameter r in the following way After substituting relationship (35) into Eq. (5), the procedure analogous to the one discussed in Sect. 4 is carried out. Two variants of MSM with two and three time scales have been applied. The modulation equation obtained using two time scales are The approximate analytical form of the function describing the body displacement using two time scales is as follows whereas, the calculations performed with three time scales give modulation equations in the following form and the approximate analytical form of the function describing the body displacement for the resonant motion is as follows where a s ð Þ and w s ð Þ occurring in (38) and (41) denote the solutions to the modulation Eqs. (36)-(37) or (39)-(40), respectively. Unfortunately, the modulation equations cannot be solved analytically, therefore a hybrid method must be used to solve the original problem i.e. the numerical solution of the modulation equations should be introduced into the analytical description of the body displacement (38) or (41).
The absolute body displacements are presented in Figs. 3 and 4. The graphs are obtained for the following fixed values of the parameters: a 1 ¼ 1:025; a 2 ¼ 1:1; f 0 ¼ 0:00205; k ¼ 0:9; r ¼ 0:03685; c ¼ 0:001; x 0 ¼ 0:1; v 0 ¼ 0: In Fig. 3, the transient vibrations are presented. In the graph on the left side, the approximate solution obtained using the variant of MSM with two time scales is compared with the numerical solution obtained using standard NDSolve procedure of  Similarly as in the non-resonant case, the asymptotic solutions regarding to the main resonance are highly compatible with the numerical solution, irrespectively of the applied variant of MSM. The error calculated according to the definition (34) is now d 1 ¼ 0:000015; d 2 ¼ 7:145389 Â 10 À8 for three time scales and d 1 ¼ 0:0000181; d 2 ¼ 2:8904498 Â 10 À6 for two time scales, assuming s s ¼ 0 and s e ¼ 800, which validates the very high accuracy of the used analytical approach. Analyzing the error values one can state that the solution that utilizing three time scales is visibly more accurate than the one with two time scales in the case of the main resonance. Thus, we decide to base the detailed study of the main resonance on the results obtained using the variant of MSM with three time scales.
The results presented above confirm the usefulness of the employed approximate analytical method.

Steady-state resonant responses
In further analysis, it is convenient to introduce into modulation Eqs. (39)-(40) the modified phase h s ð Þ ¼ rs À w s ð Þ. That allows to transform the equations into the autonomous form more suitable for analysis of the stationary and nonstationary motion.
The amplitude and modified phase that are constant in the steady-state motion correspond to the fixed point of the modulation equations. To avoid introducing additional symbols we denote them as a and h (without the argument s). According to the assumption about steady state, zeroing the time derivatives in modulation Eqs. (39)-(40) yields the following set of the algebraic equations that correspond to the steady-state conditions: Using trigonometric identity one can obtain the amplitude-frequency relationship in the form The resonance curves presenting the amplitude and the modified phase, obtained as solutions of Eqs. (42)-(43) are shown in Figs. 5, 6 for the following fixed parameters: a 1 ¼ 1:025; a 2 ¼ 1:1; k ¼ 0:9; f 0 ¼ 0:00205; c ¼ 0:001. The same parameters have been adopted to obtain time histories presented in Fig. 4. Studying the amplitude resonance response depicted in Fig. 5, one can state that for r = 0.03685 there are two stable values of the amplitude, namely 0.623505 and 0.0278721. The same values are drawn in Fig. 4 in red line as the steady-state amplitudes.
In Figs. 5, 6 the stable branches of the resonace response curves are depicted in red color, whereas the unstable ones in blue color.

Stability of the resonance curves
To examine the stability of the steady-state solution in the sense of Lyapunov, we analyze the non-stationary solutions of Eqs.
for which the components of the characteristic matrix are as follows: The fixed point ða s ; h s Þ relating to the steady-state solution is asymptotically stable in the sense of Lyapunov if the real parts of all eigenvalues of the matrix A are negative.
The Eqs. (42)-(43), which determine the resonance response curves, predict the system behavior that depends on a few parameters. The influence of the external excitation amplitude f 0 on the shape of the response curves is presented in Figs. 7, 8 for the following fixed parameters: a 1 ¼ 0:05; a 2 ¼ 0:05; k ¼ 0:9; c ¼ 0:001.   The influence of the other parameters k, a 1 , a 2 or c, can be also easily investigated using Eqs. (42)-(43).

Non-stationary vibration
Usually, when investigating an anharmonic oscillator, high emphasis is placed on steady-state vibrations as the most simple and practically important process. However, in many mechanical systems relaxation to steady-state continues for a long time. The reasons can be a strong resonance and small damping. Therefore analysis of the non-stationary motion could be also important. Usually, this problem is studied numerically due to the occurred essential mathematical difficulties.
The MSM with two time scales gives a significantly simpler form of the solution than using three scales. However, the results obtained in this way are still very good to describe the behavior of the system. Therefore, further analysis is performed employing a 2-scale variant of MSM.
A good way to qualitatively examine the system is to analyze the non-damped vibration, that can be studied analytically. When ignoring the energy dissipation (c ¼ 0), the modulation Eqs. (36)-(37) take the form The exact differential equation corresponding to the set (49)-(50) of the form has the first integral The Eq. (52) governs phase trajectories on the plane (a; h). The constant on the right side depends on the initial conditions of the system (49)-(50).
Let us analyze the case when the maximum of the energy transfer between the system and the external loading appears. This situation takes place for const ¼ 0 at Eq. (52), for which the first integral has the form where aH 0 ¼ H 0 : The trajectory in the plane (a; h) corresponding to the maximum energy transfer H 0 ¼ 0 is called the Limiting Phase Trajectory (LPT) [8]. From the total differential dH 0 ¼ 0; we obtain and hence the maximum amplitudes are achieved at h ¼ p=2 þ kp, for k 2 {. Let us consider h ¼ p=2. The number of the real roots of Eq. (53) for r [ 0 and f 0 [ 0 depends on the sign of the expression The discriminant D depends on five parameters. A qualitative change in the system behavior occurs when sign of D changes itself. The critical value of any of the five parameters causes D ¼ 0. For example, the critical value a 1cr of the nonlinearity parameter a 1 while other parameters are fixed is The important qualitative transition in the dynamical behavior of the system appears when the value of the parameter a 1 crosses its critical value a 1cr . Three limiting phase trajectories, for a 1 ¼ a 1cr , a 1 \a 1cr and a 1 [ a 1cr , are presented in Fig. 9. The graph is drawn for the data: f 0 ¼ 0:001; k ¼ 0:7; r ¼ 0:004; a 2 ¼ 0:1, for which a 1cr % 0:4328.
The comparison presented in Fig. 9 shows the qualitative change in the trajectory corresponding to the maximum energy exchange, i.e. for H ¼ 0. For a 1 ¼ 0:30\a 1cr there are two trajectories. One of them encircles quasilinear center at h ¼ Àp=2 or at h ¼ 3p=2, and the maximum value of the amplitude reaches 0.3 (point ''4'' in Fig. 9). The second one is open and does not describe vibrations. At a 1 ¼ a 1cr ¼ 0:4328 the qualitative transition appears and the trajectory changes its shape and location. Then, for a 1 ¼ 0:56 [ a 1cr the trajectory encircles strongly nonlinear center at h ¼ p=2, and the maximum amplitude of the vibration equals now 0.69 (point ''2'' in Fig. 9). The above results indicate strong modulation of the amplitude in the case of the resonance. One can observe the long term history of the non-steady resonance vibration of the undamped system. The result obtained according to the analytical solution (38) for a 1 ¼ 0:43\a 1cr is presented in Fig. 10 for the following fixed parameters: f 0 ¼ 0:001; k ¼ 0:7; r ¼ 0:004; a 2 ¼ 0:1.
Similarly, the time history of x 2 s ð Þ for the slightly greater nonlinearity a 1 ¼ 0:48 [ a 1cr is presented in Fig. 11. The values of other parameters are the same as above.
The time history of the displacement of the body obtained numerically for the original problem (5)- (7) is qualitatively and quantitatively very similar to the one obtained analytically which is seen in Figs. 10, 11. That validates the correctness of the analytical approach. The results presented in these figures show that the long term amplitude modulation is qualitatively different in the case of the quasilinear vibration for a 1 \a 1cr than in the case of strongly nonlinear vibration i.e. for a 1 [ a 1cr .

Non-linear damping
In the analysis carried out in the previous sections, the elasticity was assumed as the only source of the nonlinearities in the system. However, among many possible reasons of the nonlinearity in structural dynamics, also the energy dissipation mechanism can cause nonlinear behaviors that often play a dominant role in many real-world engineering systems. The problem of the nonlinear damping and its modeling is described among others by Elliot et al. [9].
Let us analyze the system presented in Fig. 12. Roughly, this is the same system that was the object of our considerations in the previous sections but now the nonlinearity is caused by the elasticity and the damping.
The system description and all denotations introduced in Sects. 2, 3, 4, 5 are still valid. It is assumed, however, that the damping force according to the Rayleigh model depends on the 3 rd power of the velocity. So, the external excitation and the damping force considered as the generalized force are as follows The following initial value problem written in the dimensionless form is obtained using the Lagrange formalism where c ¼ C mx ,b ¼ BD 2 mx , and the other parameters are defined the same as in Sect. 2.
The non-linear damping model adopted in the article may lead to self-excited vibrations. Such phenomena were investigated by Rayleigh [10] for vibration maintained by wind, heat, or friction.
Governing Eqs. (58)-(59) are solved using the procedure presented in Sects. 3 and 4. Let us assume the smallness of the following parameters: The variant of MSM with three time scales is applied in solving the problem governed by Eqs. (58)-(60). The amplitude and frequency modulations are given by the following modulation equations: The asymptotic analysis leads to the following form of the solution to the initial-value problem given by Eqs. (58)-(60): Fig. 11 Time history of the function x 2 s ð Þ for a 1 ¼ 0:48 (a 1 above a 1cr ) Fig. 12 The mechanical system with nonlinear elasticity and damping that could indicate the bifurcation point, which is not typical for the Duffing equation. The unexpected behavior in the resonance may indicate a limited scope of the application of the adopted non-linear damping model, which demands experimental validation.

Conclusions
The lumped system containing two serially connected nonlinear springs has been investigated. The governing equations contain both differential and algebraic equations. The dynamical problem of the forced vibration of that system has been solved in an asymptotic manner using the MSM, which had to be appropriately modified to deal with the set of the differential-algebraic system. The initial value problem has been solved by two variants of the MSM, i.e. using two and three time scales. It turns out that for the problem investigated for two time scales were enough to obtain a very good solution, however more time scales could be needed to deal with more sophisticated systems in which geometric nonlinearities appear. The analysis of vibration far from the main resonance as well as near the resonance has been presented. In the case of the resonance, both the steady-state and the non-stationary vibrations were analyzed.
The accurate solution has been obtained also for the system with smooth nonlinear damping according to the Rayleigh model.
Satisfaction of the governing equations by the approximate analytical solution has been verified. Small values of the root-mean-square error that have been recorded prove the high accuracy of the method. The correctness of the asymptotic approach has been also confirmed in the confrontation with the numerical simulation results.
It is well known that the analytical form of the solution has a great advantage over the numerical one because it allows for the qualitative analysis of the system behavior.
Finally, we would like to emphasize the obtained surprising novel results with regard to vibration with nonlinear Rayleigh type damping at resonance. Namely, we have detected that there is interval for the parameter r for which there is no any stable periodic solution (the direct numerical tests indicate that the solution tends to infinity, i.e. it is unbounded). This result is obtained for the first time to the best of our knowledge.
In addition, we have found the triple bifurcation point shown in Fig. 15 which is not typical to the known phenomena exhibited by the Duffing type oscillators.
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 need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.