Comparison of viscoelastic finite element models for laminated glass beams

TL;DR

This paper compares four viscoelastic finite element models for laminated glass beams, finding that simpler models are often sufficient but may be unsafe under rapid loading changes.

Contribution

It introduces and validates four layerwise beam theories with different non-linear and viscoelastic assumptions for laminated glass.

Findings

All four models predict similar responses.

Simpler elastic models are generally adequate.

Complex models are necessary for rapid load changes.

Abstract

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

Tables

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

Abbreviation	FSK	FSν	VKK	VKν	LINK	LINν
Kinematics	finite strains		large deflections		geometric linearity
Constant	$K^{\left(2\right)}$	${\nu }^{\left(2\right)}$	$K^{\left(2\right)}$	${\nu }^{\left(2\right)}$	$K^{\left(2\right)}$	${\nu }^{\left(2\right)}$

Table 2. Table 2: Overview of geometry of laminated glass beams

Geometry	Supports	Length [m]	Width [m]	Thicknesses [mm]
				glass/PVB/glass
I	fixed-end	3	0.15	3/0.76/3
II	simply-supported	3	0.15	6/0.76/6

Equations

$$u^{(i)}(x,z^{(i)})$$

$$w^{(i)}(x,z^{(i)})$$

$$u^{(i)}(x,\frac{h^{(i)}}{2})-u^{(i+1)}(x,-\frac{h^{(i+1)}}{2})$$

$$w^{(i)}(x,z^{(i)})-w^{(i+1)}(x,-\frac{h^{(i+1)}}{2})$$

$$u_0^{(i)}(x)-u_0^{(i+1)}(x)+\frac{h^{(i)}}{2}\sin \left({\phi }_y^{(i)}(x)\right)+\frac{h^{(i+1)}}{2}\sin \left({\phi }_y^{(i+1)}(x)\right)$$

$$w_0^{(i)}(x)-w_0^{(i+1)}(x)+\frac{h^{(i)}}{2}\cos \left({\phi }_y^{(i)}(x)\right)+\frac{h^{(i+1)}}{2}\cos \left({\phi }_y^{(i+1)}(x)\right)-\frac{h^{(i)}}{2}-\frac{h^{(i+1)}}{2}$$

$${\epsilon }_x^{(i)}(x,z^{(i)})$$

$${\gamma }_{xz}^{(i)}(x)$$

$${\epsilon }_{x,0}^{(i)}(x)$$

$${\kappa }_y^{(i)}(x)$$

$$u^{(i)}(x,z^{(i)})$$

$$w^{(i)}(x)$$

$$u_0^{(i)}(x)-u_0^{(i+1)}(x)+\frac{h^{(i)}}{2}{\phi }_y^{(i)}(x)+\frac{h^{(i+1)}}{2}{\phi }_y^{(i+1)}(x)$$

$$w_0^{(i)}(x)-w_0^{(i+1)}(x)$$

$${\epsilon }_x^{(i)}(x,z^{(i)})$$

$${\gamma }_{xz}^{(i)}(x)$$

$$\displaystyle\varepsilon_{x,0}^{(i)}(x)=\mbox{$\displaystyle\frac{\,{\mathrm{d}}u_{0}^{(i)}}{\,{\mathrm{d}}x}$}(x)+\tfrac{1}{2}\left(\mbox{$\displaystyle\frac{\,{\mathrm{d}}w_{0}^{(i)}}{\,{\mathrm{d}}x}$}(x)\right)^{2},$$

$$\log a_T=-\frac{C_1(T-T_0)}{C_2+T-T_0},$$

$$\mathbf{\epsilon }=\left[\begin{array}{cccccc}{\epsilon }_x& {\epsilon }_y& {\epsilon }_z& {\gamma }_{xy}& {\gamma }_{xz}& {\gamma }_{yz}\end{array}\right]{}^{\mathsf{T}},$$

$$\mathbf{\epsilon }(t)=\frac{1}{3}{\epsilon }_{\mathrm{V}}(t)\mathbf{i}+\mathbf{e}(t),$$

$$\mathbf{\sigma }(t)={\sigma }_{\mathrm{m}}(t)\mathbf{i}+\mathbf{s}(t),$$

$$\mathbf{e}=\left[\begin{array}{cccccc}{e}_x& e_y& e_z& {\gamma }_{xy}& {\gamma }_{xz}& {\gamma }_{yz}\end{array}\right]{}^{\mathsf{T}},$$

$${\boldsymbol{\sigma}}(t)=\sigma_{\mathrm{m}}(t){\boldsymbol{i}}+{\boldsymbol{s}}(t)={\boldsymbol{i}}\int_{0}^{t}K(t-t^{\prime})\mbox{$\displaystyle\frac{\,{\mathrm{d}}\varepsilon_{V}}{\,{\mathrm{d}}t^{\prime}}$}(t^{\prime})\,{\mathrm{d}}t^{\prime}+2{\boldsymbol{Q}}\int_{0}^{t}G(t-t^{\prime})\mbox{$\displaystyle\frac{\,{\mathrm{d}}{\boldsymbol{e}}}{\,{\mathrm{d}}t^{\prime}}$}(t^{\prime})\,{\mathrm{d}}t^{\prime},$$

$$G(t-t^{\prime })=G_{\infty }+\sum _{p=1}^P{G}_p{\exp }^{-\frac{t-t^{\prime }}{{\theta }_p}}=G_0-\sum _{p=1}^P{G}_p(1-{\exp }^{-\frac{t-t^{\prime }}{{\theta }_p}}),$$

$${\boldsymbol{\sigma}}(t)=K{\boldsymbol{i}}\varepsilon_{\mathrm{V}}(t)+2\int\limits_{0}^{t}G(t-t^{\prime}){\boldsymbol{Q}}\mbox{$\displaystyle\frac{\,{\mathrm{d}}{\boldsymbol{e}}}{\,{\mathrm{d}}t^{\prime}}$}(t^{\prime})\,{\mathrm{d}}t^{\prime},$$

