Acta Polytechnica


DOI:10.14311/AP.2020.60.0502
Acta Polytechnica 60(6):502–511, 2020 © Czech Technical University in Prague, 2020

available online at https://ojs.cvut.cz/ojs/index.php/ap

NEWMARK ALGORITHM FOR DYNAMIC ANALYSIS WITH
MAXWELL CHAIN MODEL

Jaroslav Schmidt∗, Tomáš Janda, Alena Zemanová, Jan Zeman,
Michal Šejnoha

Czech Technical University in Prague, Faculty of Civil Engineering, Department of Mechanics, Thákurova 7,
166 29 Praha, Czech Republic

∗ corresponding author: Jaroslav.Schmidt@cvut.cz

Abstract. This paper investigates a time-stepping procedure of the Newmark type for dynamic
analyses of viscoelastic structures characterized by a generalized Maxwell model. We depart from a
scheme developed for a three-parameter model by Hatada et al. [1], which we extend to a generic
Maxwell chain and demonstrate that the resulting algorithm can be derived from a suitably discretized
Hamilton variational principle. This variational structure manifests itself in an excellent stability and
a low artificial damping of the integrator, as we confirm with a mass-spring-dashpot example. After
a straightforward generalization to distributed systems, the integrator may find use in, e.g., fracture
simulations of laminated glass units, once combined with variationally-based fracture models.

Keywords: Newmark method, Maxwell chain model, variational integrators.

1. Introduction
The motivation of this work comes from the field of dy-
namics of laminated glass structures. These sandwich
structures consist of multiple glass layers connected
with transparent polymer interlayers. Combining stiff,
brittle glass with compliant viscoelastic polymers im-
proves structural safety, but the through-thickness
heterogeneity renders mechanics of laminated glass
structures intricate, e.g. [2]. In particular, time- and
temperature-dependent interlayer properties must be
accounted for even in quasi-static analyses, e.g., [3–6]
and references therein.

Earlier studies [7–10] have shown that the response
of commonly used interlayer materials can be captured
well with the Maxwell chain model combined with the
time-temperature superposition principle. Because
the viscoelastic model concurrently predicts material
damping, vibrations of laminated glass structures can
be described more accurately than in conventional
structural analyses that mostly employ the Rayleigh
damping, e.g. [11, Section 12.5]. This added value has
been addressed in detail for free vibrations of lami-
nated glass units, e.g. [12–14]; an extension towards
the response under general dynamic loads requires the
development of dedicated time-stepping schemes that
are in the focus of the current work.

Related work. Dynamics of viscoelastic solids de-
scribed by the Maxwell chain model leads to the sys-
tem of initial value problems coupling the equation of
motion with the local evolution of constitutive vari-
ables, see Section 2.1 for illustration. Because numer-
ically integrating the full system would be costly, we
will follow an alternative approach in which only the
equations of motion are solved approximately, whereas
the evolutionary constitutive equations are resolved in

the closed form, leading to an inexpensive update for-
mulas for internal variables entering the equations of
motion. This approach has been pioneered for quasi-
static problems by Zienkiewicz et al. [15]; see also [16,
Section 5.2] for a comprehensive review. To the best
of our knowledge, Hatada et al. [1] were the only ones
who used this strategy in dynamics, although no refer-
ence to the original work [15] was made. In particular,
they developed a Newmark-type [17] algorithm for
the three-parameter Maxwell model and used it to
predict the response of planar frames to earthquake
loading.

Novelty. Our work further develops the contribu-
tion [1] in three aspects. First, in Section 2.2, we
present a compact derivation of the Newmark scheme
for a generic Maxwell chain, closely following the orig-
inal exposition [15]. Second, in Section 3, we show
that the algorithm can be interpreted as a variational
integrator [18], in the sense that it can be derived from
the Hamilton variational principle combined with a
suitable time discretization. The variational structure
provides the integrator with a good numerical stability
and low numerical dissipation, as demonstrated in Sec-
tion 4.1 with selected examples. Moreover, the scheme
can easily be combined with variational approaches to
fracture, e.g., [19–21], which is of independent inter-
est when simulating the behaviour of laminated glass
under impact loads. Third, in Section 4.2, we outline
how to extend the algorithm to a continuum formu-
lation and complement the theoretical considerations
with an illustrative 3D finite element simulation.

Notation. We employ the conventional notation
through the text, in which scalar quantities are de-
noted by a plain font, whereas bold-face letters in-

502

https://doi.org/10.14311/AP.2020.60.0502
https://ojs.cvut.cz/ojs/index.php/ap


vol. 60 no. 6/2020 Newmark algorithm for dynamic analysis of viscoelastic materials

k∞

ηPη3η2η1

k1 k2 k3 kP

m

F(t)

r(t)
re,1(t)

rv,1(t)

Figure 1. Scheme of the single-degree-of-freedom
viscoelastic dynamic problem.

dicate vectors or higher-order tensors. Additional
nomenclature is introduced when needed.

2. Newmark method
In this section, we analyse a single degree of freedom
(SDOF) model of a mass supported with a Maxwell
chain, consisting of the parallel connection of an elastic
spring and multiple spring-dashpot cells, see Figure 1
for illustration and, e.g., [16, Section A] for further
details. In particular, in Section 2.1, we review the
equations of motion, which we subsequently discretize
with the average acceleration version of the Newmark
method [17] in Section 2.2.

