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 through-thickness heterogeneity renders mechanics of laminated glass structures intricate, e.g. . In particular, time- and temperature-dependent interlayer properties must be accounted for even in quasi-static 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 time-temperature 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 time-stepping schemes that are in the focus of the current 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 quasi-static problems by Zienkiewicz et al. ; see also [16, Section 5.2] for a comprehensive review. To the best of our knowledge, Hatada et al.  were the only ones who used this strategy in dynamics, although no reference to the original work  was made. In particular, they developed a Newmark-type  algorithm for the three-parameter Maxwell model and used it to predict the response of planar frames to earthquake loading.
Our work further develops the contribution  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 . Second, in Section 3, we show that the algorithm can be interpreted as a variational integrator , 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.
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 spring-dashpot cells, see Figure1 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  in Section 2.2.
2.1 Governing equations
As follows from the scheme in Figure 1, the problem under consideration is specified with the time-dependent 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
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 :
recall Figure 1. The restoring force of the -th Maxwell cell satisfies
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  with an average acceleration on the time interval , the velocity and displacement within the interval varies as
where is the local time variable within the interval and abbreviates to render the notation compact.
This Cauchy problem with the initial condition has the solution
Thus, at the end of the time interval with , we have
denotes the effective relaxation time of the -th cell and the auxiliary factors are given by
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.  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 Euler-Lagrange (E-L) 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 Euler-Lagrange equations, e.g., [22, Section 1.3]
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
For the problem from Figure 1, the state variable
collects the total displacement and the displacements of both components of each Maxwell cell; the state space . The Lagrangian has the standard form
involving the kinetic energy , potential energy of deformation , and external forces given by
The kinematical constraints take the form
the space of the Lagrange multiplies then becomes . The last component of the general framework (14) is provided by the dissipation potential
involving solely the viscous displacements of all cells.
In this setting, the E-L equation (14a) represents the system of optimality conditions. The first one, corresponding to the total displacement , attains the form
while the remaining conditions read as
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.
Recall that the incremental algorithm of Section 2.2 relies on the discretization of the total displacements, from which the evolution of cell-related variables , , and follows in the closed form. To mimic this structure, only the total displacements will be determined from the discrete (non-dissipative) E-L equations, whereas the remaining quantities are determined from the non-discretized 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
with the discrete Lagrangian given by [18, Eq. (2)]
The stationarity conditions at time , with read as
which delivers the governing equations of the variational integrator in the form111Notice 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.
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 discrete-in-time quantities satisfy the equations of motion (1).
Following [18, Section 2.2], we start from introducing auxiliary accelerations
for . Summing with comparing the result with (3.2) provides
The discrete linear momenta follow standardly from
and, using from Eq. (27), they can be evaluated as
Expressing according to the previous relation and employing (28) provides
from which we obtain
Further, by expressing the difference using (27), we find that
which can be reduced to the final form
3.4 Energy balance
with the internal energy , dissipated energy , and the work done by external forces given by
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 single-degree-of-freedom system from Figure 1. The follow-up example in Section 4.2 outlines an extension of the scheme towards continuum models.
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 . All results presented in this section are reproducible with Python-based scripts available at .
4.1 Discrete problem
We consider the following two types of loading:
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.
of the Newmark algorithm is checked by comparing its trajectories with the reference ones, obtained with the adaptive solver
lsoda  — available through
odeint function of Scipy library  — 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 Newmark-family 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.
As follows from the energy equality (39
), the additional dissipation induced by the integrator can be estimated as
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 .
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 time-independent material constant. This assumption considerably simplifies the multi-dimensional constitutive law, e.g., [16, Section 2.4], and provides the same results as the conventional volumetric-deviatoric split for our target applications .
Under these assumptions, the weak form of the equations of motion attains the form, e.g. [11, Part III]:
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 , long-term 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
The evolution of the state variables and is specified with the initial conditions on the displacements, velocities, and cell stresses:
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 Newmark-type scheme follows exactly the steps as in Section 2.2. As a result, the following variational problem needs to be solved at time :
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  for full details on the simulation.
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 structure-preserving time discretization. This variational structure is then reflected in the long-term 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 Newmark-type 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. 19-15326S.
-  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)1096-9845(200002)29:2<159::AID-EQE895>3.0.CO;2-1.
-  M. Haldimann, A. Luible, M. Overend. Structural Use of Glass, vol. 10 of Structural Engineering Documents. IABSE, Zürich, Switzerland, 2008.
-  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)0733-9399(1999)125:4(435).
-  L. Galuppi, G. Royer-Carfagni. The design of laminated glass under time-dependent loading. International Journal of Mechanical Sciences 68:67–75, 2013. doi:10.1016/j.ijmecsci.2012.12.019.
-  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.
-  A. Zemanová, J. Zeman, T. Janda, M. Šejnoha. Layer-wise 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.
-  L. Andreozzi, S. B. Bati, M. Fagone, et al. Dynamic torsion tests to characterize the thermo-viscoelastic properties of polymeric interlayers for laminated glass. Construction and Building Materials 65:1–13, 2014. doi:10.1016/j.conbuildmat.2014.04.003.
-  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.
-  I. Mohagheghian, Y. Wang, L. Jiang, et al. Quasi-static 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.
-  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.
-  R. W. Clough, J. Penzien. Dynamics of Structures. Computers & Structures, Inc, Berkeley, 3rd edn., 2003.
-  Y. Koutsawa, et al. Static and free vibration analysis of laminated glass beam on viscoelastic supports. International Journal of Solids and Structures 44(25-26):8735–8750, 2007. doi:10.1016/j.ijsolstr.2007.07.009.
-  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.
-  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.
-  O. C. Zienkiewicz, M. Watson, I. P. King. A numerical method of visco-elastic stress analysis. International Journal of Mechanical Sciences 10(10):807–827, 1968. doi:10.1016/0020-7403(68)90022-2.
-  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/978-94-024-1138-6.
-  N. M. Newmark. A Method of Computation for Structural Dynamics. Journal of the Engineering Mechanics Division 85(3):67–94, 1959.
-  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/1097-0207(20001210)49:10<1295::AID-NME993>3.0.CO;2-W.
-  B. Bourdin, G. A. Francfort, J.-J. Marigo. The variational approach to fracture. Journal of Elasticity 91(1-3):5–148, 2008. doi:10.1007/s10659-007-9107-3.
-  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.
-  B. Bourdin, C. J. Larsen, C. L. Richardson. A time-discrete model for dynamic fracture based on crack regularization. International Journal of Fracture 168(2):133–143, 2011. doi:10.1007/s10704-010-9562-x.
-  A. Bedford. Hamilton’s principle in continuum mechanics. Pitman Publishing, Boston, 1985.
-  A. Mielke, T. Roubíček. Rate-independent systems: theory and application, vol. 655 of Applied Mathematical Sciences. Springer, New York, 2015. doi:10.1007/978-1-4939-2706-7.
-  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.
-  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.
E. Jones, T. Oliphant, P. Peterson, et al. SciPy: Open source scientific tools for Python, 2001–.http://www.scipy.org, [accessed November 11, 2019].
-  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. North-Holland, Amsterdam, 1983.
-  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. Springer-Verlag, Berlin Heidelberg, 2012. doi:10.1007/978-3-642-23099-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.
-  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/s40323-016-0080-x.
-  T. Roubíček. Models of dynamic damage and phase-field fracture, and their various time discretisations, 2019. 1906.04110.