Engineering, Technology & Applied Science Research Vol. 4, No.

6, 2014, 714-723 714

Simulating Nonlinear Oscillations of Viscoelastically

Damped Mechanical Systems
M. D. Monsia Y. J. F. Kpomahou
Department of Physics/FAST Department of Physics/FAST
University of Abomey-Calavi University of Abomey-Calavi
Abomey-Calavi, Benin Abomey-Calavi, Benin
monsiadelphin@yahoo.fr kpomaf@yahoo.fr

Abstract—The aim of this work is to propose a mathematical non-allowable or excessive deformations, self-excited
model in terms of an exact analytical solution that may be used in deformations, material fatigue and failure [1].
numerical simulation and prediction of oscillatory dynamics of a
one-dimensional viscoelastic system experiencing large The one-dimensional dynamics of continuous viscoelastic
deformations response. The model is represented with the use of media idealized as a bar is described formally by the Cauchy's
a mechanical oscillator consisting of an inertial body attached to wave equation. However, when only homogenous deformation
a nonlinear viscoelastic spring. As a result, a second-order first- is considered in longitudinal forced vibration experiment for a
degree Painlevé equation has been obtained as a law, governing bar carrying a tip mass at the free end, under the condition that
the nonlinear oscillatory dynamics of the viscoelastic system. the mass of the bar may be neglected compared to that of the
Analytical resolution of the evolution equation predicts the attached mass at the end point, the bar may be assumed to
existence of three solutions and hence three damping modes of behave as a simple viscoelastic spring. So, the Cauchy's
free vibration well known in dynamics of viscoelastically damped equation may be reduced, based on the second Newton's law,
oscillating systems. Following the specific values of damping if u denotes the time history of displacement response, then in
strength, over-damped, critically-damped and under-damped the concise form [2, 3]:
solutions have been obtained. It is observed that the rate of decay
is not only governed by the damping degree but, also by the
2
magnitude of the stiffness nonlinearity controlling parameter. d u (t )
Computational simulations demonstrated that numerical m 2
= Fext - Fint (1)
solutions match analytical results very well. It is found that the dt
developed mathematical model includes a nonlinear extension of 2
the classical damped linear harmonic oscillator and incorporates d
the Lambert nonlinear oscillatory equation with well-known where 2
denotes the ordinary second time derivative
dt
solutions as special case. Finally, the three damped responses of
the current mathematical model devoted for representing applied to the attached mass displacement u (t ) depending only
mechanical systems undergoing large deformations and on the time t, m is the mass of the attached body, Fext designates
viscoelastic behavior are found to be asymptotically stable. an applied exciting force, and Fint notes the internal axial force
due to the viscoelastic stress induced in the considered
Keywords-mathematical modeling; nonlinear oscillations;
viscoelastic oscillator; Painlevé equation; exact solution; numerical
mechanical system. In this context the question of viscoelastic
simulation systems vibration transforms into that of finding an appropriate
constitutive equation for structural materials.
I. INTRODUCTION Viscoelastic behavior is widely mathematically analyzed
through the use of rheological approach which discretizes a
In this work the dynamics of mechanical systems
viscoelastic body in its elementary elastic and viscous
undergoing large deformations and viscoelastic response is
components represented in terms of mechanical analogs. So,
investigated. A major topic in the dynamics of viscoelastic
the rheodynamical approach models the motion of a
systems is the problem of vibration. Vibration phenomenon
viscoelastic system in terms of ordinary differential equations
arises in all rigid or deformable systems, such as machines and
[3, 4]. According to Laroze [5], internal damping can be
engineering structures subjected to dynamic loading. So, the
schematized for most structural materials in practice as viscous
vibration problem is of vital importance for many fields of
damping. In this perspective, if it is assumed that the
science and technology. Vibration experiments are widely used
viscoelastic spring behaves as a Kelvin-Voigt medium, that is,
in the characterization of dynamical mechanical properties of
if:
engineering materials. Vibration is also desired for machines
under working conditions. However, for most structures in d f (u )
mechanical, biomechanical, civil, aeronautical and automotive Fint = k f (u ) + b (2)
engineering, oscillatory events prediction and control is dt
intensively required in order to reduce noise, and to prevent

www.etasr.com Monsia and Kpomahou: Simulating Nonlinear Oscillations of Viscoelastically Damped Mechanical Systems
 Engineering, Technology & Applied Science Research Vol. 4, No. 6, 2014, 714-723 715