2.1. Governing equations
As follows from the scheme in Figure 1, the prob-
lem under consideration is specified with the time-
dependent load F(t), the particle mass m and the
Maxwell chain model parameters: stiffness of the elas-
tic spring k∞, spring stiffness kp and damper viscosity
ηp of the p-th Maxwell cell; P stands for the number
of Maxwell cells.

Equilibrium of the forces acting on the mass requires

mr̈(t) + k∞r(t) +
P∑
p=1

fp(t) = F(t), (1)

where r denotes the displacement of the mass, r̈ its
acceleration, and fp the restoring force of the p-th
cell.

For the p-th Maxwell cell, the displacement r splits
into an elastic part of the spring re,p and a viscous
part of the damper rv,p:

r(t) = re,p(t) + rv,p(t); (2)

recall Figure 1. The restoring force of the p-th Maxwell
cell satisfies

fp(t) = kpre,p(t) = ηpṙv,p(t), (3)

because of the serial arrangement of the spring and the
damper in the cell. Differentiating (2) with respect to
time and using (3), we obtain

ḟp(t)
kp

+
fp(t)
ηp

= ṙ(t). (4)

In summary, the motion of the SDOF model is de-
scribed with the coupled system of (P + 1) ordinary
differential equations (ODEs) (1) and (4), comple-
mented with the initial conditions

r(0) = r0, ṙ(0) = v0, fp(0) = fp,0, (5)

where r0 and v0 stand for the initial mass displacement
and velocity, respectively, fp,0 is the initial force in
the p-th Maxwell cell, and p = 1, . . . ,P .

2.2. Discretization
The time interval of interest 〈0,T〉 is divided into
(N + 1) time instants 0 = t0 < t1 < t2 < · · · <
tN−1 < tN = tmax; for notational simplicity, we as-
sume equidistant partitioning of the constant time
step ∆t = ti+1 − ti.
Considering the Newmark integration scheme [17]

with an average acceleration (r̈i + r̈i+1)/2 on the time
interval 〈ti, ti+1〉, the velocity and displacement within
the interval vary as

ṙ(ti + τ) = ṙi + 12 (r̈i + r̈i+1) τ, (6a)
r(ti + τ) = ri + ṙiτ + 14 (r̈ + r̈i+1) τ

2, (6b)

where τ ∈ 〈0, ∆t〉 is the local time variable within
the interval and •i abbreviates •(ti) to render the
notation compact.

Substituting (6a) into (4) reveals that the evolution
of the restoring force fp satisfies

ḟp(ti + τ)
kp

+
fp(ti + τ)

ηp
= ṙi + 12 (r̈i + r̈i+1)τ. (7)

This Cauchy problem with the initial condition
fp(ti) = fp,i has the solution

fp(ti + τ) =(
fp,i −ηpṙi +

η2p
2kp

(r̈i + r̈i+1)
)

exp
(
−
kp
ηp
τ

)
+ηpṙi +

ηpτkp −η2p
2kp

(r̈i + r̈i+1). (8)

Thus, at the end of the time interval with τ = ∆t, we
have

fp,i+1 = Apfp,i + kpθ̂pṙi + Bp (r̈i + r̈i+1) , (9)

where

θ̂p =
ηp
kp

(
1 − exp

(
−
kp
ηp

∆t
))

(10)

denotes the effective relaxation time of the p-th cell
and the auxiliary factors are given by

Ap =
(

1 −
kp
ηp
θ̂p

)
, Bp = 12ηp

(
∆t− θ̂p

)
. (11)

503



J. Schmidt, T. Janda, A. Zemanová et al. Acta Polytechnica

Finally, substituting equations (6b) and (9) into (1)
expressed at ti+1 = ti + ∆t and rearranging the terms
yields

(
m + 14k∞∆t

2 +
P∑
p=1

Bp

)
r̈i+1 = Fi+1 −

P∑
p=1

Apfp,i

−k∞ri −
(
k∞∆t +

P∑
p=1

kpθ̂p

)
ṙi

−
(

1
4k∞∆t

2 +
P∑
p=1

Bp

)
r̈i. (12)

After solving Eq. (12) for the acceleration r̈i+1, we
update velocity ṙi+1, displacement ri+1, and restoring
forces fp,i according to equations (6a), (6b), and (9),
respectively, and proceed to the next time interval.

Note that the initial acceleration r̈0, needed in the
first step of the algorithm, is set to

r̈0 =
1
m

(
F(0) −k∞r0 −

P∑
p=1

fp,0
)
, (13)

according to the equilibrium (1) and initial (5) condi-
tions.

3. Variational integrators
Having derived the Newmark viscoelastic algorithm
by conventional means, we now demonstrate its varia-
tional structure, by adapting the general arguments
on variational integrators by Kane et al. [18] to the
current setting.
We will proceed in four steps. In Section 3.1, we

show that the governing equations from Section 2.1 fol-
low from the Euler-Lagrange (E-L) equations of a suit-
ably defined energy functional. Its discretization then
provides the governing equations of the corresponding
variational integrator introduced in Section 3.2. In
Section 3.3, we demonstrate the equivalence of the
integrator to the Newmark algorithm from Section 2.2.
In the last step, Section 3.4, we comment on the en-
ergy conservation properties of the time integration
scheme.

3.1. Variational framework
We postulate that the trajectory q : (0,T) → Q of a
constrained dissipative mechanical system in the state
space Q is given by the generalised Euler-Lagrange
equations1, e.g., [22, Section 1.3]

∂R(q̇(t))
∂q̇

