AP08_6.vp


1 Introduction
Glass is the most frequently used transparent material in

the building envelopes. It is a fragile material, which fails in
a brittle manner. For this reason, safety glasses are used in a
situation when there is a possibility of human impact or where
the glass could fall if shattered.

Laminated glass is a multi-layer material produced by
bonding two or more layers of glass together with a plastic
interlayer, typically made of polyvinyl butyral (PVB). The
interlayer keeps the layers of glass bonded even when broken,
and its high strength prevents the glass from breaking up
into large sharp pieces. This produces a characteristic “spi-
der web” cracking pattern when the impact is not powerful
enough to completely pierce the glass. Multiple laminae and
thicker glass decrease the stress level, thereby also increasing
the load-carrying capacity of the structural member.

The focus of this study is on the establishing a simple and
versatile framework for an analysis of the mechanical behav-
ior of laminated glass units. To keep the discussion compact,
we restrict our attention to the linearly elastic response of
layered glass beams in the small strain regime. The rest of the
paper is organized as follows. Methods for an analysis of lami-
nated glass beams are introduced in Section 2, together with a
brief characterization of the proposed numerical model. The
principles of the method are described in detail in Sections 3
and 4. In particular, the mechanical formulation of the model
is shown in Section 3. Finite Element discretization is pre-
sented in Section 4. In Section 5, the proposed numerical
technique is verified and validated against a reference analyti-
cal solution and publicly available experimental data. Finally,
Section 6 concludes the paper and discusses future extensions
of the method.

2 Brief overview of available methods
The most frequent approach to the analysis of glass

structural elements was, for a long time, based on empirical
knowledge. Such relations are sufficient for the design of tra-
ditional windows glasses. However, in modern architecture
there has been a steadily growing demand in recent decades
for transparent materials for large external walls and roof
systems. Therefore, a detailed analysis of layered glass units is
becoming increasingly important in order to ensure reliable
and efficient design.

In general, the complex behavior of laminated glass
can be considered as an intermediate state of two limiting

cases [1]. In the first case, the structure is treated as an assem-
bly of two independent glass plates without any interlayer (the
lower bound on stiffness and strength of a member), while
in the second case, corresponding to the upper estimate of
strength and stiffness, the glass unit is modeled as monolithic
glass (one glass plate equal in thickness to the total thickness
of the glass plates). Both elementary cases, however, fail to
correctly capture the complex interaction among the indi-
vidual layers, leading to non-optimal layer thickness designs.
Therefore, several alternative approaches to the analysis of
layered glass structures have been proposed in the literature.
These methods can be categorized into three basic groups:

� methods calibrated with respect to experimental measure-
ments [2],

� analytical approaches [3, 4, 5],

� numerical models typically based on detailed Finite Ele-
ment simulations [6, 7].

The applicability of analytical approaches to practical
(usually large-scale) structures is far from being straightfor-
ward. In particular, the closed-form solution of the resulting
equations is possible only for very specific boundary con-
ditions, and the equations therefore have to be solved by
an appropriate numerical method. Moreover, the analytical
approaches are rather difficult to generalize to beams with
multiple layers. Therefore, it appears to be advantageous
to formulate the problem directly in the discretized form,
typically based on the Finite Element Method (FEM). Never-
theless, we would like to avoid fully resolved two- or three-di-
mensional simulations, cf. [6, 7], which lead to unnecessarily
expensive calculations.

In this paper, we propose a simple FEM model inspired by
a specific class of refined plate theories [8, 9, 10]. In this
framework, each layer is treated as a Timoshenko beam with
independent kinematics. The interaction between the indi-
vidual layers is captured by the Lagrange multipliers (with a
physical meaning of shear stresses), which result from the con-
ditions of the compatibility of the inter-layer displacements.
Such a refined approach circumvents the limitation of similar
models available in typical commercial FEM systems, which
employ a single set of kinematic variables and average the
mechanical response through the thickness of the beam,
e.g. [11]. Unlike the proposed approach, the averaging
operation is too coarse to correctly represent the inter-layer
interactions, see Section 5 for a concrete example.

22 ©  Czech Technical University Publishing House http://ctn.cvut.cz/ap/

Acta Polytechnica Vol. 48  No. 6/2008

Simple Numerical Model of Laminated
Glass Beams
A. Zemanová, J. Zeman, M. Šejnoha