where φ(u)=u, then the basic evolution equation governing the [7, 8]. The above show that an oscillatory viscoelastic system is
forced viscoelastically damped linear vibrating mechanical intrinsically characterized at least by its stiffness, damping and
system may be written in the form: inertia nonlinearities. Therefore, a reliable representation of the
nonlinear oscillatory dynamics of a viscoelastic system should
2 have the ability to mathematically incorporate in the governing
d f (u ) d f (u )
m +b + k f (u ) = Fext (3) equation these basic nonlinearity principles [7]. From the above
dt
2
dt analysis, the most general second order ordinary differential
equation that can model the nonlinear oscillatory dynamics of a
or in terms of only u (t ) under free vibration: single degree of freedom viscoelastic system under unforced
conditions may take the form:
2
d u du
2
+l
2
+ w o u = 0 (4) h(u, u )u  g (u, u )u  f (u )u  0 (5)
dt dt
b k where the dot over a symbol designates the time derivative,
h(u , u ) and g (u , u ), are function of the displacement u (t ) and
2
in which l = , and w o = .
m m
The quantity ωο is defined as the angular or undamped first time derivative u (t ), and f(u) is function of the
natural frequency of the physical system, and λ is a damping displacement u (t ) . These functions lead to the inertia, damping
factor [5, 6]. Equation (4) is a classical prototype of second- and stiffness nonlinearity properties, respectively, so that for
order ordinary differential equations used to describe the the function h  0 , and g / h   and f / h  o2 , (5) gives (4)
damped linear oscillation of a single degree of freedom
oscillator. It represents the dynamics of the motion of a one- modeling a damped linear harmonic oscillator. It is interesting
dimensional viscoelastic system in terms of Kelvin-Voigt to note that in the majority of existing evolution models of
element with the addition of a mass in the range of linear mechanical systems, only one of these nonlinearities is often
deformation. Equation (4) is the simplest differential equation considered (see for example [16] for more details). This shows
that may reproduce all different types of behavior exhibited by that, due to the mathematical complication arising quickly in
a damped second order oscillation system in the range of linear governing equations, the enhancement, for example, of
deformation. The type of response given by (4) representing a stiffness nonlinearity is often performed in models to the
damping linear oscillation system depends, in effect, on the detriment of system damping nonlinearity and, inversely, the
strength of damping degree. In other words, the response of a improvement of damping nonlinearity is made in a prejudicial
damped oscillatory system is very sensitive to changes in fashion to that of stiffness nonlinearity [1]. Moreover, very few
specific value of damping. of proposed mathematical models for studying the dynamics of
viscoelastic systems with a single degree of freedom have been
It has been observed that the oscillatory response of a real performed to enclose on the one hand the inertia nonlinearity,
system is nonlinearly damped and geometrically nonlinear, and simultaneously combined inertia, damping and stiffness
leading to the dependence of the stiffness on the induced nonlinearities governing the viscoelastic response as indicated
deformation or stress. This phenomenon is known as stiffening by (5) on the other hand [7].
or softening and terminates often by failure, following the
stiffness increases or decreases. Such a phenomenon could not It is well known again that nonlinear problems having
be predicted and explained by any linear model equation [6- explicit exact solutions in terms of elementary standard
11]. Many systems in engineering applications are designed to functions are very limited in physical and engineering fields.
behave not only viscoelastically but also nonlinearly, to say, to So, a large part of nonlinear analysis has only been performed
undergo large deformations exceeding the limiting value on the basis of qualitative theory or particular solutions derived
predicted by linear theory in loading environment [2, 6, 7, 11]. from analytical approximate or numerical integration methods.
Viscoelastic systems have the specific ability to experience In particular, homotopy perturbation analysis is intensively
large deformations even for moderate force levels. Hence, the used to investigate nonlinear vibration problems in mechanical
linear theory becomes quite unable to explain the dynamics of structures [17]. In mechanical system design calculations, for
viscoelastic systems experiencing finite deformations [6, 7, 12- example, the accurate determination of dynamical
15] . characteristics from observations requires mathematical models
having appropriate analytical solutions. In this context, the
Geometrical, damping and material nonlinearities are design of mathematical models capable of representing the
fundamental causes for the viscoelastic behavior of systems to nonlinear dynamics of viscoelastic systems satisfactorily in
be nonlinear [5]. If the stiffness or geometrical nonlinearities terms of analytical exact solutions, capturing also
are relatively well captured in terms of polynomial expansion, simultaneously the combined nonlinear phenomena, becomes a
the damping nonlinearity properties of mechanical systems are major necessity. In that, the Bauer’s rheological-dynamical
very difficult to be known [4, 7, 8]. According to [7, 8], another theory [7] consists of a notable progress in the field of the
non-negligible source contributing to nonlinear response of a dynamics of continuous viscoelastic media, since it meets the
mechanical system is the inertia properties. Inertia produces, to essential of preceding criteria about the necessity to handle
a certain extent, a force which is essentially nonlinear. simultaneously and in combined fashion, in mathematical
Nonlinear inertia forces are generally proportional to higher modeling of viscoelastic systems, the nonlinearity properties.
powers of velocity and acceleration of the mechanical system

www.etasr.com Monsia and Kpomahou: Simulating Nonlinear Oscillations of Viscoelastically Damped Mechanical Systems
 Engineering, Technology & Applied Science Research Vol. 4, No. 6, 2014, 714-723 716

Although the Bauer's method [7] seems to be simple in II. THE ONE-DIMENSIONAL NONLINEAR THEORY OF
formulation, it is in reality a powerful approach for solving VISCOELASTIC SYSTEMS
dynamical nonlinear problems arising in viscoelastic
continuum mechanics successfully. This theory [7] entails a A. Continuum mechanical nonlinear evolution equation
significant mathematical modification and extension of the In this part, the Monsia’s formulation [12, 13] of the
classical Kelvin-Voigt viscoelastic solid law for capturing Bauer’s theory [7] will be briefly reviewed to make the
simultaneously the combined inertia, damping and stiffness transition from its continuum version to a constitutive
nonlinearities characterizing the dynamics of real viscoelastic expression in terms of displacement and force. The advanced
systems. The Bauer’s approach [7] was originally developed procedure for modeling the dynamics of real mechanical
for modeling the dynamics of nonlinear viscoelastic response systems, to say, viscoelastic systems, developed by Bauer [7]
of soft biological tissues, arterial walls in particular. The basic was expressed within the framework of continuous mechanics,
idea underlying this theory consists in developing the total in other words, in terms of stresses and strains experienced by
stress within a material system as the sum of three basic the mechanical system. The Monsia’s formulation [12, 13] of
stresses: elastic, viscous and inertial stresses, operating in the Bauer’s theory [7] enables us to describe this theory, for a
parallel. Very recently, the Bauer’s rheological-dynamical given stiffness nonlinearity function f , by a single second
theory was formulated in a simple mathematical expression
that may be described by a single second order evolution order evolution equation, with only few system parameters to
equation within the framework of continuum mechanics for be determined through fitting procedure, of the form:
investigating the dynamics of viscoelastic material systems [12,
13]. This formulation has successfully been applied in several f ¢ ( e ) e + f ¢¢ ( e ) e + lf ¢ ( e ) e + w o f ( e ) = (1 / c ) s (t ) (6)
2 2