$${\boldsymbol{\sigma}}(t)={\boldsymbol{D}}_{\nu}\int\limits_{0}^{t}G(t-t^{\prime})\mbox{$\displaystyle\frac{\,{\mathrm{d}}{\boldsymbol{\varepsilon}}}{\,{\mathrm{d}}t^{\prime}}$}(t^{\prime})\,{\mathrm{d}}t^{\prime},$$

$$\displaystyle{\boldsymbol{D}}_{\nu}=\frac{1}{1-2\nu}\left[\begin{array}[]{cccccc}{2(1-\nu)}&{2\nu}&{2\nu}&0&0&0\\ {2\nu}&{2(1-\nu)}&{2\nu}&0&0&0\\ {2\nu}&{2\nu}&{2(1-\nu)}&0&0&0\\ 0&0&0&1&0&0\\ 0&0&0&0&1&0\\ 0&0&0&0&0&1\end{array}\right].$$

$$\mathbf{e}(t_{n+1})=\mathbf{e}(t_n)+\Delta \mathbf{e},$$

$$\mathbf{s}(t_{n+1})=\mathbf{s}(t_n)+\Delta \mathbf{s},$$

$$\Delta \mathbf{s}=2\hat{G}\mathbf{Q}\Delta \mathbf{e}+\Delta \hat{\mathbf{s}},$$

$$\hat{G}=G_{\infty }+\sum _{p=1}^P{G}_p\frac{{\theta }_p}{\Delta t}(1-{\exp }^{-\frac{\Delta t}{{\theta }_p}}),$$

$$\Delta {\sigma }_{\mathrm{m}}=K\Delta {\epsilon }_{\mathrm{V}}.$$

$$\Delta {\epsilon }_x$$

$$\Delta {\gamma }_{xz}$$

$$\Delta {\sigma }_x=\hat{E}\Delta {\epsilon }_x+(1+\hat{\nu })\Delta {\hat{s}}_x,$$

Full text

††journal: arXiv

1 Introduction

Laminated glass units consist of multiple glass sheets bonded together with polymer foils. The main reason to combine brittle glass sheets with soft polymer layers is to enhance the post-fracture behavior, so that a unit retains the load-bearing capacity after fracture of one or multiple glass plies [1]. Of course, the interlayer foil also influences the overall response of the unfractured elements, which is in the focus of this study, by mediating the shear interaction among the glass layers.

The glass layers can be accurately modeled as linear elastic, but the polymer foil exhibits a strongly time/temperature-dependent response, e.g. [2, 3, 4, and references therein]. Several types of interlayers have found their use in practice. From the first application of laminated glass members with polyvinyl butyral foil (PVB) in car windscreens, laminated glass has expanded into the building constructions, and new types of interlayer foils became available, such e.g. ionoplast polymer foil (IP) that is stiffer and less sensitive to viscoelastic phenomena than PVB at common temperatures.

As discussed in our previous contribution [5], the behavior of polymeric interlayer is often approximated as linear elastic, assuming a known temperature profile and duration of loading. The following methods have been developed for elastic laminated glass beams under this simplification:

The basic approach to approximate the mechanical response of laminated glass units involves two limiting cases: the layered case representing an assembly of independent glass layers and the monolithic case with thickness equal to the combined thickness of glass layers and interlayers, under the assumption of geometric linearity [6].

• 2.

Single-layered approaches, such as effective thickness methods that approximate laminated units by monolithic, geometrically linear, systems that achieve the equivalent maximum deflection or stress [7, 8, 9]; see also [10, 11] for an overview.

• 3.

Multi-layered approaches, explicitly resolving all the layers. These include analytical formulations, which are solved in a closed form for specific boundary conditions [12, 13, 14] and geometrically linear behavior, and numerical formulations based on the finite-difference [15] or finite-element methods [16, 5, 17] that can treat general boundary conditions and geometrically non-linear effects.

Solely the multi-layered approaches have been extended to account for viscoelastic effects, because they offer sufficient accuracy in stresses and strain distributions at the ply level. An analytical model of geometrically linear, simply supported, beams with viscoelastic interlayer represented by a Maxwell chain was developed by Galuppi and Royer-Carfagni [18, 19], who approximated the viscoelastic solution by the Fourier series in the spatial coordinate and solved the ensuing integral equations in the closed form to account for the time dependent effects. The same approach was later successfully used to study the stability of laminated columns [20] or time-dependent response of cold-bent laminated beams [21]. The most recent analytical study has been performed by Wu et al. [22], who develop a closed-form solutions for two-dimensional stress and strain fields in two-layer simply-supported beam with interlayer response described by the standard linear model. The only detailed numerical studies we are aware of include geometrically non-linear three-dimensional finite element simulations of viscoelastic plates under transverse loading by Duser et al. [23] and Bennison et al. [24], or laminated beams under lateral-torsional buckling by Bedon et al. [25].

The goal of this paper is to propose a viscoelastic geometrically non-linear finite element formulation for laminated glass beams based on layerwise refined laminate theories initiated by Mau [26]. The basis of our formulation is the elastic finite strain model developed earlier by Zemanová et al. [5], where we demonstrated that it offers accuracy comparable to detailed two-dimensional simulations at much lower computational cost. Here, this formulation is extended to account for the time/temperature-dependent behavior of the interlayer, which is incorporated into the numerical model by the incremental exponential algorithm by Zienkiewicz et al. [27]. In particular, we will compare the formulations arising from the two modeling aspects:

beam kinematics, involving the finite strain shear-deformable Reissner beam formulation and a simpler variant based on the Timoshenko beam theory complemented with the von Kármán assumptions of large deflections,

• 2.

constitutive assumptions for viscoelastic models, where we will consider models based on time/temperature-dependent shear modulus, complemented with either constant (time-independent) bulk modulus or constant Poisson ratio.

The differences in the mechanical results of the corresponding models will be studied in detail and verified against detailed finite element simulations, with the goal to determine the simplest, but still sufficiently accurate, formulation that is suitable for the generalization of elastic models of laminated glass plates [28] to the viscoelastic regime.

The rest of the paper is organized as follows. In Section 2, we introduce the different variants of laminated glass beam models, starting from kinematics assumptions in Section 2.1, proceeding to interlayer constitutive relations in Section 2.2, up to the finite element discretization in Section 2.3. Section 3 is devoted to the comparison of different formulations, along with an additional verification, validation, and parametric studies. The most important findings are summarized in Section 4. Finally, in Appendix A we gather technical details regarding the finite element technology, in order to make the paper self-contained.

