Simulating Nonlinear Oscillations of Viscoelastically Damped Mechanical Systems

The aim of this work is to propose a mathematical model in terms of an exact analytical solution that may be used in numerical simulation and prediction of oscillatory dynamics of a one-dimensional viscoelastic system experiencing large deformations response. The model is represented with the use of a mechanical oscillator consisting of an inertial body attached to a nonlinear viscoelastic spring. As a result, a second-order firstdegree Painlevé equation has been obtained as a law, governing the nonlinear oscillatory dynamics of the viscoelastic system. Analytical resolution of the evolution equation predicts the existence of three solutions and hence three damping modes of free vibration well known in dynamics of viscoelastically damped oscillating systems. Following the specific values of damping strength, over-damped, critically-damped and under-damped solutions have been obtained. It is observed that the rate of decay is not only governed by the damping degree but, also by the magnitude of the stiffness nonlinearity controlling parameter. Computational simulations demonstrated that numerical solutions match analytical results very well. It is found that the developed mathematical model includes a nonlinear extension of the classical damped linear harmonic oscillator and incorporates the Lambert nonlinear oscillatory equation with well-known solutions as special case. Finally, the three damped responses of the current mathematical model devoted for representing mechanical systems undergoing large deformations and viscoelastic behavior are found to be asymptotically stable. Keywords-mathematical modeling; nonlinear oscillations; viscoelastic oscillator; Painlevé equation; exact solution; numerical simulation