This paper presents a simple Finite Element model aimed at efficient simulation of layered glass units. The approach is based on considering
the independent kinematics of each layer, tied together via Lagrange multipliers. Validation and verification of the resulting model against
independent data demonstrate its accuracy, showing its potential for generalization towards more complex problems.

Keywords: laminated glass beams, finite element method, Lagrange multipliers.



3 Mechanical model of laminated
beams
As illustrated in Fig. 1, laminated glasses consist mostly of

three layers. A local coordinate system is attached to each
layer to allow for an efficient treatment of independent kine-
matics. In the following text, a quantity a expressed in the lo-
cal coordinate system associated with the i-th layer is denoted
as a i( ), whereas a variable without an index represents a glob-
ally defined quantity, cf. Fig. 1.

Each layer is modeled using the Timoshenko beam theory
supplemented with membrane effects. Hence, the following
kinematic assumptions are adopted
� the cross sections remain planar but not necessarily per-

pendicular to the deformed beam axis,
� the vertical displacement does not vary along the height

of the beam (when compared to the magnitude of the
displacement).

Under these assumptions, the non-zero displacement
components in each layer are parametrized as:

u x z u x x zi i i i i( ) ( ) ( ) ( ) ( )( , ) ( , ) ( )� �0 � ,
w x z w xi i( ) ( )( , ) ( )� ,

where i � 1 2 3, , and z i( ) is measured in the local coordinate
system from the middle plane of the i-th layer. The inter-layer
interaction is ensured via the continuity conditions specified
on the interfaces between the layers in the form (i � 1 2, )

u x
h

u x
hi

i
i

i
( )

( )
( )

( )
( , ) ( , )

2 2
01

1
� � ��

�
. (1)

The strain field in the i-th layer follows from the strain-dis-
placement relations [12, 11]

�
�

�

�
x
i i

i
i

i i
x z

u
x

x z
u

x
x

x
( ) ( )

( )
( )

( ) ( )
( , ) ( , ) ( , )� � �

d
d

d
d

0 ( ) ( )x z i ,

�
�

�

�

�
�xz

i
i

i
i ix

u
z

x z
w
x

x x
w
x

( )
( )

( )
( ) ( )( ) ( , ) ( ) ( ) (� � � �

d
d

x) , (2)

which, when combined with the constitutive equations of each
layer expressed in terms of Young’s modulus E and the shear
modulus G:

� �x
i i i

x
i ix z E x z( ) ( ) ( ) ( ) ( )( , ) ( , )� and � �xz

i i
xz
ix G x( ) ( ) ( )( ) ( )� ,

yield the expressions for the internal forces as

N x E A
u

x
xx

i i i
i

( ) ( ) ( )
( )

( ) ( , )�
d

d
0 ,

V x kG A x
w
x

xz
i i i i( ) ( ) ( ) ( )( ) ( ) ( )� �

�
�

�
	

�

d
d

,

M x E I
x

xy
i i i

i
( ) ( ) ( )

( )
( ) ( )�

d
d
�

,

where b and h i( ) are the width and height of the i-th layer,
recall Fig. 1, and k � 56, A bh

i i( ) ( )� and I b hi i( ) ( )( )� 112
3 stand

for the shear correction factor, the cross-section area and the
moment of inertia of the i-th layer, respectively.

To proceed, consider the weak form of the governing
equations, written for the i-th layer (the subscripts �x and �z
related to internal forces and kinematics-related quantities
are omitted from now on, for the sake of brevity)

d
d

d
d

d(i) (i) (i) (i)

(i)

x
u x E A

x
u x x

u x f

L

x

�	

	

( )) ( ( ))

( ) (

0
� �

� 
i
L

L

x x u x N x) ( ) ( ) ,( )d (i)�� 	 0
0

d
d

d

( )d

(i) (i) (i)

x
w x k G A x x

w x f x x w

L

z
i

�	 �

	 	

( )) ( )

( ) ( )

0
� �

� � 
( ) ( ) ,x V x
L

L
(i)

0
0
�

d
d

d
d

d(i) (i) (i) (i) (i) (i
x

x E I
x

x x x M
L

�	� � 	�( )) ( ( )) ( )
0
� �� 
)( ) ,x

L

0

	� � �(i) (i) (i) (i) (i)
d
d

d( ) ( ) ( ) ( ( ))x k G A x x
x

