SJ Mepla Manual Theory Eng

SJ MEPLA
Theory Manual
Version 5.0
October 2018
SJ Software GmbH Aachen

Dr.- Ing. Dirk Bohmann
SJ MEPLA Theory Manual, Version 5.0 Page 2
Table of Contents
1 Notations 3
2 Introduction 4
3 Elements 5
3.1 Multi-layered elements 5
3.1.1 Approach 5
3.1.2 9 node multi-layered element 7
3.2 Building of the finite element mesh 8
3.3 Point Fixing 9
3.3.1 Description 9
3.3.2 Approaches 16
3.3.3 Connection of point fixings 18
3.3.4 Forces applied at point fixings 18
3.4 Insulating glass units and the gas element 20
3.5 Spring element 22
3.6 Spacer element for insulated glass units 23
3.7 Elastic base 25
3.8 Elastic edge support element 25
3.9 Elastic line support element 26
3.10 Elastic edge beam element 27
3.11 Elastic bonding element 28
4 Basic approaches for calculation 29
4.1 Geometrically non-linear calculations 29
4.2 Contact approaches 30
4.3 Dynamic calculation approaches 32
4.3.1 Mass matrix 32
4.3.2 Newmark procedure 32
4.4 Pendulum body approach 33
4.5 Kelvin-Voigt Model 35
4.6 Principal stresses 37
4.7 Equivalent Von Mises stress 37
4.8 Loading approaches 38
4.8.1 Face loads 38
4.8.2 Line loads 39
4.8.3 Border Lineload 39
4.8.4 Concentrated loads 39
4.8.5 Dead weight 40
4.8.6 Temperature constraints 41
4.8.7 Pressure hit 41
4.8.8 Load case calculation 42
4.9 Residual load-bearing capacity 43
4.9.1 Calculation approach 43
4.9.2 Transformation of stress- and force vectors 46
4.9.3 Basic assumptions 47
4.9.4 Comparisons 48
5 Comparisons 48
5.1 Convergence study 48
5.2 Pendulum impact on heat strengthened glass 49
5.3 Pendulum impact on laminated safety glass 51
5.4 Hailstone Force Path 51
6 Notes 53
7 Literature 54

1 Notations
u midplane displacement vector
w, u, v, φ, θ degrees of freedom
ε strains
σ stresses
ε strain vector
g shear deformation
gL thermal volume expansion coefficient
D elasticity matrix
B strain matrix
t strain matrix
k shear correction factor
N, M, Q stress resultant vector
N k , Nˆ k shape function for node k
η,ζ local coordinate system
K stiffness matrix
P nodal load vector
U unbalanced residual forces
R reaction force vector
τ shear stress
E modulus of elasticity
G shear modulus
M mass matrix
m mass density matrix
H matrix of shape functions
δ,β stability parameters for Newmark-method
Δt time step length
σ internal forces and moments vector
AR tyre contact area of pendulum (foot print)
RP reaction force for pendulum
mP mass of pendulum
CR stiffness of pendulum body
σh principal stresses
ΔT temperature difference
ΔH difference of height
V, ΔV volume, volume change
pi 0 initial gas pressure
k material density

2 Introduction
The analysis of multi-layered plates is based on the following assumptions:
1. The deflections in plane except the transverse displacements w are small and the
material behaviour of each layer is isotropic.
2. According to the Mindlin plate-theory the displacement vector u = [u,v,w]T is expressed

as displacements of the mid-plane of each covering decking layer (i.e. the panes) uo, vo,
wo and independent normal rotations n and θ as:
u = u o(x, y) + zj (x, y)
v = vo(x, y) + zq (x, y) (1.1)
o
w = w (x, y)
n und θ are chosen as constant rotations related to the un-deformed system. They are not
necessarily normal to the midplane after deflection, shear deformations are considered.
3. The lateral displacements wo of each layer is set equal to w. Transverse strains εzz are
neglected, so that a plate like stress behaviour will occur, which can be expressed with
internal forces and moments N, M und Q for each layer.
4. The mid-plane displacements for the internal layers are build like (1.1) with
u = u z(x, y) + z j z (x, y)
v = v z(x, y) + z q z (x, y) (1.2)
o
w = w (x, y)
whereby the still to calculate functions for uz, vz, φz and θz result from the continuity of
deformations of the two enclosing layers.
The layer numbering counts from the bottom and increases to the top layers.
u 2 ( h2 /2) = u 3 (- h3 /2) ; u 2 (- h 2 /2) = u 1 ( h1 /2)

(1.3)
v 2 ( h 2 /2) = v 3 (- h3 /2) ; v 2 (- h2 /2) = v1 ( h1 /2)
A sharp bend in the displacements, which will lead to a skip in the layer wise constant shear
deformation, is thereby enabled.

3 Elements
3.1 Multi-layered elements
3.1.1 Approach
The displacements for each layer of a multi-layered element are defined as follows:
z, w y v0
2
n
u0
x
Degree of freedoms for a multi-layered plate
Each node of a n-layered element has 4n+1 degrees of freedom.