=
∂L(t,q(t), q̇(t),λ(t))

∂q

−
∂

∂t

∂L(t,q(t), q̇(t),λ(t))
∂q̇

, (14a)

0 = α(q(t)). (14b)
1The term generalised corresponds to the appearance of the

non-conservative forces in Eq. (14a).

Here, R stands for the dissipation potential, L for
the Lagrangian of the problem, and λ : (0,T) →
Λ denotes the Lagrange multipliers associated with
the kinematic constraint function α. Besides, these
equations correspond to the stationarity conditions of
the action functional

S(q̂,λ̂) =
∫ T

0
L
(
t, q̂(t), ̂̇q(t),λ̂(t)) dt (15)

perturbed by the dissipative forces ∂q̇R. Note that
the hat symbol in (15) now distinguishes the test
quantities from the true trajectories defined with (14).
For the problem from Figure 1, the state variable

q̂(t) =
[
r̂(t), {r̂(t)e,p}Pp=1, {r̂(t)v,p}Pp=1

] T (16)
collects the total displacement and the displacements
of both components of each Maxwell cell; the state
space Q = R2P+1. The Lagrangian has the standard
form

L(t, q̂, ̂̇q,λ̂) = K(̂̇q) −E(q̂) + fext(t)Tq̂
+ λ̂Tα(q̂) (17)

involving the kinetic energy K, potential energy of
deformation E, and external forces fext given by

K(̂̇q) = 12m(̂̇r)2, (18a)
E(q̂) = 12k∞r̂

2 + 12
P∑
p=1

kp(r̂e,p)2, (18b)

fext(t) =
[
F(t), 01×2P

] T. (18c)
The kinematical constraints take the form

αp(q̂) = r̂e,p + r̂v,p − r̂, p = 1, . . . ,P ; (19)

the space of the Lagrange multipliers Λ then becomes
RP . The last component of the general framework (14)
is provided by the dissipation potential

R(̂̇q) = 12 P∑
p=1

ηp(̂̇rv,p)2 (20)
involving solely the viscous displacements of all cells.
In this setting, the E-L equation (14a) represents

the system of (1 + 2P) optimality conditions. The
first one, corresponding to the total displacement r,
takes the form

mr̈(t) + k∞r(t) +
P∑
p=1

λp(t) = F(t), (21)

while the remaining 2P conditions read as

λp(t) = kpre,p(t), λp(t) = ηpṙv,p(t), (22)

with p = 1, . . . ,P . It is thus evident that the multipli-
ers λp play a role of the viscous force fp and, because
the optimality (14b) and compatibility (4) conditions
coincide, the current setting is equivalent to the one
of Section 2.1.

504



vol. 60 no. 6/2020 Newmark algorithm for dynamic analysis of viscoelastic materials

3.2. Discretization
Recall that the incremental algorithm of Section 2.2
relies on the discretization of the total displacements,
from which the evolution of cell-related variables re,p,
rv,p, and fp follows in the closed form. To mimic
this structure, only the total displacements r will
be determined from the discrete (non-dissipative) E-
L equations, whereas the remaining quantities are
determined from the non-discretized optimality condi-
tions (22) and (14b).

To this goal, we consider the same discretization of
the time interval 〈0,T〉 as in Section 2.2 and introduce
the discretized action functional

S(r̂,λ̂) ≈Sd
(
{r̂i}Ni=0,{λ̂i}

N
i=0
)

= ∆t
N−1∑
i=0
Ld
(
r̂i, r̂i+1,λ̂i,λ̂i+1

)
, (23)

with the discrete Lagrangian given by [18, Eq. (2)]

Ld
(
r̂i, r̂i+1,λ̂i,λ̂i+1

)
= 12m

( r̂i+1 − r̂i
∆t

)2 (24)
− 12k∞

( r̂i+1 + r̂i
2

)2
+ 14 (r̂i + r̂i+1)(Fi + Fi+1)

− 14
P∑
p=1

(r̂i + r̂i+1)(λ̂p,i + λ̂p,i+1).

The stationarity conditions at time ti, ∂Sd/∂ri = 0,
with i = 1, . . . ,N − 1,2 read as

0 =
∂Ld(ri−1,ri,λi−1,λi)

∂ri

+
∂Ld(ri,ri+1,λi,λi+1)

∂ri
(25)

which delivers the governing equations of the varia-
tional integrator in the form3

m
ri+1 − 2ri + ri−1

∆t2
+ 14k∞(ri+1 + 2ri + ri−1)

+
P∑
p=1

1
4 (λp,i+1 + 2λp,i + λp,i−1)

= 14 (Fi−1 + 2Fi + Fi+1). (26)

3.3. Equivalence to Newmark
We will proceed with two additional steps to show
that the optimality conditions (26) correspond to the
Newmark integration scheme from Section 2. First,

2Notice that the optimality conditions exclude the initial
and final times t0 and tN . This mirrors the standard treatment
of the time-continuous case, in which the stationarity conditions
of the action functional are derived under the assumption that
the state variables remain fixed at both times, e.g. [22].

3Notice that we assume the Lagrange multipliers λp,i to be
given arbitrary quantities, similarly to the forcing terms Fi.
Once we establish the equivalence to the Newmark algorithm,
their values follow from the update formula (9) from Section 2.

we demonstrate that the displacements {ri}Ni=0 pro-
vide definitions of velocities {ṙi}Ni=0 and accelerations
{r̈i}Ni=0 consistent with the kinematic assumptions in
Eq. (6). Second, we show that the discrete-in-time
quantities satisfy the equations of motion (1).