2 Finite element models

The models introduced in this work share the following assumptions:

they assume planar cross sections of individual layers but not of the whole laminated glass unit, because the stiffness of glass layer and the effective stiffness of the interlayer differs by at least three orders of magnitude, e.g. [15],

• 2.

being based on Mau’s refined plate theory [26], each layer is treated independently and the inter-layer compatibility is enforced by the Lagrange multipliers. The main advantages of this approach is that the different constitutive models for each layer can be easily combined together and delamination phenomena can be efficiently accounted for, if needed, e.g., [29]111Note, however, that the perfect inter-layer adhesion is assumed in this paper, because the application of high pressures and temperatures during the production process result in high gluing forces of chemical nature [30]..

As noted above, the formulations differ in the beam kinematics (Finite-Strain (FS) Reissner vs. Large-Deflection von Kármán (VK) models) and in the assumption of viscoelastic constitutive models (constant bulk modulus $K$ vs. the Poisson ratio $\nu$). These considerations lead to four models, which are briefly introduced in this section.

Our notation is as follows. Scalar quantities are denoted by lightface letters, e.g., $a$, and matrices are denoted in bold, e.g., ${\boldsymbol{a}}$ or ${\boldsymbol{A}}$. In addition, ${\boldsymbol{A}}{\sf{}^{T}}$ stands for the matrix transpose, ${\boldsymbol{A}}^{-1}$ for the matrix inverse, and $a^{(i)}$ denotes that the quantity $a$ is related to the $i$-th layer. Also note that in order to avoid a profusion of notation, in Section 2.1 we omit the dependence of quantities of interest on the time variable $t$, while in Section 2.2 we omit the layer index $\bullet^{(i)}$, because the constitutive description holds only for the interlayer.

2.1 Kinematics

2.1.1 Reissner finite strain model

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

[TABLE]

where $u_{0}^{(i)}$ and $w_{0}^{(i)}$ are centerline displacements, $\varphi_{y}^{(i)}$ is the cross-section rotation, $x$ is the coordinate measured along the centerline, $z^{(i)}$ is the coordinate measured along the cross-section, and $i=1,2,3$ refers to individual layers.

The inter-layer compatibility is ensured via the geometric continuity conditions at the interfaces between the layers (with $i=1,2$)

[TABLE]

where the perfect horizontal and vertical adhesion is supposed. From (1), the compatibility condition become

[TABLE]

The non-zero normal, $\varepsilon_{x}^{(i)}$, and shear, $\gamma_{xz}^{(i)}$, strain components follow from

[TABLE]

where $\varepsilon_{x,0}^{(i)}$ and $\kappa_{y}^{(i)}$ denote the generalized axial strain and pseudo-curvature of the layer reference axis, introduced by Reissner [31] in the form

[TABLE]

2.1.2 Von Kármán model

A simpler model is obtained when adopting the von-Kármán assumptions of moderate deflections $w_{0}^{(i)}$ and small rotations $\varphi_{y}^{(i)}$, e.g. [32, Section 4.3]. Under these assumptions, the displacements of the $i$-th layer (1) simplify to

[TABLE]

and the inter-layer compatibility conditions (1) become linear in $\varphi_{y}^{(i)}$,

[TABLE]

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

[TABLE]

where the axial strains $\varepsilon_{x,0}^{(i)}$ must be measured by the Green-Lagrange strain tensor because of large deflections $w_{0}^{(i)}$,

[TABLE]

while the layer pseudo-curvature $\kappa_{y}^{(i)}$ still follows from (4b).

2.2 Constitutive relations

As indicated earlier, accurate description of time/temperature-dependent response of polymeric interlayer is decisive for the development of predictive models of laminated glass structures. In Section 2.2.1 we introduce such a model in the framework of linear viscoelasticity with accelerated/retarded time to account for the temperature dependence. We will also clarify the impact of the constitutive assumptions on the structure of stress-strain relations. The integral form of constitutive equations used in Section 2.2.1 is not very convenient for the numerical implementation, because all history variables during the whole loading process must be stored in order to evaluate the integral. Therefore, in the following sections, we will develop an incremental approach based on the well-established exponential algorithm by Zienkiewicz et al. [27]. The case of the constant bulk modulus is treated in detail in Section 2.2.2, leading to a rather involved procedure because of assumptions of beam theories. Section 2.2.3 then explains how the procedure simplifies under the assumption of time-independent Poisson ratio.

2.2.1 Time/temperature-dependent behavior of polymer foil

The main engineering property relevant to the composite behavior of the units is the shear-stress versus shear-strain characteristics of the soft interlayer. The shear modulus of the polymer interlayer $G$ is experimentally determined as a function of duration of loading and temperature, see Figure 2 for data for two PVB foils.

The temperature dependence is taken into account by the time-temperature superposition principle, in which the true time $t_{\text{true}}$ is replaced by the adjusted value $t=t_{\text{true}}/a_{T}$ modified by the temperature-dependent shift factor $a_{T}$. We employ the Williams-Landel-Ferry (WLF) equation [34] for this purpose:

[TABLE]

in which $C_{1}$ and $C_{2}$ are material constants, and $T$ and $T_{0}$ are the current and reference temperatures, respectively.

For further discussion, it will be convenient to decompose the vectors of strain ${\boldsymbol{\varepsilon}}$ and stress ${\boldsymbol{\sigma}}$,

[TABLE]

into the volumetric and deviatoric parts

[TABLE]

where $t$ denotes the time instant, ${\boldsymbol{i}}=\begin{bmatrix}1&1&1&0&0&0\end{bmatrix}{\sf{}^{T}}$. The volumetric parts follow from $\varepsilon_{\mathrm{V}}={\boldsymbol{i}}{\sf{}^{T}}{\boldsymbol{\varepsilon}}$ and $\sigma_{\mathrm{m}}=\frac{1}{3}{\boldsymbol{i}}{\sf{}^{T}}{\boldsymbol{\sigma}}$, and the deviatoric parts have the components

