1 Introduction
The motivation of this work comes from the field of dynamics of laminated glass structures. These sandwich structures consist of multiple glass layers connected with transparent polymer interlayers. Combining stiff, brittle glass with compliant viscoelastic polymers enhances structural safety, but the throughthickness heterogeneity renders mechanics of laminated glass structures intricate, e.g. [2]. In particular, time and temperaturedependent interlayer properties must be accounted for even in quasistatic analyses, e.g., [3, 4, 5, 6] and references therein.
Earlier studies [7, 8, 9, 10] have shown that the response of commonly used interlayer materials can be captured well with the Maxwell chain model combined with the timetemperature superposition principle. Because the viscoelastic model concurrently predicts material damping, vibrations of laminated glass structures can be described more accurately than in conventional structural analyses that mostly employ the Rayleigh damping, e.g. [11, Section 12.5]. This added value has been addressed in detail for free vibrations of laminated glass units, e.g. [12, 13, 14]; an extension towards the response under general dynamic loads requires the development of dedicated timestepping schemes that are in the focus of the current work.
Related work.
Dynamics of viscoelastic solids described by the Maxwell chain model leads to the system of initial value problems coupling the equation of motion with the local evolution of constitutive variables, see Section 2.1 for illustration. Because numerically integrating the full system would be costly, we will follow an alternative approach in which only the equations of motion are solved approximately, whereas the evolutionary constitutive equations are resolved in the closed form, leading to an inexpensive update formulas for internal variables entering the equations of motion. This approach has been pioneered for quasistatic problems by Zienkiewicz et al. [15]; see also [16, Section 5.2] for a comprehensive review. To the best of our knowledge, Hatada et al. [1] were the only ones who used this strategy in dynamics, although no reference to the original work [15] was made. In particular, they developed a Newmarktype [17] algorithm for the threeparameter Maxwell model and used it to predict the response of planar frames to earthquake loading.
Novelty.
Our work further develops the contribution [1] in three aspects. First, in Section 2.2, we present a compact derivation of the Newmark scheme for a generic Maxwell chain, closely following the original exposition [15]. Second, in Section 3, we show that the algorithm can be interpreted as a variational integrator [18], in the sense that it can be derived from the Hamilton variational principle combined with suitable time discretization. The variational structure endows the integrator with good numerical stability and low numerical dissipation, as demonstrated in Section 4.1 with selected examples. Moreover, the scheme can be easily combined with variational approaches to fracture, e.g., [19, 20, 21], which is of independent interest when simulating the behavior of laminated glass under impact loads. Third, in Section 4.2, we outline how to extend the algorithm to a continuum formulation and complement the theoretical considerations with an illustrative 3D finite element simulation.
Notation.
2 Newmark method
In this section, we analyze a single degree of freedom (SDOF) model of a mass supported with a Maxwell chain, consisting of the parallel connection of an elastic spring and multiple springdashpot cells, see Figure
1 for illustration and, e.g., [16, Section A] for further details. In particular, in Section 2.1 we review the equations of motion, which we subsequently discretize with the average acceleration version of the Newmark method [17] in Section 2.2.2.1 Governing equations
As follows from the scheme in Figure 1, the problem under consideration is specified with the timedependent load , the particle mass and the Maxwell chain model parameters: stiffness of the elastic spring , spring stiffness and damper viscosity of the th Maxwell cell; stands for the number of Maxwell cells.
Equilibrium of the forces acting on the mass requires
(1) 
where denotes the displacement of the mass, its acceleration, and the restoring force of the th cell.
For the th Maxwell cell, the displacement splits into an elastic part of the spring and a viscous part of the damper :
(2) 
recall Figure 1. The restoring force of the th Maxwell cell satisfies
(3) 
because of the serial arrangement of the spring and damper in the cell. Differentiating (2) with respect to time and using (3), we obtain
(4) 
In summary, the motion of SDOF model is described with the coupled system ordinary differential equations (ODEs) (1) and (4), complemented with the initial conditions
(5) 
where and stand for the initial mass displacement and velocity, is the initial force in the th Maxwell cell, and .
2.2 Discretization
The time interval of interest is divided into time instants ; for notational simplicity, we assume equidistant partitioning of the constant time step .
Considering the Newmark integration scheme [17] with an average acceleration on the time interval , the velocity and displacement within the interval varies as
(6a)  
(6b) 
where is the local time variable within the interval and abbreviates to render the notation compact.
Substituting (6a) into (4) reveals that the evolution of the restoring force satisfies
(7) 
This Cauchy problem with the initial condition has the solution
(8) 
Thus, at the end of the time interval with , we have
(9) 
where
(10) 
denotes the effective relaxation time of the th cell and the auxiliary factors are given by
(11) 
Finally, substituting equations (6b) and (9) into (1) expressed at and rearranging the terms yields
(12) 
3 Variational integrators
Having derived the Newmark viscoelastic algorithm by conventional means, we now demonstrate its variational structure, by adapting the general arguments on variational integrators by Kane et al. [18] to the current setting.
We will proceed in four steps. In Section 3.1, we show that the governing equations from Section 2.1 follow from the EulerLagrange (EL) equations of a suitably defined energy functional. Its discretization then provides the governing equations of the corresponding variational integrator introduced in Section 3.2. In Section 3.3, we demonstrate the equivalence of the integrator to the Newmark algorithm from Section 2.2. In the last step, Section 3.4, we comment on the energy conservation properties of the time integration scheme.
3.1 Variational framework
We postulate that the trajectory of a constrained dissipative mechanical system in the state space is given by the EulerLagrange equations, e.g., [22, Section 1.3]
(14a)  
(14b) 
Here, stands for the dissipation potential, for the Lagrangian of the problem, and denotes the Lagrange multipliers associated with the kinematic constraint function . Besides, these equations correspond to the stationarity conditions of the action functional
(15) 
perturbed by the dissipative forces . Note that the hat symbol in (15) now distinguishes the test quantities from the true trajectories defined with (14).
For the problem from Figure 1, the state variable
(16) 
collects the total displacement and the displacements of both components of each Maxwell cell; the state space . The Lagrangian has the standard form
(17) 
involving the kinetic energy , potential energy of deformation , and external forces given by
(18a)  
(18b)  
(18c) 
The kinematical constraints take the form
(19) 
the space of the Lagrange multiplies then becomes . The last component of the general framework (14) is provided by the dissipation potential
(20) 
involving solely the viscous displacements of all cells.
In this setting, the EL equation (14a) represents the system of optimality conditions. The first one, corresponding to the total displacement , attains the form
(21) 
while the remaining conditions read as
(22) 
with . It is thus evident that the multipliers play role of the viscous force and, because the optimality (14b) and compatibility (4) conditions coincide, the current setting is equivalent to the one of Section 2.1.
3.2 Discretization
Recall that the incremental algorithm of Section 2.2 relies on the discretization of the total displacements, from which the evolution of cellrelated variables , , and follows in the closed form. To mimic this structure, only the total displacements will be determined from the discrete (nondissipative) EL equations, whereas the remaining quantities are determined from the nondiscretized optimality conditions (22) and (14b).
To this goal, we consider the same discretization of the time interval as in Section 2.2 and introduce the discretized action functional
(23) 
with the discrete Lagrangian given by [18, Eq. (2)]
(24)  
The stationarity conditions at time , with read as
(25) 
which delivers the governing equations of the variational integrator in the form^{1}^{1}1Notice that we assume the Lagrange multipliers to be given abritrary quantities, similarly to the forcing terms . Once we establish the equivalance to the Newmark algorithm, their values follow from the update formula (9) from Section 2.
(26) 
3.3 Equivalence to Newmark
We will proceed with additional two steps to show that the optimality conditions (3.2) correspond to the Newmark integration scheme from Section 2. First, we demonstrate that the displacements provide definitions of velocities and accelerations consistent with the kinematic assumptions in Eq. (6). Second, we show that the discreteintime quantities satisfy the equations of motion (1).
Kinematics.
Following [18, Section 2.2], we start from introducing auxiliary accelerations
(27) 
for . Summing with comparing the result with (3.2) provides
(28) 
with .
Equilibrium.
Employing the nodal accelerations and from (34) in the identity (28) reveals that
(35) 
Further, by expressing the difference using (27), we find that
(36) 
Now, after inserting the identity (35) into the discrete EulerLagrange equations (3.2) and subtracting (3.3) from the result, we infer that
(37) 
which can be reduced to the final form
(38) 
Indeed, the equivalence between (38) and (3.3) for holds because of the choice of the initial acceleration (13), and for it follows by induction.
3.4 Energy balance
4 Examples
In this section, we demonstrate the performance of the developed Newmark algorithm with two examples. The first one in Section 4.1 addresses the accuracy and numerical energy dissipation of the integrator for the singledegreeoffreedom system from Figure 1. The followup example in Section 4.2 outlines an extension of the scheme towards continuum models.
[kNm]  [s]  [kNm]  [s] 