Kinematics. Following [18, Section 2.2], we start
by introducing auxiliary accelerations

mr̈i+1/2 = −
k∞
2

(ri + ri+1) + 12 (Fi + Fi+1)

−
P∑
p=1

1
2 (λp,i + λp,i+1), (27)

for i = 0, 1, . . . ,N − 1. Summing mr̈i−1/2 with
mr̈i+1/2 and comparing the result with (26) provides

ri+1 − 2ri + ri−1
∆t2

= 12 (r̈i+1/2 + r̈i−1/2), (28)

with i = 1, . . . ,N − 1.
The discrete linear momenta follow standardly from

pi =
∂Ld(ri−1,ri,λi−1,λi)

∂ri
∆t, (29)

and, using r̈i−1/2 from Eq. (27), they can be evaluated
as

pi = mṙi = m
ri −ri−1

∆t
+ m

r̈i−1/2

2
∆t. (30)

Expressing pi+1 according to the previous relation
and employing (28) provides

pi+1 = pi + mr̈i+1/2∆t, (31)

from which we obtain

ṙi+1 = ṙi + ∆tr̈i+1/2. (32)

Likewise, expressing ri from (28) and employing
the velocity ṙi from (30) provides

ri+1 = ri + ṙi∆t + 12 r̈i+1/2∆t
2. (33)

Hence, expressions (33) and (32) become identical to
the ones of the Newmark method (6) once setting

r̈i+1/2 = 12 (r̈i + r̈i+1). (34)

Equilibrium. Employing the nodal accelerations
r̈i−1/2 and r̈i+1/2 from (34) in the identity (28) reveals
that

ri+1 − 2ri + ri−1
∆t2

= 14 (r̈i−1 + 2r̈i + r̈i+1). (35)

Furthermore, by expressing the difference m(r̈i+1/2 −
r̈i−1/2) using (27), we find that

1
2m(r̈i+1 − r̈i−1) +

1
2k∞(ri+1 −ri−1)

+
P∑
p=1

1
2 (λp,i+1 −λp,i−1) =

1
2 (Fi+1 −Fi−1). (36)

505



J. Schmidt, T. Janda, A. Zemanová et al. Acta Polytechnica

Now, after inserting the identity (35) into the discrete
Euler-Lagrange equations (26) and subtracting (36)
from the result, we infer that

m(r̈i−1 + r̈i) + k∞(ri−1 + ri) +
P∑
p=1

(λp,i−1 + λp,i)

=Fi−1 + Fi, (37)

which can be reduced to the final form

mr̈i + k∞ri +
P∑
p=1

λp,i = Fi. (38)

Indeed, the equivalence between (38) and (37) for
i = 1 holds because of the choice of the initial accel-
eration (13), and for i = 2, . . . ,N − 1 it follows by
induction.

3.4. Energy balance
The variational framework (14) additionally reveals
that the trajectory q satisfies the energy balance con-
dition, e.g., [23, Section 5.1]

Eint(t) + D(t) = Eint(0) + W(t) for 0 ≤ t ≤ T, (39)

with the internal energy Eint, dissipated energy D, and
the work done by external forces W given by

Eint(t) = 12mṙ
2(t) + 12k∞r

2(t)

+
P∑
p=1

1
2kpr

2
e,p(t), (40a)

D(t) =
P∑
p=1

∫ t
0
ηpṙ

2
v,p(τ)dτ, (40b)

W(t) =
∫ t

0
F(τ)ṙ(τ)dτ. (40c)

To later quantify the artificial dissipation induced by
time discretization, we also consider the time-discrete
quantities

Eint(ti) = 12mṙ
2
i +

1
2k∞r

2
i +

P∑
p=1

λ2p,i
2kp

, (41a)

Dd(ti) =
i−1∑
k=0

P∑
p=1

1
2ηp
(
λ2p,k + λ

2
p,k+1

)
∆t, (41b)

Wd(ti) =
i−1∑
k=0

1
2
(
Fkṙk + Fk+1ṙk+1

)
∆t; (41c)

the last two expressions correspond to the approxima-
tions of integrals in (40) with the trapezoidal rule and
employing the identities (22).

4. Examples
In this section, we demonstrate the performance of
the developed Newmark algorithm on two examples.

kp [kNm−1] θp [s] kp [kNm−1] θp [s]
6933.9 10−9 445.1 102
3898.6 10−8 300.1 103
2289.2 10−7 401.60 104
1672.7 10−6 348.1 105
761.60 10−5 111.6 106
2401.0 10−4 127.2 107
65.200 10−3 137.8 108
248.00 10−2 50.5 109
575.60 10−1 322.9 1010
56.30 100 100.0 1011
188.6 101 199.9 1012

Table 1. Parameters of Maxwell chain model [10],
with θp = ηp/kp and k∞ = 682.18 kNm−1. Note that
in Section 4.2, the stiffnesses k• correspond to shear
moduli G• [MPa].

The first one in Section 4.1 addresses the accuracy and
numerical energy dissipation of the integrator for the
single-degree-of-freedom system from Figure 1. The
follow-up example in Section 4.2 outlines an extension
of the scheme towards continuum models.
Data of the generalized Maxwell chain used in