w x x� ��
��

�
��

0

0
L

� � ,

to be satisfied for arbitrary admissible test fields 	u(i), 	�(i) and
	w. In particular, the first three equations correspond to equi-
librium conditions written for normal and shear forces and
bending moments, respectively. The last identity enforces the
geometrical relation (2) in the integral form, thereby allowing
to treat the shear strain as an independent field in the
discretization procedure to be discussed next. Further note
that the continuity conditions (1) will be introduced directly
into the discretized formulation, as explained in the following
Section.

4 Finite element discretization
To keep the discretization procedure transparent, it is as-

sumed that each layer of the laminated beam is divided into

©  Czech Technical University Publishing House http://ctn.cvut.cz/ap/ 23

Acta Polytechnica Vol. 48  No. 6/2008

y
x

z, w

1

2

3

x, u
(1) �

(1)

h
(1)

z , w
(1)

h
(2)

h
(3)

�
(2)

�
(3)

x, u
(2)

x, u
(3)z , w

(2)

z , w
(3)

glass
PVB
glass

Fig. 1: Kinematics of a laminated beam



an identical number of elements, leading to the discretization
scheme illustrated in Fig. 2.

Following the standard conforming Finite Element ma-
chinery, e.g. [12, 11], we express the searched and test dis-
placement fields at the element level in the form

u x x u x xe
i

e u
i

e u
i

e
i

e u
i

e
( )

,
( )

,
( ) ( )

,
( )( ) ( ) , ( ) ( )� �N r N r	 	 ,

( )

, , , ,

( )

,

( ) ( ) , ( ) ( ) ,

(

u
i

e e w e w e e w e w

e
i

w x x w x x� �N r N r	 	

� x x x xe
i

e
i

e
i

e
i

e
i) ( ) , ( ) ( ) ,,

( )
,

( ) ( )
,

( )
,

( )� �N r N r� � � �	� 	

� 	� 	� � �e
i

e
i

e
i

e
i

e
i

ex x x x
( )

,
( )

,
( ) ( )

,
( )( ) ( ) , ( ) ( )� �N r N r ,

( ) ,�
i

where e is used to denote the element number, �e and 	�e
denote a relevant searched and test field restricted to the e-th
element, Ne

i
,•

( ) is the associated matrix of the basis functions
and re

i
,•

( ) is the matrix of nodal unknowns. In the actual imple-
mentation, the fields u(i), we and �e

i( ), as well as the corre-
sponding test quantities, are assumed to be piecewise linear.
To obtain a locking-free element, the shear strain �e

i( ) is taken
as constant and is eliminated using static condensation, see
[12, 11] for additional details.

To simplify the further treatment, we consider the follow-
ing partitioning of the stiffness matrix K and the right hand
side matrix R related to the e-th element and the i-th layer af-
ter static condensation:

K K
K K

r
r

Re
i

ew
i

we
i

w
i

e
i

e w

e
( ) ( )

( ) ( )

( )

,

�

�
�
�

�

�
�
�

�

�
�
�

�

�
�
�

�
( )

,
( )

i

e w
iR

�

�
�
�

�

�
�
�
,

where K Kew
i

we
i( ) ( )( )� T and

� 
r re i e ai e bi e ai e bi e w e au u w w( ) ,( ) ,( ) ,( ) ,( ) , ,, , , , ,� �� �
T

� 
e b, .
T

Considering all three layers in Fig. 2 together gives the re-
sulting system of governing equations in the form

K K
K K E

K K
K K

e ew

e ew e

e ew

we w

( ) ( )

( ) ( )

( ) ( )

( )

1 1

2 2

3 3

1

0 0
0 0
0 0

T

e we w w w

e

e

( ) ( ) ( ) ( ) ( )2 3 1 2 3 0
0 0

K K K K
E

r

� �

�

�

�
�
�
�
�
�

�

�

�
�
�
�
�
�

( )

( )

( )

,

( )

( )

1

2

3

4 1

1

2

r
r
r

R
R

e

e

e w

( )

e

e


 �

�

�

�
�
�
�
�
�

�

�

�
�
�
�
�
�

� R
R R R

e

e w e w e w

( )

,
( )

,
( )

,
( )

,3
1 2 3

0
� �

�

�

�
�
�
�
�
�

�

�

�
�
�
�
�
�

where the matrix 
 stores the nodal values of the Lagrange