[TABLE]

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

[TABLE]

where {\boldsymbol{Q}}=\text{diag}\left[\begin{array}[]{cccccc}1&1&1&0.5&0.5&0.5\end{array}\right], and $K(t-t^{\prime})$ and $G(t-t^{\prime})$ denotes the bulk and shear relaxation moduli that completely characterize the interlayer response.

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

[TABLE]

where $G_{\infty}$ is the long-term shear modulus, $P$ stands for the number of viscoelastic units, $G_{p}$ denotes the shear modulus of the ${p}$-th unit, $\theta_{p}=\eta_{p}/G_{p}$ is its relaxation time related to the viscosity $\eta_{p}$, and $G_{0}=G_{\infty}+\sum_{p=1}^{P}G_{p}$ is the elastic shear modulus of the chain.

In order to close the constitutive description for the interlayer, two assumptions on the bulk relaxation functions will be considered. First, the constant bulk modulus approximation, e.g. [23, 24], under which the material behaves as a linear elastic solid for volumetric loading. The constitutive equation

[TABLE]

thus additionally involves the constant elastic bulk modulus. Second, the Poisson ratio $\nu$ is taken as constant, e.g. [7, 22], which simplifies the constitutive relation (11) into, e.g., [36, Eq. (29.63)]:

[TABLE]

with the constant matrix ${\boldsymbol{D}}_{\nu}$ provided by

[TABLE]

2.2.2 Incremental formulation for constant bulk

modulus

Material point

We start by dividing the time interval of interest $[0;t_{\max}]$ into non-equidistant time instants $0=t_{0}<t_{1}<t_{2}<\cdots<t_{N-1}<t_{N}=t_{\max}$ (recall that these correspond to temperature-adjusted values according to (8)). The strain history over the time interval $[t_{n};t_{n+1}]$ is assumed to be known, so that

[TABLE]

where $\Delta{\boldsymbol{e}}$ denotes the increment of deviatoric strain between the time instants $t_{n}$ and $t_{n+1}$. An analogous relation holds for the deviatoric stresses,

[TABLE]

so that, assuming that ${\boldsymbol{s}}(t_{n})$ is known, the goal is to determine the increment of the deviatoric stress $\Delta{\boldsymbol{s}}$.

For a linear strain variation over the time interval $[t_{n};t_{n+1}]$, this increment can be expressed as [27]

[TABLE]

where $\widehat{G}$ stands for the effective shear modulus over the time interval and $\Delta\widehat{{\boldsymbol{s}}}$ for the stress relaxation under the zero strain increment. For the Maxwell chain model, Figure 3, these quantities follow from

[TABLE]

where $\Delta t=t_{n+1}-t_{n}$ is the temperature-adjusted time step (we assume for simplicity that the temperature remains constant over the interval), and ${\boldsymbol{s}}_{p}$ and $\Delta\widehat{{\boldsymbol{s}}}_{p}$ stand for the values of the deviatoric stresses and its relaxation in the $p$-th unit.

Because the volumetric response is assumed to be elastic, the increment of the mean strain, $\Delta\sigma_{\mathrm{m}}$, can be obtained from the given increment of the volumetric strain $\Delta\varepsilon_{\mathrm{V}}$ directly as

[TABLE]

Beams

As follows from the standard assumptions of the conventional beam theory, e.g. [37, Section 29], the only non-zero stress components are the increments of normal and shear stresses $\Delta\sigma_{x}$ and $\Delta\tau_{xz}$; the corresponding volumetric and deviatoric stress components entering (24) and (22) thus become $\Delta\sigma_{\mathrm{m}}=1/3\Delta\sigma_{x}$, $\Delta s_{x}=2/3\Delta\sigma_{x}$, and $\Delta s_{xz}=\Delta\tau_{xz}$.

In order to obtain the corresponding strain increments $\Delta\varepsilon_{x}$ and $\Delta\gamma_{xz}$ needed in the beam formulations, we employ the volumetric-deviatoric split (9) and invert the relations (22) and (24) to obtain

[TABLE]

from which we directly obtain the incremental stress-strain relations in the form

[TABLE]

The symbols $\widehat{E}$ and $\widehat{\nu}$ in Eqs. (25) and (26) stand for the effective Young modulus and Poisson ratio, obtained as $\widehat{E}=(9K\widehat{G})/(\widehat{G}+3K)$ and $\widehat{\nu}=(3K-2\widehat{G})/(2(3K+\widehat{G}))$, e.g. [36, Table D.1]. It will also prove useful to relate the deviatoric normal strain to the full normal strain via

[TABLE]

Cross-section

In order to obtain the constitutive equations for a cross-section located at $x$, recall from (4) or (6) that the normal strain $\varepsilon_{x}$ is an affine function in the through-the-thickness coordinate $z$, whereas the shear strain $\gamma_{xz}$ is independent of $z$. Upon integrating the stress increments (26) over the cross-section, we obtain

[TABLE]

where $\Delta N_{x}$, $\Delta V_{z}$, and $\Delta M_{y}$ stand for the increments of the normal and shear forces, and the bending moment (the corresponding relaxation forces are distinguished by a hat); $\Delta\varepsilon_{0,x}$, $\Delta\gamma_{xz}$, and $\Delta\kappa_{y}$ are the corresponding increments of the centerline strain, the transverse shear, and the pseudo-curvature determined from (4) or (6); and $A$, $I_{y}$, and $A_{\mathrm{s}}$ denote the cross-section area, the second moment of area, and the shear-effective area.222For three-layered units (glass/foil/glass) we set the value of the shear-effective area to $A_{\mathrm{s}}^{(i)}=5/6A^{(i)}$ for glass layers with $i=1,3$ and $A_{\mathrm{s}}^{(2)}=A^{(2)}$ for the foil. For these values, which correspond to a quadratic shear stress variation in glass and to constant shear in the interlayer, e.g. [38], we achieved the best agreement with the response of the detailed two-dimensional model from the ADINA solver (not shown). Note that in (28), we utilized the fact that the cross-section and the effective material constants $\widehat{E}$, $\widehat{G}$, and $\widehat{\nu}$ are identical for the whole layer.