INTRODUCTION
In this work the dynamics of mechanical systems undergoing large deformations and viscoelastic response is investigated.A major topic in the dynamics of viscoelastic systems is the problem of vibration.Vibration phenomenon arises in all rigid or deformable systems, such as machines and engineering structures subjected to dynamic loading.So, the vibration problem is of vital importance for many fields of science and technology.Vibration experiments are widely used in the characterization of dynamical mechanical properties of engineering materials.Vibration is also desired for machines under working conditions.However, for most structures in mechanical, biomechanical, civil, aeronautical and automotive engineering, oscillatory events prediction and control is intensively required in order to reduce noise, and to prevent non-allowable or excessive deformations, self-excited deformations, material fatigue and failure [1].
The one-dimensional dynamics of continuous viscoelastic media idealized as a bar is described formally by the Cauchy's wave equation.However, when only homogenous deformation is considered in longitudinal forced vibration experiment for a bar carrying a tip mass at the free end, under the condition that the mass of the bar may be neglected compared to that of the attached mass at the end point, the bar may be assumed to behave as a simple viscoelastic spring.So, the Cauchy's equation may be reduced, based on the second Newton's law, if u denotes the time history of displacement response, then in the concise form [2,3]: on the time t, m is the mass of the attached body, F ext designates an applied exciting force, and int F notes the internal axial force due to the viscoelastic stress induced in the considered mechanical system.In this context the question of viscoelastic systems vibration transforms into that of finding an appropriate constitutive equation for structural materials.
Viscoelastic behavior is widely mathematically analyzed through the use of rheological approach which discretizes a viscoelastic body in its elementary elastic and viscous components represented in terms of mechanical analogs.So, the rheodynamical approach models the motion of a viscoelastic system in terms of ordinary differential equations [3,4].According to Laroze [5], internal damping can be schematized for most structural materials in practice as viscous damping.In this perspective, if it is assumed that the viscoelastic spring behaves as a Kelvin-Voigt medium, that is, if: where φ(u)=u, then the basic evolution equation governing the forced viscoelastically damped linear vibrating mechanical system may be written in the form: or in terms of only ) (t u under free vibration: The quantity ω ο is defined as the angular or undamped natural frequency of the physical system, and λ is a damping factor [5,6].Equation ( 4) is a classical prototype of secondorder ordinary differential equations used to describe the damped linear oscillation of a single degree of freedom oscillator.It represents the dynamics of the motion of a onedimensional viscoelastic system in terms of Kelvin-Voigt element with the addition of a mass in the range of linear deformation.Equation ( 4) is the simplest differential equation that may reproduce all different types of behavior exhibited by a damped second order oscillation system in the range of linear deformation.The type of response given by (4) representing a damping linear oscillation system depends, in effect, on the strength of damping degree.In other words, the response of a damped oscillatory system is very sensitive to changes in specific value of damping.
It has been observed that the oscillatory response of a real system is nonlinearly damped and geometrically nonlinear, leading to the dependence of the stiffness on the induced deformation or stress.This phenomenon is known as stiffening or softening and terminates often by failure, following the stiffness increases or decreases.Such a phenomenon could not be predicted and explained by any linear model equation [6][7][8][9][10][11].Many systems in engineering applications are designed to behave not only viscoelastically but also nonlinearly, to say, to undergo large deformations exceeding the limiting value predicted by linear theory in loading environment [2,6,7,11].Viscoelastic systems have the specific ability to experience large deformations even for moderate force levels.Hence, the linear theory becomes quite unable to explain the dynamics of viscoelastic systems experiencing finite deformations [6,7,[12][13][14][15] .Geometrical, damping and material nonlinearities are fundamental causes for the viscoelastic behavior of systems to be nonlinear [5].If the stiffness or geometrical nonlinearities are relatively well captured in terms of polynomial expansion, the damping nonlinearity properties of mechanical systems are very difficult to be known [4,7,8].According to [7,8], another non-negligible source contributing to nonlinear response of a mechanical system is the inertia properties.Inertia produces, to a certain extent, a force which is essentially nonlinear.Nonlinear inertia forces are generally proportional to higher powers of velocity and acceleration of the mechanical system [7,8].The above show that an oscillatory viscoelastic system is intrinsically characterized at least by its stiffness, damping and inertia nonlinearities.Therefore, a reliable representation of the nonlinear oscillatory dynamics of a viscoelastic system should have the ability to mathematically incorporate in the governing equation these basic nonlinearity principles [7].From the above analysis, the most general second order ordinary differential equation that can model the nonlinear oscillatory dynamics of a single degree of freedom viscoelastic system under unforced conditions may take the form: where the dot over a symbol designates the time  5) gives (4) modeling a damped linear harmonic oscillator.It is interesting to note that in the majority of existing evolution models of mechanical systems, only one of these nonlinearities is often considered (see for example [16] for more details).This shows that, due to the mathematical complication arising quickly in governing equations, the enhancement, for example, of stiffness nonlinearity is often performed in models to the detriment of system damping nonlinearity and, inversely, the improvement of damping nonlinearity is made in a prejudicial fashion to that of stiffness nonlinearity [1].Moreover, very few of proposed mathematical models for studying the dynamics of viscoelastic systems with a single degree of freedom have been performed to enclose on the one hand the inertia nonlinearity, and simultaneously combined inertia, damping and stiffness nonlinearities governing the viscoelastic response as indicated by (5) on the other hand [7].
It is well known again that nonlinear problems having explicit exact solutions in terms of elementary standard functions are very limited in physical and engineering fields.So, a large part of nonlinear analysis has only been performed on the basis of qualitative theory or particular solutions derived from analytical approximate or numerical integration methods.In particular, homotopy perturbation analysis is intensively used to investigate nonlinear vibration problems in mechanical structures [17].In mechanical system design calculations, for example, the accurate determination of dynamical characteristics from observations requires mathematical models having appropriate analytical solutions.In this context, the design of mathematical models capable of representing the nonlinear dynamics of viscoelastic systems satisfactorily in terms of analytical exact solutions, capturing also simultaneously the combined nonlinear phenomena, becomes a major necessity.In that, the Bauer's rheological-dynamical theory [7] consists of a notable progress in the field of the dynamics of continuous viscoelastic media, since it meets the essential of preceding criteria about the necessity to handle simultaneously and in combined fashion, in mathematical modeling of viscoelastic systems, the nonlinearity properties.
Although the Bauer's method [7] seems to be simple in formulation, it is in reality a powerful approach for solving dynamical nonlinear problems arising in viscoelastic continuum mechanics successfully.This theory [7] entails a significant mathematical modification and extension of the classical Kelvin-Voigt viscoelastic solid law for capturing simultaneously the combined inertia, damping and stiffness nonlinearities characterizing the dynamics of real viscoelastic systems.The Bauer's approach [7] was originally developed for modeling the dynamics of nonlinear viscoelastic response of soft biological tissues, arterial walls in particular.The basic idea underlying this theory consists in developing the total stress within a material system as the sum of three basic stresses: elastic, viscous and inertial stresses, operating in parallel.Very recently, the Bauer's rheological-dynamical theory was formulated in a simple mathematical expression that may be described by a single second order evolution equation within the framework of continuum mechanics for investigating the dynamics of viscoelastic material systems [12,13].This formulation has successfully been applied in several papers to model creep deformation [13], creep relaxation [6] and deformation restoration process under stress relaxation conditions [14,15] of a variety of viscoelastic solid bodies.
The objective of this research work is to develop a mathematical model expressed in terms of an exact analytical solution that may be useful in numerical simulations of the oscillatory dynamics of a one-dimensional viscoelastic system exhibiting large deformations.The physical properties of the mathematical model are represented by a single degree of freedom oscillator consisting of a mass attached to a viscoelastic spring that behaves as a nonlinear Kelvin-Voigt continuum medium.More precisely, given a stiffness nonlinearity function, it is proposed to develop the evolution equation of the time dependent displacement response of a mechanical system experiencing viscoelastic response and large deformations by applying the Bauer's theory as formulated mathematically by Monsia [12,13].In this sense, exact analytical solutions in terms of elementary standard functions have been determined, following various types of response able to exhibit the evolution equation under unforced regime.Numerical applications are carried out for illustrating the ability of the mathematical model to be used for numerical simulations.Graphs of numerical and exact analytical results are compared to prove the validity of the current model.
It is also shown that the developed model incorporates as special cases some oscillatory equations with well known solutions by passing to special limiting values of the model parameters.Time-asymptotic system response is investigated to check the stability character of the long time behavior of the proposed model.So, in the next section, the evolution equation governing the nonlinear oscillatory dynamics of the investigated viscoelastic system will be developed.Section 3 is devoted to solve the obtained second-order first-degree Painlevé evolution equation in closed form exact solution following various types of damping modes of oscillation.Numerical applications of the model are presented in Section 4, and a discussion of the results is performed in Section 5.The final section gives a brief conclusion of the work.

A. Continuum mechanical nonlinear evolution equation
In this part, the Monsia's formulation [12,13] of the Bauer's theory [7] will be briefly reviewed to make the transition from its continuum version to a constitutive expression in terms of displacement and force.The advanced procedure for modeling the dynamics of real mechanical systems, to say, viscoelastic systems, developed by Bauer [7] was expressed within the framework of continuous mechanics, in other words, in terms of stresses and strains experienced by the mechanical system.The Monsia's formulation [12,13] of the Bauer's theory [7] enables us to describe this theory, for a given stiffness nonlinearity function f , by a single second order evolution equation, with only few system parameters to be determined through fitting procedure, of the form:

t f e e f e e lf e e w f e s
where ( ) t s is the scalar stress function and ( ) t e the scalar strain function , the prime denotes a differentiation with respect to response variable, that is, here, the strain ( ) t e , and 0 c ¹ , is the inertia coefficient.The constants λ and o  continue to have the preceding definitions of damping factor and natural frequency, respectively.For some convenience, ( 6) may be written, by using a differential operator B as: At the present time, the function  that captures the system stiffness nonlinearity is not explicitly known, but it must be specified before any further use of the theory.In the following section, an explicit expression for  will be given as function of the deformation response.

B. Nonlinear constitutive force-displacement equation
In this section the constitutive equation in terms of force and displacement governing the dynamics of the onedimensional viscoelastic system under investigation is formulated.For a time history of displacement response ) (t u and an external exciting function , the operator B allows the equation of motion to be represented by: or by the equation written in the developed form: where the prime denotes the differentiation with respect to displacement variable u(t).The constants m, λ and ω ο are time independent and continue to have the preceding definitions.Equation ( 10) is a second-order first-degree Painlevé equation with a forcing function [18][19][20] known to be subject of many studies in mathematics.According to Roth [21], Painlevé equations have found several applications in physics, particularly in the field of circuit oscillations.So, the perspective to model the mechanical behavior of viscoelastic systems using the second-order first-degree Painlevé equation may provide a great insight on their dynamics.Equation (10) gives the differential relationship between the exciting function ) , , , ( and the resulting time displacement response u(t) for a given function φ(u) capturing the purely nonlinear elastic response of the system under consideration.This equation enters in the perspective of the general second order evolution equation of a mechanical system represented by ( 5) and consists of a general formulation requiring the specification of the nonlinear elastic spring force function φ(u) for any subsequent development.Mathematically, all nonlinear function φ(u) that tends towards u, when the displacement u tends to zero, may be selected as a possible restoring force function.In other words, at small deformations, the nonlinear model under question should reduce to the associated linear law.Therefore, an infinite number of functions φ may be designed.In doing so, various types of mathematical expressions for the function φ have been proposed in recent papers [6,[12][13][14][15]22].In general, as performed by Bauer [7], the stiffness nonlinearity function φ(u) may be expanded in a power series as: ... ... ) ( In this perception, for mathematical reasons of simplicity, the nonlinear stiffness function φ will be expressed, in this work, basically as a power law in the form: [22] ( ) where the hardening exponent, that is the stiffness nonlinearity controlling parameter 0 l ¹ .Following this point of view, the equation of motion takes the form: with the mass 0 m ¹ .Equation ( 12

III. EXACT ANALYTICAL SOLUTION FOR THE NONLINEAR OSCILLATIONS EQUATION
This part determines for the unforced oscillatory regime, that is to say, 0 ) , , , ( , the exact analytical solution for the equation of motion following the under-damped, critical damped, and over-damped responses of the mechanical system under question.In doing so, (12) becomes for unforced motion: As can be seen easily, the preceding Painlevé equation has 0 ) (  t u as trivial solution.In the usual standard form, this equation becomes: This second order autonomous differential equation can be viewed as a nonlinearly damped equation.The damping force is a nonlinear function in both the time history of response ( ) u t and the velocity ( ) u t  .The dependence on the variable response ( ) u t introduces a nonlinear Newtonian singularity at the origin as noted in central forces.The dependence on the velocity includes a combination of linear and quadratic damping terms.In this perspective (13) consists of an attractive problem to be investigated to learn more about the qualitative oscillatory dynamics of the system around the singular point origin from a viewpoint of phase plane analysis which is not the topic of the following research work.So, only a quantitative analysis of ( 13) will be performed concerned in this study.The graph in Figure 1 illustrates the oscillatory behavior of ( 13) subject to arbitrary initial conditions u(t=0)=5, and , from t=0 to t=20, obtained by numerical integration.The simulation is run by using Matlab's routine ode45 which exploits Runge-Kutta methods.The parameter values are fixed at A typical behavior of ( 13) illustrating the oscillatory nature of the mathematical model.
In the following section, the ability of the developed theoretical model to mathematically exhibit various responses of the damped nonlinear oscillatory dynamics of the system under investigation is proved.To do so, (13) will be solved in closed form exact solution by using suitable initial conditions that satisfy the dynamics of the nonlinear viscoelastic system of interest.

A. Reduction of the nonlinear oscillatory equation to the damped linear harmonic equation
It is well known in mathematics that the range of nonlinear evolution equations that can be analytically solved in the form of elementary standard functions is very limited [23].So, the objective to find exact solutions of nonlinear evolution equation leads often to the question about the explicit integration of this equation [23][24][25].In this perspective, the Painlevé analysis may be performed in order to conclude on the integration of this equation [23][24][25].However, following Kudryashov [25], the determination of exact solutions of nonlinear evolution equation can be investigated without having to apply successfully the Painlevé test.Here, the governing equation ( 13) under question is a typical equation of Painlevé which has been, according to Keckic [18,19], integrated by quadratures .Hence, the problem of integrability of ( 13) in terms of elementary standard functions is unnecessary to be considered again in this paper.Some equations similar to (13) in terms of elementary functions have recently been solved on the basis of mathematical transformations leading to a Riccati equation or a damped linear harmonic equation [12][13][14][15]22].In the present work, (13) will be turned in the damped linear harmonic oscillatory equation.Hence, to solve (13) in the form of standard functions, a change of variable is required for the expected oscillatory solutions.Performing the following suitable substitution: Equation ( 14) is the well known classical linear second order evolution equation previously noted as (4).Therefore, the general solution of (13) becomes: [19] 1/ ( ) ( ) l u t y t = where ) (t y designates the general solution of ( 14) A mathematical analysis of (14) indicates that the nature of (15) solution depends on the relative magnitudes of the damping factor λ and the stiffness module, to say, the natural frequency ω ο .In other words, the nature of solution depends on the roots of characteristic equation:  14), that can be real or complex numbers.So, three damping modes of vibration should be distinguished.

