Simple Numerical Model of Laminated Glass Beams

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


Introduction
The most frequently used transparent material in the building envelopes is glass. It is a fragile material, which fails in a brittle manner. This is the reason for using safety glasses in a situation when there is a possibility of human impact or where the glass could fall if shattered.
Laminated glass is a multi-layer material produced by bonding two or more layers of glass together with a plastic interlayer, typically made of polyvinyl butyral (PVB). The interlayer keeps the layers of glass bonded even when broken, and its high strength prevents the glass from breaking up into large sharp pieces. This produces a characteristic "spider web" cracking pattern when the impact is not powerful enough to completely pierce the glass. Multiple laminae and thicker glass decrease stress level, thereby increasing the load-carrying capacity of a structural member, too.
The focus of this study is on the establishing a simple and versatile framework for the analysis of mechanical behavior of laminated glass units. To keep the discussion compact, we restrict our attention to the linearly elastic response of layered glass beams in the small strain regime. The rest of the paper is organized as follows. Methods of analysis of laminated glass beams are introduced in Section 2, together with a brief characterization of the proposed numerical model. The principles of the method are described in detail in Sections 3 and 4. In particular, the mechanical formulation of the model is shown in Section 3. In the next section, the Finite Element discretization is presented. In Section 5, the proposed numerical technique is verified and validated against a reference analytical solution and publicly available experimental data. Finally, Section 6 concludes the paper and discusses future extensions of the method.

Brief overview of available methods
The most frequent approach to the analysis of glass structural elements was, for a long time, based on empirical knowledge. Such relations are sufficient for the design of traditional windows glasses. In modern architecture there has been a steadily growing demand on the use of transparent materials for large external walls and roof systems in the recent decades. Therefore, the detailed analysis of layered glass units is becoming increasingly important in order to ensure a reliable and efficient design.
In general, the complex behavior of laminated glass can be considered as an intermediate state of two limiting cases [1]. In the first case, the structure is treated as an assembly of two independent glass plates without any interlayer (the lower bound on stiffness and strength of a member), while in the second case, corresponding to the upper estimate of strength and stiffness, the glass unit is modeled as a monolithic glass (one glass plate with thickness equal to the total thickness of the glass plates). Both elementary cases, however, fail to correctly capture complex interaction among individual layers, leading to non-optimal layer thickness designs. Therefore, several alternative approaches to the analysis of layered glass structures have been proposed in the literature. These methods can be categorized into three basic groups: • methods calibrated with respect to experimental measurements [2], • analytical approaches [3,4,5], • numerical models typically based on detailed Finite Element simulations [6,7].
Applicability of analytical approaches to practical (usually large-scale) structures is far from being straightforward. In particular, the closed-form solution of the resulting equations is possible only for very specific boundary conditions and therefore have to be solved by an appropriate numerical method. Moreover, the analytical approaches are rather difficult to be generalized to beams with multiple layers. Therefore, it appears to be advantageous to directly formulate the problem in the discretized form, typically based on the Finite Element Method (FEM). Nevertheless, we would like to avoid fully resolved two-or three-dimensional simulations, cf. [6,7], which lead to unnecessarily expensive calculations.
In this paper, we propose a simple FEM model inspired by a specific class of refined plate theories [8,9,10]. In this framework, each layer is treated as the Timoshenko beam with independent kinematics. Interaction between individual layers is captured by the Lagrange multipliers (with a physical meaning of shear stresses), which result from the conditions of inter-layer displacements compatibility. Such a refined approach circumvents the limitation of similar models available in typical commercial FEM systems, which employ a single set of kinematic variables and average the mechanical response through the thickness of the beam, e.g. [11]. Unlike the proposed approach, the averaging operation is too coarse to correctly represent the inter-layer interactions, see Section 5 for a concrete example.