As follows from Eq. (28), the increments of internal forces involve the relaxation contributions $\widehat{N}_{x}$, $\widehat{M}_{y}$, and $\widehat{V}_{z}$, which in turn depend on the values of the deviatoric stresses in elements of the Maxwell chain (22, 23) across the cross-section. Assuming that ${\boldsymbol{s}}_{p}(t_{0})={\boldsymbol{0}}$, a closer inspection of these relations together with Eq. (27) reveals that the normal and shear deviatoric stress components exhibit the affine-constant distribution across the cross-section at $x$ at time $t_{n}$:

[TABLE]

therefore, the three quantities $s_{x,p}^{N}$, $s_{x,p}^{M}$, and $\tau^{V}_{xz,p}$ completely determine the stress distribution in the $p$-th elements of the Maxwell chain over the cross-section. They also allow for a compact expression of the incremental relaxation forces in the form,

[TABLE]

where, e.g.,

[TABLE]

recall (23). The remaining quantities $\Delta\widehat{s}^{M}_{x}$ and $\Delta\widehat{\tau}^{V}_{xz}$ follow by analogy.

At the end of the time interval, the stresses in the $p$-th unit must be updated on the basis of increments of the generalized strains $\Delta\varepsilon_{0,x}$, $\Delta\kappa_{y}$, and $\Delta\gamma_{xz}$ for a cross-section at $x$. Combining (27) and (22) with (3) or (6), one obtains

[TABLE]

2.2.3 Incremental formulation for constant Poisson ratio

Material point

Assuming the same time discretization as in Section 2.2.2, the strain and stress increments over the interval $[t_{n};t_{n+1}]$ are expressed as

[TABLE]

which involve the full strain and stress tensors. Comparing (11) with (14), we see that the incremental formula (22) from the previous section becomes

[TABLE]

where matrix ${\boldsymbol{D}}_{\nu}$ depends on the constant Poisson ratio $\nu$ via (21) and $\widehat{G}$ is defined with (23).

Beams

Increments of the non-zero components of stress vector, $\Delta\sigma_{x}$ and $\Delta\tau_{xz}$, follow directly from (32) in the form

[TABLE]

with the effective Young modulus now given by $\widehat{E}=2(1+\nu)\widehat{G}$, cf. (26)

Cross-section

Proceeding in the same way as in the previous section, the incremental constitutive relations at the cross-section level, Eq. (28), attain a similar form:

[TABLE]

Again, the increments of relaxation forces $\Delta\widehat{N}_{x}$, $\Delta\widehat{V}_{z}$, and $\Delta\widehat{M}_{y}$ depend on the distribution of stresses in the Maxwell units over the cross-section. By analogy to (29), the normal and shear stresses in the $p$-th unit admit the decomposition:

[TABLE]

from which we obtain

[TABLE]

with, e.g.,

[TABLE]

Finally, we need to relate the increments of the history variables $\Delta\sigma_{x,p}^{N}$, $\Delta\sigma_{x,p}^{M}$, and $\Delta\tau^{V}_{xz,p}$ to the increments of the generalized strains $\Delta\varepsilon_{x,0}$, $\Delta\kappa_{y}$, and $\Delta\gamma_{xz}$. The relations for the normal components follow directly from (32) as

[TABLE]

The shear stress increment $\Delta\tau^{V}_{xz,p}$ is obtained from (31c).

2.3 Finite element model

In this section, we extend the geometrically nonlinear Reissner elastic solver from [5] to the incremental viscoelastic formulation and discuss the simplifications due to the von-Kármán assumptions.

2.3.1 Variational setting

Similarly to [5], we derive the model using variational arguments. Because the glass layers are assumed to be elastic, their response is governed by a quadratic stored energy density introduced in [5, Eq. (13)]:

[TABLE]

with $i={1,3}$ and with the column matrix ${\boldsymbol{\varepsilon}}^{(i)}=\begin{bmatrix}\varepsilon_{x,0}^{(i)}&\kappa_{y}^{(i)}&\gamma_{xz}^{(i)}\end{bmatrix}{\sf{}^{T}}$ collecting the generalized strain components. For the interlayer, $i=2$, the energy density is history dependent to account for the incremental form of the cross-section constitutive relations (28) or (33). For the time interval $[t_{n};t_{n+1}]$, it reads

[TABLE]

with $W^{(2)}_{0}=0$ and

[TABLE]

Here, e.g., $\Delta\varepsilon_{x,0}^{(2)}=\varepsilon_{x,0}^{(2)}-\varepsilon_{x,0}^{(2)}(t_{n})$, and the remaining terms follow by analogy.

The energy functional at the structural level collects the contributions of layers

[TABLE]

where ${\boldsymbol{u}}=\begin{bmatrix}u_{0}^{(i)},w_{0}^{(i)},\varphi_{y}^{(i)}\end{bmatrix}_{i=1}^{3}$ stands for the centerline displacements and cross-section rotations and ${\boldsymbol{\varepsilon}}(\bullet)$ abbreviates the strain-displacement relations for the Reissner or the von-Kármán kinematics introduced in Sections 2.1.1 and 2.1.2, respectively. Complementing the internal energy functional with the energy of external forces provides

[TABLE]

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

[TABLE]

The generalized displacements at time $t_{n+1}$, ${\boldsymbol{u}}(t_{n+1})$, then follow by the minimization of the total energy functional (36) under the equality constraints (37).

2.3.2 Discretization

Discretization of the model is accomplished using the standard finite element procedures. To that purpose, each layer is discretized into the identical number of $n^{\mathrm{{}_{e}}}$ elements $\Omega^{(i)}_{e}$ and the generalized displacements are approximated as

[TABLE]

where ${\boldsymbol{N}}_{e}$ stands for the $3\times 6$ matrix of linear basis functions and ${\boldsymbol{d}}_{e}^{(i)}$ is the $6\times 1$ column matrix of element degrees of freedom; see Appendix A for details. After inserting the approximation (38) into (36), the internal energy becomes the function of the nodal degrees of freedom

[TABLE]

where ${\boldsymbol{d}}=\begin{bmatrix}{\boldsymbol{d}}_{j}^{(i)}\end{bmatrix}_{i=1,j=1}^{3,n^{\mathrm{{}_{e}}}+1}$ collects the layer nodal displacements and, e.g.,