B. Exact analytical solution for over-damped nonlinear response
This case corresponds to a strong damping [6,26], that is, a relative large damping where , when 0 t = satisfy the past history of the displacement, and taking into consideration (15), the time history of the displacement can be written in the explicit analytical expression:

C. Exact analytical solution for critically-damped nonlinear response
In this case 2 , o l w = the two solutions of the characteristic equation coincide and then, the solution of ( 14) is given by:

IV. NUMERICAL APPLICATIONS OF THE MODEL
The aim of this part is to develop numerical solutions in order to demonstrate the ability of the mathematical model to be used in numerical simulations.In doing so, analytical solutions are validated against numerical results.Numerical integration solutions have been obtained by using the ordinary differential equation solver ode45 and ode15s available in Matlab Package.Matlab routine ode45 uses the fourth-order Runge-Kutta explicit procedure.The integration function uses a variable integration step size to provide the needed accuracy as minimizing computation time.The Matlab ordinary differential equation solver ode15s is based on an implicit method.It is a multistep solver that needs the solutions at numerous previous time points in order to estimate the current solution.In this perspective, (13) should be written in state-space representation, to say, represented in terms of a set of firstorder differential equations.Therefore, defining the state variables: the following system of two first-order differential equations may be obtained: In this form, numerical integration of ( 13) becomes more appropriate to be implemented as m-files in Matlab.Figures below show the results of numerical simulation of (13).Numerical solution is plotted on the same graph as the analytical solution.Whereas the exact analytical solution is plotted in solid line, the numerical result is graphed in circles.The mean squared error is calculated with Matlab mse function to quantify the discrepancy between model predictions and numerical results.