papers to model creep deformation [13], creep relaxation [6]

and deformation restoration process under stress relaxation where s (t ) is the scalar stress function and e (t ) the scalar
conditions [14, 15] of a variety of viscoelastic solid bodies. strain function , the prime denotes a differentiation with respect
The objective of this research work is to develop a to response variable, that is, here, the strain e (t ) , and c ¹ 0 , is
mathematical model expressed in terms of an exact analytical the inertia coefficient. The constants λ and  o continue to have
solution that may be useful in numerical simulations of the the preceding definitions of damping factor and natural
oscillatory dynamics of a one-dimensional viscoelastic system frequency, respectively. For some convenience, (6) may be
exhibiting large deformations. The physical properties of the written, by using a differential operator B as:
mathematical model are represented by a single degree of
freedom oscillator consisting of a mass attached to a B ( e (t )) = (1 / c )s (t ) (7)
viscoelastic spring that behaves as a nonlinear Kelvin-Voigt with:
continuum medium. More precisely, given a stiffness
nonlinearity function, it is proposed to develop the evolution
B = f ¢ ( d / dt ) + f ¢¢ ( d / dt ) + lf ¢ ( d / dt ) + w o f (8)
2 2 2 2

equation of the time dependent displacement response of a

mechanical system experiencing viscoelastic response and At the present time, the function  that captures the
large deformations by applying the Bauer's theory as
formulated mathematically by Monsia [12, 13]. In this sense, system stiffness nonlinearity is not explicitly known, but it
exact analytical solutions in terms of elementary standard must be specified before any further use of the theory. In the
functions have been determined, following various types of following section, an explicit expression for  will be given as
response able to exhibit the evolution equation under unforced function of the deformation response.
regime. Numerical applications are carried out for illustrating
the ability of the mathematical model to be used for numerical B. Nonlinear constitutive force-displacement equation
simulations. Graphs of numerical and exact analytical results In this section the constitutive equation in terms of force
are compared to prove the validity of the current model. and displacement governing the dynamics of the one-
It is also shown that the developed model incorporates as dimensional viscoelastic system under investigation is
special cases some oscillatory equations with well known formulated. For a time history of displacement response
solutions by passing to special limiting values of the model u (t ) and an external exciting function F (t , u , u , u) , the
parameters. Time-asymptotic system response is investigated to operator B allows the equation of motion to be represented
check the stability character of the long time behavior of the by:
proposed model. So, in the next section, the evolution equation
governing the nonlinear oscillatory dynamics of the B (u (t )) = (1 / m ) F (t , u , u , u) (9)
investigated viscoelastic system will be developed. Section 3 is
devoted to solve the obtained second-order first-degree or by the equation written in the developed form:
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, u (u )  u 2 (u )  u (u )   o2 (u )  (1 / m ) F (t , u , u , u) (10)
and a discussion of the results is performed in Section 5. The where the prime denotes the differentiation with respect to
final section gives a brief conclusion of the work. displacement variable u(t). The constants m, λ and ωο are time

www.etasr.com Monsia and Kpomahou: Simulating Nonlinear Oscillations of Viscoelastically Damped Mechanical Systems
 Engineering, Technology & Applied Science Research Vol. 4, No. 6, 2014, 714-723 717

independent and continue to have the preceding definitions. damped, and over-damped responses of the mechanical system
Equation (10) is a second-order first-degree Painlevé equation under question. In doing so, (12) becomes for unforced motion:
with a forcing function [18-20] known to be subject of many
studies in mathematics. According to Roth [21], Painlevé uu + ( l - 1) u + l uu + ( w o / l ) u = 0
2 2 2

equations have found several applications in physics,

particularly in the field of circuit oscillations. So, the As can be seen easily, the preceding Painlevé equation has
perspective to model the mechanical behavior of viscoelastic u (t )  0 as trivial solution. In the usual standard form, this
systems using the second-order first-degree Painlevé equation
may provide a great insight on their dynamics. Equation (10) equation becomes:
gives the differential relationship between the exciting function
F (t , u , u , u) and the resulting time displacement response u(t) u + [ (l - 1)u / u + l ] u + ( w o / l )u = 0
2
(13)
for a given function φ(u) capturing the purely nonlinear elastic
response of the system under consideration. This equation with u (t )  0 .
enters in the perspective of the general second order evolution
equation of a mechanical system represented by (5) and This second order autonomous differential equation can be
consists of a general formulation requiring the specification of viewed as a nonlinearly damped equation. The damping force
the nonlinear elastic spring force function φ(u) for any is a nonlinear function in both the time history of response
subsequent development. Mathematically, all nonlinear u (t ) and the velocity u (t ) . The dependence on the variable
function φ(u) that tends towards u, when the displacement u response u (t ) introduces a nonlinear Newtonian singularity at
tends to zero, may be selected as a possible restoring force the origin as noted in central forces. The dependence on the
function. In other words, at small deformations, the nonlinear velocity includes a combination of linear and quadratic
model under question should reduce to the associated linear damping terms. In this perspective (13) consists of an attractive
law. Therefore, an infinite number of functions φ may be problem to be investigated to learn more about the qualitative
designed. In doing so, various types of mathematical oscillatory dynamics of the system around the singular point
expressions for the function φ have been proposed in recent origin from a viewpoint of phase plane analysis which is not
papers [6, 12-15, 22]. In general, as performed by Bauer [7], the topic of the following research work. So, only a quantitative
the stiffness nonlinearity function φ(u) may be expanded in a analysis of (13) will be performed concerned in this study. The
power series as: graph in Figure 1 illustrates the oscillatory behavior of (13)
subject to arbitrary initial conditions u(t=0)=5, and
 (u )  a1u  a 2 u 2  a3u 3  ...  a n u n  ... u (t  0)  1 , from t=0 to t=20, obtained by numerical