[TABLE]

In analogy to the previous section, the column matrix of nodal displacements ${\boldsymbol{d}}(t_{n+1})$ follows from minimization of the total energy function

[TABLE]

in which ${\boldsymbol{f}}^{(i)}_{\mathrm{ext},e}$ denotes the generalized nodal forces due to external loading of the $i$-th layer. The kinematically admissible displacements must satisfy the constraint

[TABLE]

with the meaning of the inter-layer continuity conditions (37) enforced at the nodes.

2.3.3 Solution procedure

Following the standard procedures of equality-constrained optimization, e.g., [39, Chapter 14], we consider the Lagrange function

[TABLE]

where the column matrix of Lagrange multipliers ${\boldsymbol{\lambda}}$ has the meaning of interface nodal forces that ensure the inter-layer compatibility, (39). The nodal displacements ${\boldsymbol{d}}(t_{n+1})$ and the corresponding Lagrange multipliers ${\boldsymbol{\lambda}}(t_{n+1})$ then arise as the saddle point of the Lagrange function.

The optimality conditions lead to a system of equations non-linear in ${\boldsymbol{d}}(t_{n+1})$ which is solved iteratively using the Newton method. To this purpose, we express the $(k+1)$-th approximation to ${\boldsymbol{d}}(t_{n+1})$ as

[TABLE]

where ${{}^{k}}{\boldsymbol{d}}$ denotes a known previous iterate; for $k=0$, we set ${{}^{0}}{\boldsymbol{d}}={\boldsymbol{d}}(t_{n})$. The displacement increment ${{}^{k+1}}\delta{\boldsymbol{d}}$ and the vector of Lagrange multipliers ${{}^{k+1}}{\boldsymbol{\lambda}}$ are determined from the linearized system

[TABLE]

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

[TABLE]

Because the unknown displacements are approximated independently at each layer, recall Eq. (38), the internal forces ${\boldsymbol{f}}_{\mathrm{int}}$ and tangent stiffness ${\boldsymbol{K}}_{\mathrm{t}}$ matrices possess the block-diagonal structure:

[TABLE]

The sub-matrices for the interlayer must be determined using the effective values $\widehat{G}$ and $\widehat{E}$ (as emphasized by the hat), and additionally contain the history-dependent contributions $\Delta\widehat{{\boldsymbol{K}}}_{\mathrm{t}}$ and $\Delta\widehat{{\boldsymbol{f}}}_{\mathrm{int}}$ that result from the incremental energy split (2.3.1). As usual, all sub-matrices in (40) are obtained by the assembly of element contributions, e.g. [32, Section 2.3.6]. For the readers’ convenience, these are discussed in Appendix A and for the Reissner model particularly in Appendix A.2.

To the same purpose, in Algorithm 1 we present a pseudo-code of one-step geometrically non-linear viscoelastic problem that summarizes the developments of the current section. Notice that the algorithm is terminated by two residuals, cf. [5, Eq. (34)],

[TABLE]

related to the nodal equilibrium and continuity conditions, respectively.

A similar procedure applies also to the von Kármán model with a few simplifications. Indeed, because the continuity constraint (39) is linear in ${\boldsymbol{d}}$, recall (5), we have ${{}^{k}}{\boldsymbol{c}}={\boldsymbol{0}}$, matrix ${\boldsymbol{C}}$ remains constant during iterations, and matrix ${\boldsymbol{K}}_{\lambda}={\boldsymbol{0}}$. The remaining terms needed in (40) are provided in Appendix A.1.

3 Comparison, verification, and validation

The developments described in the previous section resulted in four geometrically non-linear finite element models summarized in Table 1; recall that they emerge as the combinations of Reissner (FS) or von Kármán (VK) kinematics and the assumptions of constant bulk modulus (K) or Poisson ratio ($\nu$) of the interlayer. We also include the results of geometrically linear formulation (LIN) [40], in order to highlight the importance of geometric non-linearity where relevant.

In addition, the following data are considered in the studies reported in this section:

detailed two-dimensional viscoelastic finite element simulations (2D), performed in ADINA 9.0.2 (ADINA R&D, Inc.) finite element system,

• 2.

experimental data, obtained by López-Aenlle et al. [10] for three simply-supported beams and one two-span continuous beam with different layer thicknesses;

• 3.

results of our earlier finite strain model with elastic interlayer [5], under the assumption of the constant bulk modulus $K^{(2)}$; and

• 4.

responses of the two limiting cases, corresponding to the monolithic model, which assumes the full layer interaction, and the layered model with no inter-layer interaction, both in the small strain regime.

The remainder of this section is organized as follows. In Section 3.1, we provide dimensions and boundary conditions for beams used in the analysis and specify material data, loading procedures, and details on finite element models. Section 3.2 is dedicated to the comparison of the four finite element models from Table 1. Based on these results, the model combining the von Kármán kinematics with the time-independent Poisson ratio is verified against a two-dimensional finite element model and partially validated against experiments in Section 3.3. We proceed with a parametric study concerning the effects of geometric nonlinearity and temperature, Section 3.4 and conclude our investigations by highlighting the difference in full viscoelastic and simplified elastic analyses, Section 3.5.

3.1 Overview of examples

3.1.1 Geometry and finite element model

The most common laminated beams with three layers (glass/PVB/glass) are considered, either simply supported or clamped. Note that in practice, the laminated glass elements are not perfectly fixed or simply supported, but these two cases represent limits of the real support conditions concerning the effects of geometric non-linearity, e.g. [41, 7]. The layer thicknesses and spans were selected to ensure that the maximum deflection does not exceed $1/200$ of the span, which we estimated as a suitable serviceability limit for load-bearing glass structures.

Similarly to [5], each layer was discretized with $500$ elements and linear stress smoothing was adopted, in order to reach the four-digit accuracy of displacements and stress values reported below. Both tolerances in Algorithm 1 were set $\epsilon_{1}=\epsilon_{2}=10^{-5}$.

3.1.2 Materials