A. Numerical results for over-damped nonlinear response
The graph of the numerical solution for over-damped nonlinear regime is compared to that of the exact analytical solution ( 16) in Figure 2

www.etasr.com Monsia and Kpomahou: Simulating Nonlinear Oscillations of Viscoelastically Damped Mechanical Systems
C. Numerical results for under-damped nonlinear oscillations Figure 4 compares the numerical solution to the exact analytical solution ( 18) of ( 13

A. Analysis of model predictions
Vibrating material systems in real operating situation experience, as it is well mentioned previously, large deformations and viscoelastic behavior.So, their dynamics is characterized by a damped nonlinear oscillation.In this context, a reliable theory devoted for modeling the damped nonlinear oscillatory dynamics of these systems should, at least, have the ability to handle some fundamental nonlinearity problems, and predict all of three damping modes of oscillation known for a real damped oscillatory system.Consequently, the following subsections investigate the aptitude of the proposed mathematical model to satisfactorily reflect these three damped oscillation responses.

1)
Over-damped Response Analysis The graphs in Figure 2 illustrate the over-damped behavior of the proposed model.It can be seen, as expected, that after attaining its maximum value, the time history of the displacement u(t) declines gradually without oscillations to approach asymptotically equilibrium zero value with time.The time history of response u(t) of the system presents only, as expected again, one maximum value.Analytical result, to say, solution (16) shows that the system response has a hyperbolic behavior modulated by decay exponential, similarly to the over-damped linear harmonic oscillator response, but raised to the l / 1 th power.Equation (16) shows that the rate of decay is not only proportional to the damping factor λ, but depends also on the stiffness nonlinearity rising parameter l .So, it may be possible to use this parameter to control the oscillatory response amplitude of the system concurrently with the damping factor λ.