integration. The simulation is run by using Matlab's routine
In this perception, for mathematical reasons of simplicity, ode45 which exploits Runge-Kutta methods. The parameter
the nonlinear stiffness function φ will be expressed, in this
work, basically as a power law in the form: [22] values are fixed at l  0.987 ,   0.01 , and o  1 .

φ (u ) = u l (11)
where the hardening exponent, that is the stiffness nonlinearity
controlling parameter l ¹ 0 . Following this point of view, the
equation of motion takes the form:

u l -1u + (l - 1)u l -2 u 2 + lu l -1u + ( w o2 / l )u l = (1 / lm ) F (t , u , u , u) (12)

with the mass m ¹ 0 . Equation (12) models then the nonlinear

dynamics of the viscoelastic system with single degree of
freedom for l ¹ 1 under an external exciting function
F (t , u , u , u) . It 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 Fig. 1. A typical behavior of (13) illustrating the oscillatory nature of the
mathematical model.
the case of unforced regime, to say, for the external exciting
function F (t , u , u , u) equal to zero. In the following section, the ability of the developed
theoretical model to mathematically exhibit various responses
III. EXACT ANALYTICAL SOLUTION FOR THE NONLINEAR of the damped nonlinear oscillatory dynamics of the system
OSCILLATIONS EQUATION under investigation is proved. To do so, (13) will be solved in
closed form exact solution by using suitable initial conditions
This part determines for the unforced oscillatory regime, that satisfy the dynamics of the nonlinear viscoelastic system of
that is to say, F (t , u , u , u)  0 , the exact analytical solution for interest.
the equation of motion following the under-damped, critical

www.etasr.com Monsia and Kpomahou: Simulating Nonlinear Oscillations of Viscoelastically Damped Mechanical Systems
 Engineering, Technology & Applied Science Research Vol. 4, No. 6, 2014, 714-723 718

A. Reduction of the nonlinear oscillatory equation to the the characteristic equation are real, distinct and negative. The
damped linear harmonic equation solution of (14) becomes:
It is well known in mathematics that the range of nonlinear
evolution equations that can be analytically solved in the form y (t ) = A1 exp( r1t ) + A2 exp( r2 t )
of elementary standard functions is very limited [23]. So, the
objective to find exact solutions of nonlinear evolution where:
equation leads often to the question about the explicit
integration of this equation [23-25]. In this perspective, the r1 = -l (1 + d ) / 2 and r2 = -l (1 - d ) / 2
Painlevé analysis may be performed in order to conclude on the
integration of this equation [23-25]. However, following with:
Kudryashov [25], the determination of exact solutions of
nonlinear evolution equation can be investigated without 2 2
having to apply successfully the Painlevé test. Here, the d= 1 - 4 wo l
governing equation (13) under question is a typical equation of
Painlevé which has been, according to Keckic [18, 19], A1 and A2 are two constants of integration defined by the
integrated by quadratures . Hence, the problem of integrability initial conditions. Assuming that the conditions
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 u (t ) = uo , and u (t ) = vo , when t = 0
recently been solved on the basis of mathematical
transformations leading to a Riccati equation or a damped satisfy the past history of the displacement, and taking into
linear harmonic equation [12-15, 22]. In the present work, (13) consideration (15), the time history of the displacement can be
will be turned in the damped linear harmonic oscillatory written in the explicit analytical expression:
equation. Hence, to solve (13) in the form of standard
functions, a change of variable is required for the expected u (t ) = u o exp(-lt / 2l ) ⋅
oscillatory solutions. Performing the following suitable
é é 2lv ùù
1/ l
substitution: 2 2
ê ê ( o + l ) sinh( l - 4wo t ) +úú
ê 1
ê u 2
úú
ê( ) êê o úú (16)
y = ul ê úú
ê l - 4w o ê 2 l - 4w o t úú
2 2
2 2
Equation (13) transforms after a few algebraic ê ê l - 4wo cosh(
2
) úú
manipulations, into: êë êë 2 úúûû

y + l y + w o2 y = 0 (14)
C. Exact analytical solution for critically-damped nonlinear
Equation (14) is the well known classical linear second
response
order evolution equation previously noted as (4). Therefore, the
general solution of (13) becomes: [19] In this case l = 2w , the two solutions of the characteristic
o

equation coincide and then, the solution of (14) is given by:

u (t ) = y (t )1/ l (15)
y (t ) = ( B1t + B2 ) exp(-l t / 2)
where y (t ) designates the general solution of (14)
A mathematical analysis of (14) indicates that the nature of where B1 and B 2 are real constants determined by initial
(15) solution depends on the relative magnitudes of the conditions. Proceeding in the same way for the initial
damping factor λ and the stiffness module, to say, the natural conditions:
frequency ωο. In other words, the nature of solution depends on
the roots of characteristic equation: u (t ) = uo , and u (t ) = vo , when t = 0
2 2
r + lr + wo = 0 and using (15), the solution of the equation of motion becomes
associated with (14), that can be real or complex numbers. So,
u (t ) = uo exp(-lt / 2l ) [1 + (lvo / uo + l / 2)t ]
1/ l
three damping modes of vibration should be distinguished. (17)

B. Exact analytical solution for over-damped nonlinear D. Exact analytical solution for under-damped nonlinear
response oscillations
This case corresponds to a strong damping [6, 26], that is, a This case assumes a relative weak damping, to say,
relative large damping where l > 2w . The roots r1 and r2 of
o l < 2w . Here, the roots of the characteristic equation are
o

complex numbers. Therefore, the time history of the response

www.etasr.com Monsia and Kpomahou: Simulating Nonlinear Oscillations of Viscoelastically Damped Mechanical Systems
 Engineering, Technology & Applied Science Research Vol. 4, No. 6, 2014, 714-723 719

of the vibrating viscoelastic system, under the same initial A. Numerical results for over- damped nonlinear response
conditions taking into account the past history of the The graph of the numerical solution for over-damped
displacement: nonlinear regime is compared to that of the exact analytical
solution (16) in Figure 2 over the time range from t  0 to
u (t ) = uo , and u (t ) = vo , when t = 0 t  20 under the initial conditions u (0) = 5 and u (0)  1 .
Numerical results are obtained by using the Matlab function
can, for: ode45 for reasonable values of model parameters fixed at
l = 1.999,   3 , and o  0.9998 to produce the expected
y (t ) = (C1 cos w t + C 2 sin w t ) exp(-l t / 2) response. The calculated value of the mean squared error
provided by Matlab mse function is mse=6.5189e-010 .
C1 and C 2 are two constants of
2 2
where w = w o - l / 4 , and,
integration fixed by the initial conditions, take the explicit
analytical form:
u (t ) = uo exp(-lt / 2l ) ⋅
1/ l
é ù
ê 1 lv l
ê( )( o + )sin( wo2 - l 2 / 4t ) + cos( wo2 - l 2 / 4t )úú (18)
ê w o - l / 4 uo
2 2 2 ú
ë û

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 Fig. 2. Comparison between analytical (solid line) and numerical (circles)
Runge-Kutta explicit procedure. The integration function uses a solutions for over-damped nonlinear response
variable integration step size to provide the needed accuracy as
minimizing computation time. The Matlab ordinary differential B. Numerical results for critically- damped nonlinear
equation solver ode15s is based on an implicit method. It is a response
multistep solver that needs the solutions at numerous previous
time points in order to estimate the current solution. In this
Figure 3 shows the comparison of the numerical solution
perspective, (13) should be written in state-space
with the exact analytical solution (17) of (13) subject to the
representation, to say, represented in terms of a set of first- initial conditions u (0)  0.1 and u (0)  vo  0.5 , over the time
order differential equations. Therefore, defining the state range from t  0 , to t  20 . Reasonable model parameter
variables: values for this simulation are l  1 .5 ,   2 , and then o  1 .
Matlab routine ode45 has been used for the simulation. The
x1 (t ) = u (t ) , x2 (t ) = u (t ) mean squared error calculated is mse=7.3336e-010 .

the following system of two first-order differential

equations may be obtained:

x1 (t ) = x (t )
2
(19)
x 2 (t ) = -(l -1) x22 / x1 - l x2 - (wo2 / l ) x1
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.
Fig. 3. Simulation results of the numerical solution compared with the
exact analytical solution for critically-damped regime.

www.etasr.com Monsia and Kpomahou: Simulating Nonlinear Oscillations of Viscoelastically Damped Mechanical Systems
 Engineering, Technology & Applied Science Research Vol. 4, No. 6, 2014, 714-723 720