u = [w, u1o , v1o , j1 , q 1 , uo3 , vo3 , j 3 , q 3 ,..,q n ]T (2.1)
t5 n
uo5 5
t4
t3 n
uo3 3
t2
t1 n
uo1 1
z
x w
n freedom for x-direction
Degrees of
The defections w for all layers i=1, 3, 5, ... are set identical, which is possible as shear
deformations are considered.
n
3
~n
(3
z
n
2
x Mw
Mx
n (2
n
1 (1
Shear deformation for the decking layers and the intermediate layer
The displacements uiz, viz and rotations niz , θiz result from the existing degrees of freedom for
the respectively two enclosing layers.

1 æ 0 t i +1 t ö
z
ui = ç ui +1 - j i +1 + ui0-1 + i -1 j i -1 ÷ (2.2)
2è 2 2 ø
1æ t t ö
j i = ç ui0+1 - i +1 j i +1 - ui0-1 - i -1 j i -1 ÷
z
(2.3)
ti è 2 2 ø
By considering large deflections the strains will result by derivation the equation (1.1) to
2
i ¶ u io ¶ j i 1 æ ¶w ö
e =
xx +z + ç ÷ (2.4)
¶x ¶x 2 è ¶x ø
2
¶ vio ¶ 1 æ ¶w ö
e =
i
yy + z q i + çç ÷÷ (2.5)
¶y ¶y 2 è ¶y ø
and the in plane shear deformation from
¶ u i ¶ vi æ ¶w ¶w ö
g ixy = + +ç × ÷
¶y ¶x çè ¶x ¶y ÷ø
æ ¶ u io ¶ voi ö æ ¶ j i ¶ q i ö æ ¶w ¶w ö (2.6)
= çç + ÷ + z çç + ÷+ç × ÷
è ¶y ¶x ÷ø è ¶y ¶x ÷ø çè ¶x ¶y ÷ø
and the transverse shear strains
i ¶w
g xz = +j i (2.7)
¶x
¶w
g iyz = +q i (2.8)
¶y
with i = 1,3,5... for decking plates and zi as co-ordinate between -ti/2 < zi < ti/2.
The internal layer strains i = 2,4,6,... are coupled with the decking plates to build the total
system.
Shear strains of the internal layer (i.e. for x-direction)
¶w ¶w 1 æ 0 ti +1 t ö
g ixz = + j iz = + ç u i +1 - j i +1 - ui0-1 - i -1 j i -1 ÷ (2.9)
¶x ¶x t i è 2 2 ø
The in-plane and bending stiffness of the internal layers (laminate material) are expressed by
the following equations (i.e. x-direction):
2
¶ u iz 1 æ ¶ u oi+1 t i+1 ¶ j i+1 ¶ u io-1 t i -1 ¶ j i-1 ö 1 æ ¶w ö
= ç - + + ÷+ ç ÷
¶x 2 çè ¶x 2 ¶x ÷ø 2 è ¶x ø
(2.10)
2 ¶x ¶x

¶ j iz 1 æ ¶ uio+1 t i+1 ¶ j i+1 ¶ uio-1 t i -1 ¶ j i -1 ö

= ç - - - ÷
¶x t i çè ¶x 2 ¶x ÷ø
(2.11)
2 ¶x ¶x
The internal layers with indexes i = 2,4,6,.. are coupled by the degrees of freedom from the
deck plates with the underlying and overlying layer to build the over all stiffness matrix.
The plate element uses a modified shear interpolation function. These are necessary for
compatibility due to the shear deformations. Followed by the use of deflection derivatives,
non according polynomials would be used. Shear deformation results of the sum MN/Mx and
N, which results in element locking.
Therefore shear deformations will be build with an approach from [7]. These modified shape
functions can be easily calculated for each element and can be taken for each element
integration. The continuity of shear deformation is ensured and locking will not occur.
According to the degrees of freedom, the load vector P allows the following imposed forces:
kn T
P = { Fkz , Fk1 k1 k1 k1 k3 k3 k3 k3
x , F y , M x , M y , F x , F y , M x , M y ,..., M y } (2.12)
The resulting internal forces and moments are referred to the layers. So the multi-layered
element is a combination of a plate element with in-plane degrees and a 3-dimensional
element, where transverse strains εzz are not regarded and the deflections (u, v) are set
constant or linear (zAn, zAθ) for each layer.
3.1.2 9 node multi-layered element

The here used isoparametric element is build with 9 nodes, using the Lagrange - shape
interpolation functions:
4
7
3
>
8 9
06
1
5
z y 2
x definition of nodes
9-node element,
Integration order:
3x3 = 9 Gaussian points, (full integration)

Use of modified shear interpolation functions
Degrees of freedom:
The number of degrees of freedom for each node of one element with n-layers (one
package) is 4n+1:
T
u = { w, u , v , j , q , u , v , j , q ,..,q }
o
1
o
1 1 1
o
3 3
o
3 3 n

Load vector:
Load vector for each node of each package:

kn T
P = { Fkz , Fk1 k1 k1 k1 k3 k3 k3 k3
x , Fy , M x , M y , Fx , Fy , M x , M y ,..., M y }
Element building, order of layer numbering:
The number of layers must be odd, as two decking plates must enclose one intermediate
layer. The count ordering of layers begins with 1 and increases to the top layer 3, 5, 7,... If
more than one package (multi-layered laminated glass panes) is given, the element layers
will be continued for the next package, then including the gas volume in between (insulation
glass). Each layer is described with the young's modulus, the thickness and the poisson’s
ratio.
3
2 package 2
1
SZR intermediate space
5
z 4
3
2 package 1
1 t
Order of layers and packages
Stresses:
For each layer the stresses on the top and bottom faces σxx, σyy, σxy of each the glass panel
are calculated. These stresses are declared in SJ MEPLA with the following convention: Sxx,
Syy und Sxy (Components in x, y- direction and in plane shear stress).
Transverse shear stresses are used and calculated in this program, but are not written in the
resulting output file.
Principal stresses are calculated in the Gaussian points by the following formulae
2
s xx +s yy ( s xx - s yy )
h
s (+,-) = ± + s 2xy
2 4
and are identified in the program by Sp+ und Sp- (principal stresses +,-).
3.2 Building of the finite element mesh

The finite element mesh is formed by a free meshing algorithm. On the basis of any outline
contour, the finite element mesh is automatically produced in the inside. Beside from straight
outline edges, also circular curved edges are possible.

The formation of holes is done in that way that an additional edge inside the system is
described, in order to build the inside hole and then on the way back this additional edge will
be closed again. Such additional “stitched up” edges are not regarded as system edges.
The mesh refinement is controlled via the given element size.
Remark:
The element size used in this program will be determined in parallel from the smallest edge
length, which was set to describe the outline contour. If a system edge (e.g. a small circular
rounded corner) with a very short length (e.g. 30 mm) was set, then a very fine mesh with
elements of 30x30 mm edge length are produced, so that the circular border can be at least
formed by one element!
The following basics must be regarded to describe the system contour:

- The input of the contour lines must take place in counter clockwise direction.
- The edge lines may not over-cross. This applies also to the circular edges.
- If circular contour lines are described, the following border point must lie on this circle.
If this is not the case, then “auto-correction” can be used to generate the correct
location of this point. The actual position of this point is shifted in radial direction to
become a point lying exactly on the circle.
- The building of cut-outs or any holes on the inside of the structure is formed by
additional contour lines, which must be closed on the way back with a second contour
line having identical coordinates however in opposite direction. These double
however parallel in the opposite direction running edges, are the only exception for
attaching contour lines.
Remark:
Sometimes it can occur that the finite element net cannot be formed correctly. In these cases
the element size must be changed to accomplish a renewed meshing attempt.
3.3 Point Fixing

3.3.1 Description
The glass fixing element is a combined element between the multi-layered element and the
further use of one new node with 5 degrees of freedom. This node is introduced at the base
point (or rotation center) of the point fixing (point where the springs are applied). The core of
the fixing made from steel is set undeformable (blue region). So core deformations are not
considered due to the rigid body approach. But the elastic properties of the imbedded circular
disk pad (shim) (hs, red region) and the bush or the protective layer at the contour edge
(green region) to prevent direct glass-steel contact, are considered within this point fixing
element.

rk
th
ri
hk
z
y
x hs
zh
ra Cx
Cy
C : Drehsteifigkeit um y-Achse
C : Drehsteifigkeit um x-Achse Cz
Definition of point fixing element with rotational stiffness Cφ and Cθ
The supporting stiffness of the fixings base point is mapped by 5 springs. Cx, Cy, Cz for
stiffness of translation in x, y and z-direction and the rotational rigidities Cφ and Cθ for rotation
about the y- and x- axis!
By use of the distance Zh (down, against z-direction a negative value) an eccentric loading of
the supports may be considered, as the springs are applied to this base point and the
reaction forces will there be transferred. Stiff or fixed boundary conditions may be set by use
of high spring rigidities like 1.e6 N/mm and 1.e10 Nmm/rad for the rotation springs. Shall only
translation be forbidden and rotation be allowed, the rotational rigidity must be set to Cj = Cq
= 0 Nmm/rad. The base point of the support must therefore be set by the z-ordinate to the
position of the rotation center, while no rotational rigidity is set.
Using this springs the elastic behaviour of the underlying structure or special mechanisms
like a ball-shaped head can be considered.
Within this program 10 different kinds of fittings can be used:
Type Description
1 Countersunk fixing
2 Disk fixing
3 Circular clamp fixing
4 Angular clamp fixing
5 Circular downholder
6 Angular downholder
7 Bonded fixing without building of holes
8 Countersunk fixing within layer 1, laminated glass (LG)
9 Countersunk fixing over all layers in package 1, insulated glass unit
10 Disk fixing over all panes from package 1, insulated glass unit
1 3 5 7 9
2 4 6 8 10

Type 1: Countersunk fixing

rk
ri th
Eh hk
+z
ts
Es
Zh
ra
Cx
Cz
The counter sunk head always lies aligned with the highest glass panel. In the program it can
be selected only the conical range or also the cylindrical range of the bushing is lying in
contact to the glass and thus can introduce forces into the pane.
Type 2: Disk fixing
ri
ts
Eh
+z
ts
Es
Zh
ra th
Cx
Cz
The disk fixing exhibits two separating disk pads (red) on both sides of the surface. The disk
layer (green) transfers thereby in-plane the forces into the bore hole rim; the plate layer (red)
transfers the forces, which only act perpendicularly to the glass surface.
It can be chosen, which glass layer shall lie in direct contact to the borehole rim. For
boreholes filled with glue (e.g. plastics filling material, green region) can be set, that for all
rims within this borehole loads shall be transferred. This possibility enables to consider the
manufacturing method, where usually only one borehole is drilled with the correct diameter
and additional above lying holes will have greater diameters. Such it can be ensured that
after lamination the bush can be inserted. The in-plane load then only can be transferred via
the defined panel.
Type 3: Circular clamp (balustrade clip)
th
ts
Eh
+z r
ts
Es
Zh Cx
r
Cx
Cy
Cz

This clamp support exhibits a round, circular design. The load transfer takes place in the
same way as the disk fixings, via the flexible intermediate layers within the clamping, which
separate the glass edges and the surface from the metallic parts.
Eccentricities from the distance of the base point Zh are considered. When a position for the
clamping is set, which is not exactly lying on the border lines, the resulting eccentricity in x,y-
direction (distance to the border) is regarded too.
Type 4: Angular clamp (balustrade clip)
th
ts
Eh
+z b
ts
Es
Zh Cx a
b
Cx
Cy
Cz
Instead of a circular clamping, an angular design can be chosen. This form is defined by “a”
and “b” value, which describes the half width and the depth. The radius for the corner
rounding is done automatically and can not be changed.
The shape of the finite element mesh is formed from diophantical curves, thus ellipses of 4th
order. The mesh is build with two dense rows of elements inside and two outside the
clamping. The inside of this clamp is meshed automatically.
Type 5: Circular downholder
Es
ts
+z r
r Zh
Cx
Cx
Cz Cy
The circular downholder corresponds to the circular clamp support, with the restriction that
only forces from above (thus in negative z-direction) can be introduced and no forces into the
edges in plane direction can be applied. This downholder can thus be used as suction
protection. That the downholder acts only in one direction for compression, the contact for
the interfaces must be marked. When contact isn’t marked, the downholder acts for
compression and suction like bonded!

Type 6: Angular downholder
Es
ts
+z b
b Zh a
Cx
Cx
Cz Cy
The angular downholder is equivalent in function to the circular downholder.

For insulated glass it can be given for type 5 and 6 which glass package shall be clamped.
Type 7: Bonded disk fixing
+z
ts
E s ;G s
Zh
ra
Cx C,C
Cz
The bonded disk fixing can only be set from the bottom side onto layer 1 of package 1.
The characteristic of this fixing is a glued disk without a generated hole. In contrast to the
other point fittings, here the elastic module and the shear modulus must be set for the
interface (adhesive film).
Type 8: Countersunk fixing within layer 1, (laminated safety glass) LSG

rk
ri th
Eh hk
+z
ts
Es
Z
h
ra
Cx
Cz
This type of fixing shows the same properties than type 1; but is only fixed within the first
layer of package 1. All other additional layers or as well packages are generated without
holes. According to the element approach used in this program here as well all layers of one
package are coupled using the one degree of freedom w. Components of forces in z-
direction are transferred directly to all layers, what will lead to a small inaccuracy at the
borehole rim. This is first of all depending upon the shell element approach, where the
stresses and strains distribution in z-direction is neglected.
Type 9: Countersunk fixing for package 1, insulated glass unit

rk
ri th
Eh hk
+z
ts
Es
Z
h
ra
Cx
Cz
In contrast to type 8 this point fixing will always clamp all layers from package 1. All further
packages are build without a hole. It’s to take care, that all larger packages (>1) must be
fixed in a static determined way. This can be done by giving a small shear stiffness G to the
spacer. Otherwise the higher packages are free to move in plane direction which may cause
an error due to a statically undetermined situation.
Type 10: Disk fixing for package 1, insulated glass unit
ts
ri th
+z Eh
ts
Es
Z
h
ra
Cx
Cz
Type 10 is very similar to type 9 but defines a double sided disk fitting clamping the first glass
package.
Building the mesh:
The glass point fixings need for integration a special finite element mesh, where the rubber
disks and the edge- and borehole separating layers can be fixed. This mesh is generated
automatically depending on the mesh density.
The fixings of type 1, 2, 7-10 can’t be set too near the border (distance greater than 20 mm).
Otherwise a mesh can not be build.
The clamps 3 to 6 can only be set at the borders. When they are located farther than 70 mm
from edges, they can’t be generated and an error message will appear. If they are positioned
within this region, a perpendicular will be dropped to the border line, in order to specify the
position at the edge:

Clamping with offset to the edges
The arising offset (see picture) is considered in the calculation.
When a clamping or downholder shall be set onto corners, this will only be done if
a.) a circular shape is set (type 3 or 5), and
b.) the given distance to the corner is less than 30 mm.
In any case the given position is modified to lie exactly at the corner position. An offset can
not be regarded! That means, the point fixings base point (the position of the axis, where the
springs are connected) is modified in the program and will always be equal to the given
corner co-ordinates! A corner position is formal equal to a contour point, also when no real
corner has been build (if and additional point along a straight line has been generated)!
Remark:
Against the simplified illustration in the “Graphic Surface” the clamps of type 3 do not only
attach the top and bottom faces of the glass:
The two glass edges in the corner are as well supported over the glass thickness with
separators (green indicated).
The number of elements, which will be build around and within the point fixings, depends
upon the chosen element size. If a small element size is used, the density and the amount of
elements will be very high, accompanied with higher efforts for solving the equation system.

3.3.2 Approaches
All elastic properties of the elastic separators are considered when assembling the stiffness
matrix. According to the young's modulus, the thickness of the disk pads (shim) and the
properties of the bush and edge separators, the stiffness matrix is assembled by integrating
over the volume.
Hereby special properties are set for the separating layers:
Only forces perpendicular to the glass surfaces can be transmitted to the supports (except
type 7). There is no shear transfer between the glass and the separating layers. These layers
only can slide over the glass faces where the materials are touching. This is the most
accurate approach without taking contact-algorithms additionally into account. The bush and
shim may only be stressed through perpendicular forces.
E bush
loaded area
steel E =
glass panels
z
E shim
steel E =
eccentricity (-)
point fixing node
(5 degrees of freedom)
Force transmission for countersunk fixing
steel E =
E shim
E bush
steel E =
z glass panels
E shim
steel E =
eccentricity (-)
point fixing node
(5 degrees of freedom)
Force transmission for the disk fixing
This approach will give most accurate stresses at the drilling rim, without to pursue a great
effort with interface-elements. In this case (without contact approaches for the touching
faces) no iterations are necessary and the equation system will be solved in one step; but
then perpendicular tensile stresses may arise also.
In case of calculations with contact approaches (section 4.2) the contact problem will be
simplified, as only one stress direction must be regarded in this approach. This describes the
most realistically calculation, as tensile stresses may not occur at the surfaces, the bush and
edges. Only compression forces will arise due to the possibility of separation.

The bonded disk fixing forms an exception:

The adhesive film can transfer normal stresses as well as shear stresses.
Additionally this (as well type 3 and 5) can be set to free rotation.
=0
For that purpose a degree of freedom for rotation around the z-axis has been added. This
property can be set to “free” or “stiff” behaviour for the z-rotation. But, if such a fixing is
connected using a bar and not with 5 springs, it is always set to free rotation.
The default (unmarked) is the rotational stiff setting!
Special note:
With this multi-layered element approach, the definition of the countersunk fixing is
underlying a simplification, which will lead to good results, when the chamfer at the drilling
rim is small compared with the thickness of the pane. The loads acting onto the chamfer are
considered at the correct place, but bending and in-plane stresses in the pane are not
influenced by the chamfer. This is founded by the fact, that this element approach may not
consider the real chamfer like volume-elements can do. The chamfer really doesn’t exist but
will be treated as a stress free region.
The double sided disk fixing is not underlying this restriction.
Point fixings in insulation glass:
When point fixings are used in insulating glass units, a spacer will automatically be inserted.
Its properties are taken form the spacer definition. The spacer is aligned in between the gap
using the predefined gap height (see section 3.6).
Spacer at the borehole rim
The elastic and mechanical property of the spacer is inserted directly at the rim of the first
inner element row of the borehole to couple the glass packages.
Load transferring at the borehole rim:
The definition which layer shall lie in direct contact to the bush of a disk fixing (type 2, 10),
can only be selected for the lowest glass package 1. When here a layer is set to transfer
loads, other packages then don’t have touching faces between the glass and the bush.
When the conic face of a countersunk fixing (type 1, 8 or 9) is defined as the only contact
region, then no other layers (glass plies) are regarded with touching faces to the bush! The

conic face (except the underlying disk pad) then is the only force transmission area between
the glass and the fixing.
3.3.3 Connection of point fixings

The point fixing element is formed from one additional node with 5 degrees of freedom. The
deflection is described with 3 displacements and 2 rotational degrees of freedom at the
fixings base point. A rotation about the transversal z-axis of the point fixings is only
considered for type 3, 4 or 7. The steel body, build from the bolt, the disks and the taper
shank, is set as a rigid body. The eccentricity between the base and the lower bottom face of
the panes will be considered and is defined by the distance Zh (negative value, if the base
point lies under the pane).
The displacements and rotations rigidities are described with springs (or by definition of a bar
or cable), fixed at this additional node. These springs represent the stiffness of the sub-
construction or the design of the supports (e.g. balled shaped head). The 5 springs are fixed
to the ground.
High rigidities for the springs represent a fixed boundary condition, which is nearly
undeformable (e.g. Cx = 1.e6 N/mm). Note that such a high stiffness can’t really arise (1 MN
or 100 tonnes will lead to a deformation of only 1 mm!). A stiff rotation spring may be defined
using 1.e10 Nmm/rad, but not more! Low rigidities of C = 0.0 will form a free moving or
rotation possibility. Setting values exactly to zero, may leed to problems when contact
conditions shall be used in parallel. In this case after releasing the contact, the point fixing
may become kinematically, what can’t be solved. In such cases it’s recommended to use low
stiffness values (e.g. 10 N/mm or Nmm/rad), that prevent moving or rotation of a de-attached
point fixing!
Using only point fixings as supports, a static determined system must also be observed.
In case of a tension rod only normal forces and no bending can be transmitted. The two
endings of the bar are fixed using hinges.
Zo
Z
Y X
Xo
Definition of the tie bar position
The rigidity of the bar will be described by the young’s modulus E and the cross section area
A. These values will be used to calculate the bar normal stiffness. An eccentricity within the
point fixings definition will also induce bending forces into the panel, when the basic point
(defined by Zh) doesn’t lie exactly within the midplane of the glass package. Due to the load
directions, the pane must be fixed as well in pane direction to prevent moving. In the above
example, face loads will lead to y- reaction forces within the pane, so that the edge along the
x-axis must get a boundary conditions acting in y-direction.
A possible buckling analysis – nor in the case of a non-linear calculation - of the bar is not
considered and must be checked by hand calculation.
3.3.4 Forces applied at point fixings

At all point fixings loads can be applied directly. Instead of defining spring stiffness or bar
cross sections rigidities now as well 3 forces Fx, Fy, Fz and moments the Mφ und Mθ can be
set. All combination of such different bearing and loading approaches are simultaneously
possible.
Degrees of freedom:
5 degrees at the base point, (u, v, w, j, q)
Input:
- point fixing dimension and position in the glass panel

- type of fixing (type 1 or 7)
- material properties for the disk pad, the bush and edge separators
0: Springs:
- 5 spring rigidities [N/mm] and [Nmm/rad]
1: Bar:
- Connecting a bar to the ground with definition of the young’s modulus E and the cross
section A.
2: Forces:
- Applying forces and moments
Output:
- 5 displacements and rotations u, v, w and φ, θ

- 5 spring forces (Fx, Fy, Fz and Mφ, Mθ)
or
- FStab normal force in the bar
or
- 5 displacements and rotations u, v, w and φ, θ

3.4 Insulating glass units and the gas element

The approach for the multi-layered elements allows the calculation of changing gas pressure
of insulated glass units by use of the gas-equations. With this element description with
varying layers, the behaviour of the gas in-between the panes can be correctly considered
for any design. Due to the same topology of the mesh for all layers, this allows the correct
integration of the gas stiffness by the Gaussian points.
pe
intermediate space
volume Vo
w
gas pi
w
insulation element pe
Insulating glass unit element
Every pane package has its own degrees of freedom in z-direction. The change in volume
due to deflection of each pane can therefore be correctly integrated.
In consequence to the high non-linear gas law, the equilibrium solution is solved through
iterations. The re-arrangement of the gas volume, caused by the deflection of the panes are
considered. With this approach, pressure non-conforming loads (point loads and bearings)
and geometrically non-linear behaviour of the panes can be calculated all together.
Here are to name:

- large displacement approaches
- insulation glass with laminated glass panes
- any geometry design
- any combination of pressure loads, point loads, dead weight,...
- pendulum impact on insulating glass units
- climate loads like temperature change, barometric changes,...
- loads due to temperature changes in the gap
- point fixed insulating glass units
From the gas law
pi × V = konst
and the volumetric change of the initial volume V0 by temperature changes
V = V 0 × e g L DT » V 0 (1 + g L × DT)
with
(L : thermal volume expansion coefficient
ΔT : temperature change of the filling gas
the resulting gas pressure of a deformed insulating glass unit

0 V 0 (1 + g L × DT)
pi = pi
V 0 + DV
with
pi 0 : initial pressure during production

ΔV : volume change
can be calculated.
The change of barometric pressure (e.g. on sea level: normal 10133 mbar, high: 1030, low
970 mbar) due to the differences of assembling- and production height might be considered
by the barometric formula [54]
5.255
æ 0.0065 × DH ö
p ( h 1) = p ( h 0) çç1 - ÷÷
è T ( h 0) ø
with
p (h0): barometric pressure on reference height h0 (production height)

ΔH: difference of height [m], (installation height - production height)
T(h0): temperature in K related to the reference height h0 (273,15 + t(h0) C°)
up to a height of 11km, when the real barometric outside pressure is unknown. To simplify
this formula a middle temperature for the reference height (sea level) of 15° C is considered
(273.15 + 15 = 288.15 K). If a different altitude with a lower air pressure is used, a
temperature change from 0.65°C for 1000m altitude change is automatically applied
according to the standard atmosphere. P(h0) is thus directly linked to T(h0).
The iterations are performed by use of the Newton-Raphson iteration procedure. If necessary
the tangential stiffness matrix (also the gas stiffness) is build, to enhance the solution
process.
If geometrically linear calculations are performed, the stiffness matrix is only once assembled
and used for all further iterations. Although with this procedure more iteration must be carried
out, the solution takes less time, as the equation system only once has to be triangulated
(section 4.1). If convergence can not be achieved, a new assembly of the stiffness matrix
and equilibrium solution is performed.
The maximum number of glass packages is 4. (3 internal intermediate spaces with gas)
Input:
- external (outside) pressure pa [N/mm²], barometrically changed pressure on production

height p(h0)
- internal gas pressure pi [N/mm²] for each intermediate space (worst case values)
- temperature difference in intermediate gas space from the difference between installation
and production temperature ΔT [K°]
- height difference between assembly and installation ΔH [m]
(by this the new pressure on installation height p(h1) will be calculated)
- tolerance values for internal pane contact [mm], (section 4.2)
Output:
- internal pressure pi in each intermediate space
- stresses in the spacers sealing
- formed volume

Special note:
The edges of insulated glass panels must be coupled internally or symmetry conditions may
be used. Therefore 2 methods are possible:
- The use of the boundary condition with type 0, which will simply support all packages at
the given edges and so will keep in parallel a constant distance.
- The edges may be connected by spacers with elastic/rigid behaviour. Then distance
changes in dependency of the defined material properties are possible. In addition as
well beams may be selected at the same edges, so that no edge supports (type 0 - 7) are
necessary.
- Is possible too, to fix double glass only by glass fixings and to seal all other unsupported
borders only with spacers.
- The difference of height should only be regarded, when the insulation glass is not filled
with gas at site. In this case ΔH can be set to zero.
3.5 Spring element

Any desired springs may be set at any position within the pane area. These springs operate
in 5 the dimensions for displacement x, y, z and rotation j and q. Spring rigidities in x- or y-
direction only acts onto the lowest layer of each package. Rigidities in z-direction take effect
on all layers due the element convention (each layer of a laminated glass has the same
transversal degree of freedom w). Rotational rigidities are also applied only at the first layer.
All of these springs operate between the plate and the ground.
In SJ MEPLA they are used to establish a statically determined system. A pane simply
supported (bearing type 0) must be prevented for undefined moving in plane direction.
Therefore a proposal is done by the program for each pane package:
- The first corner point receives 2 springs in x- and y- direction.

- The second corner point receives 1 spring in y- direction.
Cx
1 2
z
y
Cy Cy
This proposal can of course be removed when other kind of bearings are necessary or will
take over this fixation.

Cz
Supports at the corners can be applied as well by single springs.
The precision of the stress calculation depends as always by finite element analysis on the
mesh (element) discretisation. Local supports always generate a singularity, so that getting a
real stress peak will not be possible. For a more detailed stress investigation a local mesh
refinement is needed, what will be done automatically by using glass point fittings.
Degrees of freedom:
No own degrees of freedom (coupled between plate and ground)
Input:
- position of springs in x-y coordinate system

- rigidities in x,y, and z- direction Cx, Cy und Cz [N/mm]
- rotational rigidities for Cj , Cq [Nmm/rad]
(Rigidities of 1.e6 N/mm (1 MN/mm) denote nearly no movement possibility; rigidities of zero
will remove a spring)
Output:
- displacement u, v, w [mm] and φ, θ [rad]
- spring forces Fx, Fy, Fz [N] and Mφ , Mθ [Nmm]
3.6 Spacer element for insulated glass units

For the calculation of point fixed glass or unsupported edges in insulation glass units
spacers in between the glass packages are necessary. With these spacers the panels are
coupled at these edges. When the pane borders are supported by the boundary condition
type 0 (simply supported), this approach will act onto all packages and therefore the distance
in-between is fixed too. For free edges in insulated glass units a spacer must be used to
combine these edges and to prevent a free opening.
For this purpose a spacer should be set, so that neither opening nor closing of the free edge
is possible. To simulate the overall behaviour of the spacer (build by the aluminium profile,
the sealing and the butyl) the modulus of elasticity E and a width b is set by default - but can
be changed too.
E, G
b

In addition to that, the shear stiffness may be controlled by the shear modulus G. Now the
panes are coupled also in plane direction and the panes are combined locally by shear
effects. Using the shear modulus G will such imply a local sandwich structure. It should be
checked if such behaviour is according to used Standards. Normally in Germany it’s not
allowed to consider this.
The height of the spacer depends on the air gap height and is set automatically.
Mechanically these properties will be added at the edges of the 2 involved 9- node elements.
The stiffness matrix of the spacer is produced by numerical integration over 3 Gaussian
points along this element side and is added to the total stiffness matrix so that both edges
are coupled. When 2 or 3 air gaps are present, the same is done there.
In extension the spacer may act non-linear:
When non-linear behaviour is chosen, the distance change is non-linear controlled. When
enlarging the distance the behaviour is linear due to the given elasticity. When the distance is
reduced the panes will get in stiffer contact with the aluminium profile. This behaviour then is
simulated by stiffening the material properties to prevent a negative distance. The maximum
rigidity is reached when the deformation due to compression is larger than the tolerance
value. This represents a rigid body for the spacer (intermediate frame), which will not
transversal deform.
A constant distance may also be set by a higher rigidity. This distance will then be kept if
larger opening forces arise. A non-linear spacer approach is not a standard calculation and
should only be used for scientific investigations.
When the young’s modulus is changed, the preset shear modulus G will be modified too. The
standard value of zero for no shear transition (which is not allowed to set into account in
Germany) is set new to G = E/3. This value must again be set to zero, when shear effects
shall not taken into account!
The selected spacer property settings will automatically be taken for the borehole rim, when
point fixings in insulated glass units are used. A very high rigidity may result in local problems
for the stress calculation. So real stiffness values should be used.
Degrees of freedom:
No own degrees of freedom (only coupled in between the packages of an insulation glass)
Input:
E [N/mm²]: modulus of elasticity (in z-direction)
G [N/mm²]: shear modulus G (in x, y- direction)
b [mm]: width of sealing (representing the combination of the
aluminium profile and the silicone sealing)
edge number: panel border, which will get a spacer
contact tolerance: non linear behaviour for tension and compression (compression: no
change in distance, tension: deformable due to the given properties
e.g. silicone)
Output:
- maximum stresses within the spacers sealing material (uniformly and related to the
used width b)
- symbols for the spacer in the graphics surface
Special note:

1. Bending rigidity is never regarded for the spacers (no bending stiffness of the
aluminium profile or the sealing material is used).
2. For point fixed insulated glass, these properties will automatically be used at all bore
hole rims.
3. As a first suggestion to consider all effects at a spacer a young’s modulus of 100
N/mm² using a width of 5mm is used and no shear stiffness is set (according to
German standards).
3.7 Elastic base

The elastic base element acts on the bottom side of the lowest layer (layer 1 and package 1).
This behaviour is integrated over all elements (9 Gaussian points) resulting from the stiffness
out of C = EA/h and is added into the global stiffness matrix. For these elastic base elements
the effects from the lateral contraction is neglected (< = 0) and so the calculation can take
place using springs, which acts only in transversal direction. With this approach no
rearrangements in the elastic base can be considered (e.g. when the material is compressed
at one position, this will have no effect to the strains at places aside).
h
E
Contact condition can be chosen as well:
When the lateral displacement w exceeds the given tolerance, the panel will detach from the
elastic base, so that no combination will exist any longer. When the material is under
compression, forces are transmitted.
Degrees of freedom:
No own degrees of freedom, (coupled between ground and bottom face of the panels)
Input:
E [N/mm²]: Modulus of elasticity (in z-direction)
h [mm]: Height of elastic base
contact tolerance [mm]: Detachment from the base material
Output:
No output, only spring symbols are shown in the graphics surface
3.8 Elastic edge support element

When the edge of the plate shall be elastically supported, this can be done by these “elastic
edge support”. In dependency of the defined underlying supporting structure (e.g. a rubber
bearing) the stiffness is calculated and integrated at 3 Gaussian points in z-direction due to
the young's modulus, width and given height. When the shear modulus G is set as well, this
will act as a bonded edge (structural glazing) where now also shear forces are transferred.
The width b is only used to calculate the stiffness of such a rubber strip. The pane is only

supported along this (small) line which generates no clamping stiffness laterally. Only if
several elastic and spaced lines are used, a clamping effect can considered.
E,G
h
Related to the defined stiffness E and G the corresponding reaction stresses are written into
the protocol.
In addition, contact approaches can be considered, so that lifting corners may occur. In this
case the shear modulus is set internally to zero!
Degrees of freedom:
No own degrees of freedom (coupled between layer 1 of package 1 and the ground)
Input:
E [N/mm²]: Modulus of elasticity (in z-direction)
G [N/mm²]: Shear modulus (action in x- and y-direction)
b [mm]: Width of the elastic profile
h [mm]: Height of the elastic profile
contact tolerance [mm]: Distance change, when separation shall take place
Outputs:
- Reaction forces Fz, Fx and Fy for each line
- stresses in x-, y- and z-direction
- symbols in the „graphics surface“
3.9 Elastic line support element

The elastic line supports analogously operate like the elastic edge supports, but can be set
anywhere within the pane and acts only in z-direction. The edge supports from 3.8 are
integrated along the border of the elements using 3 Gaussian points. The line supports are
generated in a different way: The supporting line is set by two co-ordinates, defining the
beginning- and ending point. Along this line equally spaced springs are applied, which acts
between the pane and the ground. The spring rigidity is calculated by C = E·A/h, whereby the
area A is build form the spacing length for each spring and the width b of the underlying
elastic profile.
b E
The stiffness matrix is assembled by adding each spring rigidity, so that the integration in this
case is done by summation of these segments at defined positions. This spring rigidity only
acts in z-direction. In plane movement is not suppressed and the plate may slide over this
type of bearing.

The starting and ending position of this supporting line should lie within the pane. If co-
ordinates exactly at the border are given, it could be possible that this position can’t be
evaluated as it’s not definitely clear if such a point is now lying within or outside of the pane
area. Such co-ordinates should always be given to lie some millimetres within the pane, so
that no error message is produced.
Hole and cut-offs are also areas outside of the pane!
Degrees of freedom:
No own degrees of freedom (coupled between layer 1 of package 1 and the ground)
Input:
Point 1: x, y- coordinate from the start point
Point 2: x, y- coordinate from the end point
E [N/mm²] Modulus of elasticity (in z-direction)
b [mm] Width of the elastic profile
h [mm] Height of the elastic profile
contact tolerance [mm]: Distance change, when separation shall take place
Output:
- reaction forces Fz for each line
- maximum stresses σzz
- symbols in the „graphics surface“
3.10 Elastic edge beam element

The elastic edge beam element is a reinforcing element, which will be overlaid at the borders
of the 9-node elements. The displacement approaches for this beam element corresponds to
a 3-node beam element with shear effects. As no new degrees of freedom are necessary,
this beam element is build only with a bending rigidity and has no torsion stiffness. In this
way the pane can freely rotate transversal to the beam axis. The behaviour is equivalent to a
simply on a beam supported pane. The length of beam is related to the length of the selected
edge, where the beam is located. The support conditions of the beam beginning (A) and
ending (B) depends on the degrees of freedom of the nodes located there. If a system border
got a simply supported condition, the beam will use this support at it’s ending too. It’s only
possible to remove additional degrees of freedom, so here to set clamped beam endings.
The boundary condition is set using the same type of bearings used for the edge supports
(type 0 to 7). The direction from (A) to (B) is defined counter-clockwise such as the system
contour has been given. Elastic beam supports may be combined with springs or clamps at
any position.
B
E, I, A
According to system borders which shall be reinforced by such a beam, the related element
edges are determined. Along these sides the beam element is build using the 3 nodes of this
side and integrating over 3 Gaussian points. These edges may not be curved, due to the
restriction, that torsion is not included in the beam approach. This additional stiffness acts
only on the lowest package 1 and layer 1. Similarly, normal and torsion forces can't be
transmitted, as only a bending and shear stiffness exists in the beam. Because the density is

given too, the dynamic response of the swinging beam and the loads due to the dead weight
are considered for calculation.
Degrees of freedom:
No own degrees of freedom (reinforcing element)
Input:
E [N/mm²]: Modulus of elasticity
I [mm4]: Moment of inertia
A [mm²]: Cross section area
ρ [t/mm³]: Density
Outputs:
- location and values for maximum and minimum bending moments
- visualisation in the graphics surface
Important note:
This beam approach can only be used for straight beams. When the edge beam would be
used for curved borders, the interaction between bending and torsion wouldn’t be possible!
Due to this setting, the usage of such a beam at curved borders is not allowed and an error
message will appear.
3.11 Elastic bonding element

The elastic edge bonding is mathematically used in the same way like the elastic edge
support. This elastic bearing is set at the edges with the thickness of each glass panel. The
width of the bonding material is given by “b”. The young’s modulus E and the shear modulus
G describe the elastic behaviour normal and transverse to the edge.
E ,G
When the shear modulus is set to zero, the shear behaviour will be ignored. Such, only
normal forces perpendicular to the glass edge can be transferred. The resulting stresses
exhibit thereby a characteristic:
Fnz
Fnt
Fnn
The stress components and their associated forces are related to the local co-ordinate
system along this edge! Such, the stresses σnt and σnn are components tangential and
normal to the edge. Theses stresses can only be compared with the global stresses σxx and

σxy when the global co-ordinate system is parallel to this local system. The component σnz is
the shear stress acting in z-direction and will result in forces Fz.
If de-attachment is set using the contact approach, the shear modulus G is automatically set
to zero, as shear effects can’t arise in such situations.
Degrees of freedom:
No additional degrees of freedom
Input:
E [N/mm²]: Young’s modulus of elasticity

G [N/mm²]: Shear modulus of elasticity
b [mm]: Width of the bonding
Contact tolerance [mm]: distance change, when separation shall take place
Output:
- Reaction forces Fnz, Fnt und Fnn for each border and layer
- maximum and minimum stresses σnz, σnt and σnn
- Symbols within the graphics surface
4 Basic approaches for calculation

4.1 Geometrically non-linear calculations
The layer wise approach leads to membrane forces (in plane normal forces), which are
necessary for equilibrium. For large transversal displacements the plate stiffness will
increase due to internal membrane effects, now also considered.
Regarding the displacement of a section dx, which is rotated by w',
w
du
dx +
w´
w w + w´ dx
x
dx
the following formulae
2
dx + du æ ¶w ö
= 1+ ç ÷ (3.1)
dx è ¶x ø
results. In second approximation the strains then result

2
¶u 1 æ ¶w ö
e xx = + ç ÷ (3.2)
¶x 2 è ¶x ø
and analogously
2
¶u 1 æ ¶w ö
e yy = + çç ÷÷ (3.3)
¶y 2 è ¶y ø
¶u ¶v ¶w ¶w
e xy = + + (3.4)
¶y ¶x ¶x ¶y
Solving these non-linear equations must take place with iterative procedures. The initial
stiffness matrix Ko, the tangential matrix Kt or combinations of both may therefore be used.
The first procedure with Ko has the advantage, that the stiffness matrix only once must be
assembled and can be used for further iterations. This is very fast, as the matrix triangulation
is the most time consuming process. But poor or no convergence may occur using only this
way for non-linear problems. Using the second solution of building a tangential stiffness
matrix is accompanied by good convergence but higher computational effort, as for every
iteration the stiffness matrix should be build anew.
In SJ MEPLA a combination of both methods is used: If a given convergence ratio will not be
achieved by use of a constant stiffness matrix, the over all tangential stiffness matrix is build
anew and is used for the next iterations. This procedure may be turned off by the button
<disable automatic>. Then the stiffness matrix is build and solved for each iteration, what
takes much more time. Another possibility is to apply the loads in several steps to the
system. Any loads like pressure, point, gravity and climatic load are subdivided and applied
to the system in x steps using the button <apply loads in x steps>. The solution of the
system in this case is the last step. This option can’t be used for the dynamic pendulum
impact calculation!
The tangential stiffness matrix is build out of 3 components. The initial stiffness Ko for small
deflections, the non-linear matrix Knl for large displacements and the stress matrix Kσ which
depends on the membrane forces.
K t = K o + K nl + Ks (3.5)
with
T
K o = ò A Bo D Bo da
T T T
K nl = ò A Bo D Bnl da + ò A Bnl D Bo da + ò A Bnl D Bnl da (3.6)
é N xx N xy ù
Ks = ò A G ê ú G da
T
êë N xy N yy úû
4.2 Contact approaches

In this program SJ MEPLA five possibilities are given to calculate contact conditions:
Point fixings:
- contact between the circular disks (shim) and the glass surface
- contact between the bush and the bore hole rim
- contact between the separators and the glass edges

Insulation glass:
- contact between the glass packages
Elastic edge and line supports:

- contact between the glass and the elastic line or edge supports
- contact separation at the edges
Elastic base:
- contact between the elastic base and the plate
Spacers in insulation glass:

- non-linear behaviour for tension and compression
These contact conditions are also possible for dynamic calculations like the pendulum impact
calculation.
Because the contact condition within insulation panes may cause solving problems, this
option may only be used in exception. The solving of such situations within static calculations
when the panes are overlapping is very difficult to perform.
For solving the contact problem between glass and the separators, the following function is
used, to describe a consistent separation:
When the distance between the glass and the separators is larger than the given tolerance
value, the stiffness of these layers is reduced to zero. Are the contact faces lying in between
the tolerance value, than the young's modulus is continuously changed to a factor of 1 (full
elastic modulus). So under compression (negative strains) always the given modulus is used.
The solution process is controlled and improved by the given tolerance value and cannot skip
undefined between two states.
The contact problem for insulation glass is much more complex by two stiff bodies under high
contacting velocity. The multi-layered elements here used, do not have any elasticity in
transverse direction and interface elements must be used. Using the glass material the
impact calculation would be possible, but will take place in very short time steps of a millionth
second. For practical use this way is herein not followed and instead an approximation is
done which will allow a small overlapping of the panes.
Therefore non-linear springs will be used. When contact condition is achieved (distance is
less than the tolerance value) springs are provided which rapidly increase their stiffness,

when the panes start overlapping. In this way the glass faces will return during iteration to
form the approximate contact of the faces.
The tolerance entry helps to introduce the contact smoothly, in order to increase it
continuously. This contact behaviour is comparable with an intermediate layer with the
thickness of the used tolerance and the possibility of a slight overlapping of the panes. Such
a system can be solved very much faster and therefore, is favoured here; but it represents an
approximation.
The solution of the contact problem requires from time to time a new creation and solving of
the stiffness matrix, if many relocation occur through large contact area which have to be
included into the stiffness matrix.
Note:
With increasing stiffness of separating layers the contact tolerance shall correspondingly be
reduced, what increases the requirements on the convergence of the solution process and
therefore the computation time can extend extremely. In some cases it can be possible that
no convergence can be achieved (the force error increases instead of the fact that in
converges to zero). Then the contact tolerance must set to higher values. Best solutions will
be achieved, when the pre-settings (0.001 to 0.0001 mm for rigidities of about 50 to 500
N/mm²) are used.
4.3 Dynamic calculation approaches

4.3.1 Mass matrix
The mass matrix is build through numerical integration of
n
M = å ò A H Ti mi Hi da
i =1
and summation over all layers i, with mi for the multi-layered element:
é1 0 0 0 0ù
ê ú
ê0 1 0 0 0ú
ê0 0 1 0 0ú
ê ú
ê ú
mi = r i t i ê 0 0 0 2
ti 0ú
ê ú
ê 12 ú
ê ú
ê0 0 0 0 2 ú
ti
ê ú
ë 12 û
4.3.2 Newmark procedure
Solving the differential equations of motion is done by use of Newmark´s time stepping
procedure.
In contrast to explicit procedures, here 3 equations of motion in the unknown state at time
t+Δt are denoted:
M &u&t + Dt + R t + Dt = P t + Dt (3.7)
u t + Dt = u t + Dt + D t b &u&t + Dt
pred 2
(3.8)

u& t + Dt = u& t + Dt + Dt d &u& t + Dt

pred
(3.9)
with
D t2
u
pred
t + Dt = u t + Dt u& t + (1 - 2 b ) &u&t (3.10)
2
u& t + Dt = u& t + Dt (1 - d ) u&& t
pred
(3.11)
With use of (3.10) and (3.11) a first approximation for displacement and velocity are
predicted for t+Δt (predictor-phase). After solving the equation system
K × Du = U
*
in the estimated state with the residual forces
U = P - M&u& - R
and the stiffness matrix
* M
K = + Kt
D t2 b
the displacement is build through summation of all incremental displacements in time step
Δt,
u t + Dt = u t + å D u
Using equation (3.8), the acceleration can be solved, which when inserted in (3.9) results to
a corrected velocity (corrector-phase). This iteration is then repeated during each time step
Δt until convergence is achieved. The constants δ and β controls the solution process. The
use of
1
d ³ 0.5 b = ( d + 0.5 )2
4
will lead to absolute convergence, but with numeric conditioned damping attenuation, which
can be avoided with δ = 0.5 (stability limit).
4.4 Pendulum body approach

If panes of glass are planned with drop securing function, e.g. for balustrades, a pendulum
experiment has to be carried out to simulate an impacting person. In this case, the pendulum
may not break through the pane of glass neither larger pieces of broken glass may separate
themselves and fall to ground.

glass panel
pendulum
drop height
Pendulum test for a parapet fixed at two opposite edges
The standard pendulum according DIN EN 12600 consists of two tires with 3.5 to 4.0 bar air
pressure and a mass of 50 kg.
The computational simulation of this experiment is a classical contact problem of two elastic
interacting bodies.
AR CR mp
Rp Rp Rp Rp
ws wp
Model of pendulum with deformations approaches
Since the contact occurs between the tires and the glass panel, the system is split at this
point. On the side of the glass plate, the pendulum force Rp acts with the corresponding
deflection of the glass plate ws. On the right side, it equilibrium exists within the tire spring of
rigidity CR and further equilibrium with the acceleration of the pendulum mass. This leads to
the 3 equations:
é ... ù
ê ú
ê ... ú
Ku + M&u& = P + ê ú (3.12)
ê R p
ú
ê ú
ë ... û
R p = C R ( w p - ws ) (3.13)
&& p = - R p
mp w (3.14)
Contact is achieved, when (wp - ws) < 0 and so the reaction force Rp = CR (wp - ws) is also
negative. The calculation of impact is carried out from positive z-direction with negative
velocity, which is computed from the given drop height. The impact force is distributed over

the tire foot prints, which depends on the deflection ΔwR = (wp - ws) of the tires. With the radii
of tires R1 = 51 mm and R2 = 194.5 mm
R2
R1 AR
Geometry of pendulum tires
the pressurised area AR can be calculated by use of ellipse equations due to the tire
deformation
AR = p (R 2
1 )(
- ( R1 - | D wR | )2 × R 22 - ( R 2 - | D wR | )2 )
The entire tire force RP is distributed uniformly over the tire foot print area by dividing in n x n
squares of the same size, where in each centre 0.5·RP/(n·n) is applied.
For calculation of equivalent nodal forces from single loads, acting anywhere within an
element and doesn't applied any node directly, a special procedure is used to calculate the
local position η,ξ out of global coordinates x and y.
The coupling of the equations of motion for the multi-layered element and the pendulum is
done be use of a sub-system. So the routines of the finite element program do not change
and the calculation of the pendulum impact may be switched on or off.
Therefore a single approach is used for the pendulum using the Newmark´s method. In
dependency from the relative deflection ΔwR, the interacting force RP (<0: contact or $0: no
contact with RP = 0) is coupled into the equations of motion for the plate and the pendulum.
The rigidity CR for the tire is determined from experiments which were carried out in [38, 47].
From acceleration measurements of the pendulum that impacts from different drop heights a
rigid wall has been used to get the non-linear rigidity of the spring, which was determined by
integrating the acceleration-time curve to
C R = 300 + 2 | D wR | [N/mm]
The impact velocity is derived by the drop height h:
v = 2 gh / 1000 [m / s ]
with : h [ mm ]
4.5 Kelvin-Voigt Model

A more enhanced impact model is described by this Kelvin-Voigt model. While the pendulum
impactor from 4.4 is using fixed values resulting from tests for representing the pendulum
behaviour according to DIN EN 12600, this new model is kept open for much more
possibilities:

It consists from 2 masses m1 and m2 which are internally connected by a spring – damper
system. The contact spring is described by a non-linear behaviour. The radius contacting the
plate surface can freely be defined.
Movement differential equations:
m2 × &z&2 = m2 ×g - c d ( z 2 - z1 ) - m d ( z&2 - z&1 )

m1 × &z&1 = m1 ×g + c d ( z 2 - z1 ) + m d ( z& 2 - z&1 ) - c k ( z1 - w )
with
R p = c k ( z1 - w)
Contact will arise when Rp <= 0 so when z1 < w whereas the contact body separates, when
z1 > w is.
By the deformation of the contacting spring Ck the sphere of radius R will enlarge the
contacting area A, which is described directly form the given radius R.
The contact spring ck is described by a non-linear stiffness condition as follows:
Ck = ck1 + ck2 * |w|
So two mechanical values can be adjusted to best simulate the impactor behaviour.
The considered dead weight of this 2 mass body and its actions onto the plate depends on
the gravity direction. For a vertical plate the gravity will not give any weight effects lateral to
the plate – but for a horizontal plate the full weight will be active.
The velocity is again given by the drop height (see section 4.4).

Input:
m1, m2 Mass [kg]

Ck Spring stiffness [N/mm]
R Radius of the contact spring [mm]
µd Damping [Ns/mm]
Cd Spring stiffness [N/mm]
Output:
- Reaction forces, acceleration and deflection at the impact position

- Deformations of the mass positions, the impact point and the reaction forces within
damper and springs are shown in the curve plot
- Symbols in the graphic surface show the model during impact
4.6 Principal stresses

The principal stresses are evaluated at the outer faces of a layer using the following
formulae:
2
s + s yy æ s - s yy ö
s = xx
h
+ çç xx ÷÷ + s xy2 major principal stress
2 è 2 ø
2
s + s yy æ s - s yy ö
s = xx
h
- çç xx ÷÷ + s xy2 minor principal stress
2 è 2 ø
The maximum principal stress will be used for structural analysis of glass plates.
4.7 Equivalent Von Mises stress

The equivalent Von Mises stress is used for metallic materials and evaluated by the
formulae:
sv = s xx2 + s yy2 - s xx ×s yy + 3 × s xy2
The stress components of transversal shear σxz and σyz are normally zero at the surfaces
and so are not considered. These mechanical properties are often used for structural
analyses of metal plates. Only in cases of sandwich structures a shear component at the
inner surfaces can appear. These stress components are here not considered. For laminated
glass plates the evaluation of VonMieses stresses doesn’t make any sense, as these values
can’t be used for glass.

4.8 Loading approaches

4.8.1 Face loads
Each glass package (a laminated compound of glass panels) can be loaded by face loads. A
face load affects the entire plate in the z-direction, so always perpendicular to the area and
without any holes or cut-outs.
V
p
z
y
The load vector is gained from numerical integration over all elements.
Beside the constant applied load from above as well a linear increasing load can be given.
For this purpose a first and second reference height in y-direction y0 and y1 with their related
pressures p0 and p1 must be defined. y1 – y0 must therefore be positive, so y1 shall lie over
y0.
p1
z y
p0
y0
x y1
This input describes the loading situation in between both reference lines The internal load is
again calculated by integrating every element which is lying completely within both lines, to
calculate the total nodal load vector. For elements which will be cut, a special iteration is
used.
The shown element is subdivided into several squares to approximate the loaded area. Due
to this approach, a little error may result for the totally applied force, which normally is less
than 0.5 Newton.
With this tool water pressure or culminations of snow can be simulated anywhere on the
plate.

4.8.2 Line loads

Defining the starting and ending point a line load can be defined. This line load has the unit
of [N/mm] or without conversion [kN/m].
qz
x0 ,y0
z
y
x1 ,y1
x
Along this line only the part of the line lying within the panel is considered. The load will be
separated into several not less than 20 segments and in a second condition equally spaced
with a minimum of 20 mm length. Along this line, the line load can be defined in 3 directions
(qx, q y, qz). For laminated glass each in-plane load is distributed according to the thickness of
each layer.
4.8.3 Border Lineload

Along any formed border line loads may be set. These lineloads are directly related to the
shape of the elements lying there. It’s possible to define 2 load directions qz and qn in units
N/mm. Qz is standing perpendicular onto the pane area and so is acting as a transversal
force. The component qn is acting perpendicular onto the cross section area of the edge. A
positive value acts away from the edge and so will induce tension stresses.
z
y
x qz
qz
qn
qn
The resulting nodal forces are integrated along the element border.
4.8.4 Concentrated loads

Any concentrated loads (or point loads) can be applied in 3 directions at any place. Applying
of moments is not possible.
The force Fz acts simultaneously onto all layers in accordance with the degree of freedom
approach for the multi-layered element (see section 3.1.1).

z Fz
y x ,y
Ly
Lx
x
Distribution of concentrated loads over the area Lx · Ly
Reference point is the middle of the area. The extension of the load area is defined in x
direction with the length Lx and in y direction with Ly [mm]. The given total loads Fx, Fy or Fz
[N] is distributed uniformly over the specified area. The forces Fx and Fy are distributed
according to the thickness of each layer, when laminated glass is used. In this way, it’s
possible at the edges of the panes, e.g. for glass fins, to place local in-plane forces (see as
well section 4.7.2 line loads and 4.7.3 border loads).
Ly
Lx
x
z
Concentrated loads in in-plane direction
In above picture, a load is shown in negative y-direction. If the reference point is placed close
to the edge (middle of the face) and the length of the load area Ly will be chosen very small a
line load is achieved; the load Fy is applied uniformly into the glass cross section.
The position of the reference point (middle of the load distribution area) is independent from
the element nodes. In the program, such related elements are automatically determined,
which are influenced from this loading.
This is achieved mathematically through a specific procedure, which first determines the
affected element and secondly the local coordinates which will receive a load part. The
numeric integration is performed by subdividing the applied loads into many partial loads,
which are distributed uniformly via the distribution area. The load distribution area is not
based on the element faces and in this way, is independent of the existing element mesh.
But, the accuracy of any finite element calculation depends (as always) on the mesh quality
and density. If a high accuracy for stresses and deflections in the region of a loaded area
should be achieved, this may need a finer mesh in order to represent correctly the local
effects. The number of elements is to balance between the necessary accuracy and the
tolerable computation time.
4.8.5 Dead weight

The calculation of the dead loads is carried out by integration the mass density and by
multiplication with the earth's gravity acceleration (9.81 m/s²), so that a volume load vector
results which is applied to the system.
For dead weight calculations the density of the layers material must be given. Glass has for
example a density of 2550 kg/m³, which must be given in MEPLA in units [tonnes/mm³], so
that 2.55e-9 [to/mm³] must be set for the density. This must be done in predefined units, to
avoid later conflicts when dynamic calculations are carried out.

Example:
v ·9.81
z y
o
45
x
Input of v = (0.0, -1.0, -1.0) which is converted within the program to (0.0, -0.707, -0.707 )
The pane is lying in x-y plane, but the dead weight (gravity vector) is acting as if the pane
has been rotated by 45° against the horizontal.
Alternatively it’s possible to define the angle of rotation about the x-axis in degree. This will
calculate the vector of acceleration automatically. When applying dead weight onto rotated
panes, not only perpendicular forces will arise! As well forces in plane will occur and must be
supported such, that a statically determined system is set!
4.8.6 Temperature constraints

A temperature difference in the layers of a laminated glass pane (or a general sandwich)
forms a constraining force state and secondary bending of the pane is caused. This
behaviour may be considered by taking a constant temperature per layer into account. A
temperature gradient within one layer (difference between top and bottom side of a layer)
can not be applied. A temperature gradient in a laminated pane can only be represented by a
layer wise constant temperature gradient.
T3
T2 = 0
T1
Constraints due to temperature differences in the layers
Each temperature difference will lead to an elongation or a shortening of the panel. When
they are glued together this change in length will lead to bending, depending on the
intermediate glue stiffness.
Necessary for this calculation is the coefficient of thermal expansion "T for the layer material.
4.8.7 Pressure hit

The use of pressure hit allows varying face loads and point loads to be considered over time.
This calculation is carried out dynamically. Measured or free data points can be given. The
value of the y- co-ordinate is taken as a load factor to be multiplied with previously given
loads, before they are applied to the system within the regarded actual time step length.
When a face load is set to 1.0, the y co-ordinate can be given in the unit of N/mm². When
measured values are given in kN/m² the values can be converted to N/mm² using a face load
value of 0.001. As well the loading co-ordinate can for example be regarded as load factor in
between the values 0.0 and 1.0 (100%).
The calculation shall always start with zero forces (factors), that no sudden load applying can
occur. As in reality it’s not possible to apply loads within no time frame!

The curve peek smoothing is a value which describes the rounding of sharp curve edges by
an inserted arc of given radius.
The time stepping will take place in the same way as shown in the card <pendulum
impact>. With the given time step length the values of the curve are interpolated and applied
onto the system until no more values can be read or the given duration has been reached.
With this card it’s possible to simulate wind blasts or detonations and their effects onto glass
panels.
This load option can be used with all other possibilities within the program except the
pendulum calculation. If both is set, the pressure hit calculation will get the priority.
4.8.8 Load case calculation

The calculation of load cases allows the combination of different loads, which all are applied
automatically within one calculation run. Calculations which are neither non-linear nor needs
iterations due to contact conditions, can be solved very fast, as the stiffness matrix must be
assembled only once and after a first solutions process only back substitutions are needed.
When non-linear calculations are set, the effort for solution will increase, as the system can’t
be solved by using constant stiffness matrices.
Additional to new setting of loads like dead weight, wind, snow, point loads or climatic loads
these values can be combined with safety factors, valid for each of those loads within the
load case set.
The setting of safety factor for climatic loads needs some additional explanations:
A climatic load consists of temperature changes within the gap filled with gas, the arising
pressure difference in the gap, which will result from changes of the barometric ambient air
pressure and a possible difference of height between assembly and installation of this
insulation glass, which will also change the ambient air pressure.
Climatic behaviour of insulation glass
The change of temperature within the gap can be multiplied directly by the safety factor, as
it’s the same way possible for the difference of height. As in MEPLA no approximations are
used and the real gas pressure law is regarded so that any system can be solved, the effect
of a safety factor g set onto pressure difference must be recalculated from the acting outside
and inside pressures, which are only used.
Therefore the following calculation for g to get a new factor g* is performed:
pa + ( pi - pa ) × g
g *=
pi
pa: barometric outside pressure
pi : inside pressure within the gap (the pressure when closing/sealing the gap)

This new factor g* can now be used to act directly onto the inside pressure pi to describe
originally the same effect for the factor g acting onto pressure difference.
pi * = pi × g *
The outside pressure so will be hold constant and the inside pressure is chosen such, that
the same effect is achieved.
When the factor g is set to zero, the pressure difference equals to zero and so the outside
pressure is equal to the inside pressure of the insulation glass.
When the factor is set to 1.0 the pressure difference between the preset load cases for pi and
pa (winter/summer condition) will be kept with original unchanged values.
4.9 Residual load-bearing capacity
With this approach to residual load capacity, broken panes, their crack formation and the
change in load-bearing and deformation behaviour can be approximately investigated.
4.9.1 Calculation approach

The calculation of the fracture behaviour is based on the principal stresses evaluated in the
Gaussian points.
At each Gaussian point of the 9-node element, the stress field is rotated in the respective
principal axis and the angle of rotation is determined.
In this rotated principal stress state, the principal stress components on the upper and lower
side are examined and compared.

If these principal stresses exceed the pre-set limit value, the breakage is activated at this
point (Gaussian point), where
a.) the break angle α for this section is determined, and
b.) only compressive forces at this crack surface (green) are further on transferred.
If tensile forces are detected at this Gaussian point, this section is used neither for the
transmission of normal force nor for the transmission of bending moments. The reaction
forces are corrected accordingly and used to form the imbalance vector.
However, if the calculation of the stresses at this current integration point shows a
compressive stress, this is still considered effective. However, bending stiffness is only
applied as long as no fracture state has been detected. In fracture state, the pressure point
will shift to the outside of the splinter according to the deformation (curvature).
This effect, which is only activated under pressure, in the cracked state and only in the
principal axes, is regarded here as additional elongation in the neutral fibre of the cracked
layer.
In the neutral axis, an additional change in length

