Acta Polytechnica CTU Proceedings doi:10.14311/APP.2020.26.0086 Acta Polytechnica CTU Proceedings 26:86–93, 2020 © Czech Technical University in Prague, 2020 available online at https://ojs.cvut.cz/ojs/index.php/app EULER AND EXPONENTIAL ALGORITHM IN VISCOELASTIC ANALYSES OF LAMINATED GLASS Jaroslav Schmidt∗, Alena Zemanová Czech Technical University in Prague, Faculty of Civil Engineering, Department of Mechanics, Thákurova 7, 166 29 Prague 6, Czech Republic ∗ corresponding author: jaroslav.schmidt@fsv.cvut.cz Abstract. Laminated glass combines two remarkable materials: glass and a polymer ply. While glass is stiff and brittle, the polymer ply is a rate-dependent compliant material. Together, they form a material which keeps the aesthetic value of glass, and due to the polymer, no fragile collapse appears. The polymer ply exhibits time- and temperature-dependency, whereas glass suffers from brittle fracture, which makes the analysis difficult. In this article, a 2D sectional plane-stress model for the viscoelastic analysis of laminated glass is presented. This study presents one step in the development of a phase-field-based damage solver for laminated glass to select the optimal time-integration scheme for a quasistatic-analysis and later also for dynamics. The validation against experimental data is provided, and the model reduction is also discussed. Keywords: Laminated glass, viscoelasticity, generalised Maxwell model, exponential algorithm, Euler method. 1. Introduction Nowadays laminated glass is a popular composite multi-layer material which consists of several glass plates and polymer interlayers. The laminated pro- cess improves post-breakable behaviour compared to conventional glass but retains its aesthetic value [1]. This makes laminated glass a suitable material for transparent structures in architecture and other fields. Unfortunately, despite its advantages, laminated glass is still a very brittle material which makes the sim- ulations of the initiation and propagation of cracks difficult. This paper partially follows the project objectives specified in our research project focused on the design and advanced modelling of forced-entry and bullet resistant glass structures. One of the goals of the project is to simulate crack development in laminated glass using phase-field fracture method [2] to better understand and explain the pre- and post-fracture behaviour of laminated glass under impact loading. Several dynamic and quasi-static experiments were performed, and the next step is the validation of the collected experimental data against numerical models for laminated glass. At this moment, the numerical solver is imple- mented for laminated beams and plates made of sev- eral elastic layers under quasi-static loading. This study presents another step in the development of the phase-field-based damage solver with the focus on the time- and temperature-dependent behaviour of the polymer interlayer. At this stage, the analysis deals with a quasi-static regime only. In particular, two questions concerning interlayer modelling are dis- cussed: • What is the optimal time-integration scheme with respect to the computational time and implementa- tion demands? • How is the global response of the model influenced by involving only certain cells of the generalised Maxwell chain into the numerical analysis? Moreover, this study provides also validation of mate- rial parameters of two interlayer foils obtained in the previous step of the project [3]. The paper is organised as follows: Section 2 in- troduces the 2D plane-stress model and presents the constitutive relationships. Next, in Section 3, two time-integrators are recalled, and finally, Section 4 provides numerical simulations, validation, and para- metric study, concluded with the discussion of results in Section 5. 2. Laminated glass model In this section, the mechanical model for the numerical analysis is presented. The constitutive assumptions for glass and interlayer are as follows: • Glass has an almost purely elastic response. There- fore, it is modelled as an elastic isotropic material described by two constant: Young’s modulus E and Poisson’s ratio ν. • The interlayer polymer is strongly time- and temperature-dependent material effectively mod- elled by the generalised Maxwell model [4, Ap- pendix A], see Figure 1, where the material con- stants are the stiffnesses of springs Gp and the viscosities ηp and the long-term shear modulus G∞ corresponding to the single stiffness. The constitu- tive relationship is expressed in the shear variables 86 https://doi.org/10.14311/APP.2020.26.0086 https://ojs.cvut.cz/ojs/index.php/app vol. 26/2020 Viscoelastic analysis of laminated glass instead of in the normal components as the inter- layer material is almost incompressible, and shear in the interlayer is dominant in bending regime. Then, τ and γ represent shear stress and strain respectively. γ(t) γe,1(t) γv,1(t) τ(t) G1 G2 G3 GP G∞ η1 η2 η3 ηP Figure 1. Generalised Maxwell constitutive model. A practical choice for the mechanical model of lami- nated glass is a layer-wise model consisting of several independent beams or plates connected together at the interfaces. Despite its advantages, we validate experiments against the plane-stress 2D model, where the longitudinal cross-section is discretizted. One of the reasons for this approach is to check if the through- the-thickness compression of the interlayer is really negligible and can be neglected. 3. Time integrators The brief overview of the exponential and backward- Euler integrators is presented in this section. 3.1. Governing equations The following model assumptions are considered: • A special two-dimensional state of stress called plane stress is assumed to reduce the originally three-dimensional task and subsequently the com- putational cost. • Because all of the laminated glass samples were simply-supported, the analysis assumes only small deformations and deflections. The comparison with experiments shows that this assumption is sufficient as no stabilisation due to membrane forces occurred. • We assume that the Poisson ratio of the polymer interlayer is constant. For further comparison with a constant bulk modulus assumption see also [5]. As mentioned before, the polymer material point is described by the generalised Maxwell chain which consists of one single linear spring with stiffness G∞ and several Maxwell cells. Each cell is characterised by a spring with stiffness Gp and a damper with viscosity ηp. The cells are indexed by p ∈ P, where P = 〈0,P〉 ⊂ N. The time integrators are used only for the viscoelastic interlayer and glass is considered as rate-independent. Under these assumptions, the behaviour of the in- terlayer can be described by a weak form of the equi- librium equation ∫ Ωfoil δε :  G∞Dν : ε + ∑ p∈P σp   dΩ − δFext = 0, (1) complemented by ordinary differential equations (ODE) σ̇p Gp + σp ηp = Dν : ε̇, p ∈P. (2) In the former equations δFext stands for the virtual work done by the external loads on the virtual displace- ment δu, and the small strain tensor ε is obtained as the symmetric displacement gradient ε = ∇su. The dimensionless tensor Dν corresponds to the stiffness tensor of an isotropic material expressed by a unit shear modulus and the Poisson ratio of ν. The stress carried by p-th cell is denoted by σp. On the other hand, the governing equation for the deformation of glass does not contain viscous stress tensors σp, and the corresponding weak form becomes∫ Ωglass δε : GglassDν : ε dΩ − δFext = 0. (3) The implementation for laminated glass, which com- bines the elastic glass domain and the viscoelastic interlayer, is described in Section 4. Equations (1)–(2) contain both space and time gradients. The space discretization is captured by using conventional finite element methods (FEM), whereas the type of the time discretization is a part of this analysis. Two discretiza- tion options are recalled in the following subsection. 3.2. Numerical integrators The common feature for all presented time integra- tors is the discretization into a set of time instants T = {0, . . . , tmax}, where tmax is the largest time of interest. Backward Euler method For the backward Euler method, we approximate the time gradient by the first- order backward finite difference •̇(ti+1) ≈ •(ti+1) −•(ti) ∆ti . (4) To shorten the notation, the time increment ∆ti is defined as ti+1 − ti, and the explicit time dependency is omitted and replaced by a lower index, therefore •(ti) ≡ •i. If we use the above-mentioned gradient approximation of the ODE evaluated in time instant ti+1, we get σp,i+1 −σp,i Gp∆ti + 1 ηp σp,i+1 = 1∆ti Dν : (εi+1 −εi) . (5) Then, we can express explicitly σi+1 as σp,i+1 = 1 ζp,i (Dν : (εi+1 −εi)Gp + σp,i) , (6) 87 Jaroslav Schmidt, Alena Zemanová Acta Polytechnica CTU Proceedings where ζp,i = 1 + ∆ti/τp is an auxiliary dimensionless parameter. The stress tensor can be directly substi- tuted into the equilibrium equation to get the weak form for εi+1 and the resulting equation ∫ Ωfoil δε :  G∞ + ∑ p∈P Gp ζp,i   Dν : εi+1 dΩ − ∫ Ωfoil δε : ∑ p∈P Gp ζp,i Dν : εi dΩ (7) + ∫ Ωfoil δε : ∑ p∈P σp,i ζp,i dΩ −δFext = 0 governs the FE analysis for εi+1 determination. Then, the stress tensors are updated using (6). This proce- dure is repeated for all time instants. Exponential algorithm The main difference for the exponential algorithm compared with the Euler method is the discretization of the stress tensor. As we suppose that the strain is distributed linearly within the time interval 〈ti, ti+1〉, the viscous stress tensor is solved exactly using the corresponding ODE. The linearly distributed strain induces a constant strain rate, and the governing equation for the Maxwell unit becomes σ̇p(ti + si) Gp + σp(ti + si) ηp = 1 ∆ti Dν : (εi+1 −εi) , (8) where si is an auxiliary time variable starting from ti and going forward in time. The relation between the original time t and the variable si is si = t− ti. Considering the boundary condition σp(si = 0) = σp,i, the stress evaluated in t = ti+1 reads as σp,i+1 = σp,ie −∆ti τp + ηp ∆ti ( 1 − e− ∆ti τp ) Dν : (εi+1 −εi). (9) The backward substitution into the equilibrium con- dition (1) in ti+1 gives ∫ Ωfoil δε :  G∞ + ∑ p∈P ηp ∆ti ξp,i   Dν : εi+1 dΩ − ∫ Ωfoil δε : ∑ p∈P ηp ∆ti ξp,iDν : εi dΩ (10) + ∫ Ωfoil δε : ∑ p∈P e− ∆ti τp σp,i dΩ −δFext = 0, where another auxiliary variable ξp,i = 1−e −∆ti τp is introduced. 4. Numerical results 4.1. Example All calculations were performed on a longitudinal cross- section of a laminated glass beam assuming a plane- stress state, see Figure 2. The dimensions of the beam are: the thicknesses h1 = h3 = 10 mm, h2 = 0.76 mm, the overhanging edges l1 = l3 = 30.0 mm and their span l2 = 1040.0 mm. h1 h2 h3 l1 l2 l3 Figure 2. Scheme of laminated glass sample. The laminated glass samples were simply-supported and loaded by prescribed uniformed pressure in a vacuum chamber. The time evolution of pressure was almost the same for the two tested samples with EVA and PVB interlayer, see Figure 3. At the beginning of loading, the pressure increase is almost quadratic, whereas later the slope becomes linear. 0 5 10 15 20 25 30 35 40 Time [s] 5 0 5 10 15 20 25 30 35 40 Pr es su re lo ad [k Pa ] Pressure load Figure 3. Time evolution of pressure recorded during the experimental testing. 4.2. Material parameters Glass is treated as a perfectly elastic material de- scribed by the Young modulus Eglass = 76.6 GPa obtained from indentation tests and by the Poisson ratio νglass = 0.22. Whereas the interlayer polymer is effectively described by the generalised Maxwell model to capture its viscoelastic time/temperature depen- dent nature. In the tested samples, an ethylene-vinyl acetate (EVA) and a polyvinyl butyral (PVB) foils were used as the interlayer material, with the material constants taken from [3] and given in Table 1. The Poisson ratio was assumed to be constant, νfoil = 0.49. 88 vol. 26/2020 Viscoelastic analysis of laminated glass EVA PVB τp [s] Gp [kPa] Gp [kPa] 10−9 6933.9 – 10−8 3898.6 – 10−7 2289.2 – 10−6 1672.7 – 10−5 761.6 1,782,124.2 10−4 2401.0 519,208.7 10−3 65.2 546,176.8 10−2 248.0 216,893.2 10−1 575.6 13,618.3 100 56.3 4988.3 101 188.6 1663.8 102 445.1 587.2 103 300.1 258.0 104 401.6 63.8 105 348.1 168.4 106 111.6 – 107 127.2 – 108 137.8 – 109 50.5 – 1010 322.9 – 1011 100.0 – 1012 199.9 – Table 1. Prony series for generalised Maxwell model complemented with the long-term moduli G∞ = 682.18 kPa for EVALAM 80-120 (EVA-based foil) and G∞ = 232.26 kPa for TROSIFOL BG R20 (PVB- based foil) for T0 = 20 ◦C from [3]. The room temperature during the experimental test- ing was 24◦C. For the polymer ply, the temperature dependency is incorporated into the material model using the time-temperature superposition principle [6, Chapter 11]. Time and temperature are tied together via a multiplicative parameter aT as τp(T) = aT(T) ηp Gp . (11) The parameter aT can be evaluated through the William–Landel–Ferry equation (WLF), [7], log10 aT = − C1(T −T0) C2 + T −T0 , (12) where the constants C1 and C2 are additional material parameters, and T0 is a reference temperature. For the reference temperature T0 = 20 ◦C, the calibration process described in [3] provides the following values: C1 = 339.1 and C2 = 1185.8◦C, or C1 = 8.635 and C2 = 339.102, for EVA or PVB foil respectively. 4.3. Solver For a given time instant ti+1, the unknown displace- ment field u(ti+1) is calculated numerically by the conventional finite element method. The part of the discretized domain is illustrated in Figure 4, where the red colour represents glass layers and the blue colour the interlayer. Due to the disproportionality between the glass and interlayer thicknesses, the mesh is re- fined near the centerline with at least two elements through the ply thickness. Moreover, the quadratic basis functions were employed. The computational FEM library FEniCS was used to quickly implement the presented models for all numerical experiments. This solver with a high-level Python interface pro- vides an efficient tool for solving partial differential equations [8]. Figure 4. Finite element mesh: glass domain in red and the interlayer in blue. For modelling of the glass/polymer coupling, an indicator function I(x) was introduced as follows I(x) = { 0, if x ∈ Ωglass, 1, if x ∈ Ωfoil. (13) All viscous members of the governing equations, Sec- tion 3.1, are multiplied by this function. Therefore, the rate-dependent expressions appear only in the interlayer domain Ωfoil. Let us highlight that no delamination, glass fracture, or damage of polymer is considered in this study. 4.4. Validation and comparison The main advantage of the exponential algorithm is that the time interval can be covered by a series of logarithmically increasing steps without significant loss of accuracy. Therefore, the computational cost and the number of time increments can be significantly reduced, for example for long-term stress relaxation or creep quasi-static tests. The results for 20 logarithmically or uniformly dis- tributed time steps are presented in Figure 5. The reference solution correspondes to the backward Euler method with the time step ∆t = 0.05 s. The detailed comparison of responses on the right side shows the advantage of the exponential algorithm for logarith- mically increasing time steps, as was supposed. On the other hand, the performance of both methods is the same for uniformly distributed time steps. For uniformly distributed time steps, three different time increments were tested, i.e., 0.05 s, 0.5 s and 2.0 s. The numerical responses are presented together with the experimental data in Figure 6. This vali- dation shows that even the numerical response with 89 Jaroslav Schmidt, Alena Zemanová Acta Polytechnica CTU Proceedings 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Time [s] 0.00000 0.00005 0.00010 0.00015 0.00020 0.00025 0.00030 0.00035 0.00040 0.00045 M ax im um d ef le ct io n [m m ] Reference Exp. log BE log 0.700 0.705 0.710 0.715 0.720 Time [s] 0 2 4 6 8 10 M ax im um d ef le ct io n [m m ] 1e 7+3.4e 5 DETAIL (a) logarithmically distributed time steps 0.0 0.5 1.0 1.5 2.0 2.5 Time [s] 0.00000 0.00005 0.00010 0.00015 0.00020 0.00025 0.00030 0.00035 0.00040 0.00045 M ax im um d ef le ct io n [m m ] Reference Exp. BE 0.7880 0.7885 0.7890 0.7895 0.7900 Time [s] 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 M ax im um d ef le ct io n [m m ] 1e 7+4.2e 5 DETAIL (b) uniformly distributed time steps with ∆t = 0.125 s Figure 5. Time evolution of deflection of laminated glass sample under uniformly distributed pressure with quadratically increasing intensity according Figure 3 obtained by: Reference – Backward Euler method with ∆t = 0.05 s, Exp – Exponential algorithm, and BE – Backward Euler method. 15 20 25 30 Time [s] 0 5 10 15 M ax im um d ef le ct io n [m m ] Exponential Backward Euler Experiment ∆t = 0.05 s 15 20 25 30 Time [s] 0 5 10 15 Exponential Backward Euler Experiment ∆t = 0.5 s 15 20 25 30 Time [s] 0 5 10 15 Exponential Backward Euler Experiment ∆t = 2.0 s (a) response for samples with EVA foil 14 16 18 20 22 24 26 28 Time [s] 0 5 10 15 M ax im um d ef le ct io n [m m ] Exponential Backward Euler Experiment ∆t = 0.05 s 14 16 18 20 22 24 26 28 Time [s] 0 5 10 15 Exponential Backward Euler Experiment ∆t = 0.5 s 14 16 18 20 22 24 26 28 Time [s] 0 5 10 15 Exponential Backward Euler Experiment ∆t = 2.0 s (b) response for samples with PVB foil Figure 6. Validation of numerical responses provided by the Exponential and Backward Euler method against the experimental measurement for three different time steps: ∆t = 0.05 s, ∆t = 0.5 s and ∆t = 2.0 s. 90 vol. 26/2020 Viscoelastic analysis of laminated glass 15 20 25 30 Time [s] 0 5 10 15 M ax im um d ef le ct io n [m m ] dt=0.05s dt=0.5s dt=2.0s dt=5.0s dt=10.0s B-Euler Figure 7. Numerical response of backward Euler method for different time steps. the coarser time steps is in good agreement with the measured deflections. The dashed line representing the exponential algorithm almost coincides with the solid one for the backward Euler method in all cases. For larger time steps, the differences of numeri- cal predictions against the measured experimental response would increase, see Figure 7. On the other hand, we aim to describe the crack initiation and prop- agation in glass layers, which requires small time steps. It results from this comparison that both methods are able to capture the time/temperature-dependent response of polymer interlayers well. Further, the validation against experimental data proved that the parameters of the generalised Maxwell model for both foils, obtained from material tests and summarised in [3], provides very good predictions for quasi-static behaviour of laminated glass samples with EVA-based and PVB-based foils of the same type. The response is just slightly stiffer at the beginning of loading for the sample with EVA foil and a slightly softer before the fracture of glass for the sample with PVB. However, the largest error is under 5%. 4.5. Sufficient number of parameters The comparison of the viscoelastic response of lami- nated glass samples with three elastic simulations in Figure 8 demonstrates the influence of shear coupling of glass plies due to the interlayer. The behaviour of the polymer was approximated as an elastic one with three different values of shear moduli. The response of the laminated glass sample is often bounded by the behaviour of a monolithic glass plate with the same overall thickness and by the response of two indepen- dent plies representing the two glass plates without any coupling. We narrowed this range, assuming a constant long-term or initial shear modulus, G∞ or G0. For the sample with EVA foil, a reasonable ap- proximation of the viscoelastic response was found using the shear modulus of the interlayer correspond- ing to a half of the loading time G(t/2). For PVB, this assumption gives still a softer prediction of deflections compared with the measured data. In the next stage, we tried to reduce the number of spring-dashpot cells as each of these Maxwell cells is represented by a second-order viscous stress tensor. Therefore, it is desirable to balance the number of cells and save the computational time. The exponential and backward Euler methods provided the same response for the fine time step. Thus, only the backward Euler method was used for the analysis. The influence of the individual Maxwell cell can be assessed using the multiplicator ζp,i. It can be seen from equation (7) that the significant viscous stresses are those where the expression 1/ζp,i is not close to zero. Recalling the definition of ζp,i, 1 ζp,i = 1 1 + ∆ti τp = τp τp + ∆ti , (14) the whole fraction is close to one if the relaxation time is much more greater than the given time step τp � ∆ti, and the viscous stress is activated. On the other hand if τp � ∆ti, the denominator is large and the whole fraction is close to zero. To verify this observation we perform two numerical tests. Fig- ure 9 shows that if we remove the low-relaxation-time cells, the response is not affected (left), and conversely, removing of cells with large relaxation time have signif- icant influence on the response (right). The reference solution plotted by the red line corresponds to the original generalised Maxwell model with all cells, Ta- ble 1. The number of cells was gradually decreased. Figure 9 (left) illustrates the case, for which the cells with the smallest relaxation times were neglected. It can be seen that neglecting even first seven cells for EVA or three for PVB with the smallest relaxation times did not affect the solution as the time evolution of deflection is indistinguishable from the response of the full generalised Maxwell model. On the other hand, removing a few Maxwell cells with the largest relaxation times significantly affected the deflection of the sample, Figure 9 (right). This is caused by the fact that we start in the proposed formulation with the same long-term modulus. 5. Conclusions The numerical analysis of laminated glass samples under quasi-static loading is presented in this pa- per and validated against experimental measurements. The laminated glass plate is discretized as a 2D lon- gitudinal cross-section with two elastic glass layers and the viscoelastic interlayer. The polymer ply was characterised by the generalised Maxwell model with parameters derived from material tests in [3]. The conclusions from this study are: • Both tested methods for time integration, i.e., the exponential algorithm and the backward Euler method, can be used for the subsequent finite ele- ment analysis of glass fracture. The difference in results is negligible for values of time steps needed in proposed phase-field formulation. 91 Jaroslav Schmidt, Alena Zemanová Acta Polytechnica CTU Proceedings 15 20 25 30 Time [s] 0 5 10 15 M ax im um d ef le ct io n [m m ] Linear G_inf Linear G_0 Linear G(t/2) Exponential ∆t = 0.05 s 15 20 25 30 Time [s] 0 5 10 15 M ax im um d ef le ct io n [m m ] Linear G_inf Linear G_0 Linear G(t/2) Exponential ∆t = 0.05 s (a) sample with EVA foil (b) sample with PVB foil Figure 8. Comparison of the viscoelastic response of laminated glass samples obtained by the exponential algorithm with the elastics ones with three different shear moduli of interlayer assumed. 15 20 25 30 Time [s] 0 5 10 15 M ax im um d ef le ct io n [m m ] n=0 n=5 n=10 n=15 n=22 (Ref.) 15 20 25 30 Time [s] 0 5 10 15 M ax im um d ef le ct io n [m m ] n=0 n=5 n=10 n=15 n=22 (Ref.) (a) samples with EVA foil 15 20 25 30 Time [s] 0 5 10 15 M ax im um d ef le ct io n [m m ] n=0 n=4 n=8 n=11 (Ref.) 15 20 25 30 Time [s] 0 5 10 15 M ax im um d ef le ct io n [m m ] n=0 n=4 n=8 n=11 (Ref.) (b) samples with PVB foil Figure 9. Time evolution of deflections predicted by the backward Euler method with time step ∆t = 0.5 s: Reference – numerical simulation with all Maxwell (spring-dashpot) cells or Reduced – with n Maxwell cells if the cells with the smallest relaxation times (left) or the largest relaxation times (right) are removed. 92 vol. 26/2020 Viscoelastic analysis of laminated glass • For laminated glass samples with EVA-based or PVB-based foils, the formerly fitted material pa- rameters provides reliable numerical predictions of deflections. The errors in deflections were under 5%. • The computational cost can be reduced without the loss of accuracy using a constant shear modulus for the interlayer G(t/2) for EVA-based samples. For PVB-based samples, three Maxwell cells with the smallest relaxation times can be neglected without an observable change in response of the sample. Finally, let us conclude that this study was fo- cused on the comparison of the response of lami- nated glass before fracture with emphasis on the time/temperature-dependent behaviour of polymer interlayers. No phase-field description of glass frac- ture is presented and discussed in this paper, and also no damage of polymer foils is assumed. The reason for this is the fact that based on the phase-field model for- mulation and setting, the linear elastic phase can also be affected by the damage evolution if the growth of damage starts immediately when the strain becomes nonzero. Then, the model does not reproduce the linear elastic behaviour of laminated glass in the first stage, and the stress-strain diagram becomes curved. Therefore, the phase-field model for brittle fracture of glass will be discussed separately to avoid the possi- ble deviation from linearity and misinterpretation of results. Acknowledgements This publication was supported by the Czech Science Foundation, the grant No. 19-15326S and by the Grant Agency of the Czech Technical University in Prague, grant No. SGS19/033/OHK1/1T/11. References [1] M. Haldimann, A. Luible, M. Overend, et al. Structural Use of Glass. Structural engineering documents. International Association for Bridge and Structural Engineering, 2008. [2] B. Bourdin, G. A. Francfort, J. J. Marigo. The variational approach to fracture. Journal of elasticity 91:5—-148, 2008. doi:10.1007/s10659-007-9107-3. [3] T. Hána, T. Janda, J. Schmidt, et al. Experimental and Numerical Study of Viscoelastic Properties of Polymeric Interlayers Used for Laminated Glass : Determination of Material Parameters. Materials 12(14), 2019. doi:10.3390/ma12142241. [4] Z. P. Bažant, M. Jirásek. Creep and hygrothermal effects in concrete structures. Springer, 2018. doi:10.1007/978-94-024-1138-6. [5] A. Zemanová, J. Zeman, M. Šejnoha. Comparison of viscoelastic finite element models for laminated glass beams. International Journal of Mechanical Sciences 131-132:380–395, 2017. doi:10.1016/j.ijmecsci.2017.05.035. [6] J. D. Ferry. Viscoelastic Properties of Polymers. John Wiley and Sons, 1980. doi:10.1149/1.2428174. [7] M. L. Williams, R. F. Landel, J. D. Ferry. The Temperature Dependence of Relaxation Mechanisms in Amorphous Polymers and Other Glass-forming Liquids. Journal of the American Chemical society 77(14):3701–3707, 1955. doi:10.1021/ja01619a008. [8] M. S. Alnæs, J. Blechta, J. Hake, et al. The fenics project version 1.5. Archive of Numerical Software 3(100), 2015. doi:10.11588/ans.2015.100.20553. 93 https://doi.org/10.1007/s10659-007-9107-3 https://doi.org/10.3390/ma12142241 https://doi.org/10.1007/978-94-024-1138-6 https://doi.org/10.1016/j.ijmecsci.2017.05.035 https://doi.org/10.1149/1.2428174 https://doi.org/10.1021/ja01619a008 https://doi.org/10.11588/ans.2015.100.20553 Acta Polytechnica CTU Proceedings 26:86–93, 2020 1 Introduction 2 Laminated glass model 3 Time integrators 3.1 Governing equations 3.2 Numerical integrators 4 Numerical results 4.1 Example 4.2 Material parameters 4.3 Solver 4.4 Validation and comparison 4.5 Sufficient number of parameters 5 Conclusions Acknowledgements References