The elastic constants of glass, the Young modulus $E^{(1)}=E^{(2)}=72$ GPa and the Poisson ratio $\nu^{(1)}=\nu^{(3)}=0.23$, correspond to the values determined by López-Aenlle et al. [10] from a static bending test. The experimental viscoelastic characterization of PVB was performed by Pelayo et al. [33] on specimens of thickness $0.38$ mm by tensile relaxation tests with the duration of 10 min at different temperatures ranging from $-15$ to $50^{\circ}$ C. The resulting parameters of the generalized Maxwell chain are summarized in Table 3, the parameters of the Williams-Landel-Ferry (WLF) equation (8) were set to $C_{1}=12.6$, $C_{2}=74.46$, and $T_{0}=20\,^{\circ}$C. The corresponding shear relaxation functions $G^{(2)}$ are plotted in Figure 2(b). According to the adopted constitutive assumption, we further set the bulk modulus $K^{(2)}=2$ GPa, e.g. [23, 24], or the Poisson ratio $\nu^{(2)}=0.49$, e.g. [1, Section 1.1.3].

3.1.3 Loading

We considered loading scenarios with spatially uniform distributed loads of intensities shown in Figure 4. Scenario A, Figure 4 (a), approximates the instantaneous loading reported in [33] by a loading branch, in which the intensity linearly increases from [math] to $10$ Nm*-1* in $10^{-5}$ s, followed by the constant load of intensity $10$ Nm*-1* up to $10^{5}$ s. The time steps in the incremental algorithms were distributed uniformly in the logarithmic scale, in particular $6$ time steps were used to discretize the interval $\langle 10^{-6};10^{-5}\rangle$ s, whereas the interval $\langle 10^{-5};10^{5}\rangle$ s was discretized into $24$ time steps.

Histories B and C represent loading-unloading scenarios shown in Figure 4 (b), where the steps in intensity were approximated by linearly varying loading within $10^{-5}$ s. The overall loading program was discretized into $90$ time steps uniformly distributed in the logarithmic scale and refined in the vicinity of jumps.

3.2 Comparison

Comparison of the proposed models and of the effects of geometric nonlinearity are discussed in this section. Three examples are considered, namely fixed-end beam (geometry I) under one-load history (loading A), fixed-end beam (geometry I) under two-load history (loading B), and simply-supported beam (geometry II) under two-load history (loading C). The temperature of $25\,^{\circ}$C was kept constant in all examples, as the temperature effects are in the focus of Section 3.4.

In Figures 7–7, we plot the deflections and the maximum of stresses predicted by four geometrically nonlinear (FSK, FSν, VKK, VKν) and two geometrically linear (LINK, LINν) approaches. Their comparison reveals that:

The viscoelastic approaches based on the assumption of constant bulk modulus $K$ (FSK, VKK, and LINK), or of the Poisson ratio $\nu$ (FSν, VKν, resp. LINν) provide the same results for all examples because the interlayer mostly transfers shear loading (compare cf. Eq. (31) and Eq. (34)). The errors in deflections and stresses are much smaller than 0.1%.

• 2.

The geometrically nonlinear approaches based on the assumptions of large deflections (VKK, VKν), or finite strains (FSK, FSν) give again almost identical results for tested examples; the errors in deflections and stresses do not exceed $0.1\%$.

• 3.

For statically determinate example, the geometrically linear and nonlinear solver provide the same results, whereas for statically indeterminate examples, the error of linear approach can reach 100% up to 300%; cf. [5].

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

In conclusion, it follows from the results presented in this section that all four geometrically non-linear models provide results indistinguishable from the practical design viewpoint. Therefore, we will report only the results of VKν model in what follows, because its formulation is the most transparent from the four options.

3.3 Verification and validation

Validation of the numerical model was performed using data from experimental campaign by López-Aenlle et al. [10], utilizing three types of simply-supported one- or two-span laminated glass beams (exact dimensions are specified in the caption of Figure 8). The beams were subjected to seven concentrated loads that approximated uniformly distributed loading and the evolution of mid-span deflections and of strains was measured with laser sensors and tensometers, respectively. In addition, the average temperature during the experiment was recorded. In order to generate verification data in ADINA system, we discretized each beam with a very fine mesh of quadratic finite strain elements of approximately square shape with edge size of $0.38$ mm, and determined its response under a ramp loading (analogous to the one shown in Figure 4(a)) while assuming the constant value of the bulk modulus $K^{(2)}=2$ GPa.

Results of the simulations and experiments are summarized in Figure 8 and involve the values of deflections, Figure 8(a), and stresses, Figure 8(b), as functions of loading time. Notice that the experimental results for stresses appear only in Figure 8(3b), in order to demonstrate that they are oscillatory and even do not fit in the layered-monolithic bounds. This is probably because the temperature was not kept constant during the experiment and the polymer interlayer shows very high sensitivity to temperature variations, as can be deduced from data in Table 3. Therefore, the experimental values of stresses are excluded from consideration.

We conclude from the collected graphs that our computational model is successfully verified and only partially validated. Indeed, the results of the detailed two-dimensional finite element model (2D) are indistinguishable from the results of VKν model, both for displacements and stresses. The best match between the measured displacements (EXP) is achieved for non-symmetric 4/0.38/8 beam, Figure 8(1a); the other experiments show visible deviations, which we attribute again to only partial temperature control.

In order to make these claims more quantitative, in Tables 6 and 6 we provide detailed results of all tested beams at 10 hours. Note that, for example, the VKν-EXP error is understood as $(w_{\mathrm{VK}_{\nu}}-w_{\mathrm{EXP}})/w_{\mathrm{EXP}}$. The results confirm an excellent match between detailed FE model and VKν model, the errors do not exceed $0.3\%$ in deflections and $0.2\%$ in stresses. The difference between experimental data and VKν formulation is somewhat less satisfactory, and ranges from about $5$–$20$%. However, even this lower accuracy is superior to the accuracy of monolithic (LIMMON) and layered (LIMLAY) limits: the monolithic bound predicts values that can be up to $50\%$ smaller than for the more refined models, whereas the overestimation by the layered limit can reach up to $200\%$.

3.4 Geometric nonlinearity and temperature