2)
Critically-damped Response Analysis The curves in Figure 3 show that the time dependent displacement ) (t u increases during the first time period to reach a single maximum, followed by an exponential decaying to asymptotically approach equilibrium zero value with time.
Here again, the rate of decay depends not only on the damping coefficient λ, but also on the stiffness nonlinearity parameter l , and the same preceding remark is again valid .

3)
Under-damped Response Analysis Equation (18) shows that the under-damped response predicted by the mathematical model consists of a product of sinusoidal oscillations with a decaying exponential behavior, raised to the l / 1 th power.This prediction is illustrated by the curves in Figure 4 showing qualitatively the dynamic behavior of under-damped mechanical systems.The curves display a decaying sinusoidal response in which the amplitude of successive peaks decreases to finally stabilize asymptotically to equilibrium zero value with time.It may be noted that following the parameter values selected the oscillations are only occurred over a brief time range.As mentioned previously, the decaying rate of the displacement is not only proportional to the damping factor λ, but depends also on the stiffness nonlinearity parameter l .This demonstrates that the damped oscillatory dynamics of the mechanical system under investigation can be again controlled from the magnitude of the stiffness rising parameter .l This again means that an assessment of the time constant characterizing the time displacement exponential decay may provide more info not only about the damping degree, but also on the stiffness strength of the investigated mechanical system.Moreover, the formula of t shows that the stiffness nonlinearity increases or decreases as the damping factor λ respectively increases or decreases.So, the current model allows simultaneously the enhancement of stiffness nonlinearity and damping of the mechanical system under question.