t
Du = Dj ×
2
arises which, in relation to the length Δx of the regarded piece

t
Dj ×
Du 2
=
Dx Dx
is the strain in relation to the curvature φ’:

t
u¢ = j ¢ ×
2
This strain is applied analogue to the temperature strain as an additive term per Gaussian
point and thus forces the splinter into a position, as if it had to rotate around the now outer
pressure point. The position of the neutral axis of the fragment under normal stress is thus
maintained. However, the normal force is increased, resulting in a back drilling moment.
This approach represents an approximation of this effect, as the inner lever arm for the
pressure transmission (moment formation) is not put outward. If this approach were followed,
the normal force rearrangement within a splitter would also have to be considered, which
would then also depend on its size and a reduced pressure surface and thus changed
normal force stiffness would have to be taken into account.
In this approach, even within one element a "quasi" continuum mechanical behaviour can be
assumed and integrated via the element as usual. A back and forth bouncing neutral axis
would not allow this.
It should be borne in mind, however, that a loss of stiffness is also taken into account here if
a fracture condition is present and tensile stresses in the splinters are detected. In addition,
the entire fracture state is always viewed in the principal axes and from there always turned
back into the global axis system x-y. This also applies to the curvature φ’ and θ° in the plate,
from which the additional elongation for a broken piece is derived.
If the condition of a completely broken pane, e.g. FT glass in a crumbly fracture, is examined,
the above approach has no influence on the size of the broken fragments themselves.
However, this is different if a single crack formation (from a limit stress) is investigated. Here,
the curvature is located significantly on the broken section and thus increases the local angle
of rotation Δφ so that here the forced strains and their effect on other possibly still intact
layers is increased.
This position can be specified as a percentage dp for the calculation. With dp = 1.0 the
pressure point is on the upper or lower side; with dp = 0. the reset moment is switched off.
Mentally this corresponds to a splinter with a spherical fracture surface.