both examples appear in Table 1; they represent real
PolyVinyl Butyral (PVB) material with sufficiently
short and long relaxation times for testing the algo-
rithm robustness. For more information on experimen-
tal procedures to determine these parameters, see [10].
All results presented in this section are reproducible
with Python-based scripts available at [24].

4.1. Discrete problem
We consider the following two types of loading:

F(t) = F for t ≥ 0, (42a)
F(t) = F sin t for t ≥ 0, (42b)

corresponding to step and harmonic loads, respectively.
In both cases, we set the amplitude F = 1 MN and the
mass m = 106 kg to scale the displacement amplitude
to ≈ 1 m. Initial displacement, velocity, and forces in
Maxwell cells were set to zero; recall (5). As for the
Newmark algorithm, we set the time steps ∆t to 1.0,
0.5, and 0.2 s.

Accuracy of the Newmark algorithm is checked by
comparing its trajectories with the reference ones,
obtained with the adaptive solver lsoda [25] — avail-
able through odeint function of Scipy library [26]
— applied to the full initial value problem (1), (4),
and (5).

The results appear in Figure 3 and demonstrate that
the Newmark algorithm is stable even for coarse time
steps. The errors behave consistently with findings for
Newmark-family methods applied to linearly damped
systems, e.g. [27, Section B.II.5]. In particular, the
numerical dispersion (understood as the error in peri-
ods) and dissipation (error in amplitudes) decays as

506



vol. 60 no. 6/2020 Newmark algorithm for dynamic analysis of viscoelastic materials

10-2 10-1 100
Time step ∆t [s]

10-6

10-5

10-4

10-3

10-2

10-1

100
Re

la
tiv

e 
er

ro
r [

-]

Numerical errors
Quadratic order

Figure 2. Convergence diagram for the step loading.
Crosses represent the relative numerical errors for dif-
ferent time steps, the solid line a quadratic function.

O((ω∆t)2), where ω stands for the eigenfrequency of
the undamped system. For ∆t = 0.2 s, the trajectories
predicted by the Newmark scheme closely match the
reference ones.

To demonstrate the second-order accuracy qualita-
tively, in Figure 2, we plot the decay of a relative error
with a decreasing time step ∆t. The relative errors
were calculated according to ||r∆t − rref∆t||2/||r

ref
∆t||2,

where r∆t stores displacements of all time instants for
time step ∆t, rref∆t the reference solution of odeint at
the same times, and ‖•‖ denotes the Euclidean norm.
A comparison with the quadratic function (solid line)
confirms the second-order accuracy. Note that the
deviation visible for short time steps results from a
prescribed tolerance of the odeint solver.

Numerical dissipation. As follows from the en-
ergy equality (39), the additional dissipation induced
by the integrator can be estimated as

∆d(ti) = |Eint(0) + Wd(ti) −Eint(ti) −Dd(ti)| , (43)

with the individual terms provided by Eq. (41). The
evolution of these quantities for the step and harmonic
loading appears in Figure 4, considering the time
interval 〈0, 300〉 s.
For both loads, we observe that the work done

by external forces eventually distributes between the
internal and dissipated energies; the ratio Dd/Wd
stabilizes at 0.4 for the step load and for the harmonic
loading, the ratio reaches about 0.9. The artificial
dissipation is only significant for the coarsest step of
∆t = 1.0 s; for ∆t = 0.1 s it reaches the value of
about 1 ‰ and further deteriorates with a decreasing
time step. This confirms excellent energy conservation
properties of the scheme, especially when taking into
account that the error introduced by the trapezoidal
rule in (41) is of order O(∆t2).

4.2. Generalization
Additional constitutive assumptions must be adopted
to extend the SDOF models into a continuum for-
mulation. Here, we assume that the Maxwell model
applies when modelling the response under shear. The
spring stiffnesses k• thus become shear moduli G•,
and the Poisson ratio ν is a time-independent mate-
rial constant. This assumption considerably simpli-
fies the multi-dimensional constitutive law, e.g., [16,
Section 2.4], and provides the same results as the
conventional volumetric-deviatoric split for our target
applications [5].
Under these assumptions, the weak form of the

equations of motion take the form, e.g. [11, Part III]:

δFext(t,δu) =
∫

Ω
δu ·ρü(t) dΩ

+
∫

Ω
δε :

(
G∞Dν : ε(t) +

P∑
p=1

σp(t)
)

dΩ, (44)

where δFext stands for the virtual work done by exter-
nal loads on a virtual displacement δu, ü denotes the
acceleration, and the small strain tensor ε is obtained
as the symmetric part of the displacement gradient,
ε = ∇su; virtual strain δε is defined in the same
way. The material is characterized by its density ρ,
long-term shear modulus of the Maxwell chain G∞,
and the dimensionless tensor Dν corresponding to
the stiffness tensor of an isotropic material of the unit
shear modulus and the Poisson ratio of ν. The stresses
σp carried by individual cells follow as the solution of
initial value problems

σ̇p(t)
Gp

+
σp(t)
ηp

= Dνε̇(t), p = 1, 2, . . . ,P. (45)

The evolution of the state variables u and σp is speci-
fied with the initial conditions on the displacements,
velocities, and cell stresses:

u(0) = u0, u̇(0) = v0, σp(0) = σp,0. (46)

The comparison of the initial value problems spec-
ified with Eqs. (1), (4), and (5) and Eqs. (44), (45),
and (46) reveals that the derivation of the Newmark-
type scheme follows exactly the steps as in Section 2.2.
As a result, the following variational problem needs
to be solved at a time ti+1:∫

