Comparison of viscoelastic finite element models for laminated glass beams
Alena Zemanová∗, Jan Zeman, Michal Šejnoha
Department of Mechanics, Faculty of Civil Engineering, Czech Technical University in Prague,
Thákurova 7, 166 29 Prague 6, Czech Republic

arXiv:1701.03636v2 [cs.CE] 15 Jun 2017

Abstract
Laminated glass elements, which consist of stiff elastic glass layers connected with a compliant viscoelastic polymer foil, exhibit geometrically non-linear and time/temperature-sensitive behavior.
In computational modeling, the viscoelastic effects are often neglected or a detailed continuum
formulation typically based on the volumetric-deviatoric elastic-viscoelastic split is used for the interlayer. Four layerwise beam theories are introduced in this paper, which differ in the non-linear
beam formulation at the layer level (von Kármán/Reissner) and in constitutive assumptions for
the interlayer (a viscoelastic solid with the time-independent bulk modulus/Poisson ratio). We
perform detailed verification and validation studies at different temperatures and compare the accuracy of the selected formulation with simplified elastic solutions used in practice. We show that
all the four formulations predict very similar responses. Therefore, our suggestion is to use the most
straightforward formulation that combines the von Kármán model with the assumption of timeindependent Poisson ratio. The simplified elastic model mostly provides response in satisfactory
agreement with full viscoelastic solutions. However, it can lead to unsafe or inaccurate predictions for rapid changes of loading. These findings provide a suitable basis for extensions towards
laminated plates and glass layer fracture, owing to the modular format of layerwise theories.
Keywords: Laminated glass, finite element method, finite strain Reissner model, von Kármán
assumptions, generalized Maxwell model, Williams-Landel-Ferry equation

1. Introduction
Laminated glass units consist of multiple glass sheets bonded together with polymer foils. The
main reason to combine brittle glass sheets with soft polymer layers is to enhance the post-fracture
behavior, so that a unit retains the load-bearing capacity after fracture of one or multiple glass
plies [1]. Of course, the interlayer foil also influences the overall response of the unfractured
elements, which is in the focus of this study, by mediating the shear interaction among the glass
layers.
The glass layers can be accurately modeled as linear elastic, but the polymer foil exhibits a
strongly time/temperature-dependent response, e.g. [2–4, and references therein]. Several types of
interlayers have found their use in practice. From the first application of laminated glass members with polyvinyl butyral foil (PVB) in car windscreens, laminated glass has expanded into the
∗

Corresponding author
Email addresses: alena.zeman@fsv.cvut.cz (Alena Zemanová), jan.zeman@fsv.cvut.cz (Jan Zeman),
sejnom@fsv.cvut.cz (Michal Šejnoha)
Preprint submitted to arXiv

June 16, 2017

building constructions, and new types of interlayer foils became available, such e.g. ionoplast polymer foil (IP) that is stiffer and less sensitive to viscoelastic phenomena than PVB at common
temperatures.
As discussed in our previous contribution [5], the behavior of polymeric interlayer is often
approximated as linear elastic, assuming a known temperature profile and duration of loading. The
following methods have been developed for elastic laminated glass beams under this simplification:
• The basic approach to approximate the mechanical response of laminated glass units involves
two limiting cases: the layered case representing an assembly of independent glass layers
and the monolithic case with thickness equal to the combined thickness of glass layers and
interlayers, under the assumption of geometric linearity [6].
• Single-layered approaches, such as effective thickness methods that approximate laminated
units by monolithic, geometrically linear, systems that achieve the equivalent maximum deflection or stress [7–9]; see also [10, 11] for an overview.
• Multi-layered approaches, explicitly resolving all the layers. These include analytical formulations, which are solved in a closed form for specific boundary conditions [12–14] and
geometrically linear behavior, and numerical formulations based on the finite-difference [15]
or finite-element methods [5, 16, 17] that can treat general boundary conditions and geometrically non-linear effects.
Solely the multi-layered approaches have been extended to account for viscoelastic effects, because they offer sufficient accuracy in stresses and strain distributions at the ply level. An analytical
model of geometrically linear, simply supported, beams with viscoelastic interlayer represented by
a Maxwell chain was developed by Galuppi and Royer-Carfagni [18, 19], who approximated the
viscoelastic solution by the Fourier series in the spatial coordinate and solved the ensuing integral
equations in the closed form to account for the time dependent effects. The same approach was later
successfully used to study the stability of laminated columns [20] or time-dependent response of
cold-bent laminated beams [21]. The most recent analytical study has been performed by Wu et al.
[22], who develop a closed-form solutions for two-dimensional stress and strain fields in two-layer
simply-supported beam with interlayer response described by the standard linear model. The only
detailed numerical studies we are aware of include geometrically non-linear three-dimensional finite
element simulations of viscoelastic plates under transverse loading by Duser et al. [23] and Bennison
et al. [24], or laminated beams under lateral-torsional buckling by Bedon et al. [25].
The goal of this paper is to propose a viscoelastic geometrically non-linear finite element formulation for laminated glass beams based on layerwise refined laminate theories initiated by Mau
[26]. The basis of our formulation is the elastic finite strain model developed earlier by Zemanová
et al. [5], where we demonstrated that it offers accuracy comparable to detailed two-dimensional
simulations at much lower computational cost. Here, this formulation is extended to account for
the time/temperature-dependent behavior of the interlayer, which is incorporated into the numerical model by the incremental exponential algorithm by Zienkiewicz et al. [27]. In particular, we
will compare the formulations arising from the two modeling aspects:
• beam kinematics, involving the finite strain shear-deformable Reissner beam formulation and
a simpler variant based on the Timoshenko beam theory complemented with the von Kármán
assumptions of large deflections,
2

• constitutive assumptions for viscoelastic models, where we will consider models based on
time/temperature-dependent shear modulus, complemented with either constant (time-independent)
bulk modulus or constant Poisson ratio.
The differences in the mechanical results of the corresponding models will be studied in detail and
verified against detailed finite element simulations, with the goal to determine the simplest, but
still sufficiently accurate, formulation that is suitable for the generalization of elastic models of
laminated glass plates [28] to the viscoelastic regime.
The rest of the paper is organized as follows. In Section 2, we introduce the different variants of
laminated glass beam models, starting from kinematics assumptions in Section 2.1, proceeding to
interlayer constitutive relations in Section 2.2, up to the finite element discretization in Section 2.3.
Section 3 is devoted to the comparison of different formulations, along with an additional verification, validation, and parametric studies. The most important findings are summarized in Section 4.
Finally, in Appendix A we gather technical details regarding the finite element technology, in order
to make the paper self-contained.
2. Finite element models
The models introduced in this work share the following assumptions:
• they assume planar cross sections of individual layers but not of the whole laminated glass
unit, because the stiffness of glass layer and the effective stiffness of the interlayer differs by
at least three orders of magnitude, e.g. [15],
• being based on Mau’s refined plate theory [26], each layer is treated independently and the
inter-layer compatibility is enforced by the Lagrange multipliers. The main advantages of
this approach is that the different constitutive models for each layer can be easily combined
together and delamination phenomena can be efficiently accounted for, if needed, e.g., [29]1 .
As noted above, the formulations differ in the beam kinematics (Finite-Strain (FS) Reissner vs.
Large-Deflection von Kármán (VK) models) and in the assumption of viscoelastic constitutive
models (constant bulk modulus K vs. the Poisson ratio ν). These considerations lead to four
models, which are briefly introduced in this section.
Our notation is as follows. Scalar quantities are denoted by lightface letters, e.g., a, and
matrices are denoted in bold, e.g., a or A. In addition, AT stands for the matrix transpose, A−1
for the matrix inverse, and a(i) denotes that the quantity a is related to the i-th layer. Also note
that in order to avoid a profusion of notation, in Section 2.1 we omit the dependence of quantities
of interest on the time variable t, while in Section 2.2 we omit the layer index •(i) , because the
constitutive description holds only for the interlayer.
1

Note, however, that the perfect inter-layer adhesion is assumed in this paper, because the application of high
pressures and temperatures during the production process result in high gluing forces of chemical nature [30].

3

2.1. Kinematics
2.1.1. Reissner finite strain model
In the Reissner finite strain beam theory [31], the non-zero displacement components of the
i-th layer can be parameterized as, Figure 1,


(i)
u(i) (x, z (i) ) = u0 (x) + sin ϕ(i)
(x)
z (i) ,
(1a)
y




(i)
(i)
w(i) (x, z (i) ) = w0 (x) + cos ϕ(i)
(1b)
y (x) − 1 z ,
(i)

(i)

(i)

where u0 and w0 are centerline displacements, ϕy is the cross-section rotation, x is the coordinate
measured along the centerline, z (i) is the coordinate measured along the cross-section, and i = 1, 2, 3
refers to individual layers.
The inter-layer compatibility is ensured via the geometric continuity conditions at the interfaces
between the layers (with i = 1, 2)
h(i+1)
h(i)
) − u(i+1) (x, −
) = 0,
2
2
h(i+1)
w(i) (x, z (i) ) − w(i+1) (x, −
) = 0,
2
where the perfect horizontal and vertical adhesion is supposed. From (1), the compatibility condition become

 h(i+1)


h(i)
(i)
(i+1)
(i+1)
sin ϕ(i)
(x)
+
sin
ϕ
(x)
= 0,
(2a)
u0 (x) − u0
(x) +
y
y
2
2

 h(i+1)

 h(i) h(i+1)