If a break was detected in a Gaussian point, it may still be a longer fragment which can still
transmit both bending and normal stresses in the longitudinal direction.
Only if this second principal axis (under 90°) also fails because the bending stresses exceed
also here the limit value, the glass segment is applied as a broken quadratic element in all
directions. Now it can only transmit compressive stresses in the two unchanging principal
axes.
Since this state can occur at each Gaussian point, changes in direction can also be
represented within one element:
vector plot of principal stresses (red: tension; blue: compression)
Tensile stresses (red) cannot be transmitted transversely to this crack orientation (green
lines). However, tension stresses (red) can still develop (in parallel to green cracking
directions) as long as only the first principal direction is broken.
4.9.2 Transformation of stress- and force vectors
From the deformation field u
e = B ×u
the strains in the investigated fracture layers are determined. With

s = D ×e
the stress vector in the Gaussian point is determined.
From the determined angle

2 ×t
a = arctan( )/2
sx -sy
for rotation into the principal axis, the transformation matrices are then created in order to be
able to transform the global x-y components into the principal rotated axis.
cos2 a sin 2 a 2 sin a × cosa

T= sin a
2
cos a
2
- 2 sin a × cosa
cosa × sin a - sin a × cosa cos2 a - sin 2 a
s1 sx
s2 = T × s y
t 12 t xy
After modification according to the crack formation from principal tensile stress excess, the
corrected stress components are turned back into the global x-y system by the inverse T-1
matrix.
The principal stress direction can change freely as long as no fracture has been detected.
This means that forces can be rearranged in all non-cracked glass sections, which influence
the further position and alignment of the crack fronts.
4.9.3 Basic assumptions