Ω
δu ·ρüi+1 dΩ

+
∫

Ω
δε :

(
1
4G∞∆t

2 +
P∑
p=1

Bp

)
Dν : ε̈i+1 dΩ

= δFext(ti+1,δu) −
∫

Ω
δε :

( P∑
p=1

Apσp,i

)
dΩ

−
∫

Ω
δε : G∞Dν : εi dΩ

507



J. Schmidt, T. Janda, A. Zemanová et al. Acta Polytechnica

0 5 10 15 20 25 30
time t [s]

0.0

0.1

0.2

0.3

0.4

0.5

0.6
di

sp
la

ce
m

en
t r

(t
) [

m
]

Reference
Newmark

0 5 10 15 20 25 30
time t [s]

0.8

0.6

0.4

0.2

0.0

0.2

0.4

0.6

di
sp

la
ce

m
en

t r
(t

) [
m

]

Reference
Newmark

0 5 10 15 20 25 30
time t [s]

0.0

0.1

0.2

0.3

0.4

0.5

0.6

di
sp

la
ce

m
en

t r
(t

) [
m

]

Reference
Newmark

0 5 10 15 20 25 30
time t [s]

0.8

0.6

0.4

0.2

0.0

0.2

0.4

0.6

di
sp

la
ce

m
en

t r
(t

) [
m

]

Reference
Newmark

0 5 10 15 20 25 30
time t [s]

0.0

0.1

0.2

0.3

0.4

0.5

0.6

di
sp

la
ce

m
en

t r
(t

) [
m

]

Reference
Newmark

0 5 10 15 20 25 30
time t [s]

0.6

0.4

0.2

0.0

0.2

0.4

0.6

di
sp

la
ce

m
en

t r
(t

) [
m

]

Reference
Newmark

Figure 3. Accuracy of the viscous Newmark method for step (left) and harmonic (right) loadings defined by Eq. (42).
Top, center, and bottom graphs show trajectories for time steps ∆t = 1.0 s, ∆t = 0.5 s, and ∆t = 0.2 s respectively.

−
∫

Ω
δε :

(
G∞∆t +

P∑
p=1

Gpθ̂p

)
Dν : ε̇i dΩ

−
∫

Ω
δε :

(
1
4G∞∆t

2 +
P∑
p=1

Bp

)
Dν : ε̈i dΩ, (47)

with the parameters θ̂p, Ap, and Bp provided by
Eqs. (10) and (11); recall that δFext denotes the vir-
tual work done by external forces. Once the the accel-
erations üi+1 are obtained from the weak form (47),
the displacements ui+1, velocities u̇i+1, and the cell

stresses σp,i+1 are updated according to Eqs. (6)
and (9), respectively.
The outlined formulation (47) was further dis-

cretized with the finite element method and imple-
mented in FEniCS project [28, 29] version 2018.1. As
an indicative example, we consider a unit cube, see
Figure 5, fixed on the bottom surface and subjected to
a step load (42a) with the tensile traction of intensity
1.0 Nm−2 perpendicular to the top surface. The ma-
terial response is characterized by the Maxwell chain
parameters from Table 1 and the value of the Pois-
son ratio ν = 0.49. In the numerical resolution, we

508



vol. 60 no. 6/2020 Newmark algorithm for dynamic analysis of viscoelastic materials

0 50 100 150 200 250 300
time t [s]

0.0

0.2

0.4

0.6

0.8

1.0
Re

la
tiv

e 
en

er
gy

 [-
]

Eint/Wd
Dd/Wd

0 50 100 150 200 250 300
time t [s]

10-9
10-8
10-7
10-6
10-5
10-4
10-3
10-2
10-1
100

Re
la

tiv
e 

di
ss

ip
at

io
n 

∆
d
/
W

d
 [-

]

∆t = 1.00

∆t = 0.10

∆t = 0.01

0 50 100 150 200 250 300
time t [s]

0.0

0.2

0.4

0.6

0.8

1.0

Re
la

tiv
e 

en
er

gy
 [-

]

Eint/Wd
Dd/Wd

0 50 100 150 200 250 300
time t [s]

10-5

10-4

10-3

10-2

10-1

100

Re
la

tiv
e 

di
ss

ip
at

io
n 

∆
d
/
W

d
 [-

]

∆t = 1.00

∆t = 0.10

∆t = 0.01

Figure 4. Normalized energies corresponding to SDOF response to step (top) and harmonic loads (bottom). Left:
the evolution of internal energy Eint and dissipation Dd, normalized by the work done by external forces Wd for
time step ∆t = 0.5 s. Right: the evolution of numerical dissipation ∆d, normalized by the work done by external
forces Wd.

discretize the sample into identical 1,000 hexahedron
elements and consider the time step of 0.01 s.
The snapshots of the vibrations reveal a similar

behaviour to the SDOF example, recall Figure 3,
namely the attenuation of the propagating waves by
viscous damping. An interested reader is invited to
the dataset [24] for full details on the simulation.

5. Conclusions
In this contribution, we have developed a Newmark
integration scheme for viscoelastic solids character-
ized by the generalized Maxwell model. Besides the
direct derivation, we have shown the scheme can be
derived from the Hamilton variational principle com-
bined with a suitable structure-preserving time dis-
cretization. This variational structure is then reflected
in low numerical energy dissipation of the resulting
scheme, which has been confirmed with selected nu-
merical examples.
As the next step, we will combine the continuum