6933.9  10  445.1  10 
3898.6  10  300.1  10 
2289.2  10  401.60  10 
1672.7  10  348.1  10 
761.60  10  111.6  10 
2401.0  10  127.2  10 
65.200  10  137.8  10 
248.00  10  50.5  10 
575.60  10  322.9  10 
56.30  10  100.0  10 
188.6  10  199.9  10 
Data of the generalized Maxwell chain used in both examples appear in Table 1; they represent real PVB material with sufficiently short and long relaxation times for testing algorithm robustness. For more information on experimental procedures to determine these parameters, see [10]. All results presented in this section are reproducible with Pythonbased scripts available at [24].
4.1 Discrete problem
We consider the following two types of loading:
(42a)  
(42b) 
corresponding to ramp and harmonic loads, respectively. In both cases, we set the amplitude MN and the mass kg to scale the displacement amplitude to m. Initial displacement, velocity, and forces in Maxwell cells were set to zero; recall (5). As for the Newmark algorithm, we set the time steps to , , and s.
Accuracy
of the Newmark algorithm is checked by comparing its trajectories with the reference ones, obtained with the adaptive solver lsoda
[25] — available through odeint
function of Scipy library [26] — applied to the full initial value problem (1), (4), and (5).
Results appear in Figure 2 and demonstrate that the Newmark algorithm is stable even for coarse time steps, thanks to its variational structure. The errors behave consistently with findings for Newmarkfamily methods applied to linearly dampened systems, e.g. [27, Section B.II.5]. In particular, the numerical dispersion (understood as the error in periods) and dissipation (error in amplitudes) decays as , where stands for the angular frequency of the response. For s, the trajectories predicted by the Newmark scheme closely match the reference ones.
Numerical dissipation.
As follows from the energy equality (39
), the additional dissipation induced by the integrator can be estimated as
(43) 
with the individual terms provided by Eq. (41). The evolution of these quantities for the step and harmonic loading appears in Figure 3, considering the time interval s.
For both loads, we observe that the work done by external forces eventually distributes between the internal and dissipated energies; the ratio stabilizes at for the step load and for harmonic loading the ratio reaches about . The artificial dissipation is only significant for the coarsest step of s; for s it reaches the value of about ‰ and further deteriorates with a decreasing time step. This confirms excellent energy conservation properties of the scheme, especially when taking into account that the error introduced by the trapezoidal rule in (41) is of order .
4.2 Generalization
Additional constitutive assumptions must be adopted to extend the SDOF models into a continuum formulation. Here, we assume that the Maxwell model applies when modeling the response under shear. The spring stiffnesses thus become shear moduli , and that the Poisson ratio is a timeindependent material constant. This assumption considerably simplifies the multidimensional constitutive law, e.g., [16, Section 2.4], and provides the same results as the conventional volumetricdeviatoric split for our target applications [5].
Under these assumptions, the weak form of the equations of motion attains the form, e.g. [11, Part III]:
(44) 
where stands for the virtual work done by external loads on a virtual displacement , denotes the acceleration, and the small strain tensor is obtained as the symmetric part of displacement gradient, ; virtual strain is defined in the same way. The material is characterized by its density , longterm shear modulus of the Maxwell chain , and the dimensionless tensor corresponding to the stiffness tensor of an isotropic material of unit shear modulus and the Poisson ratio of . The stresses carried by individual cells follow as the solution of initial value problems
(45) 
The evolution of the state variables and is specified with the initial conditions on the displacements, velocities, and cell stresses:
(46) 
The comparison of the initial value problems specified with Eqs. (1), (4), and (5) and Eqs. (4.2), (45), and (46) reveals that the derivation of the Newmarktype scheme follows exactly the steps as in Section 2.2. As a result, the following variational problem needs to be solved at time :
(47) 
with the parameters , , and provided by Eqs. (10) and (11); recall that denotes the virtual work done by external forces. Once the the accelerations are obtained from the weak form (47), the displacements , velocities , and the cell stresses are updated according to Eqs. (6) and (9), respectively.
The outlined formulation (47) was further discretized with the finite element method and implemented in FEniCS project [28, 29] version 2018.1. As an indicative example, we consider a unit cube, see Figure 4, fixed on the bottom surface and subjected to a ramp load (42a) with the tensile traction of intensity 1.0 Nm perpendicular to the top surface. The material response is characterized by the Maxwell chain parameters from Table 1 and the value of the Poisson ratio . In the numerical resolution, we discretize the sample into identical 1,000 hexahedron elements and consider the time step of 0.01 s.
The snapshots of the vibrations reveal a similar behavior to the SDOF example, recall Figure 2, namely the attenuation of the propagating waves by viscous damping. An interested reader is invited to the dataset [24] for full details on the simulation.
5 Conclusions
In this contribution, we have developed a Newmark integration scheme for viscoelastic solids characterized by the generalized Maxwell model. Besides the direct derivation, we have shown the scheme can be derived from the Hamilton variational principle combined with a suitable structurepreserving time discretization. This variational structure is then reflected in the longterm stability and low energy dissipation of the resulting scheme, which has been confirmed with selected numerical examples.
As the next step, we will combine the continuum framework outlined in Section 4.2 with Newmarktype solvers for variational fracture models, e.g. [30, 31], to extend the currently available approaches to simulating the response of laminated glass structures under impact.
Acknowledgements.
This publication was supported by the Czech Science Foundation, the grant No. 1915326S.References
 [1] T. Hatada, T. Kobori, M. Ishida, N. Niwa. Dynamic analysis of structures with Maxwell model. Earthquake Engineering & Structural Dynamics 29(2):159–176, 2000. doi:10.1002/(SICI)10969845(200002)29:2<159::AIDEQE895>3.0.CO;21.
 [2] M. Haldimann, A. Luible, M. Overend. Structural Use of Glass, vol. 10 of Structural Engineering Documents. IABSE, Zürich, Switzerland, 2008.
 [3] A. van Duser, A. Jagota, S. J. Bennison. Analysis of glass/Polyvinyl Butyral laminates subjected to uniform pressure. Journal of Engineering Mechanics 125(4):435–442, 1999. doi:10.1061/(ASCE)07339399(1999)125:4(435).
 [4] L. Galuppi, G. RoyerCarfagni. The design of laminated glass under timedependent loading. International Journal of Mechanical Sciences 68:67–75, 2013. doi:10.1016/j.ijmecsci.2012.12.019.
 [5] A. Zemanová, J. Zeman, M. Šejnoha. Comparison of viscoelastic finite element models for laminated glass beams. International Journal of Mechanical Sciences 131132:380–395, 2017. doi:10.1016/j.ijmecsci.2017.05.035.
 [6] A. Zemanová, J. Zeman, T. Janda, M. Šejnoha. Layerwise numerical model for laminated glass plates with viscoelastic interlayer. Structural Engineering and Mechanics 65(4):369–380, 2018. doi:10.12989/sem.2018.65.4.369.
 [7] L. Andreozzi, S. B. Bati, M. Fagone, et al. Dynamic torsion tests to characterize the thermoviscoelastic properties of polymeric interlayers for laminated glass. Construction and Building Materials 65:1–13, 2014. doi:10.1016/j.conbuildmat.2014.04.003.
 [8] Y. Shitanoki, S. Bennison, Y. Koike. A practical, nondestructive method to determine the shear relaxation modulus behavior of polymeric interlayers for laminated glass. Polymer Testing 37:59–67, 2014. doi:10.1016/j.polymertesting.2014.04.011.
 [9] I. Mohagheghian, Y. Wang, L. Jiang, et al. Quasistatic bending and low velocity impact performance of monolithic and laminated glass windows employing chemically strengthened glass. European Journal of Mechanics – A/Solids 63:165–186, 2017. doi:10.1016/j.euromechsol.2017.01.006.
 [10] 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):2241, 2019. doi:10.3390/ma12142241.
 [11] R. W. Clough, J. Penzien. Dynamics of Structures. Computers & Structures, Inc, Berkeley, 3rd edn., 2003.
 [12] Y. Koutsawa, et al. Static and free vibration analysis of laminated glass beam on viscoelastic supports. International Journal of Solids and Structures 44(2526):8735–8750, 2007. doi:10.1016/j.ijsolstr.2007.07.009.
 [13] M. L. Aenlle, F. Pelayo. Frequency Response of Laminated Glass Elements: Analytical Modeling and Effective Thickness. Applied Mechanics Reviews 65(2):020802–020802–13, 2013. doi:10.1115/1.4023929.
 [14] A. Zemanová, J. Zeman, T. Janda, et al. On modal analysis of laminated glass: Usability of simplified methods and Enhanced Effective Thickness. Composites Part B: Engineering 151:92–105, 2018. doi:10.1016/j.compositesb.2018.05.032.
 [15] O. C. Zienkiewicz, M. Watson, I. P. King. A numerical method of viscoelastic stress analysis. International Journal of Mechanical Sciences 10(10):807–827, 1968. doi:10.1016/00207403(68)900222.
 [16] Z. P. Bažant, M. Jirásek. Creep and Hygrothermal Effects in Concrete Structures, vol. 225 of Solid Mechanics and Its Applications. Springer, Dordrecht, 2018. doi:10.1007/9789402411386.
 [17] N. M. Newmark. A Method of Computation for Structural Dynamics. Journal of the Engineering Mechanics Division 85(3):67–94, 1959.
 [18] C. Kane, J. E. Marsden, M. Ortiz, M. West. Variational integrators and the Newmark algorithm for conservative and dissipative mechanical systems. International Journal for Numerical Methods in Engineering 49(10):1295–1325, 2000. doi:10.1002/10970207(20001210)49:10<1295::AIDNME993>3.0.CO;2W.
 [19] B. Bourdin, G. A. Francfort, J.J. Marigo. The variational approach to fracture. Journal of Elasticity 91(13):5–148, 2008. doi:10.1007/s1065900791073.
 [20] M. Buliga. Hamiltonian inclusions with convex dissipation with a view towards applications. Annals of the Academy of Romanian Scientists Series on Mathematics and its Applications 1(2):228–251, 2009.
 [21] B. Bourdin, C. J. Larsen, C. L. Richardson. A timediscrete model for dynamic fracture based on crack regularization. International Journal of Fracture 168(2):133–143, 2011. doi:10.1007/s107040109562x.
 [22] A. Bedford. Hamilton’s principle in continuum mechanics. Pitman Publishing, Boston, 1985.
 [23] A. Mielke, T. Roubíček. Rateindependent systems: theory and application, vol. 655 of Applied Mathematical Sciences. Springer, New York, 2015. doi:10.1007/9781493927067.
 [24] J. Schdmit, T. Janda, A. Zemanová, et al. Source codes for preprint Newmark algorithm for dynamic analysis with Maxwell chain model, 2019. doi:10.5281/zenodo.3531802.
 [25] L. Petzold. Automatic selection of methods for solving stiff and nonstiff systems of ordinary differential equations. SIAM Journal on Scientific and Statistical Computing 4(1):136–148, 1983. doi:10.1137/0904010.

[26]
E. Jones, T. Oliphant, P. Peterson, et al. SciPy: Open source scientific tools for Python, 2001–.
http://www.scipy.org, [accessed November 11, 2019].  [27] T. J. R. Hughes. Analysis of transient algorithms with particular reference to stability behavior. In T. Belytschko, T. J. R. Hughes (eds.), Computational Methods for Transient Analysis, chap. 2, pp. 67–155. NorthHolland, Amsterdam, 1983.
 [28] A. Logg, K.A. Mardal, G. Wells (eds.). Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book. Lecture Notes in Computational Science and Engineering. SpringerVerlag, Berlin Heidelberg, 2012. doi:10.1007/9783642230998.
 [29] 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.
 [30] T. Li, J.J. Marigo, D. Guilbaud, S. Potapov. Numerical investigation of dynamic brittle fracture via gradient damage models. Advanced Modeling and Simulation in Engineering Sciences 3(1):26, 2016. doi:10.1186/s403230160080x.
 [31] T. Roubíček. Models of dynamic damage and phasefield fracture, and their various time discretisations, 2019. 1906.04110.