The following simplifications are used to calculate this post breakage behaviour:
- The breakage of a layer at the Gaussian point occurs when the applied break stress limit
in its principal stress component is exceeded.
- With the first break in the first direction of the principal stress, the break line is fixed at
this point and can no longer change.
- Along the break line, a glass splinter can still transmit bending until the limit stress is also
exceeded in the 2nd main axis.
- The calculation does not take place with real time influence (thus no dynamic procedure)
but only the sequence of the iterations is pursued, which can lead to the state of
equilibrium.
- The resulting stress redistribution due to a fracture condition prevents a convergent
solution of this ongoing crack condition.
- Only constant material behaviour is assumed - no viscous-elastic calculation (which
would only be possible considering again time effects). The result obtained thus only
represents the condition that can be found shortly after the crack formation. If several
layers have been applied as failing, the film is the only supporting element that in real
would continue to deform viscously over time.
- It is also not possible to perform a non-linear calculation, as the ongoing cracking will
destroy the formation of a stable tangential stiffness.
- To achieve useful results, a fine mesh (element size approx. 20 - 30mm) should be used.
- A convergent solution can only be achieved by tracking the deformations. Due to the
recurring imbalance forces caused by the progressive formation of cracks, the error force
of the imbalance forces cannot be used to follow the crack propagation.