C. Numerical results for under- damped nonlinear attaining its maximum value, the time history of the
oscillations displacement u(t) declines gradually without oscillations to
Figure 4 compares the numerical solution to the exact approach asymptotically equilibrium zero value with time. The
analytical solution (18) of (13) under the initial conditions that time history of response u(t) of the system presents only, as
expected again, one maximum value. Analytical result, to say,
u (0)  1.01 , and u (0)  0.94 , over the time interval [ 0, 20 ] . solution (16) shows that the system response has a hyperbolic
Here, Matlab routine ode15s has been used to run the behavior modulated by decay exponential, similarly to the
numerical simulation. Reasonable system parameters that over-damped linear harmonic oscillator response, but raised to
generate the expected system response are: l  1/ 3 ,   0.5 , the 1/ l th power. Equation (16) shows that the rate of decay is
and o  0.7 . The mean squared error computed using the not only proportional to the damping factor λ, but depends also
Matlab mse function is mse=8.6863e-008 . Matlab intrinsic on the stiffness nonlinearity rising parameter l . So, it may be
function ode15s has been used here since ode45 becomes very possible to use this parameter to control the oscillatory
inefficient. It may be suspected that, for some values of the response amplitude of the system concurrently with the
stiffness rising parameter l , as the damping factor  becomes damping factor λ.
more small compared with the frequency 2o , (13) becomes 2) Critically-damped Response Analysis
more stiff. The curves in Figure 3 show that the time dependent
displacement u(t ) 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 1/ l 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
Fig. 4. Comparison of analytical and numerical solutions for under- successive peaks decreases to finally stabilize asymptotically to
damped oscillation equilibrium zero value with time. It may be noted that
following the parameter values selected the oscillations are
V. DISCUSSIONS only occurred over a brief time range. As mentioned
This section is devoted for analyzing model predictions and previously, the decaying rate of the displacement is not only
demonstrating the validity of the current mathematical model. proportional to the damping factor λ, but depends also on the
The stability character of solutions predicted by the model is stiffness nonlinearity parameter l . This demonstrates that the
also analyzed. The model predictions are discussed on the basis damped oscillatory dynamics of the mechanical system under
of numerical results presented in the preceding section. investigation can be again controlled from the magnitude of the
stiffness rising parameter l. This again means that an
A. Analysis of model predictions assessment of the time constant t = 2l l characterizing the
Vibrating material systems in real operating situation time displacement exponential decay may provide more info
experience, as it is well mentioned previously, large not only about the damping degree, but also on the stiffness
deformations and viscoelastic behavior. So, their dynamics is strength of the investigated mechanical system. Moreover, the
characterized by a damped nonlinear oscillation. In this formula of t shows that the stiffness nonlinearity increases or
context, a reliable theory devoted for modeling the damped decreases as the damping factor λ respectively increases or
nonlinear oscillatory dynamics of these systems should, at decreases. So, the current model allows simultaneously the
least, have the ability to handle some fundamental nonlinearity enhancement of stiffness nonlinearity and damping of the
problems, and predict all of three damping modes of mechanical system under question.
oscillation known for a real damped oscillatory system. B. Stability analysis of damped responses
Consequently, the following subsections investigate the
The above numerical simulations have shown that the
aptitude of the proposed mathematical model to satisfactorily
system behavior relative to all of three damping modes of
reflect these three damped oscillation responses.
oscillation previously studied converges in the case where the
1) Over-damped Response Analysis oscillatory solution (18) is globally defined for t  0 ,
The graphs in Figure 2 illustrate the over-damped behavior asymptotically to equilibrium zero value at an exponential rate
of the proposed model. It can be seen, as expected, that after

www.etasr.com Monsia and Kpomahou: Simulating Nonlinear Oscillations of Viscoelastically Damped Mechanical Systems
 Engineering, Technology & Applied Science Research Vol. 4, No. 6, 2014, 714-723 721

of decay. Exact analytical results, that is to say (16), (17) and b

(18), determining the system response for over-damped,   0 means that the system mass m is large relative to
m
critically-damped and under-damped regimes respectively,
indicate moreover that the system response tends to zero when the damping coefficient b. Substituting   0 directly into
the time t   . Therefore, it may be concluded that the (18), modeling the under-damped nonlinear oscillatory system
damped nonlinear dynamics response of the system under response, it becomes possible to obtain as exact analytical
study has an asymptotic stability character. This exponential solution for the Lambert oscillatory equation the following
asymptotic stability is not significantly sensitive to the initial explicit expression:
conditions. As a consequence, for the long time system
u (t ) = uo [ cos(wo t ) + (1 / wo )(lvo / uo ) sin(wo t ) ]
1/ l
behavior, a determination of initial conditions is not necessary (22)
to be known with accuracy. It is worth to note that this stability
character depends not only on the damping factor  as or
suggested by the classical second order damped linear u (t )  C cos(ot )  D sin(ot )1 / l (23)
oscillatory equation, but also on the hardening parameter l .
The positive value of l is secured by the definition of the luol 1vo
by setting C  uol and D 
nonlinear restoring force function  (u )  u l that should obey in o
the Bauer's theory [7] to the linear Hooke's law when The formula (23) is the same analytical solution found by
u (t )  0 . Next, special limits of (13) governing the nonlinear He [27] using a variational approach. Putting moreover l  1 ,
oscillatory dynamics of the one-dimensional viscoelastic directly into (22) or (23), the well-known sinusoidal response
system of interest will be established in order to demonstrate for the second order undamped harmonic oscillator is found.
the ability of the developed mathematical model to capture Therefore, the above proves mathematically the aptitude of the
some well known linear and nonlinear oscillatory dynamics developed model to capture the Lambert oscillatory equation as
models. special case.
C. Limiting oscillatory equations It is worth noting that, compared to solutions of the
The objective of this section is to develop some special classical damped linear harmonic equation, the argument of
limits for (13) to demonstrate that the proposed mathematical exponential functions in the developed nonlinear solutions
model captures the classical damped linear harmonic differs by a factor of l which represents the stiffness rising
oscillatory equation and the Lambert nonlinear oscillatory parameter. Thus, the factor l appeared to be the fundamental
equation with well-known analytical exact solutions as special parameter controlling the nonlinearity of the oscillatory
cases. One special limit of the current mathematical model may dynamics of the system under study. From the above, it can be
be obtained by substituting the specific value l  1 of the concluded that the proposed nonlinear oscillatory viscoelastic
stiffness nonlinearity parameter into (13). In doing so, (13) model captures in mathematical fashion the classical damped
becomes: 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
u + l u + w o2 u = 0 (20)
the ability to incorporate also the Lambert nonlinear oscillatory
It is easy to note that (20) is identical to the well known equation as a special case.
classical damped linear harmonic oscillator (4) previously D. Validation of the mathematical model
noted, which is extensively studied in academic textbooks [26]
and widely used in engineering design calculations. The The validity of the current mathematical model for
damped linear harmonic oscillator (20) has a well-known exact predicting and numerically simulating the damped oscillatory
analytical solution showing all of three damping modes of response of a one-dimensional mechanical system exhibiting
vibration following the magnitude of the damping compared to large deformations and viscoelastic properties is analyzed in
that of the natural frequency [6, 26]. Another special limit of this section. According to [26, 29, 30], exact analytical
solutions found for the developed model are validated against
(13) is achieved by setting the damping factor   0 . As a
numerical results. In this regard, exact analytical solutions
result, (13) transforms into:
determined for over-damped, critically-damped and under-
damped regimes are compared in Figures (2), (3) and (4) with
u 2 o2 corresponding solutions carried out by numerical integration,
u  (l  1)  u0 (21) respectively. So, Figure (2) shows that the exact analytical
u l
solution (16) matches associated numerical solution very well.
Equation (21) is the second order Lambert oscillatory This is confirmed by the low value of the mean squared error,
equation with well-known exact analytical solution [27]. 6.5189e - 010, measuring the existing discrepancy between
Equation (21) is explored by He [27] as a strongly nonlinear these two results. Figure 3 leads to the conclusion that both
oscillatory equation from Lambert. Following Pellicer and exact analytical solution (17) and numerical result for
Solà-Morales [28], this limiting case has a significant interest critically-damped system response are in very good agreement
for the investigated physical problem, since the special limit as underlined by the mean squared error

