Non-stationary response of a beam for a new rheological model
OLGA MARTIN
Applied Sciences Faculty
University “Politehnica” of Bucharest
ROMANIA
Abstract: The paper is intended to provide the quasi-static and dynamic analysis of beam with fractional order viscoelastic
material model, which was derived from integer order description using the Boltzmann superposition principle. The results
were obtained for a fractional Zener model by the techniques of Laplace transform and binomial series. An example proves
the accuracy of the solution for a simply-supported beam subjected to a uniform distributed load. Theoretical and numerical
solutions can be easily extending to the complex structures configurations.
Key-words: viscoelastic beam, Euler-Bernoulli beam theory, fractional derivative, Laplace transform, the constitutive law
in hereditary integral form, correspondence principle
1. Introduction
Engineering design of the past thirty years are
frequently present the structures with viscoelastic
components due to their ability to dampen out the
vibrations. The metals at elevated temperatures,
rubbers, polymers that have the characteristic of both
elastic solids as well as viscous solids are examples
of viscoelastic materials. There are remarkable
theoretical studies on these materials of Cristescu
[7], Mainardi and Spada [7], Flugge [10], Freundlich
[10]-[11], Reddy [13], Kennedy [14]. These are
complemented by the works dealing with the
analysis of viscoelastic structures from both
mathematical and engineering points of views: [1],
[8], [16], [17], [20] - [22].
Using the Euler-Bernoulli beam theory, we
present in our paper the governing equation for a
simply supported viscoelastic beam under a uniform
distributed load, [19]. This equation is accompanied
by a constitutive law presented in a hereditary
integral form. Then, extending the procedures of the
classical Zener model with the denomination of
Standard Linear Solid (SLS) to a fractional Zener
model. In order to obtain the quasi-static exact
solution (i.e. the solution ignoring inertia effects),
we will use the correspondence principle for
classical Zener model [15], [22]. This principle
relates mathematically the solution of a linear,
viscoelastic boundary value problem to an analogous
problem of an elastic body of the same geometry and
under the same initial boundary conditions. Mention
that not all problems can be solved by this principle,
but only those for which the boundary conditions do
not vary with the time. Applying the principle of
d’Alembert, this structure will be analyzed in the
dynamic case with a mixed algorithm that is based
on the Galerkin’s method for the spatial domain and
then, the Laplace transform and binomial series
expansion for the time domain.
Numerical results for both quasi-static and
dynamic analysis are presented for classical model
and fractional model, [2], [4], [10], [15]. The first
rheological model will be accompanied of the
comparative studies made with the help of graphical
representations.
2. Beam governing equation
The differential equation of the transverse
oscillations of a beam, which is subjected to
uniformly distributed forces
p
, is obtained from the
dynamic equilibrium of an element having the length
dx. If all the forces acting on the beam element are
projected on the axis Oz, we find:
0
dx
x
T
TdxqTdxpi
(1)
so
x
T
qp i
(2)
where T is the shearing force and qi are the inertial
force. Accordance with the d’Alembert’s principle,
(2) becomes
2
2
t
w
Sp
x
T
(3)
PROOF
DOI: 10.37394/232020.2022.2.2
Olga Martin
E-ISSN: 2732-9941
5
Fig. 1
with the following notations: ρ - density of
material; S – cross sectional area; w - transverse
displacement of the beam and t – time.
Let us now consider the moments in relation to
section x + dx:
0
2
)( 2
dx
pTdxdx
x
M
MM
(4)
and approximating (dx)2 ≈ 0, is obtained
T
x
M
(5)
and (2) becomes the differential equation:
2
2
2
2
t
w
Sp
x
M
(6)
Next, we will study the deformation of the beam
element. The strain on the fiber cd of the deformed
element is tensile (positive) and the strain on the
fiber ab is compressive (negative). Because the
strain is continuous throughout the cross section,
will exist an axis where the strain is zero. This is
named neutral axis and is rs in the Fig. 2.
If line mn there is at a distance z from the neutral
axis, then the strain corresponding to mn is defined
as:
z
z
x
z
d
d
d
dd)(
(7)
where ρ is the curvature radius of the neutral axis
and θ is the angle subtended by the deformed
element.
Fig. 2
The sections ac and bd remain plane and normal
on the deformed axis rs of the beam after
deformation, according to Bernoulli's hypothesis.
In differential geometry, the expression of the
curvature radius is the following:
2
2
2
1
1
x
w
x
w
(8)
Practical applications show that
x
w
is very
small with respect to unity and so we can
consider
2
2
1
x
w
(9)
x
qi
x dx
M T T+dT M+dM
dx
l
z
w
a b
r s
z m n
c d
M
dx
dθ
ρ
a’ b’
r’ s’
m’ n’
c’ d’
M
PROOF
DOI: 10.37394/232020.2022.2.2
Olga Martin
E-ISSN: 2732-9941
6
Introducing (9) in (7), we get
2
2),(
),( x
txw
ztx
(10)
Let us consider a constitutive law in the hereditary
integral form, [21]:
dx
tdG
txGtx t
),(
d
)(
),()0(),(
0
(11)
where G is the relaxation modulus for the beam
material, G(0) = E (elasticity modulus) and σ the
stress corresponding to the strain ε.
Fig. 3
The bending moment M at the beam cross section
x may be expressed in terms of stress in the form
of an integral:

S
dydzztxM ),(
(12)
and using (11) this becomes:
td
x
xw
d
tdG
I
x
txw
IGM
0
2
2
2
2),()(),(
)0(
(13)
where I is the moment of inertia of the section S
with respect to the axis Oy.
Thus, after substitution of M in (6), we obtain the
following integro - differential equation:
pd
x
xw
d
tdG
I
x
txw
IG
t
txw
S
t
0
4
4
4
4
2
2
),()(
),(
)0(
),(
(14)
Solving of equation (13) will lead to the finding of
the transverse displacements for any boundary and
initial conditions given for the viscoelastic beam
subjected to a uniformly distributed loading p.
2. Rheological model
The stretching - relaxation process is analyzed using
a Poynting model that contains two spring elements
and one damping and consists in a serial connection
of a spring and a Kelvin - Voigt model (Fig. 4).
Now, formulating equilibrium and taking the
kinematics of the rheological model into account,
we get the following equalities:
222
11
k
k
21
21
(15)
where k1, k2 are the elastic modulus of the springs
and η is the coefficient of viscosity of the dashpot
Fig. 4
These equations (15) lead to the differential
equation
1
2
1
k
k
k
(16)
or
12121 kkkkk
(17)
where σ and ε depend on the time t.
Let us consider that the material of beam is in its
relaxation phase, so, under constant strain ε = ε0
S1
x
y
z
σ dS dS
z
k2
k1
σ σ
η
PROOF
DOI: 10.37394/232020.2022.2.2
Olga Martin
E-ISSN: 2732-9941
7
the stress will decrease. In this case, the equation
(17) becomes
02121
kkkk
(18)
For the condition: σ (0) = k1 ε0, the solution σ of
equation (18) will be of the form
aa
tt
eke
kk
kk
t
01
21
021 1)(
(19)
where the relaxation time is
21 kk
a
(20)
Using the definition of the relaxation modulus
G(t), we have in this case:
0
)()(
tGt
(21)
From (19), we get
aa
tt
eke
kk
kk
tG
1
21
21 1)(
(22)
that is a function that depends on the beam
material.
It will expand this result to the fractional calculus
using the Mittag - Leffler functions, [15]:
0,10
,
)1(
)/(
)1())/((
0
n
n
n
n
t
tE
(23)
that for ν = 1 reduce to exp (- t /τ).
In the case 0 < ν <1, the differential equation (16)
becomes:
td
d
kk
k
t
kk
kk
td
d
kk
t
21
1
21
21
21
)(
)(
(24)
The relaxation modulus will be now of the form:
))/((1)(
2
1
21
21 a
tE
k
k
kk
kk
tG
(25)
where E ν (0) = 1.
4. Mixed method for solving integro –
differential equation (14)
To solve the equation (14) associated with the
boundary and initial conditions, we express the
transverse deflection w by an expansion of the form
)()(),(
1
xtatxw n
jjjn
(26)
where φj(x) is the j th shape function and aj(t) is the
corresponding time-dependent amplitude. For the
spatial domain will be used the Galerkin’s method
and then, the techniques of the Laplace transform
for the time domain. The shape functions are chosen
to be linearly independent, orthonormal
ji
dxxx ji
L
ji 0
1
)()(
0
(27)
and must satisfy all boundary condition for the
convergence of Galerkin’s method.
Although wn satisfies the boundary conditions, it
generally, does not satisfy equation (14). If the
expansion (26) is substituted into (14) will result
the residual function
n
iii
n
iii
n
xtaIGxtaS
txR
1
)4(
1
)()()0()()(
),(
pdxa
d
tdG
Itn
iii
01
)4( )()(
)(
(28)
Let us now consider the shape functions of the
form
nj
l
xj
l
x
j,....,2,1,sin
2
)(
and
PROOF
DOI: 10.37394/232020.2022.2.2
Olga Martin
E-ISSN: 2732-9941
8
)()()(
4
)4( xx
l
j
xjjjj
The Galerkin’s method requires that the residual to
be orthogonal to each of the chosen shape
functions, so
nitdxdxtxR jn ,...,2,1,0)(),(

(29)
where
],0[],0[ tL
.This leads to n equations
verified by the functions aj(t):
dxxptaGI
da
d
tdG
ItaS
L
jjj
j
t
jj
)()()0(
)(
)(
)(
0
0
(30)
Using the Laplace transform techniques and the
initial conditions
,....2,1,0,)(
)0,(
0
0
kdxx
t
xw
td
ad
j
L
k
k
t
k
j
k
(31)
the functions aj(t) are determined independent of
one another. Finally, an approximate value of the
transverse deflection w(x,t) will be found by (26).
5. Numerical example
Let us consider a simply supported beam under the
uniform distributed load
p
= 4 N/m, which is
applied as a creep load at t = 0 (the load is applied
suddenly at t = 0 and then help constant). The length
of the beam is l = 4 m, width b = 0.08 m and height
h = 0.23 m. These input data lead to the moment of
inertia of the rectangular section:
53 10812/
bhI
m4. The material is taken to
have the density of 1200 kg/m3. For this example,
we employ the three – parameter solid model (SLS
model) with the relaxation modulus expressed by
(19) - (20) [15], [17], where
k1 = 9.8∙107 N/m2 , k2 = 2.45∙107 N/m2 ,
η = 2.74∙108 N-sec/ m2
Classical quasi - static analysis
For ν = 1, the relaxation modulus will be
24.2/77 1084.71096.1)( t
etG
N/m2 (32)
with t in seconds. To get the creep compliance
D(t), will be calculated first, [17]:
t
etD
1)(
1
that corresponds to
24.2,4,1
1096.1
)(
)( 7
1
t
e
tG
tG
Then
2.11
7
18.011096.1)()(
t
etDtD
where
2.111 and 8.0
1
Finally, the creep compliance that corresponds to
the relaxation modulus (32), will be of the form:
2.11
77 10408.01051.0)(
t
etD
(33)
The solving of the problem (30) - (31) will be
accompanied by the appropriate boundary
conditions for the simply supported beams:
0)()0(
0),(),0(
l,tM,tM
tlwtw
(34)
In the quasi-static case the inertial forces are
ignored. An exact solution for the creep loading
applied at t = 0 can be computed using the
correspondence principle. We get the transverse
displacements of the following form:
)()(),( tDxwtxw
(35)
PROOF
DOI: 10.37394/232020.2022.2.2
Olga Martin
E-ISSN: 2732-9941
9
where D(t) is given in (33) and
w
is the solution
for a similar elastic structure that has the modulus
of elasticity E = 1, [19]:
I
LLxxxp
xw 24
)2(
)(
323
(36)
Classical dynamic analysis
For above input data, the equations (30) become:
1)1(
7.2
)(2983
)(4.1065)(8.28
14
0
(
424.2/)
j
j
t
j
t
j
j
taj
daejta
or
1)1(
1.0
)(104
)(37)(
14
0
(
424.2/)
j
j
t
j
t
j
j
taj
daejta
(37)
We remark from the form of (37) that aj(t) are zero
for even values of j. Hence, only odd values of j
need to be considered in the finding of the
Galerkin’s solution. These equations will have
exact solutions determined with the techniques of
Laplace transform. Since the initial conditions on w
and its derivatives are zero, then, and the conditions
on aj and its derivatives are also zero. If Aj(p) is
Laplace transform (
L
) of aj(t) and we use the
Multiplication Theorem, we get
)24.2/1(
)(
)(
0
24.2
s
sA
dae j
j
tL
t
(38)
and (37) becomes:
)43.910445.0(
45.02.0
)( 4423 jjsss
s
j
pAj
(39)
Using the Newton iterative formula, we find a root
of a polynomial that there is into parenthesis: p1 =
0.091. Then, the partial fraction decomposition will
lead to the solution aj(t) of (37). Finally, the
transverse displacements w are approximated by
(26) and for n = 9.
The Fig. 5 shows how the classical dynamic results
oscillate around of the quasi-static results for the
midpoint transverse deflection (x = 2 m). It may be
noted that the amplitude of vibration decreases with
increasing time, due to the presence of the
viscoelastic damping. In general, the amplitude of
the oscillations depends upon the material density,
the rate at which the loading is applied and the
amount of the viscoelastic damping.
Fig. 5
Fractional dynamic model
Using the relaxation modulus (25) in the equation
(30), we get
]1,0(,)()()0(
)(
)(
)(
0
0
dxxptaGI
da
d
tGd
ItaS
L
jjj
j
t
jj
(40)
Appealing to the theory of Laplace transform,
equation (40) becomes:
]1,0(,)(
1
)()0(
)(
)(
)(
0
2
dxxp
s
sAGI
sA
d
tGd
LIsAsS
L
jjj
jjj
where
PROOF
DOI: 10.37394/232020.2022.2.2
Olga Martin
E-ISSN: 2732-9941
10
k
n
k
ksfsFstf
dt
dL
)0()()(
1
)1(
,
F(s) – Laplace transform of f (t) and
4
)/( lj
j
.
According to the paper [15], if the
)0(1
1
1
)]0()(
~
[
)(
2
1
1
2
1
E
k
k
s
s
s
s
s
GssGs
d
tGd
L
a
(41)
where
)(
~sG
- Laplace transformation of G(u),
u = t – τ,
sec, 24.2
sec, 2.11,1096.1
21
2
2
7
21
21
kk
kkk
kk
a
a
,
2
are the retardation time and relaxation
time.
Thus, the Laplace transformation (41)
becomes:
s
ss
s
d
tGd
L24.21
2.1152.111
96.1
)( 1
and from (32), we get
777 108.91084.71096.1)0( G
.
Finally, the equation (40) will be of the form
j
sAj
sA
s
s
jsAs
27.2
)(2989
)(
24.21
2391)(8.28
4
1
42
(42)
So
)4.4645.037104(
45.0
)(
),(
2.0
)(
143443
1
sjsjsjss
ss
sC
sC
j
sA
j
jj
(43)
We notice that for ν = 1 we have (39).
To obtain the original function that
corresponds to (43), we write Cj(s) for
β = 0.45 in the following form:
i
iii
ik
kkk
kk
kkk
kk
kkkk
j
sjs
j
i
ik
sjs
s
s
sjs
j
sjs
s
s
jsjs
sjss
s
jsjs
sjss
jsjs
s
sC
)104(
83
104
)1(
)1(
104
83
1)104(
)1(
)1(
)83104(
)104()1(
)1(
37104
)104(
1)37104(
1
)(
43
4
00
43
0
1
43
4
43
0
1443
43
443
43
443
which was used equality
r
rm
z
r
rm
zr
r
m1
and
1
)1(
0
is binomial coefficient. Still obtain
PROOF
DOI: 10.37394/232020.2022.2.2
Olga Martin
E-ISSN: 2732-9941
11
mmm
m
m
ikiiki
ik
k
ii
kikii
ik
k
j
sj
m
mi
sj
i
ik
s
sjs
sj
i
ik
ssC
24
0
334
00
12433
4
00
104)1(
83)1()1(
)1041(
83
)1()1()(
)(2)(1
)1(104
83
)1(
0323)(4
0
0
sAsA
sj
m
im
i
ki
s
jj
kmikmkimm
m
iik
i
Applying the inverse Laplace transformation
formula:
!
1
1
n
t
sn
nL
we get
)323(
)1(
104
!!
)!(
83
!!
)!(
)(1
)1041(
83
)1()(1
223
)(4
0 0 0
12433
4
00
mik
t
j
im
im
ki
ik
ta
sjs
sj
i
ik
sA
mik
mk
imm
k i m
iik
j
ii
kikii
ik
k
j
Similarly, we find for
)323)1((
)1(
104
!!
)!(
83
!!
)!(
)(2
)1041(
83
)1()(2
223)1(
)(4
0 0 0
1
12433
)1(41
00
mik
t
j
im
im
ki
ik
ta
sjs
sj
i
ik
ssA
mik
mk
imm
k i m
iik
j
ii
kikii
ik
k
j
Finally, the transverse deflection is equal to
)())(2)(1(
2.0
),(
1
xtata
j
txw n
jjjjn
6. Conclusions
In this paper is presented an efficient method for
numerical simulation of a quasi-static and dynamic
response of viscoelastic beam both classical Zener
model and for fractional Zener model. The
proposed method can be easily extended to the
complex viscoelastic structures calculus.
References
[1] J. Argyris, I. Doltsini, V.D. Silva, Constitutive
modelling and computation of non-linear
viscoelastic solids. Part I: Rheological models and
numerical integration
techniques. Comput. Meth. Appl. Mech. Eng, 88
(1991), 135 -163.
[2] R. L. Bagley, P.J. Torvik, Fractional calculus –
A different approach to the analysis of
viscoelastically damped structures, AIAA Journal
219 (1983), 741 – 748.
[3] A.Barari, H.D. Kalijib,M. Ghadimic and G.
Domairry, Non-linear vibration of Euler-Bernoulli
beams, Latin American Journal of Solids and
Structures, 8 (2011), 139 – 148.
[4] M. Caputo, Vibrations of an infinitive
viscoelastic layer with a dissipative memory,
Journal of Acoustical Society of America 56 (1974)
897 – 904.
[5] Q. Chen, Y. W.Chan, Integral finite element
method for dynamical analysis of elastic-
viscoelastic composite structures,” Computers and
Structures, 74 (2000), 51–64.
[6] R. J. Crawford, Plastics Engineering, Elsevier
Butterworth – Heinemann, Oxford, 1998.
[7] N.D. Cristescu, Dynamic plasticity, North
Holland Pub., Amsterdam, 1967.
[8] R.M. Christensen, Theory of Viscoelasticity,
Academic Press, New York, 1982.
PROOF
DOI: 10.37394/232020.2022.2.2
Olga Martin
E-ISSN: 2732-9941
12
[9] W. Flugge, Viscoelasticity, Springer, Berlin,
Heidelberg, 2nd edition, 1975.
[10] J. Freundlich, Vibrations of a simply supported
beam with a fractional viscoelastic material model-
supports movement excitation, Shock and
Vibrations, 20 (2013) 1103 – 1112.
[11] J. Freundlich, M.Pietrzakowski, FEM
simulation of the vibration girder measurement
using piezoelerctric sensors, Machine Dynamics
Research 34 (2010), 28 – 36.
[12] D.W. Van Krevelen. Properties of Polymers,
Elsevier, Amsterdam, 1990.
[13] J.N. Reddy, An Introduction to Continuum
Mechanics with Applications, Cambridge
University Press, 2008.
[14] T. C. Kennedy, Nonlinear viscoelastic analysis
of composite plates and shells. Composite
Structures, 41 (1998), 265–272.
[15] F. Mainardi, G. Spada, Creep, relaxation and
viscosity properties for basic fractional models in
rheology, The European Physical Journal, 193
(2011) 133 – 160.
[16] O. Martin, Propagation of elastic-plastic
waves in bars, Proceedings of the European
Computing Conference, Springer Science, Lecture
Notes in Electrical Engineering, 28 (2009), Part II,
cap.10, 1011-1024.
[17] O. Martin, Quasi-static and dynamic analysis
for viscoelastic beams with the constitutive
equation in a hereditary integral form, Annals of the
University of Bucharest), 5 (2014), 1–13. in a
[18] S.W. Park, Y.R. Kim, Interconversion between
relaxation modulus and creep compliance for
viscoelastic solids, J. mater. Civil Eng., 11 (1999),
76 – 82
[19] G. S. Payette, J. N. Reddy, Nonlinear quasi-
static finite element formulations for viscoelastic
Euler–Bernoulli and Timoshenko beams,
International Journal for Numerical Methods in
Biomedical Engineering, 26 (2010),1736–1755.
[20] G. Pissarenko, A. Yakovlev, V. Matveev,
Aide-mémoire de résistance des matériaux,
Editions de Moscou, 1979.
[21] B.S. Sarma and T.K. Varadan, Lagrange-type
formulation for finite element analysis of nonlinear
beam vibrations, J. Sound Vibrations., 86 (1983),
61–70.
[22] S. Shaw, M. K. Warby, J. R. Whiteman, A
comparison of hereditary integral and internal
variable approaches to numerical linear solid
viscoelasticity, In Proc. XIII Polish Conf. on
Computer Methods in Mechanics, Poznan, (1997).
Creative Commons Attribution License 4.0
(Attribution 4.0 International, CC BY 4.0)
This article is published under the terms of the Creative
Commons Attribution License 4.0
https://creativecommons.org/licenses/by/4.0/deed.en_US
PROOF
DOI: 10.37394/232020.2022.2.2
Olga Martin
E-ISSN: 2732-9941
13