- The dynamic, mass-influenced vibration behaviour due to a impact body cannot be

investigated up to now. This would require a real time step calculation in the nanosecond
range (possibly as a later extension!).
- No dynamic effects such as crack release energies or changes in the shape of the glass
splinters are applied.
- A difference between AN (float glass) and HS (heat strengthened glass) can only be
made via the limit break stress that may lead to a crack. The effect of the pre-stressing
on the crack formation cannot be taken into account.
- To simulate the behaviour of FT (fully toughened glass), the immediate and complete
break of a layer should be applied. This is done, for example, by specifying a (negative)
limit stress of - 200.0 N/mm². The batch of 0.0 N/mm² is not sufficient, since zones under
constant pressure would not break with this condition and would then continue to transmit
bending moments. However, the crumbling fracture of FT glass extends over the entire
layer and can therefore never transmit bending stresses at any position.
4.9.4 Comparisons
A comparison between a broken laminated glass and a simulation shows a quite good,
qualitative agreement:
Naturally, this simulation of fracture behaviour can only represent a qualitative process and
can never predict the exact fracture pattern.
5 Comparisons
5.1 Convergence study
The accuracy of a finite element calculation depends on the fineness of the used mesh. The
following examples shows in which way the mesh density affects the result.
The following system will be used:

- square plate with the size 1000x1000mm
- layers:
5.0 mm, E = 70000. N/mm², < = 0.23 (top layer)
5.0 mm, E = 5. N/mm², < = 0.40 (intermediate layer)
5.0 mm, E = 70000. N/mm², < = 0.23 (bottom layer)

- face load: -0.005 N/mm²

- simply supported
- non-linear calculation (large displacements)
Reference calculations with the finite element program MARC [32] by making use of 300 and
2400 volume-elements leads to the following results:
10x10x3 = 300 20x20x6 = 2400 elements

elements
displacement, centre of plate [mm]
-5.72 -5.73
2
stresses, centre top side [N/mm ]
-11.00 -11.04
2
stresses, centre bottom side [N/mm ]
15.37 15.18
Under variation of the mesh refinements, the 9-node multi-layered element in SJ MEPLA
leads to the following results:
mesh refinement 1x1 2x2 3x3 4x4 5x5 6x6 8x8 10x10
displacement, centre -4.05 -5.86 -5.69 -5.69 -5.70 -5.70 -5.70 -5.70
of plate [mm]
stresses, centre top -10.82 -11.58 -11.72 -9.95 -11.39 -10.48 -10.71 -10.82
side [N/mm2]
stresses, centre 6.51 23.02 14.10 16.21 14.78 15.50 15.30 15.22
bottom side [N/mm2]
The high accuracy of the 9-node element is shown. Using more than 25 elements, the
accuracy for the stresses is better than 2.6%. The deflections converge very much faster
onto the required value of 5.73 mm.
5.2 Pendulum impact on heat strengthened glass

Within the scope of the experiments in [47], some pendulum test for pre-stressed glass
(toughened glass) where carried out. In picture 5.1 the slope of the acceleration-time curve is
shown for a drop height of 700 mm onto a 8 mm glass panel:

300
Versuch
Beschleunigung des Pendels [m/s^2] Berechnung
240
180
120
60
0
-0,005 0,005 0,015 0,025 0,035 0,045 0,055
Zeit [s]
Fig 5.1: comparison between experimental and calculated acceleration
The comparison between the calculated acceleration and the measured curve shows the
high accordance with the not broken glass.
240
Berechnung horizontal
Berechnung vertikal
Versuch horizontal
180 Versuch vertikal
Spannungen [N/mm^2]
120
60
-60
0 0,016 0,032 0,048 0,064 0,08 0,096
Zeit [s]
Fig 5.2: experimental data versus calculated stress results
The stresses in picture 5.2 are measured in orthogonal directions at the impact point and can
be simulated very well by finite element calculation.
This tests indicates that for short time effects like a pendulum impact, a higher glass strength
can be used, as this glass wasn't broken.
Allowable stresses for the pendulum impact (according to German guideline TRAV):
Float glass: 80 N/mm²

Strengthened glass: 120 N/mm²
Heat toughened glass: 170 N/mm²

Due to the short impact duration, such higher stresses do not generally lead to failure of the
panes. The onsite pendulum test is passed, when no larger broken fragments fall to ground.
Therefore an additional safety in pendulum calculation exists, as a simulation can only be
done until a possible failure - but no further ongoing breakage and energy consumption.
5.3 Pendulum impact on laminated safety glass

In [47] also safety glass (laminated glass) supported in a standard frame where tested. They
consist of 2x3mm float glass panes with 0.38 mm butatice (PVB). The drop height has been
450 mm.
210
Beschleunigung Mittelwert
180 monolithisch t = 6.38 mm
VSG mit G = 1 N/mm^2
monolithisch t = 6.0 mm
Beschleunigung [m/s^2]
150
VSG mit G = 4 N/mm^2
120
90
60
30
-30
0 0,01 0,02 0,03 0,04 0,05 0,06 0,07 0,08 0,09
Zeit [s]
Fig 5.3: acceleration-time slope for pendulum impact on 6 mm strengthened glass
In picture 5.3 several approaches for the simulation of laminated glass panes versus
experimental data are shown. The black dashed lines show the measured acceleration at
positions over and under the double tires. The unbroken black line characterizes the mean
value of these accelerations.
In reality, laminated panes by PVB behave visco-elastically. The stiffness of butatice
depends on the load, the deformation state and the temperature. It could be shown [1,26],
that for short term load an extreme stiffening effect occurs and that the laminated panes and
then act as a monolithic one.
Therefore some calculations are represented in picture 5.3, where the attempt was made to
best approximate the test results by using a linear material law instead of the real viscous-
elastic behaviour. The best compliance has been found for a monolithic approach, using the
total thickness without the compound layer. Even an approach for a laminated glass pane
with a shear modulus of G = 4 N/mm² still brought no sufficient agreement. The shear
modulus activated is even higher so that laminated safety glass under pendulum impact can
best be simulated in a monolithic approach.
5.4 Hailstone Force Path
In a study on the effect of hail stones on air traffic [55] the force-time progression of different
hailstones are shown in a diagram.

The 9g heavy hailstone (V2_06) with a speed of 144 m/s and a diameter of 24mm was
examined in detail and the characteristics of the extended Kelvin-Voigt model were adjusted
to maximum conformity.
Determined characteristic values:
m1 [kg] m2 [kg] Cd [N/mm] µd [Ns/mm] R [mm] Ck1 [N/mm] Ck2 [N/mm]
1.35e-4 8.865e-3 1. 0.08 12. 4000. 10.
The drop speed of 144m/s was generated from the drop height of 1057 m = 1057000 mm.
The very short pulse of max. 0.8ms was mapped with a time step length of 0.000001s.
Superposition between test results [55] (dark blue line) and simulation:
The comparison shows not only a good correlation of the maximum force pulse (9.3kN) but
also of the time curve, which decays to zero up to approx. 0.4ms and then oscillates only
slightly around the zero axis in the test. The small superimposed vibrations that may have
arisen from the vibration of the load cell itself can't be shown here.

Legal Notice for literature usage [55]:
6 Notes
1. In the program MEPLA, the units to be used are given. Towards the usual free choice of
units in finite element programs, we here only use Newton and millimetres (if nothing
other is indicated) in order to avoid conversion mistakes in particular for dynamic
calculations. From this default a unit of [tonnes/mm³] results e.g. for the density.
2. No further loads shall be applied for the pendulum simulation (dynamic calculation)
except the pendulum body or the pressure hit loads itself. An additional concentrated
load for the pendulum impact would lead to a suddenly applied load which is set on at
infinitely short time and will lead to sudden vibration just before the impact will act.
3. Calculating of several glass packages with in each case several layers however, is
indeed possible, but one ought to consider that the time required for the solution can
increase extremely under contact approaches, especially in the case of a finer element
mesh or dynamic calculations.
4. For dynamic pendulum calculations full shear effect for laminated glass panels can be
used by setting only one monolithic layer with the sum of all glass layers. This will
reduce the required time for solving drastically.
5. In some extreme cases, e.g. using a lot of point fixings, small contact tolerances or a
system which is nearby static indefiniteness (may be caused by contact conditions), it
can occur that the solving process does not converge. Any calculation will end
automatically after at least 6000 unsuccessful iterations (see User´s Manual).
6. The automated generation of the point fixing elements needs a specific space in order to
produce the element mesh. The mesh for the point fixing of type 1, 2, 7 - 10 needs a
small distance (about 20 mm) between the disk and the pane border. Before an
extensive calculation, we always recommend to check the mesh via <System preview>.