B. Stability analysis of damped responses
The above numerical simulations have shown that the system behavior relative to all of three damping modes of oscillation previously studied converges in the case where the oscillatory solution (18) is globally defined for 0  t , asymptotically to equilibrium zero value at an exponential rate of decay.Exact analytical results, that is to say ( 16), ( 17) and ( 18), determining the system response for over-damped, critically-damped and under-damped regimes respectively, indicate moreover that the system response tends to zero when the time   t .Therefore, it may be concluded that the damped nonlinear dynamics response of the system under study has an asymptotic stability character.This exponential asymptotic stability is not significantly sensitive to the initial conditions.As a consequence, for the long time system behavior, a determination of initial conditions is not necessary to be known with accuracy.It is worth to note that this stability character depends not only on the damping factor  as suggested by the classical second order damped linear oscillatory equation, but also on the hardening parameter l .The positive value of l is secured by the definition of the nonlinear restoring force function  that should obey in the Bauer's theory [7] to the linear Hooke's law when 0 ) (  t u . Next, special limits of ( 13) governing the nonlinear oscillatory dynamics of the one-dimensional viscoelastic system of interest will be established in order to demonstrate the ability of the developed mathematical model to capture some well known linear and nonlinear oscillatory dynamics models.

C. Limiting oscillatory equations
The objective of this section is to develop some special limits for (13) to demonstrate that the proposed mathematical model captures the classical damped linear harmonic oscillatory equation and the Lambert nonlinear oscillatory equation with well-known analytical exact solutions as special cases.One special limit of the current mathematical model may be obtained by substituting the specific value 1  l of the stiffness nonlinearity parameter into (13).In doing so, (13) becomes: It is easy to note that ( 20) is identical to the well known classical damped linear harmonic oscillator (4) previously noted, which is extensively studied in academic textbooks [26] and widely used in engineering design calculations.The damped linear harmonic oscillator (20) has a well-known exact analytical solution showing all of three damping modes of vibration following the magnitude of the damping compared to that of the natural frequency [6,26].Another special limit of ( 13) is achieved by setting the damping factor 0   .As a result, (13) transforms into: Equation ( 21) is the second order Lambert oscillatory equation with well-known exact analytical solution [27].Equation ( 21) is explored by He [27] as a strongly nonlinear oscillatory equation from Lambert.Following Pellicer and Solà-Morales [28], this limiting case has a significant interest for the investigated physical problem, since the special limit 0   m b  means that the system mass m is large relative to the damping coefficient b.Substituting 0   directly into (18), modeling the under-damped nonlinear oscillatory system response, it becomes possible to obtain as exact analytical solution for the Lambert oscillatory equation the following explicit expression: [ ]


The formula (23) is the same analytical solution found by He [27] using a variational approach.Putting moreover 1  l , directly into ( 22) or ( 23), the well-known sinusoidal response for the second order undamped harmonic oscillator is found.Therefore, the above proves mathematically the aptitude of the developed model to capture the Lambert oscillatory equation as special case.
It is worth noting that, compared to solutions of the classical damped linear harmonic equation, the argument of exponential functions in the developed nonlinear solutions differs by a factor of l which represents the stiffness rising parameter.Thus, the factor l appeared to be the fundamental parameter controlling the nonlinearity of the oscillatory dynamics of the system under study.From the above, it can be concluded that the proposed nonlinear oscillatory viscoelastic model captures in mathematical fashion the classical damped linear harmonic oscillator as a special case.The current mathematical model captures not only the classical second order damped linear harmonic equation as a subcase, but it has the ability to incorporate also the Lambert nonlinear oscillatory equation as a special case.

D. Validation of the mathematical model
The validity of the current mathematical model for predicting and numerically simulating the damped oscillatory response of a one-dimensional mechanical system exhibiting large deformations and viscoelastic properties is analyzed in this section.According to [26,29,30], exact analytical solutions found for the developed model are validated against numerical results.In this regard, exact analytical solutions determined for over-damped, critically-damped and underdamped regimes are compared in Figures (2), ( 3) and ( 4) with corresponding solutions carried out by numerical integration, respectively.So, Figure (2) shows that the exact analytical solution (16) matches associated numerical solution very well.This is confirmed by the low value of the mean squared error, 6.5189 010, emeasuring the existing discrepancy between these two results.Figure 3 15) is the exact solution, according to Keckic [18], found by Painlevé.The above confirms the reliability and accuracy of determined exact analytical solution of ( 13) and as a consequence, the ability of the developed mathematical model to numerically simulate the results of possible experiments [26].On the other hand, the preceding study on special limits of the proposed model demonstrates its powerfulness to be a nonlinear extension of the well known classical damped linear harmonic oscillator extensively used in engineering design calculations and to incorporate the Lambert oscillatory equation as special case.In doing so, the developed model offers its mathematical ability to be applied for simulating all of three damped oscillatory responses of a one-dimensional viscoelastic system under large deformations.