www.etasr.com Monsia and Kpomahou: Simulating Nonlinear Oscillations of Viscoelastically Damped Mechanical Systems
 Engineering, Technology & Applied Science Research Vol. 4, No. 6, 2014, 714-723 722

mse = 7.3336e - 010. A very good agreement is also satisfactorily used in numerical simulations of viscoelastic
observed regarding the comparison of exact analytical solution oscillators. In this sense, the Bauer's theory significantly
(18) with numerical result for under-damped response of the contributes to a better understanding of mathematical modeling
studied system, as can be noted by the value 8.6863e  008 for of viscoelastic systems. This theory represents the nonlinear
the mean squared error. Besides, (15) is the exact solution, dynamics of viscoelastic systems in the form of Painlevé
according to Keckic [18], found by Painlevé. The above equation which is subject of intensive investigation in
confirms the reliability and accuracy of determined exact mathematics. Hence, the research work developed in this paper
analytical solution of (13) and as a consequence, the ability of confirms the nature of the Bauer's theory to be a powerful
the developed mathematical model to numerically simulate the mathematical tool to model satisfactorily the nonlinear
results of possible experiments [26]. On the other hand, the dynamics of mechanical systems. Monsia's formulation of
preceding study on special limits of the proposed model Bauer's theory may permit also to formulate Painlevé equation
demonstrates its powerfulness to be a nonlinear extension of of the third order or in general of the nth order to describe
the well known classical damped linear harmonic oscillator mechanical systems by considering, instead of the linear
extensively used in engineering design calculations and to Kelvin-Voigt constitutive law, the general setting of the linear
incorporate the Lambert oscillatory equation as special case. In viscoelastic constitutive law expressed as a single linear
doing so, the developed model offers its mathematical ability to ordinary differential equation of the nth order relating the total
be applied for simulating all of three damped oscillatory strain and its time derivatives with the total stress.
responses of a one-dimensional viscoelastic system under large REFERENCES
deformations.
[1] E. I. Rivin, “Use of Stiffness /Damping/Natural Frequency Criteria
inVibration Control”, IFToMM World Congress, Besanson, June 2-6,
V. CONCLUSIONS 2007
The viscoelastically damped linear harmonic equation is [2] F. J. Lockett, Nonlinear Viscoelastic Solids, Academic Press, 1972
well known to be only applicable for a small range of [3] Y. M. de Haan, G. M. Sluimer, “Standard linear solid model for dynamic
and time dependent behaviour of building materials”, Heron , Vol. 46,
deformation of mechanical systems. In these conditions, its use No. 1, pp.49-76, 2001
to characterize flexible mechanical systems in engineering [4] H. H. Hilton, S. Yi, “Generalized viscoelastic 1-dof deterministic
applications may be the cause of catastrophe. In real working nonlinear oscillators”, Nonlinear Dynamics, Vol. 36, No. 2-4, pp. 281-
situations, mechanical systems undergo large deformations and 298, 2004
show viscoelastic behavior. This involves the necessity to build [5] S. Laroze, Mécanique des Structures. Tome 3, Masson, 1992
reliable and satisfactory models for predicting, simulating and [6] M. D. Monsia, “A mathematical model for predicting the relaxation of
analyzing their response to an excitation. The present work was creep strains in materials”, Physical Review & Research International,
intended to develop a mathematical model that may be used in Vol.2, No. 3, pp.107-124, 2012
numerical simulations of the oscillatory dynamics of a one- [7] R. D. Bauer, “Rheological approaches of arteries”, Biorheology Suppl.
dimensional viscoelastic system experiencing finite I, Vol. 1, pp.159-167, 1984
deformations. The reliability of the proposed model is secured [8] G. Kerschen, K. Worden, A. F. Vakakis, J-.C. Golinval, “Past, present
and future of nonlinear system identification in structural dynamics”,
in the regard that it takes into consideration the nonlinearity Mechanical Systems and Signal Processing, Vol. 20, No. 3, pp.505-592,
properties of a mechanical system undergoing large 2006
deformations and viscoelastic behavior. In this perspective, a [9] H. Jrad, J. L. Dion, F. Renaud, I. Tawfiq, M. Haddar, “Experimental
second-order first-degree Painlevé equation was developed characterization, modeling and parametric identification of the non linear
from the application of Bauer's theory to model the nonlinear dynamic behavior of viscoelastic components”, European Journal of
oscillatory dynamics of the system of interest. This equation is Mechanics-A/Solids, Vol. 42 pp. 176-187, 2013
an extension of the damped linear harmonic oscillator equation [10] H. Jrad, J. L. Dion, F. Renaud, I. Tawfiq, M. Haddar, “A new approach
for nonlinear generalized Maxwell model for depicting dynamic
widely employed in engineering design calculations for the behaviour of viscoelastic elements-parameters identification and
nonlinear regime of behavior of real mechanical systems. The validation”, 18th Symposium of Vibrations, Shocks and Noise, France,
developed mathematical model captures also the Lambert July 3-5, 2012
nonlinear oscillatory equation as a special case. It is found that [11] M. D. Monsia, “A Simplified nonlinear generalized Maxwell model for
the obtained Painlevé evolution equation successfully models predicting the time dependent behavior of viscoelastic materials”, World
the dynamics of over-damped, critically-damped and under- Journal of Mechanics, Vol.1, No. 3, pp. 158-167, 2011
damped nonlinear behaviors of a mechanical system under [12] M. D. Monsia, “Modeling the nonlinear rheological behavior of
large deformations and viscoelastic behavior. The presented materials with a hyper-exponential type function”, Mechanical
Engineering Research, Vol.1, No. 1, pp. 103-109, 2011
mathematical model provides the ability to control the damped
[13] M. D. Monsia, “A nonlinear mechanical model for predicting the
oscillatory dynamics of the system under investigation dynamics response of materials under a constant loading”, Journal of
concurrently from the damping coefficient or stiffness Materials Science Research, Vol.1, No. 1, pp. 90-100, 2012
nonlinearity parameter. An estimation of the time constant [14] M. D. Monsia, Y. J. F. Kpomahou, “Predicting the dynamic behavior of
characterizing the exponential decay of the time history of the materials with a nonlinear modified Voigt model”, Journal of Materials
displacement can, as shown by the current mathematical Science Research, Vol. 1, No. 2, pp. 166-173, 2012
model, give knowledge not only on the damping strength, but [15] M. D. Monsia, Y. J. F. Kpomahou, “A theoretical characterization of
also regarding the stiffness degree of the mechanical system of time dependent materials by using a hyperlogistic-type model'”,
Mechanical Engineering Research, Vol. 2, No. 1, pp. 36-43, 2012
interest. Numerical and analytical results demonstrated that the
proposed damped nonlinear oscillatory equation may be