7 Literature
[1] Bohmann, D.; Ein numerisches Verfahren zur Berechnung von Verbundglasscheiben,
Shaker Verlag Aachen, Schriftenreihe - Stahlbau, RWTH Aachen, Heft 43 (1999)
[2] Pica A.; Wood R. D.; Hinton E. Finite element analysis of geometrically nonlinear plate
behavior using a Mindlin formulation; Computer&Structures Vol. 11 pp. 203-215 (1980)
[3] Praktische Rheologie der Kunststoffe und Elastomere / VDI-Gesellschaft Kunststoff-
technik; Düsseldorf VDI-Verlag, (1991)
[4] Zienkiewicz O. C.; Methode der finiten Elemente; München,Wien Hanser (1984)
[5] Nádai A.; Die elastischen Platten; Berlin, Springer Verlag (1925,1968)
[6] Bathe K. J.; Finite-Elemente-Methoden, Berlin, Springer Verlag (1986)
[7] Donea J.; Lamain L. G.; A modified representation of transverse shear in Co
quadrilateral plate elements; Comp. Meth. Appl. Mech. Engng., 63, 183-207 (1987)
[8] Prathap G.; Ramesh Babu C.; Field-consistent strain interpolations for the quadratic
shear flexible beam element; Int. J. Numer. Meth. Engng., Vol 23 pp. 1973-1984 (1986)
[9] Prathap G.; Naganarayana B. P.; Somashekar B. R.; Field-consistency analyses of the
isoparametric eight-node plate bending element; Computer&Structures, Vol 29 No. 5 pp.
857-873 (1988)
[10] Prathap G.; Somashekar B. R.; Field- and edge-consistency synthesis of a 4-noded
quadrilateral plate bending element; Int. J. Numer. Meth. Engng. Vol. 26, 1693-1708 (1988)
[11] Whitney J. M.; Stress analyses of thick laminated composite and sandwich plates; J.
Composite Materials, Vol. 6 pp. 426-440 (1972)
[12] Bathe K.-J.; Dvorkin E. N.; Short communication, A four-node plate bending element
based on Mindlin/Reissner plate theory and mixed interpolation; Int. J. Numer. Meth. Engng.,
Vol. 21, 367-383 (1985)
[13] Kebari H.; A one point integrated assumed strain 4-node Mindlin plate element; Eng.
Comput., Vol. 7, 284-290 (1990)
[14] Huang H. C.; Hinton E.; A 9-node Lagrangian Mindlin plate element with enhanced
shear interpolation; Eng. Comput., Vol. 1, 369-379 (1984)
[15] Hinton E.; Huang H. C.; Shear forces and twisting moments in plates using Mindlin
elements; Eng. Comput.; Vol. 3, 129-142 (1986)
[16] Timoshenko S. P.; Woinowsky-Krieger S.; Theory of Plates and Shells, McGraw-Hill,
New York (1959)
[17] Owen D. R. J.; Hinton E.; Finite Elements in Plasticity: Theory and Practice; Pineridge
Press Limited, Swansea, U.K. (1980)
[18] Mukhopadhay M.; Analyses of plates using isoparametric quadratic element-shear
reactions, patch loading and some convergence studies; Computer&Structures, Vol. 17 No.
4, 587-597 (1983)
[19] Geier B.; Composite Laminate Stiffnesses and Their Sensitivities; Forschungsbereicht,
DLR, Institut für Strukturmechanik, Braunschweig, DLR-FB 92-16
[20] Stiftinger M. A.; Semi-Analytical Finite Element Formulation for Layered Composite
Shells with Consideration of Edge Effekts; Fortschrittsberichte VDI Verlag, Reihe 18: Mecha-
nik/Bruchmechanik, Nr. 203
[21] Becker W.; Beiträge zur analytischen Behandlung ebener Laminate; Fortschrittsberich-
te VDI Verlag, Reihe 18: Mechanik/Bruchmechanik, Nr. 121
[22] Reddy J. N.; Robbins D. H. Jr.; Theories and Computational Models for Composite
Laminates; Appl. Mech. Rev., Vol. 47, No. 6, 147-169 (1994)
[23] Sarma B. S.; Varadan T. K.; Certain discussions in the finite element formulation of
nonlinear vibration analyses; Computer&Structures; Vol. 15 No. 6, pp. 643-646 (1982)
[24] Abrate S.; Sun C. T.; Dynamic analyses of geometrically nonlinear truss structures;
Computer&Structures, Vol. 17, No. 4, pp. 491-497 (1983)
[25] Chow T. S.; On the propagation of flexural waves in an orthotropic laminated plate and
its response to an impulsive load; J. Composite Materials, Vol 5 p. 306 (1971)
[26] Sobek W., Kutterer M., Messmer R. ; Rheologisches Verhalten von PVB im Schubver-
bund; Forschungsbericht 4/98, Bauen mit Glas, Institut für leichte Flächentragwerke, Stutt-
gart (1998)

[27] Gross D.; Hauger W.; Schnell W.; Wriggers P; Technische Mechanik; Band 4: Hydro-
mechanik, Elemente der höheren Mechanik, Numerische Methoden; Springer-Verlag (1995)
[28] Mindlin R. D.; Influence of rotary inertia and shear on flexural motions of isotropic
plates; J. Appl. Mech., 18, 31-38 (1951)
[29] Beyle H. P.; Prüfbericht Nr. 98 - 3022; Labor für Stahl- und Leichtmetallbau, FH Mün-
chen, Prof. Dr.-Ing. Ö. Bucak, Auftraggeber: Verroplan GmbH (1998)
[30] Huang H.-C.; Static and Dynamic Analyses of Plates and Shells; Springer Verlag
[31] ABAQUS; Version 5.7, Hibbitt, Karlsson & Sorensen, Inc. Pawtucket, USA (1998)
[32] MARC/MENTAT;Version 7.2/3.22, MARC Analysis Research Corporation, Palo Alto,
USA (1998)
[33] Maksimovic S., Komnenovic M; Improved nonlinear finite element analysis of layered
composite structures using third-order theory, CIVIL-COMP Ltd. Edinburgh, Scottland,
Advances in non-linear finite element methods, 251-261 (1994)
[34] Campos A.D., Kurban A.E.A.; Vibrations of laminated beams, CIVIL-COMP Ltd.
Edinburgh, Scottland, Advances in non-linear finite element methods, 243-250 (1994)
[35] Bisegna P., Sacco E.; A layer-wise laminate theory rationally decuded from the three-
dimensional elasticity, J. Appl. Mech., Vol 64, 538-545 (1997)
[36] Levinson M.; An accurate, simple theory of the statics and dynamics of elastic plates,
Mech. Res. Com., Vol 7, 343-350 (1980)
[37] Hinton. E; The flexural analysis of laminated composites using a parabolic isometric
plate bending element, Numer. Meth. Engng., Vol 11, No 1, 174-179 (1975)
[38] Sedlacek G., Bohmann D., Schäfer E.; Untersuchung zum Schwingverhalten von Ver-
bundglasscheiben, Lehrstuhl für Stahlbau, RWTH Aachen; (unveröffentlicht) (1998)
[39] Gawehn W.; Finite-Elemente-Methode, Lehrbuch, Vieweg, Braunschweig (1985)
[40] Gummert P., Reckling K.-A.; Mechanik, Vieweg, Braunschweig (1986)
[41] HT Troplast AG; Technische Daten zur Trosifol PVB-Folie
[42] Sedlacek G., Hortmanns M.; Wellershoff F.; Windkanaluntersuchung zur Bestimmung
der Fassadenlasten am Staatlichen Bauamt Dortmund, Lehrstuhl für Stahlbau, RWTH Aa-
chen (1998)
[43] Deutsches Ingenieurblatt, Neue Produkte, Mitteilungen von HT Troplast und DuPont,
Dezember (1998)
[44] Güsgen, J.; Bemessung tragender Bauteile aus Glas, Shaker Verlag Aachen, Schrif-
tenreihe - Stahlbau, RWTH Aachen, Heft 42 (1998)
[45] Fa. Wendker & Selders, Steinglas-Verbundsysteme, patentierte Marmor- Glas Ver-
bundanwendung (1998)
[46] Sedlacek G., Blank K., Laufs W., Güsgen J.; Glas im konstruktiven Ingenieurbau,
Ernst&Sohn (1999)
[47] Wörner J.-D., Schneider J.; Zwischenbericht zur Durchführung von Versuchen zum
weichen Stoß mit dem Pendelkörper nach DIN EN 12600, Institut für Statik, Darmstadt
(1999)
[48] Pölling, R; Grundsatzuntersuchung von Einscheiben- und Verbundglas als Tragele-
ment, Diplomarbeit, Ruhruniversität Bochum, (1996)
[49] Schutte A., Hanenkamp W.; Zum Tragverhalten von Verbund- und Verbundsicher-
heitsglas bei erhöhten Temperaturen unter Einwirkung von statischen und dynamischen Las-
ten, Ernst&Sohn, Bautechnik 76, Heft 1, (1999)
[50] Laufs, W.; Ein Bemessungskonzept zur Festigkeit thermisch vorgespannter Gläser,
Shaker Verlag Aachen, Schriftenreihe - Stahlbau, RWTH Aachen, Heft 45 (2000)
[51] Dercks, Ch.; Isolierglasscheiben unter Klimalasten, Diplomarbeit RWTH Aachen, Lehr-
stuhl für Stahlbau, (1997)
[52] Bohmann, D.; FEM – New Ways in Glass Design and Engineering, Internet Forum der
Glass Processing Days (www.glassfiles.com), Tampere, Finland, 2001
[53] Schneider, J.; Bohmann D.; Glasscheiben unter Stoßbelastung, Experimentelle und
theoretische Untersuchungen für absturzsichernde Verglasungen bei weichem Stoß, Bauin-
genieur, Springer Verlag, Band 77, Dez, 2002
[54] https://en.wikipedia.org/wiki/Barometric_formula
[55] European Aviation Safety Agency, Research Project EASA.2008/5, Hail Threat
Standardisation, FINAL report for EASA.2008.OP.25

SJ Mepla Manual Theory Eng

Uploaded by

Copyright:

Available Formats

SJ Mepla Manual Theory Eng

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

SJ Mepla Manual Theory Eng

Uploaded by

Copyright:

Available Formats

SJ MEPLA

SJ Software GmbH Aachen

SJ Software GmbH Aachen

SJ Software GmbH Aachen

2. According to the Mindlin plate-theory the displacement vector u = [u,v,w]T is expressed

u 2 ( h2 /2) = u 3 (- h3 /2) ; u 2 (- h 2 /2) = u 1 ( h1 /2)

SJ Software GmbH Aachen

Each node of a n-layered element has 4n+1 degrees of freedom.

SJ Software GmbH Aachen

and the transverse shear strains

SJ Software GmbH Aachen

¶ j iz 1 æ ¶ uio+1 t i+1 ¶ j i+1 ¶ uio-1 t i -1 ¶ j i -1 ö

3.1.2 9 node multi-layered element

3x3 = 9 Gaussian points, (full integration)

SJ Software GmbH Aachen

Load vector for each node of each package:

Element building, order of layer numbering:

3.2 Building of the finite element mesh

SJ Software GmbH Aachen

The mesh refinement is controlled via the given element size.

The following basics must be regarded to describe the system contour:

3.3 Point Fixing

SJ Software GmbH Aachen

Within this program 10 different kinds of fittings can be used:

SJ Software GmbH Aachen

Type 1: Countersunk fixing

Type 2: Disk fixing

Type 3: Circular clamp (balustrade clip)

SJ Software GmbH Aachen

Type 4: Angular clamp (balustrade clip)

Type 5: Circular downholder

SJ Software GmbH Aachen

Type 6: Angular downholder

The angular downholder is equivalent in function to the circular downholder.

Type 7: Bonded disk fixing

Type 8: Countersunk fixing within layer 1, (laminated safety glass) LSG

Type 9: Countersunk fixing for package 1, insulated glass unit

SJ Software GmbH Aachen

Type 10: Disk fixing for package 1, insulated glass unit

Building the mesh:

SJ Software GmbH Aachen

Clamping with offset to the edges

The arising offset (see picture) is considered in the calculation.

SJ Software GmbH Aachen

Force transmission for countersunk fixing

Force transmission for the disk fixing

SJ Software GmbH Aachen

The bonded disk fixing forms an exception:

The default (unmarked) is the rotational stiff setting!

The double sided disk fixing is not underlying this restriction.

Point fixings in insulation glass:

Spacer at the borehole rim

Load transferring at the borehole rim:

SJ Software GmbH Aachen

3.3.3 Connection of point fixings

Definition of the tie bar position

3.3.4 Forces applied at point fixings

SJ Software GmbH Aachen

5 degrees at the base point, (u, v, w, j, q)