V. CONCLUSIONS
The viscoelastically damped linear harmonic equation is well known to be only applicable for a small range of deformation of mechanical systems.In these conditions, its use to characterize flexible mechanical systems in engineering applications may be the cause of catastrophe.In real working situations, mechanical systems undergo large deformations and show viscoelastic behavior.This involves the necessity to build reliable and satisfactory models for predicting, simulating and analyzing their response to an excitation.The present work was intended to develop a mathematical model that may be used in numerical simulations of the oscillatory dynamics of a onedimensional viscoelastic system experiencing finite deformations.The reliability of the proposed model is secured in the regard that it takes into consideration the nonlinearity properties of a mechanical system undergoing large deformations and viscoelastic behavior.In this perspective, a second-order first-degree Painlevé equation was developed from the application of Bauer's theory to model the nonlinear oscillatory dynamics of the system of interest.This equation is an extension of the damped linear harmonic oscillator equation widely employed in engineering design calculations for the nonlinear regime of behavior of real mechanical systems.The developed mathematical model captures also the Lambert nonlinear oscillatory equation as a special case.It is found that the obtained Painlevé evolution equation successfully models the dynamics of over-damped, critically-damped and underdamped nonlinear behaviors of a mechanical system under large deformations and viscoelastic behavior.The presented mathematical model provides the ability to control the damped oscillatory dynamics of the system under investigation concurrently from the damping coefficient or stiffness nonlinearity parameter.An estimation of the time constant characterizing the exponential decay of the time history of the displacement can, as shown by the current mathematical model, give knowledge not only on the damping strength, but also regarding the stiffness degree of the mechanical system of interest.Numerical and analytical results demonstrated that the proposed damped nonlinear oscillatory equation may be satisfactorily used in numerical simulations of viscoelastic oscillators.In this sense, the Bauer's theory significantly contributes to a better understanding of mathematical modeling of viscoelastic systems.This theory represents the nonlinear dynamics of viscoelastic systems in the form of Painlevé equation which is subject of intensive investigation in mathematics.Hence, the research work developed in this paper confirms the nature of the Bauer's theory to be a powerful mathematical tool to model satisfactorily the nonlinear dynamics of mechanical systems.Monsia's formulation of Bauer's theory may permit also to formulate Painlevé equation of the third order or in general of the nth order to describe mechanical systems by considering, instead of the linear Kelvin-Voigt constitutive law, the general setting of the linear viscoelastic constitutive law expressed as a single linear ordinary differential equation of the nth order relating the total strain and its time derivatives with the total stress.
) models then the nonlinear dynamics of the viscoelastic system with single degree of freedom for 1 l ¹ under an external exciting function ) can be solved analytically and numerically to investigate various types of oscillatory response of the system under study.In the sequel of this work, (12) will be studied in the case of unforced regime, to say, for the external exciting function ) Fig. 1.A typical behavior of (13) illustrating the oscillatory nature of the mathematical model.
obtained by using the Matlab function ode45 for reasonable values of model parameters fixed at 1expected response.The calculated value of the mean squared error provided by Matlab mse function is mse=6.5189e-010 .

Fig. 3 .
Fig. 3. Simulation results of the numerical solution compared with the exact analytical solution for critically-damped regime.

.
routine ode15s has been used to run the numerical simulation.Reasonable system parameters that generate the expected system response are: The mean squared error computed using the Matlab mse function is mse=8.6863e-008 .Matlab intrinsic function ode15s has been used here since ode45 becomes very inefficient.It may be suspected that, for some values of the stiffness rising parameter l , as the damping factor  becomes more small compared with the frequency o  2 , (13) becomes more stiff.

Fig. 4 .
Fig. 4. Comparison of analytical and numerical solutions for underdamped oscillation V. DISCUSSIONS This section is devoted for analyzing model predictions and demonstrating the validity of the current mathematical model.The stability character of solutions predicted by the model is also analyzed.The model predictions are discussed on the basis of numerical results presented in the preceding section.