Mechanical model of laminated beams
As illustrated in Figure 1, laminated glasses consist mostly of three layers. A local coordinate system is attached to each layer to allow for an efficient treatment of independent kinematics. In the following text, a quantity a expressed in the local coordinate system associated with the i-th layer is denoted as a (i) , whereas a variable without an index represents a globally defined quantity, cf. Figure 1.
Each layer is modeled using the Timoshenko beam theory supplemented with membrane effects. Hence, the following kinematic assumptions are adopted • the cross sections remain planar but not necessarily perpendicular to the deformed beam axis, • vertical displacement does not vary along the height of the beam (when compared to the magnitude of the displacement).
Under these assumptions, the non-zero displacement components in each layer are parametrized as: where i = 1, 2, 3 and z (i) is measured in the local coordinate system from the middle plane of the i-th layer. The inter-layer interaction is ensured via the continuity Figure 1: Kinematics of laminated beam conditions specified on interfaces between layers in the form (i = 1, 2) The strain field in the i-th layer follows from the strain-displacement relations [12,11] which, when combined with the constitutive equations of each layer expressed in terms of Young's modulus E and the shear modulus G: yield the expressions for the internal forces as where b and h (i) are the width and height of the i-th layer, recall Figure 1, and k = 5 6 , A (i) = bh (i) a I (i) = 1 12 b(h (i) ) 3 stand for the shear correction factor, the cross-section area and the moment of inertia of the i-th layer, respectively.
To proceed, consider the weak form of the governing equations, written for the i-th layer (the subscripts • x and • z related to internal forces and kinematics-related quantities are omitted in the sequel for the sake of brevity) to be satisfied for arbitrary admissible test fields δu (i) , δϕ (i) and δw. In particular, the first three equations correspond to equilibrium conditions written for normal and shear forces and bending moments, respectively. The last identity enforces the geometrical relation (2) in the integral form, thereby allowing to treat the shear strain as an independent field in the discretization procedure discussed next. Further note that the continuity conditions (1) will be introduced directly into the discretized formulation, as explained in the following Section.

Finite element discretization
To keep the discretization procedure transparent, it is assumed that each layer of the laminated beam is divided into identical number of elements, leading to the discretization scheme illustrated in Figure 2.
Following the standard conforming Finite Element machinery, e.g. [12,11], we express the searched and test displacement fields at the element level in the form where e is used to denote the element number, • e and δ• e denote a relevant searched and test field restricted to the e-th element, N (i) e,• is the associated matrix of basis functions and r (i) e,• the matrix of nodal unknowns. In the actual implementation, the fields u (i) , w e and ϕ (i) e , as well as the corresponding test quantities, are assumed to be piecewise linear. To obtain a locking-free element, the shear strain γ (i) e is taken as constant and is eliminated using the static condensation, see [12,11] for additional details.
To simplify the further treatment, we consider the following partitioning of the stiffness matrix K and the right hand side matrix R related to the e-th element and the i-th layer after the static condensation: T , r e,w = w e,a , w e,b T .
Considering all three layers in Figure 2 together gives the resulting system of governing equations in the form where the matrix λ stores the nodal values of the Lagrange multipliers, associated with the compatibility constraint (1), and the matrix implements the tying conditions.  To verify and validate the performance of the present approach, the previously described FEM model was implemented using MATLAB R system and compared against predictions of an analytical model and experimental data for a three-point bending test on a simply supported laminated glass beam presented in [5], see also Figure 3. The width of the beam is b = 0.1 m and material data of individual components of the structure are available in Table 1. The glass modulus of elasticity is slightly lower than the conventional values of 70-73 GPa reported in the literature and is specific for the material employed in the experiment. Moreover, as the PVB layer shows viscoelastic and temperature-dependent behavior, the modulus of elasticity corresponds to an effective secant value related to load duration of 60 s and temperature of 22 • C.  Table 2 summarizes values of the mid-span deflection for a representative load level determined by FE-based discretization using 60 elements (30 when symmetry of the problem is exploited) and the corresponding experimental values. Note that the discretization is sufficient to achieve three-digit accuracy in the midspan deflection. In addition to the results obtained by an analytical method proposed by Asik and Tezcan in [5], the results of the analysis using ADINA R system and the elementary lower and upper bounds are included. In particular, the ADINA R model is based on the classical laminate theory, cf. [11], whereas the two simplified approaches assume zero or infinite compliance of the interlayer, recall also discussion in Section 2. In the following discussion, e.g. η num exp denotes the relative error between the numerical prediction and reference experimental value, while e.g. η an is used for the error of analytical solution when compared to candidate approaches. Clearly, the results of the last three methods differ substantially from experimental data as well as the analytical results. The proposed numerical model, on the other hand, shows a response almost identical to the analytical method, which deviates from experimental measurement by less then 6%. Such accuracy can be considered as sufficient from the practical point of view. To further confirm predictive capacities of the proposed numerical scheme, a response corresponding to a proportionally increasing load was investigated. The results appear in Tables 3 and 4. Again, the method seems to be sufficiently accurate in the investigated range of loads when considering the values of deflections as well as values of local stresses and strains.

Conclusions
As shown by the presented results, the proposed numerical method is wellsuited for the modeling of laminated glass beams, mainly because of its low computational cost and accurate representation of the structural member behavior.  Future improvements of the model will consider large deflections and the timedependent response of the interlayer and will be reported separately.