www.etasr.com Monsia and Kpomahou: Simulating Nonlinear Oscillations of Viscoelastically Damped Mechanical Systems
 Engineering, Technology & Applied Science Research Vol. 4, No. 6, 2014, 714-723 723
[16] L. Cveticanin, “Oscillator with strong quadratic damping force”, [23] R. Conte, “Partial integrability of the anharmonic oscillator”, Journal of
Publications de l'Institut Mathématique, Nouvelle série, Vol. 85, No. 99, Nonlinear Mathematical Physics , Vol. 14, No. 3, pp. 462-473, 2007
pp. 119-130, 2009 [24] R. Conte, M. Musette, The Painlevé Handbook, Springer, 2008
[17] M. M. Rashidi, A. Shooshtari, O. Anwar Bég, “Homotopy Perturbation [25] N. A. Kudryashov, “Nonlinear differential equations with exact
Study of Nonlinear Vibration of Von Karman Rectangular Plates”, solutions”, arXiv:nlin/0311058v1[nlin. SI]
Computer & Structures, Vol. 106-107, pp. 46-55, 2012
[26] D. W. Jordan, P. Smith, Nonlinear Ordinary Differential Equations,
[18] J. D. Keckić, “Additions to Kamke's treatise ,VI: A nonlinear second Oxford University Press, 2007
order equation”, Publications de l'Institut Mathématique, nouvelle série,
[27] J.-H. He, “Some asymptotic methods for strongly nonlinear equations”,
Vol. 19, No. 33, pp.81-82, 1975
International Journal of Modern Physics B, Vol. 20, No. 10, pp.1141-
[19] J. D. Keckić, “Additions to Kamke's treatise ,VII: Variation of 1199, 2006
parameters for nonlinear second order differential equations”, Ser. Mat.
[28] M. Pellicer, J. Solà-Morales, “Analysis of a viscoelastic spring-mass
Fiz., Vol. 544-576 , pp. 31-36, 1946, available at:
http://pefmath2.etf.bg.ac.rs/files/100/549.pdf model Journal of Mathematical Analysis and Applications, Vol. 294,
No. 2, pp. 687-698, 2004
[20] P. Painlevé, “Mémoire sur les équations différentielles dont l'intégrale
[29] L. Cveticanin, “Vibrations of the nonlinear oscillator with quadratic
générale est uniforme”, Bull. Soc. Math. France, Vol. 28, pp. 201-261,
1900 nonlinearity”, Physica A: Statistical Mechanics and its Applications,
Vol. 341, pp. 123-135, 2004
[21] L. Roth, “On the solutions of certain differential equations of the second
order”, Phil. Mag. Vol. 32, No. 211, pp.155-164, 1941 [30] R. Jauregni, F. Silva, “Numerical validation methods” in Numerical
Analysis-Theory and Application , IntechOpen, 2011, available at:
[22] M. D. Monsia, “A mathematical model for predicting the nonlinear www.intechopen.com/download/pdf/19916
deformation response of viscoelastic materials”, International Journal of
Applied Mathematics and Mechanics, Vol. 8, No. 16, pp. 22-30, 2012

www.etasr.com Monsia and Kpomahou: Simulating Nonlinear Oscillations of Viscoelastically Damped Mechanical Systems