h(i)
(i)
(i+1)
w0 (x) − w0
(x) +
cos ϕ(i)
cos ϕ(i+1)
(x) −
−
= 0,
(2b)
y (x) +
y
2
2
2
2
u(i) (x,

b
L
y (i)

x

h(i)
z (i)
(i)

ϕy

(i)

w0

z (i)

z (i)

(i)

u0

Figure 1: Kinematics of a cross section of the bottom layer of a laminated beam (i = 3).
(i)

(i)

The non-zero normal, εx , and shear, γxz , strain components follow from
(i)

(i)
(i)
(i)
ε(i)
x (x, z ) = εx,0 (x) + κy (x)z ,
!
(i)



 dw(i) (x)
du
(x)
(i)
0
0
γxz
(x) = sin ϕ(i)
1+
+ cos ϕ(i)
,
y (x)
y (x)
dx
dx

4

(3a)
(3b)

(i)

(i)

where εx,0 and κy denote the generalized axial strain and pseudo-curvature of the layer reference
axis, introduced by Reissner [31] in the form
!
(i)

 dw(i) (x)


du
(x)
(i)
(i)
0
0
1
+
−
sin
ϕ
(x)
εx,0 (x) = cos ϕ(i)
(x)
− 1,
(4a)
y
y
dx
dx
(i)

κ(i)
y (x) =

dϕy (x)
.
dx

(4b)

2.1.2. Von Kármán model
A simpler model is obtained when adopting the von-Kármán assumptions of moderate deflec(i)
(i)
tions w0 and small rotations ϕy , e.g. [32, Section 4.3]. Under these assumptions, the displacements of the i-th layer (1) simplify to
(i)

u(i) (x, z (i) ) = u0 (x) + ϕy(i) (x)z (i) ,
(i)

w(i) (x) = w0 (x),
(i)

and the inter-layer compatibility conditions (1) become linear in ϕy ,
(i)

(i+1)

u0 (x) − u0

(x) +

h(i+1) (i+1)
h(i) (i)
ϕy (x) +
ϕy
(x) = 0,
2
2
(i)
(i+1)
w0 (x) − w0
(x) = 0.

(5a)
(5b)

In addition, the normal and shear strains are expressed as, cf. (3),
(i)

(i)
(i)
(i)
ε(i)
x (x, z ) = εx,0 (x) + κy (x)z ,

(6a)

(i)

(i)
γxz
(x) = ϕ(i)
y (x) +

dw0
(x),
dx

(6b)

(i)

where the axial strains εx,0 must be measured by the Green-Lagrange strain tensor because of large
(i)

deflections w0 ,
(i)

du0
(i)
εx,0 (x) =
(x) +
dx

1
2

!2
(i)
dw0
(x) ,
dx

(7)

(i)

while the layer pseudo-curvature κy still follows from (4b).
2.2. Constitutive relations
As indicated earlier, accurate description of time/temperature-dependent response of polymeric
interlayer is decisive for the development of predictive models of laminated glass structures. In
Section 2.2.1 we introduce such a model in the framework of linear viscoelasticity with accelerated/retarded time to account for the temperature dependence. We will also clarify the impact of
the constitutive assumptions on the structure of stress-strain relations. The integral form of constitutive equations used in Section 2.2.1 is not very convenient for the numerical implementation,
because all history variables during the whole loading process must be stored in order to evaluate
5

the integral. Therefore, in the following sections, we will develop an incremental approach based
on the well-established exponential algorithm by Zienkiewicz et al. [27]. The case of the constant
bulk modulus is treated in detail in Section 2.2.2, leading to a rather involved procedure because
of assumptions of beam theories. Section 2.2.3 then explains how the procedure simplifies under
the assumption of time-independent Poisson ratio.
2.2.1. Time/temperature-dependent behavior of polymer foil
The main engineering property relevant to the composite behavior of the units is the shearstress versus shear-strain characteristics of the soft interlayer. The shear modulus of the polymer
interlayer G is experimentally determined as a function of duration of loading and temperature,
see Figure 2 for data for two PVB foils.
3

3

10
10°C
20°C
30°C
40°C
50°C
60°C
70°C

2

10

1

10

SHEAR MODULUS [MPa]

SHEAR MODULUS [MPa]

10

0

10

-1

10

-2

10 -12
10

(a)

10°C
20°C
30°C
40°C
50°C
60°C
70°C

2

10

1

10

0

10

-1

-9

10

-6

10

-3

10

0

10

3

10

6

10

9

10

10 -12
10

12

10

DURATION OF LOADING [s]

(b)

-9

10

-6

10

-3

10

0

10

3

10

6

10

9

10

DURATION OF LOADING [s]

Figure 2: Shear modulus G as a function of the duration of loading and temperature, after (a) Bennison et al. [24]
and (b) Pelayo et al. [33].

The temperature dependence is taken into account by the time-temperature superposition principle, in which the true time ttrue is replaced by the adjusted value t = ttrue /aT modified by the
temperature-dependent shift factor aT . We employ the Williams-Landel-Ferry (WLF) equation [34]
for this purpose:
C1 (T − T0 )
log aT = −
,
(8)
C2 + T − T0
in which C1 and C2 are material constants, and T and T0 are the current and reference temperatures, respectively.
For further discussion, it will be convenient to decompose the vectors of strain ε and stress σ,




ε = εx εy εz γxy γxz γyz T ,
σ = σx σy σz τxy τxz τyz T ,
into the volumetric and deviatoric parts
1
ε(t) = εV (t)i + e(t),
3

(9)

σ(t) = σm (t)i + s(t),
(10)

T
where t denotes the time instant, i = 1 1 1 0 0 0 . The volumetric parts follow from
εV = iT ε and σm = 31 iT σ, and the deviatoric parts have the components




e = ex ey ez γxy γxz γyz T ,
s = sx sy sz τxy τxz τyz T .
6

Assuming isotropic material behavior and a smooth strain history ε(t) with ε(0) = 0, the stress
at time t reads, e.g., [35, Section 1.2]
t

Z
σ(t) = σm (t)i + s(t) = i

dεV 0
(t ) dt0 + 2Q
dt0

K(t − t0 )

0

Z

t

G(t − t0 )

0

de 0
(t ) dt0 ,
dt0

(11)



where Q = diag 1 1 1 0.5 0.5 0.5 , and K(t − t0 ) and G(t − t0 ) denotes the bulk and shear
relaxation moduli that completely characterize the interlayer response.
s

G1

G2

G3

GP

ev η1

η2

η3

ηP

ee

G∞

e

s
Figure 3: Generalized Maxwell chain consisting of P viscoelastic units and one elastic spring.

In what follows, the response under shear will be represented by a generalized Maxwell chain,
Figure 3, whose relaxation function is provided by the Dirichlet-Prony series, e.g., [35, page 32]
0

G(t − t ) = G∞ +

P
X

0

− t−t
θ

Gp exp

p

= G0 −

p=1

P
X

Gp (1 − exp

− t−t
θ
p

0

),

(12)

p=1

where G∞ is the long-term shear modulus, P stands for the number of viscoelastic units, Gp denotes
the shear modulus of the p-th unit, θp = ηp /Gp is its relaxation time related to the viscosity ηp ,
P
and G0 = G∞ + Pp=1 Gp is the elastic shear modulus of the chain.
In order to close the constitutive description for the interlayer, two assumptions on the bulk relaxation functions will be considered. First, the constant bulk modulus approximation, e.g. [23, 24],
under which the material behaves as a linear elastic solid for volumetric loading. The constitutive
equation
Zt
de
σ(t) = KiεV (t) + 2 G(t − t0 )Q 0 (t0 ) dt0 ,
(13)
dt
0

thus additionally involves the constant elastic bulk modulus. Second, the Poisson ratio ν is taken
as constant, e.g. [7, 22], which simplifies the constitutive relation (11) into, e.g., [36, Eq. (29.63)]:
Zt
σ(t) = D ν

G(t − t0 )

0

7

dε 0
(t ) dt0 ,
dt0

(14)

with the constant matrix D ν provided by

2(1 − ν)
2ν
2ν

2ν
2(1
−
ν)
2ν

1 
2ν
2ν
2(1 − ν)

Dν =
0
0
0
1 − 2ν 


0
0
0
0
0
0

0
0
0
1
0
0

0
0
0
0
1
0

0
0
0
0
0
1





.




(15)

2.2.2. Incremental formulation for constant bulk modulus
Material point. We start by dividing the time interval of interest [0; tmax ] into non-equidistant time
instants 0 = t0 < t1 < t2 < · · · < tN −1 < tN = tmax (recall that these correspond to temperatureadjusted values according to (8)). The strain history over the time interval [tn ; tn+1 ] is assumed to
be known, so that
e(tn+1 ) = e(tn ) + ∆e,
where ∆e denotes the increment of deviatoric strain between the time instants tn and tn+1 . An
analogous relation holds for the deviatoric stresses,
s(tn+1 ) = s(tn ) + ∆s,
so that, assuming that s(tn ) is known, the goal is to determine the increment of the deviatoric
stress ∆s.
For a linear strain variation over the time interval [tn ; tn+1 ], this increment can be expressed
as [27]
b
∆s = 2GQ∆e
+ ∆b
s,

(16)

b stands for the effective shear modulus over the time interval and ∆b
where G
s for the stress relaxation under the zero strain increment. For the Maxwell chain model, Figure 3, these quantities
follow from
b = G∞ +
G

P
X
p=1

θp
− ∆t
Gp (1 − exp θp ),
∆t

∆b
s=

P
X

∆b
sp = −

p=1

P
X

sp (tn )(1 − exp

− ∆t
θ
p

),

(17)

p=1

where ∆t = tn+1 − tn is the temperature-adjusted time step (we assume for simplicity that the
temperature remains constant over the interval), and sp and ∆b
sp stand for the values of the
deviatoric stresses and its relaxation in the p-th unit.
Because the volumetric response is assumed to be elastic, the increment of the mean strain,
∆σm , can be obtained from the given increment of the volumetric strain ∆εV directly as
∆σm = K∆εV .

(18)

Beams. As follows from the standard assumptions of the conventional beam theory, e.g. [37, Section
29], the only non-zero stress components are the increments of normal and shear stresses ∆σx and
∆τxz ; the corresponding volumetric and deviatoric stress components entering (18) and (16) thus
become ∆σm = 1/3∆σx , ∆sx = 2/3∆σx , and ∆sxz = ∆τxz .
8

In order to obtain the corresponding strain increments ∆εx and ∆γxz needed in the beam
formulations, we employ the volumetric-deviatoric split (9) and invert the relations (16) and (18)
to obtain
1
1
1
1
1
∆σm +
∆εx = ∆εV + ∆ex =
(∆sx − ∆b
sx ) = ∆σx −
∆b
s ,
b
b
b x
3
3K
E
2G
2G
1
τxz ) ,
∆γxz = (∆τxz − ∆b
b
G

(19a)
(19b)

from which we directly obtain the incremental stress-strain relations in the form
b x + (1 + νb)∆b
∆σx = E∆ε
sx ,

b xz + ∆b
∆τxz = G∆γ
τxz .

(20)

b and νb in Eqs. (19) and (20) stand for the effective Young modulus and Poisson
The symbols E
b = (9K G)/(
b G
b + 3K) and νb = (3K − 2G)/(2(3K
b
b e.g. [36, Table D.1]. It
ratio, obtained as E
+ G)),
will also prove useful to relate the deviatoric normal strain to the full normal strain via
2
(1 + νb)
∆ex = (1 + νb)∆εx −
∆b
sx .
3
9K

(21)

Cross-section. In order to obtain the constitutive equations for a cross-section located at x, recall
from (4) or (6) that the normal strain εx is an affine function in the through-the-thickness coordinate
z, whereas the shear strain γxz is independent of z. Upon integrating the stress increments (20)
over the cross-section, we obtain
Z
b
b
b
∆Nx (x) = EA∆ε
(x)
+
(1
+
ν
b
)
∆b
sx (x, z) dA = EA∆ε
(22a)
x,0
x,0 (x) + ∆Nx (x),
A
Z
b y ∆κy (x) + ∆M
cy (x),
b
∆b
sx (x, z)z dA = EI
(22b)
∆My (x) = EIy ∆κy (x) + (1 + νb)
A
Z
b s ∆γxz (x) + ∆Vbz (x),
b s ∆γxz (x) +
∆b
τxz (x) dA = GA
(22c)
∆Vz (x) = GA
As

where ∆Nx , ∆Vz , and ∆My stand for the increments of the normal and shear forces, and the
bending moment (the corresponding relaxation forces are distinguished by a hat); ∆ε0,x , ∆γxz ,
and ∆κy are the corresponding increments of the centerline strain, the transverse shear, and the
pseudo-curvature determined from (4) or (6); and A, Iy , and As denote the cross-section area, the
second moment of area, and the shear-effective area.2 Note that in (22), we utilized the fact that
b G,
b and νb are identical for the whole layer.
the cross-section and the effective material constants E,
As follows from Eq. (22), the increments of internal forces involve the relaxation contributions
b
cy , and Vbz , which in turn depend on the values of the deviatoric stresses in elements of
Nx , M
the Maxwell chain (16, 17) across the cross-section. Assuming that sp (t0 ) = 0, a closer inspection of these relations together with Eq. (21) reveals that the normal and shear deviatoric stress
components exhibit the affine-constant distribution across the cross-section at x at time tn :
M
sx,p (tn , x, z) = sN
x,p (tn , x) + sx,p (tn , x)z,

V
τxz,p (tn , x, z) = τxz,p
(tn , x) for p = 1, 2, . . . , P ;
(i)

2

(23)

For three-layered units (glass/foil/glass) we set the value of the shear-effective area to As = 5/6A(i) for glass
(2)
layers with i = 1, 3 and As = A(2) for the foil. For these values, which correspond to a quadratic shear stress
variation in glass and to constant shear in the interlayer, e.g. [38], we achieved the best agreement with the response
of the detailed two-dimensional model from the ADINA solver (not shown).

9

M
V
therefore, the three quantities sN
x,p , sx,p , and τxz,p completely determine the stress distribution
in the p-th elements of the Maxwell chain over the cross-section. They also allow for a compact
expression of the incremental relaxation forces in the form,

bx (x) = (1 + νb)A∆b
∆N
sN
x (x),

cy (x) = (1 + νb)Iy ∆b
∆M
sM
x (x),

V
∆Vbz (x) = As ∆b
τxz
(x),

(24)

where, e.g.,
∆b
sN
x (x)

=

P
X

∆b
sN
x,p (x),

− ∆t
θ

N
∆b
sN
x,p (x) = −sx,p (tn , x)(1 − exp

p

),

p=1
V follow by analogy.
recall (17). The remaining quantities ∆b
sM
τxz
x and ∆b
At the end of the time interval, the stresses in the p-th unit must be updated on the basis of
increments of the generalized strains ∆ε0,x , ∆κy , and ∆γxz for a cross-section at x. Combining (21)
and (16) with (3) or (6), one obtains

b p (1 + νb)
4
2G
b p ∆ε0,x (x) + ∆b
∆sN
b)G
sN
∆b
sN
x,p (x) = (1 + ν
x,p (x) −
x (x),
3
9
K
b p (1 + νb)
4
2G
b p ∆κy (x) + ∆b
∆sM
b)G
sM
∆b
sM
x,p (x) = (1 + ν
x,p (x) −
x (x),
3
9
K
V
V
b p ∆γxz (x) + ∆b
(x) = G
τxz,p
(x).
∆τxz,p

(25a)
(25b)
(25c)

2.2.3. Incremental formulation for constant Poisson ratio
Material point. Assuming the same time discretization as in Section 2.2.2, the strain and stress
increments over the interval [tn ; tn+1 ] are expressed as
ε(tn+1 ) = ε(tn ) + ∆ε,

σ(tn+1 ) = σ(tn ) + ∆σ,

which involve the full strain and stress tensors. Comparing (11) with (14), we see that the incremental formula (16) from the previous section becomes
b
∆σ = D ν G∆ε
+ ∆b
σ,

∆b
σ=

P
X
p=1

∆b
σp = −

P
X

σ p (tn )(1 − exp

− ∆t
θ
p

),

(26)

p=1

b is defined with (17).
where matrix D ν depends on the constant Poisson ratio ν via (15) and G
Beams. Increments of the non-zero components of stress vector, ∆σx and ∆τxz , follow directly
from (26) in the form
b x + ∆b
∆σx = E∆ε
σx ,

b xz + ∆b
∆τxz = G∆γ
τxz ,

b = 2(1 + ν)G,
b cf. (20)
with the effective Young modulus now given by E

10

Cross-section. Proceeding in the same way as in the previous section, the incremental constitutive
relations at the cross-section level, Eq. (22), attain a similar form:
Z
b
b
b
∆b
σx (x, z) dA = EA∆ε
(27a)
∆Nx (x) = EA∆εx,0 (x) +
x,0 (x) + ∆Nx (x),
A
Z
b y ∆κy (x) + ∆M
cy (x),
b y ∆κy (x) +
∆b
σx (x, z)z dA = EI
(27b)
∆My (x) = EI
A
Z
b s ∆γxz (x) + ∆Vbz (x),
b s ∆γxz (x) +
∆b
τxz (x) dA = GA
(27c)
∆Vz (x) = GA
As

bx , ∆Vbz , and ∆M
cy depend on the distribution
Again, the increments of relaxation forces ∆N
of stresses in the Maxwell units over the cross-section. By analogy to (23), the normal and shear
stresses in the p-th unit admit the decomposition:
N
M
σx,p (tn , x, z) = σx,p
(tn , x) + σx,p
(tn , x)z,

V
τxz,p (tn , x, z) = τxz,p
(tn , x) for p = 1, 2, . . . , P,

from which we obtain
bx (x) = A∆b
∆N
σxN (x),

cy (x) = Iy ∆b
∆M
σxM (x),

V
∆Vbz (x) = As ∆b
τxz
(x),

with, e.g.,
∆b
σxN (x) =

P
X

N
∆b
σx,p
(x),

N
N
∆b
σx,p
(x) = −σx,p
(tn , x)(1 − exp

− ∆t
θ
p

).

p=1
N , ∆σ M , and ∆τ V
Finally, we need to relate the increments of the history variables ∆σx,p
x,p
xz,p to
the increments of the generalized strains ∆εx,0 , ∆κy , and ∆γxz . The relations for the normal
components follow directly from (26) as
N
N
b p ∆εx,0 (x) + ∆b
∆σx,p
(x) = 2(1 + ν)G
σx,p
(x),
M
M
b p ∆κy (x) + ∆b
∆σx,p (x) = 2(1 + ν)G
σx,p (x).

(28a)
(28b)

V
is obtained from (25c).
The shear stress increment ∆τxz,p

2.3. Finite element model
In this section, we extend the geometrically nonlinear Reissner elastic solver from [5] to the
incremental viscoelastic formulation and discuss the simplifications due to the von-Kármán assumptions.
2.3.1. Variational setting
Similarly to [5], we derive the model using variational arguments. Because the glass layers are
assumed to be elastic, their response is governed by a quadratic stored energy density introduced
in [5, Eq. (13)]:

 
 2 


 2
(i) 2
(i)
(i) (i)
(i)
W (i) ε(i) = 21 E (i) A(i) εx,0 + G(i) A(i)
,
κ
+
E
I
γxz
s
y
y
11

h
i
(i)
(i) T
collecting the generalized strain
with i = 1, 3 and with the column matrix ε(i) = ε(i)
γxz
x,0 κy
components. For the interlayer, i = 2, the energy density is history dependent to account for the
incremental form of the cross-section constitutive relations (22) or (27). For the time interval
[tn ; tn+1 ], it reads






(2)
Wn+1 ε(2) = Wn(2) ε(2) (tn ) + ∆W (2) ε(2) ,
(2)

with W0

= 0 and

∆W

(2)



ε

(2)



=

1
2



2

2

2 
(2)
(2)
(2)
(2)
(2)
(2)
(2)
(2)
(2)
b A
b As ∆κy
b Iy ∆γxz
E
∆ε
+G
+E
x,0






b (2) ∆ε(2) + V (2) (tn ) + ∆Vb (2) ∆γ (2)
+ Nx(2) (tn ) + ∆N
x
z
z
xz
x,0


cy(2) ∆κ(2)
+ My(2) (tn ) + ∆M
y .
(2)

(2)

(29)

(2)

Here, e.g., ∆εx,0 = εx,0 − εx,0 (tn ), and the remaining terms follow by analogy.
The energy functional at the structural level collects the contributions of layers
Z


XZ


(2)
Wn+1 ε u(2) (x) dx,
Πint,n+1 (u) =
W (i) ε u(i) (x) dx +
Ω(2)

(i)
i=1,3 Ω

h
i3
(i)
(i)
where u = u(i)
stands for the centerline displacements and cross-section rotations
,
w
,
ϕ
y
0
0
i=1
and ε(•) abbreviates the strain-displacement relations for the Reissner or the von-Kármán kinematics introduced in Sections 2.1.1 and 2.1.2, respectively. Complementing the internal energy
functional with the energy of external forces provides
Πn+1 (u) = Πint,n+1 (u) +

3
X

(i)

Πext (tn+1 , u(i) ).

(30)

i=1

Note that the admissible generalized displacements u must satisfy the inter-layer continuity conditions (2) or (5), written compactly as

c u(x) = 0 for x ∈ [0; L].
(31)
The generalized displacements at time tn+1 , u(tn+1 ), then follow by the minimization of the total
energy functional (30) under the equality constraints (31).
2.3.2. Discretization
Discretization of the model is accomplished using the standard finite element procedures. To
(i)
that purpose, each layer is discretized into the identical number of ne elements Ωe and the generalized displacements are approximated as
(i)
u(i) ≈ N e (x)d(i)
e for x ∈ Ωe ,

12

(32)

(i)

where N e stands for the 3 × 6 matrix of linear basis functions and de is the 6 × 1 column matrix
of element degrees of freedom; see Appendix A for details. After inserting the approximation (32)
into (30), the internal energy becomes the function of the nodal degrees of freedom
ne
XX

Πint,n+1 (d) =

ne

 X


(i)
(2)
(2)
+
Π
d
,
Πint,e d(i)
e
e
int,n+1,e

i=1,3 e=1

e=1

i3,ne +1
h
collects the layer nodal displacements and, e.g.,
where d = d(i)
j
i=1,j=1

(2)
Πint,n+1,e



d(2)
e



Z
=

(2)

Ωe



(2)
Wn+1 ε N e (x)d(i)
dx.
e

In analogy to the previous section, the column matrix of nodal displacements d(tn+1 ) follows
from minimization of the total energy function
Πn+1 (d) = Πint,n+1 (d) −

3 X
ne
X

(i)

T
d(i)
f ext,e (tn+1 ) ,
e

i=1 e=1
(i)

in which f ext,e denotes the generalized nodal forces due to external loading of the i-th layer. The
kinematically admissible displacements must satisfy the constraint
c(d) = 0,

(33)

with the meaning of the inter-layer continuity conditions (31) enforced at the nodes.
2.3.3. Solution procedure
Following the standard procedures of equality-constrained optimization, e.g., [39, Chapter 14],
we consider the Lagrange function
Ln+1 (d, λ) = Πn+1 (d) + λT c(d),
where the column matrix of Lagrange multipliers λ has the meaning of interface nodal forces that
ensure the inter-layer compatibility, (33). The nodal displacements d(tn+1 ) and the corresponding
Lagrange multipliers λ(tn+1 ) then arise as the saddle point of the Lagrange function.
The optimality conditions lead to a system of equations non-linear in d(tn+1 ) which is solved
iteratively using the Newton method. To this purpose, we express the (k + 1)-th approximation to
d(tn+1 ) as
k+1

d = k d + k+1 δd for k = 0, 1, . . .

where k d denotes a known previous iterate; for k = 0, we set 0 d = d(tn ). The displacement
increment k+1 δd and the vector of Lagrange multipliers k+1 λ are determined from the linearized
system
k


k

K k C T k+1 δd
f int − f ext (tn+1 )
=
−
,
(34)
kC
k+1 λ
kc
0
13

in which the individual terms correspond to, cf. [5, Eq. (31)]
 

k
f int = ∇Πint,n+1 k d = f int k d ,



k
K = ∇2 Πint,n+1 k d = K t k d + K λ k d, k λ ,

k


c = c kd ,

(35a)


k
C = ∇c k d .

(35b)

Because the unknown displacements are approximated independently at each layer, recall Eq. (32),
the internal forces f int and tangent stiffness K t matrices possess the block-diagonal structure:

(1)
f int (d(1) )
 b (2) (2)
b (2) (2) 
f int (d) = f
int (d ) + ∆f int (d )


(36a)

(3)
f int (d(3) )

 (1) (1)
K t (d )

c(2) (d(2) )
c(2) (d(2) ) + ∆K
K t (d) = 
K
t
t


(3)
K t (d(3) )


.

(36b)

b and E
b (as emThe sub-matrices for the interlayer must be determined using the effective values G
b
ct and ∆f
phasized by the hat), and additionally contain the history-dependent contributions ∆K
int
that result from the incremental energy split (29). As usual, all sub-matrices in (34) are obtained
by the assembly of element contributions, e.g. [32, Section 2.3.6]. For the readers’ convenience,
these are discussed in Appendix A and for the Reissner model particularly in Appendix A.2.
To the same purpose, in Algorithm 1 we present a pseudo-code of one-step geometrically nonlinear viscoelastic problem that summarizes the developments of the current section. Notice that
the algorithm is terminated by two residuals, cf. [5, Eq. (34)],
k

η1 =

kk f int − f ext + k C Tk λk2
,
max (kf ext k2 , 1)

k

η2 =

kk ck2
,
min3i=1 h(i)

related to the nodal equilibrium and continuity conditions, respectively.
A similar procedure applies also to the von Kármán model with a few simplifications. Indeed,
because the continuity constraint (33) is linear in d, recall (5), we have k c = 0, matrix C remains
constant during iterations, and matrix K λ = 0. The remaining terms needed in (34) are provided
in Appendix A.1.
3. Comparison, verification, and validation
The developments described in the previous section resulted in four geometrically non-linear
finite element models summarized in Table 1; recall that they emerge as the combinations of Reissner (FS) or von Kármán (VK) kinematics and the assumptions of constant bulk modulus (K) or
Poisson ratio (ν) of the interlayer. We also include the results of geometrically linear formulation (LIN) [40], in order to highlight the importance of geometric non-linearity where relevant.
In addition, the following data are considered in the studies reported in this section:
• detailed two-dimensional viscoelastic finite element simulations (2D), performed in ADINA 9.0.2 (ADINA R&D, Inc.) finite element system,
14

Algorithm 1: One-step exponential algorithm for geometrically nonlinear laminated beams
Data:
tolerances 1 and 2 , load f ext (tn+1 ),
displacements d(tn ), interlayer internal forces and generalized strains at tn
 (2)
ne
(2)
(2)
Nx,e (tn ), Vz,e (tn ), My,e (tn ) e=1 ,

P,ne
M
V
. constant K: Maxwell chain stresses sN
x,p,e (tn ), sx,p,e (tn ), τxz,p,e (tn ) p=1,e=1 ,
 N
P,ne
M (t ), τ V
. constant ν: Maxwell chain stresses σx,p,e
(tn ), σx,p,e
n
xz,p,e (tn ) p=1,e=1
Initialization:
k ← 0, 0 λ ← 0, 0 d ← d(tn ), assemble 0 f int , 0 c, and 0 C in (36a) and (34)
while (k η1 > 1 ) or (k η2 > 2 ) do
assemble k K from Eq. (36b)
solve for (k+1 δd, k+1 λ) from Eq. (34)
k+1 d ← k d + k+1 δd
assemble k+1 f int , k+1 c, and k+1 C in (36a) and (34)
k ←k+1
d(tn+1 ) ← k+1 d, λ(tn+1 ) ← k+1 λ
. constant K: update internal forces from (22) and Maxwell chain stresses from (25)
. constant ν: update internal forces from (27) and Maxwell chain stresses from (28)

Table 1: Overview of finite element models of laminated glass beams with viscoelastic interlayer.

Abbreviation
Kinematics
Constant

FSK
FSν
finite strains
K (2)
ν (2)

VKK
VKν
large deflections
K (2)
ν (2)

LINK
LINν
geometric linearity
K (2)
ν (2)

• experimental data, obtained by López-Aenlle et al. [10] for three simply-supported beams
and one two-span continuous beam with different layer thicknesses;
• results of our earlier finite strain model with elastic interlayer [5], under the assumption of
the constant bulk modulus K (2) ; and
• responses of the two limiting cases, corresponding to the monolithic model, which assumes
the full layer interaction, and the layered model with no inter-layer interaction, both in the
small strain regime.
The remainder of this section is organized as follows. In Section 3.1, we provide dimensions and
boundary conditions for beams used in the analysis and specify material data, loading procedures,
and details on finite element models. Section 3.2 is dedicated to the comparison of the four finite
element models from Table 1. Based on these results, the model combining the von Kármán
kinematics with the time-independent Poisson ratio is verified against a two-dimensional finite
element model and partially validated against experiments in Section 3.3. We proceed with a
parametric study concerning the effects of geometric nonlinearity and temperature, Section 3.4
and conclude our investigations by highlighting the difference in full viscoelastic and simplified
elastic analyses, Section 3.5.
15

3.1. Overview of examples
3.1.1. Geometry and finite element model
The most common laminated beams with three layers (glass/PVB/glass) are considered, either
simply supported or clamped. Note that in practice, the laminated glass elements are not perfectly
fixed or simply supported, but these two cases represent limits of the real support conditions
concerning the effects of geometric non-linearity, e.g. [7, 41]. The layer thicknesses and spans
were selected to ensure that the maximum deflection does not exceed 1/200 of the span, which we
estimated as a suitable serviceability limit for load-bearing glass structures.
Similarly to [5], each layer was discretized with 500 elements and linear stress smoothing was
adopted, in order to reach the four-digit accuracy of displacements and stress values reported
below. Both tolerances in Algorithm 1 were set 1 = 2 = 10−5 .
Table 2: Overview of geometry of laminated glass beams

Geometry

Supports

Length [m]

Width [m]

I
II

fixed-end
simply-supported

3
3

0.15
0.15

Thicknesses [mm]
glass/PVB/glass
3/0.76/3
6/0.76/6

3.1.2. Materials
The elastic constants of glass, the Young modulus E (1) = E (2) = 72 GPa and the Poisson ratio
ν (1) = ν (3) = 0.23, correspond to the values determined by López-Aenlle et al. [10] from a static
bending test. The experimental viscoelastic characterization of PVB was performed by Pelayo et al.
[33] on specimens of thickness 0.38 mm by tensile relaxation tests with the duration of 10 min at
different temperatures ranging from −15 to 50◦ C. The resulting parameters of the generalized
Maxwell chain are summarized in Table 3, the parameters of the Williams-Landel-Ferry (WLF)
equation (8) were set to C1 = 12.6, C2 = 74.46, and T0 = 20 ◦ C. The corresponding shear relaxation
functions G(2) are plotted in Figure 2(b). According to the adopted constitutive assumption, we
further set the bulk modulus K (2) = 2 GPa, e.g. [23, 24], or the Poisson ratio ν (2) = 0.49, e.g. [1,
Section 1.1.3].
Table 3:
Generalized Maxwell chain description of the shear relaxation modulus for PVB, G∞ =
1.9454×10−4 GPa (after [33])

p
1
2
3
4
5
6
7

θp [s]
2.3660×10−7
2.2643×10−6
2.1667×10−5
2.0733×10−4
1.9839×10−3
1.8984×10−2
1.8165×10−1

Gp [GPa]
9.9482×10−2
9.0802×10−2
7.4140×10−2
5.0772×10−2
5.7856×10−2
2.9055×10−2
1.7601×10−2

p
8
9
10
11
12
13

θp [s]
1.7382×100
1.6633×101
1.5916×102
1.5230×103
1.4573×104
1.3945×105

Gp [GPa]
3.0802×10−3
1.2001×10−3
1.1523×10−4
1.8237×10−4
4.1645×10−5
2.2405×10−4

3.1.3. Loading
We considered loading scenarios with spatially uniform distributed loads of intensities shown
in Figure 4. Scenario A, Figure 4 (a), approximates the instantaneous loading reported in [33] by a
16

25

50
LOADING A
40
INTENSITY [N/m]

INTENSITY [N/m]

20

15

10

5

30

20

10

0 -5
-3
-2
0
-4
-1
10 10 10 10
10
10

10

1

10

2

10

3

10

4

10

0 -5
-3
-2
0
-4
-1
10 10 10 10
10
10

5

DURATION OF LOADING [s]

(a)

LOADING B
LOADING C

(b)

10

1

10

2

10

3

10

4

10

5

DURATION OF LOADING [s]

Figure 4: Intensity of (a) one- and (b) two-load history.

loading branch, in which the intensity linearly increases from 0 to 10 Nm−1 in 10−5 s, followed by
the constant load of intensity 10 Nm−1 up to 105 s. The time steps in the incremental algorithms
were distributed uniformly in the logarithmic scale, in particular 6 time steps were used to discretize
the interval h10−6 ; 10−5 i s, whereas the interval h10−5 ; 105 i s was discretized into 24 time steps.
Histories B and C represent loading-unloading scenarios shown in Figure 4 (b), where the steps
in intensity were approximated by linearly varying loading within 10−5 s. The overall loading
program was discretized into 90 time steps uniformly distributed in the logarithmic scale and
refined in the vicinity of jumps.
3.2. Comparison
Comparison of the proposed models and of the effects of geometric nonlinearity are discussed in
this section. Three examples are considered, namely fixed-end beam (geometry I) under one-load
history (loading A), fixed-end beam (geometry I) under two-load history (loading B), and simplysupported beam (geometry II) under two-load history (loading C). The temperature of 25 ◦ C was
kept constant in all examples, as the temperature effects are in the focus of Section 3.4.
In Figures 5–7, we plot the deflections and the maximum of stresses predicted by four geometrically nonlinear (FSK , FSν , VKK , VKν ) and two geometrically linear (LINK , LINν ) approaches.
Their comparison reveals that:
• The viscoelastic approaches based on the assumption of constant bulk modulus K (FSK ,
VKK , and LINK ), or of the Poisson ratio ν (FSν , VKν , resp. LINν ) provide the same results
for all examples because the interlayer mostly transfers shear loading (compare cf. Eq. (25)
and Eq. (28)). The errors in deflections and stresses are much smaller than 0.1%.
• The geometrically nonlinear approaches based on the assumptions of large deflections (VKK ,
VKν ), or finite strains (FSK , FSν ) give again almost identical results for tested examples;
the errors in deflections and stresses do not exceed 0.1%.
• For statically determinate example, the geometrically linear and nonlinear solver provide the
same results, whereas for statically indeterminate examples, the error of linear approach can
reach 100% up to 300%; cf. [5].
17

4.5

16

12

4.0
3.5
STRESS σx [MPa]

DEFLECTION [mm]

14

FSK ~ FSν
VKK ~ VKν
LINK ~ LINν

10
8
6
4

3.0
2.5
2.0
FSK ~ FSν
VKK ~ VKν
LINK ~ LINν

1.5
1.0

2

0.5

0 -6
-5
-3
-2
0
2
3
5
-4
-1
1
4
10 10 10 10 10 10 10 10 10 10 10 10

0.0 -6
-5
-3
-2
0
2
3
5
-4
-1
1
4
10 10 10 10 10 10 10 10 10 10 10 10

DURATION OF LOADING [s]

(a)

DURATION OF LOADING [s]

(b)

Figure 5: Comparison of mid-point (a) deflections and (b) maximum normal stresses for a fixed-end beam with
geometry I under loading A, as predicted by finite element models summarized in Table 1.
16

50

14
12
STRESS σx [MPa]

DEFLECTION [mm]

40

FSK ~ FSν
VKK ~ VKν
LINK ~ LINν

30

20

FSK ~ FSν
VKK ~ VKν
LINK ~ LINν

10
8
6
4

10

2
0 -6
-5
-3
-2
0
2
3
5
-4
-1
1
4
10 10 10 10 10 10 10 10 10 10 10 10

0 -6
-5
-3
-2
0
2
3
5
-4
-1
1
4
10 10 10 10 10 10 10 10 10 10 10 10

DURATION OF LOADING [s]

(a)

DURATION OF LOADING [s]

(b)

Figure 6: Comparison of mid-point (a) deflections and (b) maximum normal stresses for a fixed-end beam with
geometry I under loading B, as predicted by finite element models summarized in Table 1.
16

12

7
6
STRESS σx [MPa]

DEFLECTION [mm]

14

8
FSK ~ FSν
VKK ~ VKν
LINK ~ LINν

10
8
6

(a)

5
4
3

4

2

2

1

0 -6
-5
-3
-2
0
2
3
5
-4
-1
1
4
10 10 10 10 10 10 10 10 10 10 10 10

FSK ~ FSν
VKK ~ VKν
LINK ~ LINν

0 -6
-5
-3
-2
0
2
3
5
-4
-1
1
4
10 10 10 10 10 10 10 10 10 10 10 10

DURATION OF LOADING [s]

(b)

DURATION OF LOADING [s]

Figure 7: Comparison of mid-point (a) deflections and (b) maximum normal stresses for a simply-supported beam
with geometry II under loading C, as predicted by finite element models summarized in Table 1.

18

Recall that in the previous examples, the load level corresponds to maximum deflections about
1/200 of the span of the beam, which might conceal the effects of geometric non-linearity. The
results collected in Table 4 demonstrate that predictions of VK and FS models remain very close
for increased load levels; the difference in deflection is smaller than 0.1% and in stresses than
2%. Notice that for the largest load intensity of 5, 000 Nm−1 in Table 4, the maximum deflection
reaches 1/50 of the span, which exceeds the limit value of 1/65 suggested in [42, Paragraph 9.1.4]
for non-load bearing glass elements.
Table 4: Comparison of the maximum values of deflections and normal stresses at time 105 s

load [N/m]
50
500
5,000

Deflections [mm]
VKK ≈ VKν FSK ≈ FSν
13.239
13.239
30.064
30.067
65.783
65.820

error
-0.00%
-0.01%
-0.06%

Stresses σx [MPa]
VKK ≈ VKν FSK ≈ FSν
22.901
22.919
105.49
105.92
471.62
480.22

error
-0.08%
-0.41%
-1.79%

In conclusion, it follows from the results presented in this section that all four geometrically
non-linear models provide results indistinguishable from the practical design viewpoint. Therefore,
we will report only the results of VKν model in what follows, because its formulation is the most
transparent from the four options.
3.3. Verification and validation
Validation of the numerical model was performed using data from experimental campaign
by López-Aenlle et al. [10], utilizing three types of simply-supported one- or two-span laminated
glass beams (exact dimensions are specified in the caption of Figure 8). The beams were subjected
to seven concentrated loads that approximated uniformly distributed loading and the evolution of
mid-span deflections and of strains was measured with laser sensors and tensometers, respectively.
In addition, the average temperature during the experiment was recorded. In order to generate
verification data in ADINA system, we discretized each beam with a very fine mesh of quadratic
finite strain elements of approximately square shape with edge size of 0.38 mm, and determined
its response under a ramp loading (analogous to the one shown in Figure 4(a)) while assuming the
constant value of the bulk modulus K (2) = 2 GPa.
Results of the simulations and experiments are summarized in Figure 8 and involve the values
of deflections, Figure 8(a), and stresses, Figure 8(b), as functions of loading time. Notice that the
experimental results for stresses appear only in Figure 8(3b), in order to demonstrate that they
are oscillatory and even do not fit in the layered-monolithic bounds. This is probably because the
temperature was not kept constant during the experiment and the polymer interlayer shows very
high sensitivity to temperature variations, as can be deduced from data in Table 3. Therefore, the
experimental values of stresses are excluded from consideration.
We conclude from the collected graphs that our computational model is successfully verified
and only partially validated. Indeed, the results of the detailed two-dimensional finite element
model (2D) are indistinguishable from the results of VKν model, both for displacements and
stresses. The best match between the measured displacements (EXP) is achieved for non-symmetric
4/0.38/8 beam, Figure 8(1a); the other experiments show visible deviations, which we attribute
again to only partial temperature control.
In order to make these claims more quantitative, in Tables 5 and 6 we provide detailed results
of all tested beams at 10 hours. Note that, for example, the VKν -EXP error is understood as
19

4.5

1.5

EXP
VKν
2D
LIMLAY
LIMMON

1.1

STRESS σx [MPa]

DEFLECTION [mm]

1.3

0.9
0.7

4

4

1×10

4

2×10

4

3×10

5×10

1.2

STRESS σx [MPa]

DEFLECTION [mm]

4

4

2×10

4

3×10

4

4×10

5×10

DURATION OF LOADING [s]

VKν
2D
LIMLAY
LIMMON

4.0

1.0
0.8
0.6

3.5
3.0
2.5
2.0

0.4
4

2×10

4

4×10

4

6×10

1.5
0

4

8×10

DURATION OF LOADING [s]

4

2×10

4

4×10

4

6×10

4

8×10

DURATION OF LOADING [s]

(2b)

3.5

6
EXP
VKν
2D
LIMLAY
LIMMON

2.5

5
STRESS σx [MPa]

3.0
DEFLECTION [mm]

4

1×10

4.5
EXP
VKν
2D
LIMLAY
LIMMON

1.4

(2a)

2.5

(1b)

1.6

0.2
0

3.0

1.5
0

4

4×10

DURATION OF LOADING [s]

(1a)

3.5

2.0

0.5
0.3
0

VKν
2D
LIMLAY
LIMMON

4.0

2.0
1.5
1.0

4

EXP
VKν
2D
LIMLAY
LIMMON

3

0.5
0.0
0

(3a)

4

2×10

4

4×10

4

6×10

DURATION OF LOADING [s]

2
0

4

8×10

(3b)

4

2×10

4

4×10

4

6×10

4

8×10

DURATION OF LOADING [s]

Figure 8: Comparison of maximum mid-span (a) deflections and (b) normal stresses for simply supported beams of
span 1 m and of width 0.1 m composed of three glass/PVB/glass layers with thicknesses of (1) 4/0.38/8 mm or (2)
4/0.76/8 mm and (3) a simply supported continuous beam of spans 2 × 0.7 m and of width 0.1 m, composed of three
4/0.38/4 mm layers. The beam is subjected to a ramp load with uniform intensity of (1) 38.25 N/m, (2) 38.25 N/m,
and (3) 94.22 N/m under average temperature of (1) 17.4 ◦ C, (2) 18.3 ◦ C, or (3) 17.8 ◦ C. Response of VKν formulation
is compared with experimental data (EXP), 2D FE model (2D), and monolithic (LIMMON ) and layered (LIMLAY )
limits.

20

Table 5: Comparison of mid-span deflections at 10 hours

thickness [mm]
glass/PVB/glass
4/0.38/8
4/0.76/8
4/0.38/4

EXP
0.8158
0.9947
1.254

deflection [mm]
VKν
2D
LIMMON
0.7839 0.7840
0.4370
0.9234 0.9237
0.4000
1.018
1.021
0.3340

LIMLAY
1.441
1.441
3.068

error [%]
VKν -EXP VKν -2D
-4
-0.01
-7
-0.03
-19
-0.29

Table 6: Comparison of mid-span maximum normal stresses at 10 hours

thickness [mm]
glass/PVB/glass
4/0.38/8
4/0.76/8
4/0.38/4

VKν
2.567
2.846
4.261

stress σx [MPa]
2D
LIMMON LIMLAY
2.567
1.872
3.984
2.847
1.762
3.984
4.266
2.465
5.410

error [%]
VKν -2D
0.00
-0.04
-0.12

(wVKν −wEXP )/wEXP . The results confirm an excellent match between detailed FE model and VKν
model, the errors do not exceed 0.3% in deflections and 0.2% in stresses. The difference between
experimental data and VKν formulation is somewhat less satisfactory, and ranges from about 5–
20%. However, even this lower accuracy is superior to the accuracy of monolithic (LIMMON ) and
layered (LIMLAY ) limits: the monolithic bound predicts values that can be up to 50% smaller than
for the more refined models, whereas the overestimation by the layered limit can reach up to 200%.
3.4. Geometric nonlinearity and temperature
Herein, we complement the verification and validation results with a study on a structure
exhibiting significant effects of geometric nonlinearity. Attention is also paid to the effect of temperature variation, accounted for by the Williams–Landel–Ferry formula (8). In particular, we
consider the fixed-end beam structure from Table 2 subjected to loading A, recall Figure 4(a),
under three constant temperatures 0 ◦ C, 25 ◦ C, and 50 ◦ C.
The resulting time evolution of deflections and stresses is plotted in Figure 9 for geometrically
linear and nonlinear models. The standard upper and lower limits, provided by the geometrically
linear monolithic and layered approximations, are plotted in Figure 9 (a) and they indeed bound
the behavior of the laminated unit. For geometrically nonlinear response, we also provide the
elastic response (FSEL ), determined by elastic model from [5] with the interlayer shear modulus
determined from (12) for the given temperature at the final simulation time of 105 s. This is
complemented with selected data of mid-point deflections, Table 7 and stresses, Table 8.
The results in Figure 9(a) are in full agreement with the outcomes presented in [43], namely
that the laminated units exposed to temperatures around 0 ◦ C effectively behave as the monolithic
ones. In addition, even at 50◦ C, the shear stiffness is sufficient to ensure the inter-layer interaction,
so that the overall response is much closer to the monolithic limit rather than to the layered
limit. Behavior of geometrically non-linear models, Figure 9(b), exhibits similar trends, but the
displacement and stress magnitudes are reduced to 40% and 60%, respectively; see also Tables 7
and 8. In addition, unlike in the geometrically linear model, the extreme stresses at the mid-span
decrease with increasing time. We consider this to be a direct consequence of the development
of membrane stresses in glass layers, which leads to the redistribution of extreme stresses from
the beam mid-section towards the supports (the same mechanism was observed in the detailed 2D
finite element model).
21

8

30

7
DEFLECTION [mm]

DEFLECTION [mm]

40

50°C: LIN ν
25°C: LIN ν ~ 2D
0°C: LIN ν ~ 2D
LIMLAY
LIMMON

20

10

6
5
4
3
2
1

0 -6
-5
-3
-2
0
2
3
5
-4
-1
1
4
10 10 10 10 10 10 10 10 10 10 10 10
DURATION OF LOADING [s]

(1a)

50°C: VK ν
25°C: VK ν ~ 2D
0°C: VK ν ~ 2D
5

FSEL 50°C 10 s
5

FSEL 0°C 10 s

0 -6
-5
-3
-2
0
2
3
5
-4
-1
1
4
10 10 10 10 10 10 10 10 10 10 10 10
DURATION OF LOADING [s]

(1b)
3.5

6

50°C: LIN ν
25°C: LIN ν ~ 2D
0°C: LIN ν ~ 2D
LIMLAY
LIMMON

3.0
STRESS σx [MPa]

STRESS σx [MPa]

8

4

2.5
2.0
1.5
1.0

2

5

0.5

0 -6
-5
-3
-2
0
2
3
5
-4
-1
1
4
10 10 10 10 10 10 10 10 10 10 10 10
DURATION OF LOADING [s]

(2a)

14
12

8
50°C: LIN ν
25°C: LIN ν ~ 2D
0°C: LIN ν ~ 2D
LIMLAY
LIMMON

6

10
8
6

4

2

4

50°C: VK ν
25°C: VK ν ~ 2D
0°C: VK ν ~ 2D
5

FSEL 50°C 10 s
5

FSEL 0°C 10 s

2
0 -6
-5
-3
-2
0
2
3
5
-4
-1
1
4
10 10 10 10 10 10 10 10 10 10 10 10

(3a)

5

FSEL 0°C 10 s

DURATION OF LOADING [s]

(2b)

STRESS σx [MPa]

STRESS σx [MPa]

16

FSEL 50°C 10 s

0.0 -6
-5
-3
-2
0
2
3
5
-4
-1
1
4
10 10 10 10 10 10 10 10 10 10 10 10

20
18

50°C: VK ν
25°C: VK ν ~ 2D
0°C: VK ν ~ 2D

DURATION OF LOADING [s]

0 -6
-5
-3
-2
0
2
3
5
-4
-1
1
4
10 10 10 10 10 10 10 10 10 10 10 10

(3b)

DURATION OF LOADING [s]

Figure 9: Fixed-end beam with geometry I under loading A: Mid-point (1) deflections, (2) maximum normal stresses,
and (3) the largest values of maximum normal stresses at temperatures of 0, 25, and 50 ◦ C. Response of (a) geometrically linear (LINν ) and (b) nonlinear (VKν ) beam models are compared with detailed model in ADINA (2D),
monolithic (LIMMON ) and layered (LIMLAY ) limits, and elastic finite strain model (FSEL ), using the shear modulus
of the interlayer corresponding to temperatures from 0 to 50 ◦ C and duration of loading 105 s.

22

Table 7: Comparison of mid-span deflections at 105 s. The symbol – indicates divergence of ADINA solver.

temperature [ ◦ C]
0
25
50

linear response [mm]
LINν
2D(LIN)
8.192
8.191
16.15
16.15
16.63
–

nonlinear response [mm]
VKν
2D
FSEL
5.596 5.595
5.701
6.838 6.838
6.857
6.863
–
6.863

error [%]
VKν -2D FSEL -VKν
0.02
1.88
0.00
0.28
–
0.00

Table 8: Comparison of maximum mid-span normal stresses σx at 105 s. The symbol – indicates divergence of
ADINA solver.

temperature [ ◦ C]
0
25
50

linear response [MPa]
LINν
2D(LIN)
3.332
3.332
4.170
4.170
4.237
–

nonlinear response [MPa]
VKν
2D
FSEL
2.724 2.724
2.706
2.437 2.438
2.433
2.431
–
2.431

error [%]
VKν -2D FSEL -VKν
0.00
-0.66
-0.04
-0.21
–
0.00

Furthermore, the results in this section serve as an additional verification of VKν model, because
the errors in displacements and stresses do not exceed 0.5% for temperatures of 0 and 25◦ C with
respect to the results of detailed 2D finite element model, consult again Tables 7 and 8. Note that
simulations in ADINA did not convergence at 50◦ C, because of poor conditioning of the stiffness
matrix caused by a very high contrast in the interlayer and glass properties. In addition, the results
reveal that the elastic model provides very good estimates at the final time of simulations: the
errors stay below 2% for all quantities of interest. This aspect will be studied in more detail in the
next section.
3.5. Comparison with simplified elastic solution
Being inspired by earlier studies by Galuppi and Royer-Carfagni [19], we close this section by
comparing the full viscoelastic solution with a simplified ”secant” solution, in which we determine
the structural response at the relevant times by the geometrically non-linear model [44] with
shear modulus determined from (12). Similarly to [19], we consider the two-load history B, see
Figure 4(b), consisting of a sudden loading, sudden unloading, and holding parts. In addition,
both simply supported and fixed-end beams are considered under constant temperatures of 25 and
50◦ C.
The simulation results collected in Figures 10 and 11 reveal that the correspondence between the
simplified model and full viscoelastic solution is rather satisfactory for holding phases of the loading,
but differences appear in the vicinity of displacement jumps, where the simplified solution can lead
to unsafe predictions. This effects is the most pronounced for stresses in the fixed-end beam and it
is attributed to a different stress redistribution in the viscoelastic and elastic beams as discussed in
the previous section. In addition, the magnitude seems to be highly temperature-dependent: the
differences are negligible at 0◦ C (not shown) and increase with increasing temperature.
More quantitatively, for the fixed-end beam, Figure 10, the differences remain below 5%; at
25◦ C they do not exceed 2% for both stresses and strains, whereas for 50◦ C they reach 4% for
displacements and 5% for stresses. For the simply supported beam, the differences between elastic
and viscoelastic solution increase and reach up to 27% for displacements and 14% for stresses.
These results again highlight the complex interplay between structural heterogeneity, geometric
nonlinearity, and viscoelastic effects.
23

13

6.0
STRESS σx [MPa]

DEFLECTION [mm]

12

6.5
FSEL
VKν

11
10
9
8

1

10

2

10

3

10

4

10

5.0
4.5

3.5 -5
-3
-2
0
-4
-1
10 10 10 10
10
10

5

10

DURATION OF LOADING [s]

(1a)

1

10

2

10

3

10

4

10

5

10

DURATION OF LOADING [s]

(2a)

13

6.0
FSEL
VKν

5.5
STRESS σx [MPa]

DEFLECTION [mm]

5.5

4.0

7 -5
-3
-2
0
-4
-1
10 10 10 10
10
10

12

FSEL
VKν

11
10
9
8

5.0
4.5
4.0
3.5

7 -5
-3
-2
0
-4
-1
10 10 10 10
10
10

(1b)

FSEL
VKν

1

10

2

10

DURATION OF LOADING [s]

3

10

4

10

3.0 -5
-3
-2
0
-4
-1
10 10 10 10
10
10

5

10

(2b)

1

10

2

10

3

10

4

10

5

10

DURATION OF LOADING [s]

Figure 10: Fixed-end beam with geometry I under loading B: Mid-point (1) deflections or (2) maximum normal
stresses at temperatures of (a) 25 ◦ C and (b) 50 ◦ C. Response of geometrically nonlinear VKν models is compared
with simplified elastic model (FSEL ).

Note in conclusion that all these findings agree well with the analytical study by Galuppi and
Royer-Carfagni [19] for a simply supported beam with an interlayer modeled using simplified oneand three-unit Maxwell chains, disregarding the temperature effects.
4. Conclusions
In this paper, we have introduced and thoroughly tested four layer-wise finite element formulations to determine displacements and stresses in laminated glass beams with interlayer(s)
exhibiting time/temperature-dependent material response. The models differ in the adopted beam
kinematics — large-deflection von Karmán formulation vs. finite strain Reissner formulation — and
the assumptions adopted in the constitutive model — constant bulk modulus vs. constant Poisson ratio. By comparing their response with detailed two-dimensional finite element simulations,
experimental data, and simplified approaches, we have reached the following conclusions:
1. All four considered models provide almost identical response under monotone and nonmonotone loading for deflections up to 1/50 of span. Therefore, the most straightforward
24

6.5

16

6.0
5.5
STRESS σx [MPa]

DEFLECTION [mm]

14

FSEL
VKν

12
10
8

FSEL
VKν

5.0
4.5
4.0
3.5
3.0

6
2.5
4 -5
-3
-2
0
-4
-1
10 10 10 10
10
10

1

10

2

10

3

10

4

10

2.0 -5
-3
-2
0
-4
-1
10 10 10 10
10
10

5

10

DURATION OF LOADING [s]

(1a)

3

10

4

10

5

10

7.0
FSEL
VKν

6.0

16
STRESS σx [MPa]

DEFLECTION [mm]

2

10

DURATION OF LOADING [s]

(2a)

20
18

1

10

14
12
10
8

FSEL
VKν

5.0

4.0

3.0

6
4 -5
-3
-2
0
-4
-1
10 10 10 10
10
10

(1b)

1

10

2

10

DURATION OF LOADING [s]

3

10

4

10

2.0 -5
-3
-2
0
-4
-1
10 10 10 10
10
10

5

10

(2b)

1

10

2

10

3

10

4

10

5

10

DURATION OF LOADING [s]

Figure 11: Simply-supported beam with geometry II under loading C: Mid-point (1) deflections or (2) maximum
normal stresses at temperatures of (a) 25 ◦ C and (b) 50 ◦ C. Response of geometrically nonlinear VKν models is
compared with simplified elastic model (FSEL ).

formulation, combining the assumptions of the von Kármán kinematics and of constant Poisson ratio, can be used to simulate more complex structures.
2. This formulation has been successfully verified against several detailed finite element simulations performed in ADINA system and partially validated against experimental data collected
by López-Aenlle et al. [10]. The model has also consistently reproduced typical features of
mechanical response of laminated glass structures at different temperatures.
3. A simplified model, based on an elastic formulation with the interlayer shear modulus adjusted
to duration of loading and temperature, provides satisfactory agreement with full viscoelastic
solutions at low temperatures or at holding periods of loading. At room and elevated temperatures, differences have been observed for rapid changes of loading and unloading that
may reach the order of 10% and/or may be on the unsafe side.
As the next step, we will enhance our elastic finite element solver for laminated glass plates [28]
with viscoelastic effects, closely following the formulations developed in this paper. Subsequently,
we would like to extend the formulation to post-critical state, by incorporating glass fracture into
the layer-wise formulation.
25

Acknowledgments. This work was initiated when Alena Zemanová was a post-doctoral researcher,
supported by the project No. CZ.1.07/2.3.00/30.0034: Support of inter-sectoral mobility and quality enhancement of research teams at Czech Technical University in Prague (06/2014–06/2015).
All authors acknowledge the support by the Czech Science Foundation, project No. 16-14770S.
References
References
[1] M. Haldimann, A. Luible, M. Overend, Structural Use of Glass, vol. 10 of Structural Engineering Documents,
IABSE, Zürich, Switzerland, 2008.
[2] T. Serafinavičius, J. P. Lebet, C. Louter, T. Lenkimas, A. Kuranovas, Long-term laminated glass four point
bending test with PVB, EVA and SG interlayers at different temperatures, Procedia Engineering 57 (2013)
996–1004, doi:10.1016/j.proeng.2013.04.126.
[3] X. Huang, G. Liu, Q. Liu, S. Bennison, An experimental study on the flexural performance of laminated glass,
Structural Engineering and Mechanics 49 (2) (2014) 261–271, doi:10.12989/sem.2014.49.2.261.
[4] 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 (2014) 59–67, doi:10.1016/
j.polymertesting.2014.04.011.
[5] A. Zemanová, J. Zeman, M. Šejnoha, Numerical model of elastic laminated glass beams under finite strain,
Archives of Civil and Mechanical Engineering 14 (4) (2014) 734–744, doi:10.1016/j.acme.2014.03.005.
[6] R. A. Behr, J. E. Minor, M. P. Linden, C. V. G. Vallabhan, Laminated glass units under uniform lateral pressure,
Journal of Structural Engineering 111 (5) (1985) 1037–1050, doi:10.1061/(ASCE)0733-9445(1985)111:5(1037).
[7] Y. Koutsawa, E. Daya, Static and free vibration analysis of laminated glass beam on viscoelastic supports,
International Journal of Solids and Structures 44 (2007) 8735–8750, doi:10.1016/j.ijsolstr.2007.07.009.
[8] S. J. Benninson, M. H. X. Qin, P. S. Davis, High-performance laminated glass for structurally efficient glazing,
Innovative light-weight structures and sustainable facades, Hong Kong .
[9] L. Galuppi, G. F. Royer-Carfagni, Effective thickness of laminated glass beams: New expression via a variational
approach, Engineering Structures 38 (2012) 53–67, doi:10.1016/j.engstruct.2011.12.039.
[10] M. López-Aenlle, F. Pelayo, A. Fernández-Canteli, M. A. Garcı́a Prieto, The effective-thickness concept in
laminated-glass elements under static loading, Engineering Structures 56 (2013) 1092–1102, doi:10.1016/j.
engstruct.2013.06.018.
[11] J. Eisentraeger, K. Naumenko, H. Altenbach, H. Koeppe, Application of the first-order shear deformation theory
to the analysis of laminated glasses and photovoltaic panels, International Journal of Mechanical Sciences 96–97
(2015) 163–171, doi:10.1016/j.ijmecsci.2015.03.012.
[12] J. Hooper, On the bending of architectural laminated glass, International Journal of Mechanical Sciences 15 (4)
(1973) 309–323, doi:10.1016/0020-7403(73)90012-X.
[13] H. S. Norville, K. W. King, J. L. Swofford, Behavior and strength of laminated glass, Journal of Engineering
Mechanics 124 (1) (1998) 46–53, doi:10.1061/(ASCE)0733-9399(1998)124:1(46).
[14] I. V. Ivanov, Analysis, modelling, and optimization of laminated glasses as plane beam, International Journal
of Solids and Structures 43 (22–23) (2006) 6887–6907, doi:10.1016/j.ijsolstr.2006.02.014.
[15] M. Z. Aşık, S. Tezcan, A mathematical model for the behavior of laminated glass beams, Computers & Structures
83 (21–22) (2005) 1742–1753, doi:10.1016/j.compstruc.2005.02.020.
[16] S.-H. Schulze, M. Pander, K. Naumenko, H. Altenbach, Analysis of laminated glass beams for photovoltaic applications, International Journal of Solids and Structures 49 (15–16) (2012) 2027–2036, doi:10.1016/j.ijsolstr.
2012.03.028.
[17] J. Jaśkowiec, P. Pluciński, J. Pamin, Thermo-mechanical XFEM-type modeling of laminated structure with
thin inner layer, Engineering Structures 100 (2015) 511–521, doi:10.1016/j.engstruct.2015.06.035.
[18] L. Galuppi, G. F. Royer-Carfagni, Laminated beams with viscoelastic interlayer, International Journal of Solids
and Structures 49 (18) (2012) 2637–2645, doi:10.1016/j.ijsolstr.2012.05.028.
[19] L. Galuppi, G. Royer-Carfagni, The design of laminated glass under time-dependent loading, International
Journal of Mechanical Sciences 68 (2013) 67–75, doi:10.1016/j.ijmecsci.2012.12.019.
[20] L. Galuppi, G. Royer-Carfagni, Buckling of three-layered composite beams with viscoelastic interaction, Composite Structures 107 (2014) 512–521, doi:10.1016/j.compstruct.2013.08.006.

26

[21] L. Galuppi, G. Royer-Carfagni, Localized contacts, stress concentrations and transient states in bent-lamination
with viscoelastic adhesion. An analytical study, International Journal of Mechanical Sciences 103 (2015) 275–287,
doi:10.1016/j.ijmecsci.2015.08.016.
[22] P. Wu, D. Zhou, W. Liu, L. Wan, D. Liu, Elasticity solution of two-layer beam with a viscoelastic interlayer
considering memory effect, International Journal of Solids and Structures 94–95 (2016) 76–86, doi:10.1016/j.
ijsolstr.2016.05.007.
[23] A. V. Duser, A. Jagota, S. J. Bennison, Analysis of glass/polyvinyl butyral laminates subjected to uniform
pressure, Journal of Engineering Mechanics 125 (4) (1999) 435–442, doi:10.1061/(ASCE)0733-9399(1999)125:
4(435).
[24] S. J. Bennison, A. Jagota, C. A. Smith, Fracture of Glass/Poly(vinyl butyral) (Butacite R ) Laminates in Biaxial
Flexure, Journal of the American Ceramic Society 82 (7) (1999) 1761–1770, doi:10.1111/j.1151-2916.1999.
tb01997.x.
[25] C. Bedon, J. Belis, A. Luible, Assessment of existing analytical models for the lateral torsional buckling analysis
of PVB and SG laminated glass beams via viscoelastic simulations and experiments, Engineering Structures 60
(2014) 52–67, doi:10.1016/j.engstruct.2013.12.012.
[26] S. T. Mau, A refined laminated plate theory, Journal of Applied Mechanics–Transactions of the ASME 40 (2)
(1973) 606–607, doi:10.1115/1.3423032.
[27] O. C. Zienkiewicz, M. Watson, I. P. King, A numerical method of visco-elastic stress analysis, International
Journal of Mechanical Sciences 10 (10) (1968) 807–827, doi:10.1016/0020-7403(68)90022-2.
[28] A. Zemanová, J. Zeman, M. Šejnoha, Finite element model based on refined plate theories for laminated glass
units, Latin American Journal of Solids and Structures 15 (6) (2015) 1158–1180, doi:10.1590/1679-78251676.
[29] J. Kruis, J. Zeman, P. Gruber, Model of Imperfect Interfaces in Composite Materials and Its Numerical Solution
by FETI Method, in: R. Bank, M. Holst, O. Widlund, J. Xu (Eds.), Domain Decomposition Methods in Science
and Engineering XX, no. 91 in Lecture Notes in Computational Science and Engineering, Springer Berlin
Heidelberg, 337–344, doi:10.1007/978-3-642-35275-1_39, 2013.
[30] M. Larcher, G. Solomos, F. Casadei, N. Gebbeken, Experimental and numerical investigations of laminated
glass subjected to blast loading, International Journal of Impact Engineering 39 (1) (2012) 42–50, doi:10.1016/
j.ijimpeng.2011.09.006.
[31] E. Reissner, On One-Dimensional Finite-Strain Beam Theory: the Plane Problem, Journal of Applied Mathematics and Physics 23 (5) (1972) 795–804, doi:10.1007/BF01602645.
[32] J. N. Reddy, An Introduction to Nonlinear Finite Element Analysis, Oxford University Press, doi:10.1093/
acprof:oso/9780198525295.001.0001, 2004.
[33] F. Pelayo, M. López-Aenlle, L. Hermans, A. Fraile, Modal Scaling of a Laminated Glass Plate, 5th International
Operational Modal Analysis Conference (2013) 1–10URL http://iomac.eu/iomac/2013/IOMAC_Guimaraes/
files/papers/Model%20Updating/Paper-175-Fernandes.pdf.
[34] M. Williams, R. Landel, J. Ferry, The Temperature Dependence of Relaxation Mechanisms in Amorphous
Polymers and Other Glass-forming Liquids, Journal of the American Chemical Society 77 (14) (1955) 3701–
3707, doi:10.1021/ja01619a008.
[35] R. Christensen, Theory of Viscoelasticity: An Introduction, Elsevier, second edn., 1982.
[36] M. Jirásek, Z. P. Bažant, Inelastic analysis of structures, John Wiley & Sons, Ltd., 2002.
[37] I. S. Sokolnikoff, Mathematical theory of elasticity, McGraw-Hill Book Company, Inc., New York, Toronto,
London, second edn., 1956.
[38] V. Birman, C. Bert, On the Choice of Shear Correction Factor in Sandwich Structures, Journal of Sandwich
Structures & Materials 4 (1) (2002) 83–95, doi:10.1177/1099636202004001180.
[39] J. F. Bonnans, J. C. Gilbert, C. Lemaréchal, C. A. Sagastizábal, Numerical Optimization: Theoretical and
Practical Aspects, Springer, doi:10.1007/978-3-540-35447-5, 2003.
[40] A. Zemanová, J. Zeman, M. Šejnoha, Simple numerical model of laminated glass beams, Acta Polytechnica
48 (6) (2008) 22–26.
[41] M. Z. Aşık, Laminated glass plates: revealing of nonlinear behavior, Computers & Structures 81 (28–29) (2003)
2659–2671, doi:10.1016/S0045-7949(03)00325-0.
[42] prEN 16612, Glass in building - Determination of the load resistance of glass panes by calculation and testing,
European Strandard, European Committee for Standardization, 2013.
[43] R. A. Behr, J. E. Minor, H. S. Norville, Structural Behavior of Architectural Laminated Glass, Journal of
Structural Engineering 119 (1) (1993) 202–222, doi:10.1061/(ASCE)0733-9445(1993)119:1(202).
[44] A. Zemanová, Numerical modeling of laminated glass structures, Ph.D. thesis, Faculty of Civil Engineering,
CTU, Prague, 2014.

27

A. Numerical aspects
The finite element discretization is accomplished with two-node elements with linear basis
functions for each layer, and all integrals are evaluated by the one-point quadrature to mitigate
locking effects. Therefore, the kinematics for the e-th element of the i-th layer is specified by the
column matrix of nodal displacements and rotations
h
i
(i)
(i)
(i)
(i)
(i)
(i)
d(i)
=
u0,e,1 w0,e,1 ϕy,e,1 u0,e,2 w0,e,2 ϕy,e,2 T ,
e
(i)

(i)

where, e.g., u0,e,1 stands for the centerline displacements of the left node and ϕy,e,2 denotes the
rotation of the cross-section at the right node. These choices enables us to express, in the explicit
form, all matrices introduced in Eqs. (35) and (36), as shown in the remainder of this section.
A.1. von Kármán beams
Nodal forces. The matrix of internal nodal forces follows by differentiating the element internal energy function with respect to the nodal displacements, recall (35). For the von Kármán formulation,
we obtain
(i)

(i)

fint,e,1 = −E (i) A(i) εx,0,e ,
(i)
fint,e,2
(i)

fint,e,3

=

(i)
w0,e,2
(i)
fint,e,1

(i)

(i)

(A.1a)

(i)

(i)

(A.1b)

fint,e,4 = −fint,e,1 ,
(i)

− w0,e,1

(i)
− G(i) A(i)
s γxz,e ,

Le
(i) Le
− E (i) Iy(i) κ(i)
= G(i) A(i)
y,e ,
s γxz,e
2

fint,e,5 = −fint,e,2 ,
(i)

(i)

(i)
fint,e,6 = fint,e,3 + 2E (i) Iy(i) κy,e
,

(A.1c)

(i)

where we have introduced the generalized element strains associated with de :
(i)

(i)

εx,0,e =
κ(i)
y,e =

(i)

u0,e,2 − u0,e,1
Le
(i)
ϕy,e,2

−
Le

(i)
ϕy,e,1

+

1
2

 w(i) − w(i) 2
0,e,2
0,e,1
Le

(i)
γxz,e
=

,

1
2



 w(i) − w(i)
0,e,2
0,e,1
(i)
(i)
,
ϕy,e,1 + ϕy,e,2 +
Le

,

b (2)
according to (6) and (7). Note that Eq. (A.1) in fact holds only for i = 1, 3, the matrix f
int,e
related to the interlayer must be obtained with effective values of the Young and shear moduli,
b (2) and G
b (2) , introduced in Section 2.2.
E
In order to derive the expression for the additional contribution to the nodal force matrix
(2)
b
∆f
, it is convenient to introduce the following history variables
int,e

(2)
(2)
b (2) − E
b (2) A(2) ε(2) (tn ),
δNx,e
(tn ) = Nx,e
(tn ) + ∆N
x,e
x,0,e
(2)
(2)
(2)
b (2) A(2) γ (2) (tn ),
δVz,e
(tn ) = Vz,e
(tn ) + ∆Vbz,e
−G
s
xz,e
(2)
(2)
(2)
cy,e
b (2) Iy(2) κ(2)
δMy,e
(tn ) = My,e
(tn ) + ∆M
−E
y,e (tn ),

28

that account for initial internal forces and generalized strains at the beginning of time step tn and
b (2) receives a form similar to Eq. (A.1):
for the relaxation effects. Using these quantities, ∆f
int,e

(2)

(2)
(2)
∆fbint,e,4 = −∆fbint,e,1 ,

(2)
∆fbint,e,1 = −δNx,e
(tn ),
(2)

(2)
(2)
∆fbint,e,2 = ∆fbint,e,1

(2)

w0,e,2 − w0,e,1

(2)

(2)
∆fbint,e,3 = δVz,e
(tn )

Le

(2)
(2)
∆fbint,e,5 = −∆fbint,e,2 ,

(2)
− δVz,e
(tn ),

Le
(2)
− δMy,e
(tn )
2

(2)
(2)
(2)
∆fbint,e,6 = ∆fbint,e,3 + 2δMy,e
(tn ).

Stiffness matrices. Upon differentiating nodal forces with respect to nodal displacements, we obtain
the stiffness matrix in the form


(i)
(i)
(i)
(i)
Kt,e,11
Kt,e,12
0
−Kt,e,11 −Kt,e,12
0


(i)
(i)
(i)
(i)
(i)
(i)
 Kt,e,12

K
K
−K
−K
K
t,e,22
t,e,23
t,e,12
t,e,22
t,e,23


(i)
(i)
(i)
(i)

0
Kt,e,23
Kt,e,33
0
−Kt,e,23 Kt,e,36 


(i)
K t,e = 
,
(i)
(i)
(i)
(i)
 −Kt,e,11

−K
0
K
K
0
t,e,12
t,e,11
t,e,12


(i)
(i)
(i)
(i)
(i)
 −K (i)

−K
−K
K
K
−K

t,e,12
t,e,22
t,e,23
t,e,12
t,e,22
t,e,23 
(i)
(i)
(i)
(i)
0
Kt,e,23
Kt,e,36
0
−Kt,e,23 Kt,e,33
where the individual entries read
Kt,e,11 =

E (i) A(i)
,
Le

(i)
Kt,e,22

(i)
(i)
Kt,e,11 εx,0,e

(i)

=

(i)

(i)

Kt,e,12 = E (i) A(i)
+

(i)
(i) w0,e,2
Kt,e,12

(i)

− w0,e,1
Le

,

(i)

Kt,e,23 = − 12 G(i) A(i)
s ,
(i)

E (i) Iy
,
Le

(i)

(i)

Kt,e,36 = Kt,e,33 −

Likewise, the history-dependent increment of the
the tangent stiffness matrix

0
0
(2)
 0 ∆K
b

t,e,22

0
0
c(2) = 
∆K

t,e
0
 0

b (2)
 0 −∆K
t,e,22
0
0
(2)

L2e

(i)

G(i) As
,
Le

+

(i)

(i)

Kt,e,33 = 41 G(i) A(i)
s Le +

(i)

w0,e,2 − w0,e,1

2E (i) Iy
.
Le

nodal forces yields an additional contribution to
0
0
0
0
0
0

0
0
(2)
b
0 −∆Kt,e,22
0
0
0
0
(2)
b
0 ∆Kt,e,22
0
0

0
0
0
0
0
0






,




(2)

b
with ∆K
t,e,22 = δNx,e (tn )/Le .
Compatibility. The inter-layer compatibility (5) conditions, written for the j-th node for the interface between the i-th and (i + 1)-th layer, receives the form
" (i)
#
(i+1)
(i)
(i+1)
u0,j − u0,j + 21 h(i) ϕy,j + 12 h(i+1) ϕy,j
(i,i+1)
(i,i+1)
cj
=
= Cj
d,
(i)
(i+1)
w0,j − w0,j
29

from which we obtain the corresponding block of the constraining matrix


1 0 12 h(i) · · · −1 0 12 h(i+1)
(i,i+1)
Cj
=
.
0 1
0
· · · 0 −1
0
Note that this expression coincides with the geometrically linear formulation developed by Zemanová et al. [40].
A.2. Reissner beams
(i)

Nodal forces. Explicit expressions for nodal forces f int,e in the elastic formulation are available
in [5, Eq. (37)]. Therefore, we provide here only the history-dependent contribution:
(2)

(2)
(2)
∆fbint,e,1 = −δNx,e
(tn ) cos βe(2) − δVz,e
(tn ) sin βe(2) ,

(2)
(2)
∆fbint,e,4 = −∆fbint,e,1 ,

(2)
(2)
(2)
(2)
(2)
∆fbint,e,2 = δNx,e
(tn ) sin βe(2) − δVz,e
(tn ) cos βe(2) , ∆fbint,e,5 = −∆fbint,e,2 ,




(2)
(2)
(2)
(2)
(2)
(2)
1
b(2)
b(2)
∆fbint,e,3 = − 21 L(2)
e + u0,e,2 − u0,e,1 ∆fint,e,2 + 2 w0,e,2 − w0,e,1 ∆fint,e,1 − δMy,e (tn ),
(2)

(2)

(2)
∆fbint,e,6 = ∆fbint,e,3 + 2δMy,e
(tn ),
(i)

(i)

(i)

with βe = 12 (ϕy,e,1 + ϕy,e,2 ).
Stiffness matrices. Again, the expressions for the tangent stiffness matrix K t and the contribution
of the Lagrange multipliers K λ are available in [5, Eqs. (38) and (40)], therefore they are omitted
here for the sake of brevity. The history-dependent contribution attains the form


b (2)
b (2)
0
0
∆K
0
0
∆K
t,e,13
t,e,13


b (2)
b (2)


0
0
∆K
0
0
∆K
t,e,23
t,e,23 

(2)
(2)
(2)
(2)
(2)
(2)
 ∆K

b
b
b
b
b
b
∆K
−∆K
∆K
t,e,13 ∆Kt,e,23
t,e,33
t,e,13 −∆Kt,e,23
t,e,33 
c(2) = 
∆K

,
t,e
b (2)
b (i) 

0
0
−∆K
0
0
−∆K
t,e,13
t,e,13 


b (2)
b (i) 
0
0
−∆K
0
0
−∆K

t,e,23
t,e,23 
(2)
(2)
(2)
(2)
(2)
(2)
b
b
b
b
b
b
∆Kt,e,13 ∆Kt,e,23 ∆Kt,e,33 −∆Kt,e,13 −∆Kt,e,23 ∆Kt,e,33
where
(2)


1  (2)
(2)
δNx,e (tn ) sin βe(2) − δVz,e
(tn ) cos βe(2) ,
2

1  (2)
(2)
=
δNx,e (tn ) cos βe(2) + δVz,e
(tn ) sin βe(2) ,
2




1   (2)
(2)
(2)
b (2) + w(2) − w(2) ∆K
b (2)
=
− Le + u0,e,2 − u0,e,1 ∆K
t,e,23
0,e,2
0,e,1
t,e,13 .
2

b
∆K
t,e,13 =
(2)

b
∆K
t,e,23
(2)

b
∆K
t,e,33

(i,i+1)

Compatibility. Expressions for the nodal block of the compatibility condition cj
dients

(i,i+1)
Cj

can be found in [5, Eqs. (23) and (39)].

30

and its gra-