multipliers, associated with the compatibility constraint (1),
and the matrix

Ee

h h

h h
�

�

�

1 0 0 1 0 0 0 0 0 0

0 1 0 0 1 0 0 0 0 0

0 0 0 0 1

1 2

1 2
2 2

2 2

( ) ( )

( ) ( )

0 0 1 0 0

0 0 0 0 0 1 0 0 1 0

2 3

2 3
2 2

2 2

h h

h h

( ) ( )

( ) ( )

�

�

�

�

�
�
�
�
�
�

�

�

�
�
�
�
�
�

implements the tying conditions.

5 Verification and validation of the
numerical model
To verify and validate the performance of the present ap-

proach, the previously described FEM model was imple-
mented using the MATLAB® system and compared with pre-
dictions of the analytical model and experimental data for a
three-point bending test on the simply supported laminated

glass beam presented in [5], see also Fig. 3. The width of the
beam is b � 01. m and the material data of individual compo-
nents of the structure are available in Table 1. The modulus of
elasticity of the glass is slightly lower than the conventional

24 ©  Czech Technical University Publishing House http://ctn.cvut.cz/ap/

Acta Polytechnica Vol. 48  No. 6/2008

y
x

z, w

i

�
( )i

e, a

we, a
�

( )i

e
a b

u
( )i

e, b

we, b
L

e
�

( )i

e, b

u
( )i

e, a

Fig. 2: Finite element discretization of the i-th layer

F glass 5 mm

PVB 0.38 mm
glass 5 mm

10 cm40 cm

strain gage

40 cm10 cm

Fig. 3: Three point bending setup for simply supported beam

Glass PVB layer

Young’s modulus, E 64.5 GPa 1.287 MPa

Poisson’s ratio, � 0.23 0.4

Table 1: Material data



values of 70–73 GPa reported in the literature, and is specific
for the material employed in the experiment. Moreover, as
the PVB layer shows viscoelastic and temperature-dependent
behavior, the modulus of elasticity corresponds to an effective
secant value related to load duration of s and temperature of
22 °C.

Table 2 summarizes the values of the mid-span deflection
for a representative load level determined by FE-based dis-
cretization using 60 elements (30 when symmetry of the
problem is exploited) and the corresponding experimental
values. Note that the discretization is sufficient to achieve
three-digit accuracy in the mid-span deflection. In addition to
the results obtained by the analytical method proposed by
Asik and Tezcan in [5], the results of the analysis using
ADINA® system and the elementary lower and upper bounds
are included. In particular, the ADINA® model is based on
the classical laminate theory, cf. [11], whereas the two simpli-

fied approaches assume zero or infinite compliance of the
interlayer, recall also the discussion in Section 2. In the follow-
ing discussion, e.g. �exp

num denotes the relative error between
the numerical prediction and reference experimental value,
while e.g. �an is used for the error of the analytical solution
when compared to the candidate approaches. Clearly, the
results of the last three methods differ substantially from the
experimental data and also from the analytical results. The
proposed numerical model, on the other hand, shows a
response almost identical to the analytical method, which
deviates from the experimental measurement by less than
6 %. Such accuracy can be considered as sufficient from the
practical point of view.

To further confirm the predictive capacities of the pro-
posed numerical scheme, a response corresponding to a
proportionally increasing load was investigated. The results
appear in Tables 3 and 4. Again, the method seems to be suf-

©  Czech Technical University Publishing House http://ctn.cvut.cz/ap/ 25

Acta Polytechnica Vol. 48  No. 6/2008

Model Central deflection [mm] �exp [%] �an [%]

Laminated glass beam: thickness [mm] 5/0.38/5 (glass/PVB/glass)

Experiment 1.27 – �5.2

Analytical model 1.34 5.5 –

Numerical model 1.34 5.5 0.0

ADINA® (Multi-layered shell) 0.89 �30.2 �33.8

Monolithic glass beam: thickness [mm] 10 (glass+glass)

Analytical model 0.99 �21.8 �25.9

Two independent glass beams: thickness [mm] 5/5 (without any interlayer)

Analytical model 3.97 212.6 196.2

Table 2: Comparison of results for a simply supported beam (load 50 N)

Load [N]
Central deflection [mm]

wexp wan �exp
an [%] wnum �exp

num [%] �an
num [%]