framework outlined in Section 4.2 with Newmark-
type solvers for variational fracture models, e.g. [30,
31], to extend the currently available approaches to

simulating the response of laminated glass structures
under impact.

Acknowledgements
This publication was supported by the Czech Science
Foundation, the grant No. 19-15326S. We wish to thank
an anonymous reviewer for her/his constructive criticism
that helped us clarify the original manuscript at several
places.

References
[1] T. Hatada, T. Kobori, M. Ishida, N. Niwa. Dynamic
analysis of structures with Maxwell model. Earthquake
Engineering & Structural Dynamics 29(2):159–176,
2000. doi:10.1002/(SICI)1096-
9845(200002)29:2<159::AID-EQE895>3.0.CO;2-1.

[2] M. Haldimann, A. Luible, M. Overend. Structural Use
of Glass, vol. 10 of Structural Engineering Documents.
IABSE, Zürich, Switzerland, 2008.

[3] A. van Duser, A. Jagota, S. J. Bennison. Analysis of
glass/Polyvinyl Butyral laminates subjected to uniform
pressure. Journal of Engineering Mechanics
125(4):435–442, 1999.
doi:10.1061/(ASCE)0733-9399(1999)125:4(435).

509

http://dx.doi.org/10.1002/(SICI)1096-9845(200002)29:2<159::AID-EQE895>3.0.CO;2-1
http://dx.doi.org/10.1002/(SICI)1096-9845(200002)29:2<159::AID-EQE895>3.0.CO;2-1
http://dx.doi.org/10.1061/(ASCE)0733-9399(1999)125:4(435)


J. Schmidt, T. Janda, A. Zemanová et al. Acta Polytechnica

time t [s]

di
sp
la
ce
m
en
t
r
(t

)[
m

]

Figure 5. Snapshots of deformations of a viscoelastic cube subjected to step load.

[4] L. Galuppi, G. Royer-Carfagni. The design of
laminated glass under time-dependent loading.
International Journal of Mechanical Sciences 68:67–75,
2013. doi:10.1016/j.ijmecsci.2012.12.019.

[5] A. Zemanová, J. Zeman, M. Šejnoha. Comparison of
viscoelastic finite element models for laminated glass
beams. International Journal of Mechanical Sciences
131-132:380–395, 2017.
doi:10.1016/j.ijmecsci.2017.05.035.

[6] A. Zemanová, J. Zeman, T. Janda, M. Šejnoha.
Layer-wise numerical model for laminated glass plates
with viscoelastic interlayer. Structural Engineering and
Mechanics 65(4):369–380, 2018.
doi:10.12989/sem.2018.65.4.369.

[7] L. Andreozzi, S. B. Bati, M. Fagone, et al. Dynamic
torsion tests to characterize the thermo-viscoelastic
properties of polymeric interlayers for laminated glass.
Construction and Building Materials 65:1–13, 2014.
doi:10.1016/j.conbuildmat.2014.04.003.

[8] Y. Shitanoki, S. Bennison, Y. Koike. A practical,
nondestructive method to determine the shear
relaxation modulus behavior of polymeric interlayers for
laminated glass. Polymer Testing 37:59–67, 2014.
doi:10.1016/j.polymertesting.2014.04.011.

[9] I. Mohagheghian, Y. Wang, L. Jiang, et al.
Quasi-static bending and low velocity impact
performance of monolithic and laminated glass windows
employing chemically strengthened glass. European
Journal of Mechanics – A/Solids 63:165–186, 2017.
doi:10.1016/j.euromechsol.2017.01.006.

[10] T. Hána, T. Janda, J. Schmidt, et al. Experimental
and numerical study of viscoelastic properties of
polymeric interlayers used for laminated glass:
Determination of material parameters. Materials
12(14):2241, 2019. doi:10.3390/ma12142241.

[11] R. W. Clough, J. Penzien. Dynamics of Structures.
Computers & Structures, Inc, Berkeley, 3rd edn., 2003.

[12] Y. Koutsawa, et al. Static and free vibration analysis
of laminated glass beam on viscoelastic supports.
International Journal of Solids and Structures 44(25-
26):8735–8750, 2007. doi:10.1016/j.ijsolstr.2007.07.009.

[13] M. L. Aenlle, F. Pelayo. Frequency Response of
Laminated Glass Elements: Analytical Modeling and
Effective Thickness. Applied Mechanics Reviews
65(2):020802–020802–13, 2013. doi:10.1115/1.4023929.

[14] A. Zemanová, J. Zeman, T. Janda, et al. On modal
analysis of laminated glass: Usability of simplified
methods and Enhanced Effective Thickness. Composites
Part B: Engineering 151:92–105, 2018.
doi:10.1016/j.compositesb.2018.05.032.

[15] O. C. Zienkiewicz, M. Watson, I. P. King. A
numerical method of visco-elastic stress analysis.
International Journal of Mechanical Sciences 10(10):807–
827, 1968. doi:10.1016/0020-7403(68)90022-2.

[16] Z. P. Bažant, M. Jirásek. Creep and Hygrothermal
Effects in Concrete Structures, vol. 225 of Solid
Mechanics and Its Applications. Springer, Dordrecht,
2018. doi:10.1007/978-94-024-1138-6.

[17] N. M. Newmark. A Method of Computation for
Structural Dynamics. Journal of the Engineering
Mechanics Division 85(3):67–94, 1959.

