Log In Sign Up

Newmark algorithm for dynamic analysis with Maxwell chain model

by   Jaroslav Schmidt, et al.

This paper investigates a time-stepping procedure of the Newmark type for dynamic analyses of viscoelastic structures characterized by a generalized Maxwell model. We depart from a scheme developed for a three-parameter model by Hatada et al. in 2000, which we extend to a generic Maxwell chain and demonstrate that the resulting algorithm can be derived from a suitably discretized Hamilton variational principle. This variational structure manifests itself in an excellent stability and a low artificial damping of the integrator, as we confirm with a mass-spring-dashpot example. After a straightforward generalization to distributed systems, the integrator may find use in, e.g., fracture simulations of laminated glass units, once combined with variationally-based fracture models.


page 1

page 2

page 3

page 4


Frequentist Consistency of Generalized Variational Inference

This paper investigates Frequentist consistency properties of the poster...

A minimizing-movements approach to GENERIC systems

We present a new time discretization scheme adapted to the structure of ...

A Generalized Markov Chain Model to Capture Dynamic Preferences and Choice Overload

Assortment optimization is an important problem that arises in many prac...

Variational Osmosis for Non-linear Image Fusion

We propose a new variational model for nonlinear image fusion. Our appro...

An arbitrary order scheme on generic meshes for miscible displacements in porous media

We design, analyse and implement an arbitrary order scheme applicable to...

Connecting beams and continua: variational basis and mathematical analysis

We present a new variational principle for linking models of beams and d...

Can Subnetwork Structure be the Key to Out-of-Distribution Generalization?

Can models with particular structure avoid being biased towards spurious...

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 through-thickness heterogeneity renders mechanics of laminated glass structures intricate, e.g. [2]. 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.

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 quasi-static 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 Newmark-type [17] 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 [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.


We employ the conventional notation through the text, in which scalar quantities are denoted by a plain font, whereas bold-face letters indicate vectors or higher-order tensors. Additional nomenclature is introduced when needed.

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 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.

Figure 1: Scheme of the single-degree-of-freedom viscoelastic dynamic problem.

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


because of the serial arrangement of the spring and damper in the cell. Differentiating (2) with respect to time and using (3), we obtain


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


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


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


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


Finally, substituting equations (6b) and (9) into (1) expressed at and rearranging the terms yields


After solving Eq. (12) for the acceleration , we update velocity , displacement , and restoring forces according to equations (6a), (6b), and (9), respectively, and proceed to the next time interval.

Note that the initial acceleration , needed in the first step of the algorithm, is set to


according to the equilibrium (1) and initial (5) conditions.

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 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


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


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.

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 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


with .

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


Likewise, expressing from (28) and employing the velocity from (30) provides


Hence, expressions (33) and (32) become identical to the ones of the Newmark method (6) once setting



Employing the nodal accelerations and from (34) in the identity (28) reveals that


Further, by expressing the difference using (27), we find that


Now, after inserting the identity (35) into the discrete Euler-Lagrange equations (3.2) and subtracting (3.3) from the result, we infer that


which can be reduced to the final form


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

The variational framework (14) additionally reveals that the trajectory satisfies the energy balance condition, e.g., [23, Section 5.1]


with the internal energy , dissipated energy , and the work done by external forces given by


To later quantify the articifial dissipation induced by time discretization, we also consider the time-discrete quantities


the last two expressions correspond to the approximations of integrals in (40) with the trapezoidal rule and employing the indentities (22).

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 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.

[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
Table 1: Parameters of Maxwell chain model [10], with and kNm. Note that in Section 4.2, the stiffnesses correspond to shear moduli  [MPa].

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 Python-based scripts available at [24].

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 [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 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.

Figure 2: Accuracy of the viscous Newmark method for step (left) and harmonic (right) loadings defined by Eq. (42). Top, center, and bottom graphs show trajectories for time steps  s,  s, and  s respectively.

Numerical dissipation.

As follows from the energy equality (39

), the additional dissipation induced by the integrator can be estimated as


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.

Figure 3: Normalized energies corresponding to SDOF response to step (top) and harmonic loads (bottom). Left: the evolution of internal energy and dissipation , normalized by the work done by external forces for time step  s. Right: the evolution of numerical dissipation , normalized by the work done by external forces .

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 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 [5].

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 [24] for full details on the simulation.

time [s]
Figure 4: Snapshots of deformations of a viscoelastic cube subjected to ramp load.

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 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.

This publication was supported by the Czech Science Foundation, the grant No. 19-15326S.


  • [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)1096-9845(200002)29:2<159::AID-EQE895>3.0.CO;2-1.
  • [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)0733-9399(1999)125:4(435).
  • [4] 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.
  • [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] 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.
  • [7] 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.
  • [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. 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.
  • [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(25-26):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 visco-elastic stress analysis. International Journal of Mechanical Sciences 10(10):807–827, 1968. doi:10.1016/0020-7403(68)90022-2.
  • [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/978-94-024-1138-6.
  • [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/1097-0207(20001210)49:10<1295::AID-NME993>3.0.CO;2-W.
  • [19] 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.
  • [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 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.
  • [22] A. Bedford. Hamilton’s principle in continuum mechanics. Pitman Publishing, Boston, 1985.
  • [23] 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.
  • [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–., [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. North-Holland, 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. Springer-Verlag, Berlin Heidelberg, 2012. doi:10.1007/978-3-642-23099-8.
  • [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/s40323-016-0080-x.
  • [31] T. Roubíček. Models of dynamic damage and phase-field fracture, and their various time discretisations, 2019. 1906.04110.