50 1.27 1.34 5.51 1.34 5.51 0.00
100 2.55 2.69 5.49 2.68 5.10 �0.37
150 4.12 4.03 �2.18 4.02 �2.43 �0.25
200 5.57 5.38 �3.41 5.36 �3.77 �0.37

Table 3: Comparison of deflections for a simply supported beam

Load [N]
Maximum strain [×10-6] Maximum stress [MPa]

�an �num �an
num [%] �an �num �an

num [%]

50 112 114 1.79 7.23 7.34 1.52
100 224 228 1.79 14.45 14.68 1.59
150 336 341 1.49 21.68 22.02 1.57
200 448 455 1.56 28.9 29.36 1.59

Table 4. Comparison of stresses and strains for a simply supported beam



ficiently accurate in the investigated range of loads when
considering the values of deflections as well as the values of lo-
cal stresses and strains.

6 Conclusions
As shown by the presented results, the proposed numeri-

cal method is well-suited for modeling laminated glass beams,
mainly because of its low computational cost and its accurate
representation of the behavior of the structural member.
Future improvements of the model will consider large deflec-
tions and the time-dependent response of the interlayer and
will be reported separately.

Acknowledgment
The support provided by GAČR grant No. 106/07/1244 is

gratefully acknowledged.

References
[1] Vallabhan, C. V. G., Minor, J. E., Nagalla, S. R.: Stress in

Layered Glass Units and Monolithic Glass Plates, Journal
of Structural Engineering, Vol. 113 (1987), p. 36–43.

[2] Norville, H. S., King, K. W., Swofford, J. L.: Behavior
and Strength of Laminated Glass, Journal of Engineering
Mechanics, Vol. 124 (1998), p. 46–53.

[3] Vallabhan, C. V. G., Das, Y. C., Magdi, M., Asik, M.,
Bailey, J. R.: Analysis of Laminated Glass Units, Journal
of Structural Engineering, Vol. 113 (1993), p. 1572–1585.

[4] Asik, M. Z.: Laminated Glass Plates: Revealing of Non-
linear Behavior, Computers and Structures, Vol. 81 (2003),
p. 2659–2671.

[5] Asik, M. Z., Tezcan, S.: A Mathematical Model for the
Behavior of Laminated Glass Beams, Computers and
Structures Vol. 83 (2005), p. 1742–1753.

[6] Duser, A. V., Jagota, A., Bennison, S. J.: Analysis of
Glass/Polyvinyl Butyral Laminates Subjected to Uniform
Pressure, Journal of Engineering Mechanics, Vol. 125
(1999), p. 435–442.

[7] Ivanov, I. V.: Analysis, Modelling, and Optimization of
Laminated Glasses as Plane Beam, International Journal
of Solids and Structures, Vol. 43 (2006), p. 6887–6907.

[8] Mau, S. T.: A Refined Laminated Plate Theory, Journal of
Applied Mechanics – Transactions of the ASME, Vol. 40
(1973), No. 2, p. 606–607.

[9] Šejnoha, M.: Micromechanical Modeling of Unidirection-
al Fibrous Composite Plies and Laminates, Ph.D. thesis,
Rensselaer Polytechnic Institute, Troy, NY (1996).

[10] Matouš, K., Šejnoha, M., Šejnoha, J.: Energy Based Op-
timization of Layered Beam Structures, Acta Polytechnica,
Vol. 38 (1998), No. 2, p. 5–15.

[11] Bathe, K. J.: Finite Element Procedures, 2nd Edition, Pren-
tice Hall, 1996.

[12] Bittnar, Z., Šejnoha, J.: Numerical Methods in Structural
Mechanics, ASCE Press and Thomas Telford Ltd, New
York and London, 1996.

Ing. Alena Zemanová
e-mail: zemanova.alena@gmail.com

Doc. Ing. Jan Zeman, Ph.D.
phone: +420 224 354 482
Fax: +420 224 310 775
e-mail: zemanj@cml.fsv.cvut.cz
[url]: http://mech.fsv.cvut.cz/ zemanj mech

Department of Mechanics

Prof. Ing. Michal Šejnoha, Ph.D.
e-mail: sejnom@fsv.cvut.cz

Department of Mechanics
Centre for Integrated Design of Advances Structures

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

26 ©  Czech Technical University Publishing House http://ctn.cvut.cz/ap/

Acta Polytechnica Vol. 48  No. 6/2008