[18] C. Kane, J. E. Marsden, M. Ortiz, M. West.
Variational integrators and the Newmark algorithm for
conservative and dissipative mechanical systems.
International Journal for Numerical Methods in
Engineering 49(10):1295–1325, 2000.
doi:10.1002/1097-0207(20001210)49:10<1295::AID-
NME993>3.0.CO;2-W.

[19] B. Bourdin, G. A. Francfort, J.-J. Marigo. The
variational approach to fracture. Journal of Elasticity
91(1-3):5–148, 2008. doi:10.1007/s10659-007-9107-3.

[20] M. Buliga. Hamiltonian inclusions with convex
dissipation with a view towards applications. Annals of
the Academy of Romanian Scientists Series on
Mathematics and its Applications 1(2):228–251, 2009.

510

http://dx.doi.org/10.1016/j.ijmecsci.2012.12.019
http://dx.doi.org/10.1016/j.ijmecsci.2017.05.035
http://dx.doi.org/10.12989/sem.2018.65.4.369
http://dx.doi.org/10.1016/j.conbuildmat.2014.04.003
http://dx.doi.org/10.1016/j.polymertesting.2014.04.011
http://dx.doi.org/10.1016/j.euromechsol.2017.01.006
http://dx.doi.org/10.3390/ma12142241
http://dx.doi.org/10.1016/j.ijsolstr.2007.07.009
http://dx.doi.org/10.1115/1.4023929
http://dx.doi.org/10.1016/j.compositesb.2018.05.032
http://dx.doi.org/10.1016/0020-7403(68)90022-2
http://dx.doi.org/10.1007/978-94-024-1138-6
http://dx.doi.org/10.1002/1097-0207(20001210)49:10<1295::AID-NME993>3.0.CO;2-W
http://dx.doi.org/10.1002/1097-0207(20001210)49:10<1295::AID-NME993>3.0.CO;2-W
http://dx.doi.org/10.1007/s10659-007-9107-3


vol. 60 no. 6/2020 Newmark algorithm for dynamic analysis of viscoelastic materials

[21] B. Bourdin, C. J. Larsen, C. L. Richardson. A
time-discrete model for dynamic fracture based on crack
regularization. International Journal of Fracture
168(2):133–143, 2011. doi:10.1007/s10704-010-9562-x.

[22] A. Bedford. Hamilton’s principle in continuum
mechanics. Pitman Publishing, Boston, 1985.

[23] A. Mielke, T. Roubíček. Rate-independent systems:
theory and application, vol. 655 of Applied Mathematical
Sciences. Springer, New York, 2015.
doi:10.1007/978-1-4939-2706-7.

[24] J. Schmidt, T. Janda, A. Zemanová, et al. Source
codes for preprint Newmark algorithm for dynamic
analysis with Maxwell chain model, 2019.
doi:10.5281/zenodo.3531802.

[25] L. Petzold. Automatic selection of methods for
solving stiff and nonstiff systems of ordinary differential
equations. SIAM Journal on Scientific and Statistical
Computing 4(1):136–148, 1983. doi:10.1137/0904010.

[26] E. Jones, T. Oliphant, P. Peterson, et al. SciPy:
Open source scientific tools for Python, 2001–.
http://www.scipy.org, [accessed January 12, 2021].

[27] T. J. R. Hughes. Analysis of transient algorithms
with particular reference to stability behavior. In
T. Belytschko, T. J. R. Hughes (eds.), Computational
Methods for Transient Analysis, chap. 2, pp. 67–155.
North-Holland, Amsterdam, 1983.

[28] A. Logg, K.-A. Mardal, G. Wells (eds.). Automated
Solution of Differential Equations by the Finite Element
Method: The FEniCS Book. Lecture Notes in
Computational Science and Engineering. Springer-Verlag,
Berlin Heidelberg, 2012. doi:10.1007/978-3-642-23099-8.

[29] M. S. Alnæs, J. Blechta, J. Hake, et al. The FEniCS
Project Version 1.5. Archive of Numerical Software
3(100), 2015. doi:10.11588/ans.2015.100.20553.

[30] T. Li, J.-J. Marigo, D. Guilbaud, S. Potapov.
Numerical investigation of dynamic brittle fracture via
gradient damage models. Advanced Modeling and
Simulation in Engineering Sciences 3(1):26, 2016.
doi:10.1186/s40323-016-0080-x.

[31] T. Roubíček. Models of dynamic damage and phase-
field fracture, and their various time discretisations. In
M. Hintermüller, J. F. Rodrigues (eds.), Topics in
Applied Analysis and Optimisation, pp. 363–396.
Springer International Publishing, Cham, 2019.

511

http://dx.doi.org/10.1007/s10704-010-9562-x
http://dx.doi.org/10.1007/978-1-4939-2706-7
http://dx.doi.org/10.5281/zenodo.3531802
http://dx.doi.org/10.1137/0904010
http://www.scipy.org
http://dx.doi.org/10.1007/978-3-642-23099-8
http://dx.doi.org/10.11588/ans.2015.100.20553
http://dx.doi.org/10.1186/s40323-016-0080-x

	Acta Polytechnica 60(6):502–511, 2020
	1 Introduction
	2 Newmark method
	2.1 Governing equations
	2.2 Discretization

	3 Variational integrators
	3.1 Variational framework
	3.2 Discretization
	3.3 Equivalence to Newmark
	3.4 Energy balance

	4 Examples
	4.1 Discrete problem
	4.2 Generalization

	5 Conclusions
	Acknowledgements
	References