Herein, we complement the verification and validation results with a study on a structure exhibiting significant effects of geometric nonlinearity. Attention is also paid to the effect of temperature variation, accounted for by the Williams–Landel–Ferry formula (8). In particular, we consider the fixed-end beam structure from Table 2 subjected to loading A, recall Figure 4(a), under three constant temperatures $0\,^{\circ}$C, $25\,^{\circ}$C, and $50\,^{\circ}$C.

The resulting time evolution of deflections and stresses is plotted in Figure 9 for geometrically linear and nonlinear models. The standard upper and lower limits, provided by the geometrically linear monolithic and layered approximations, are plotted in Figure 9 (a) and they indeed bound the behavior of the laminated unit. For geometrically nonlinear response, we also provide the elastic response (FS${}_{\text{EL}}$), determined by elastic model from [5] with the interlayer shear modulus determined from (12) for the given temperature at the final simulation time of $10^{5}$ s. This is complemented with selected data of mid-point deflections, Table 8 and stresses, Table 8.

The results in Figure 9(a) are in full agreement with the outcomes presented in [43], namely that the laminated units exposed to temperatures around $0\,^{\circ}$C effectively behave as the monolithic ones. In addition, even at $50^{\circ}$C, the shear stiffness is sufficient to ensure the inter-layer interaction, so that the overall response is much closer to the monolithic limit rather than to the layered limit. Behavior of geometrically non-linear models, Figure 9(b), exhibits similar trends, but the displacement and stress magnitudes are reduced to $40\%$ and $60\%$, respectively; see also Tables 8 and 8. In addition, unlike in the geometrically linear model, the extreme stresses at the mid-span decrease with increasing time. We consider this to be a direct consequence of the development of membrane stresses in glass layers, which leads to the redistribution of extreme stresses from the beam mid-section towards the supports (the same mechanism was observed in the detailed 2D finite element model).

Furthermore, the results in this section serve as an additional verification of VKν model, because the errors in displacements and stresses do not exceed $0.5\%$ for temperatures of [math] and $25^{\circ}$C with respect to the results of detailed 2D finite element model, consult again Tables 8 and 8. Note that simulations in ADINA did not convergence at $50^{\circ}$C, because of poor conditioning of the stiffness matrix caused by a very high contrast in the interlayer and glass properties. In addition, the results reveal that the elastic model provides very good estimates at the final time of simulations: the errors stay below $2\%$ for all quantities of interest. This aspect will be studied in more detail in the next section.

3.5 Comparison with simplified elastic solution

Being inspired by earlier studies by Galuppi and Royer-Carfagni [19], we close this section by comparing the full viscoelastic solution with a simplified ”secant” solution, in which we determine the structural response at the relevant times by the geometrically non-linear model [44] with shear modulus determined from (12). Similarly to [19], we consider the two-load history B, see Figure 4(b), consisting of a sudden loading, sudden unloading, and holding parts. In addition, both simply supported and fixed-end beams are considered under constant temperatures of $25$ and $50^{\circ}$C.

The simulation results collected in Figures 10 and 11 reveal that the correspondence between the simplified model and full viscoelastic solution is rather satisfactory for holding phases of the loading, but differences appear in the vicinity of displacement jumps, where the simplified solution can lead to unsafe predictions. This effects is the most pronounced for stresses in the fixed-end beam and it is attributed to a different stress redistribution in the viscoelastic and elastic beams as discussed in the previous section. In addition, the magnitude seems to be highly temperature-dependent: the differences are negligible at $0^{\circ}$C (not shown) and increase with increasing temperature.

More quantitatively, for the fixed-end beam, Figure 10, the differences remain below $5\%$; at $25^{\circ}$C they do not exceed $2\%$ for both stresses and strains, whereas for $50^{\circ}$C they reach $4\%$ for displacements and $5\%$ for stresses. For the simply supported beam, the differences between elastic and viscoelastic solution increase and reach up to $27\%$ for displacements and $14\%$ for stresses. These results again highlight the complex interplay between structural heterogeneity, geometric nonlinearity, and viscoelastic effects.

Note in conclusion that all these findings agree well with the analytical study by Galuppi and Royer-Carfagni [19] for a simply supported beam with an interlayer modeled using simplified one- and three-unit Maxwell chains, disregarding the temperature effects.

4 Conclusions

In this paper, we have introduced and thoroughly tested four layer-wise finite element formulations to determine displacements and stresses in laminated glass beams with interlayer(s) exhibiting time/temperature-dependent material response. The models differ in the adopted beam kinematics — large-deflection von Karmán formulation vs. finite strain Reissner formulation — and the assumptions adopted in the constitutive model — constant bulk modulus vs. constant Poisson ratio. By comparing their response with detailed two-dimensional finite element simulations, experimental data, and simplified approaches, we have reached the following conclusions:

All four considered models provide almost identical response under monotone and non-monotone loading for deflections up to $1/50$ of span. Therefore, the most straightforward formulation, combining the assumptions of the von Kármán kinematics and of constant Poisson ratio, can be used to simulate more complex structures. 2. 2.

This formulation has been successfully verified against several detailed finite element simulations performed in ADINA system and partially validated against experimental data collected by López-Aenlle et al. [10]. The model has also consistently reproduced typical features of mechanical response of laminated glass structures at different temperatures. 3. 3.

A simplified model, based on an elastic formulation with the interlayer shear modulus adjusted to duration of loading and temperature, provides satisfactory agreement with full viscoelastic solutions at low temperatures or at holding periods of loading. At room and elevated temperatures, differences have been observed for rapid changes of loading and unloading that may reach the order of $10\%$ and/or may be on the unsafe side.

As the next step, we will enhance our elastic finite element solver for laminated glass plates [28] with viscoelastic effects, closely following the formulations developed in this paper. Subsequently, we would like to extend the formulation to post-critical state, by incorporating glass fracture into the layer-wise formulation.

Acknowledgments

This work was initiated when Alena Zemanová was a post-doctoral researcher, supported by the project No. CZ.1.07/2.3.00/30.0034: Support of inter-sectoral mobility and quality enhancement of research teams at Czech Technical University in Prague (06/2014–06/2015). All authors acknowledge the support by the Czech Science Foundation, project No. 16-14770S.