12. Dynamics

12.1. Semi-Discrete Approach

We now include the inertial terms in the discrete equation system and consider solving

(12.1)\[\mathbf{M} \ddot{\mathbf{d}}(t) + \mathbf{F}^{\mathrm{int}}(\mathbf{d}(t)) = \mathbf{F}^{\mathrm{ext}}(\mathbf{d}(t))\]

for \(t \in [0,\mathrm{T}]\) subject to the initial conditions

(12.2)\[\mathbf{d}(0) = \mathbf{d}_0\]

and

(12.3)\[\dot{\mathbf{d}}(0) = \mathbf{v}_0 .\]

Note that in (12.1) time remains continuous, whereas spatial discretization has already been achieved by the finite element interpolations summarized in Section 10. This type of finite element approach to transient problems is sometimes referred to as the semi-discrete finite element method, since the approximation in space is performed first, leaving a set of equations discrete in space but still continuous in time. To complete the approximation, a finite differencing procedure is generally applied in time as discussed next.

12.2. Time-Stepping Procedures

As discussed in Section 11, we subdivide the time interval of interest \([0,\mathrm{T}]\) via

(12.4)\[[0,\mathrm{T}] = \bigcup_{n=0}^{N-1} \left[ t_n , t_{n+1} \right]\]

and consider the problem:

Given algorithmic approximations for the solution vector \((\mathbf{d}_n)\), velocity \((\mathbf{v}_n)\), and acceleration \((\mathbf{a}_n)\) at time \(t_n\), find approximations \(\mathbf{d}_{n+1}\), \(\mathbf{v}_{n+1}\), and \(\mathbf{a}_{n+1}\) for these quantities at time \(t_{n+1}\). Note that, in contrast to the quasistatic problem, the variable \(t\) here does have the interpretation of actual time.

A thoroughly studied topic in dynamics is the construction of effective time integrators for application to the semi-discrete equations of motion. An ideal approach possesses minimal dispersion and dissipation. As shown in Fig. 12.1, a measure of numerical dispersion is period error \((\bar{T}-T)\), and a measure of numerical dissipation is amplitude decay \((\bar{A}-A)\). Fig. 12.1 depicts a single wave with amplitude and period \(A\) and \(T\) that generically is the exact solution to the wave equation (subject to the proper initial conditions and/or external force). Numerical dispersion by the time integrator causes a wave’s frequency to decrease, thus dispersing its energy to the lower frequencies. Numerical dissipation by the time integrator causes the wave’s energy to decrease and therefore is said to dissipate its energy.

Fig. 12.1 Simple illustration of approximation error in transient time integrators.

The time integrators we consider here can all be described by a 3-parameter method called the \(\alpha\)-method of time integrators. It is also referred to as the Hilber-Hughes-Taylor Method, or HHT method, as described in Reference [[1]], which is a generalization of the well-known and pervasive Newmark family of temporal integrators (Reference [[2]]). The Newmark algorithm can be summarized in a time step \(\left[ t_n , t_{n+1} \right]\) as follows:

(12.5)\[\begin{split}\begin{gathered} \mathbf{M} \mathbf{a}_{n+1} + \mathbf{F}^{\mathrm{int}}(\mathbf{d}_{n+1}) = \mathbf{F}^{\mathrm{ext}}(\mathbf{d}_{n+1}) \\ \mathbf{d}_{n+1} = \mathbf{d}_n + \Delta t \mathbf{v}_n + \frac{\Delta t^2}{2} \left[ (1- 2 \beta) \mathbf{a}_n + 2 \beta \mathbf{a}_{n+1} \right] \\ \mathbf{v}_{n+1} = \mathbf{v}_n + \Delta t \left[ (1-\gamma)\mathbf{a}_n + \gamma \mathbf{a}_{n+1} \right] , \end{gathered}\end{split}\]

where \(\beta\) and \(\gamma\) are algorithmic parameters that define the stability and accuracy characteristics of the method.

The extension of the Newmark family of integrators to the \(\alpha\)-method of integrators is accomplished with the addition of the parameter, \(\alpha\):

(12.6)\[\begin{split}\begin{gathered} \mathbf{M} \mathbf{a}_{n+1} + (1+\alpha)\mathbf{F}^{\mathrm{int}}(\mathbf{d}_{n+1}) - \alpha \mathbf{F}^{\mathrm{int}}(\mathbf{d}_n) = (1+\alpha)\mathbf{F}^{\mathrm{ext}}(\mathbf{d}_{n+1}) - \alpha \mathbf{F}^{\mathrm{ext}}(\mathbf{d}_n)\\ \mathbf{d}_{n+1} = \mathbf{d}_n + \Delta t \mathbf{v}_n + \frac{\Delta t^2}{2} \left[ (1- 2 \beta) \mathbf{a}_n + 2 \beta \mathbf{a}_{n+1} \right] \\ \mathbf{v}_{n+1} = \mathbf{v}_n + \Delta t \left[ (1-\gamma)\mathbf{a}_n + \gamma \mathbf{a}_{n+1} \right] , \end{gathered}\end{split}\]

where, as expected, setting \(\alpha\) to zero reduces the HHT integrator to Newmark’s method. Although a wide range of algorithms exist corresponding to the different available choices of \(\beta\) and \(\gamma\), two algorithms in particular are significant:

• Central Differences (\(\alpha=0\), \(\beta=0\), \(\gamma=1/2\)). This integrator is second-order accurate in time and only conditionally stable, meaning that the linearized stability is only retained when \(\Delta t\) is less than some critical value. This algorithm is an example of an explicit finite element integrator discussed in Section 12.3.

• Trapezoid rule (\(\alpha=0\), \(\beta=1/4\), \(\gamma=1/2\)). This integrator is also second-order accurate but unconditionally stable for linear problems, meaning that the spectral radii of the integrator remains less than one in modulus for any time step \(\Delta t\) (in linear problems). This algorithm is an example of an implicit finite element integrator discussed in Section 12.4

12.3. Explicit Finite Element Methods

Examining the central differences algorithm, we substitute \(\beta=0\), \(\gamma=1/2\) into (12.6) to obtain

(12.7)\[\begin{split}\begin{gathered} \mathbf{a}_{n+1} = \mathbf{M}^{-1} \left( \mathbf{F}^{\mathrm{ext}}(\mathbf{d}_{n+1}) - \mathbf{F}^{\mathrm{int}}(\mathbf{d}_{n+1}) \right) \\ \mathbf{d}_{n+1} = \mathbf{d}_n + \Delta t \mathbf{v}_n + \frac{\Delta t^2}{2} \mathbf{a}_n \\ \mathbf{v}_{n+1} = \mathbf{v}_n + \frac{\Delta t}{2} \left[ \mathbf{a}_n + \mathbf{a}_{n+1} \right] , \end{gathered}\end{split}\]

where the first equation has been written as solved for \(\mathbf{a}_{n+1}\).

(12.7) can be used to explain why this formulation is termed explicit. Given the three vectors \(\left\{ \mathbf{a}_n, \mathbf{v}_n, \mathbf{d}_n \right\}\), the data at \(t_{n+1}\), \(\left\{\mathbf{a}_{n+1}, \mathbf{v}_{n+1}, \mathbf{d}_{n+1} \right\}\) can be computed explicitly, i.e., without the need for solution of coupled equations provided the mass matrix \(\mathbf{M}\) is a diagonal matrix.

It is important to note approximation properties of the explicit time integrator (see Reference [[3]]). By itself, the explicit time integrator causes the period to be shortened. However, a lumped or diagonalized mass matrix as opposed to a consistent mass matrix causes the period to be elongated. For one-dimensional problems with uniform meshes the period error cancels exactly. In the words of Reference [[3]], these compensating errors generally produce a matched approach. Thus a lumped mass matrix gives rise to the fully explicit algorithm, requiring only an inverse of a diagonal matrix.

Although this form of the central difference formulation ((12.7)) is readily obtained from the Newmark formulas, it does not give insight into the source of the central difference terminology and, in fact, does not represent the (historical) manner in which the integrator is ordinarily developed or implemented. To see the usual form, one starts with the difference formulas for acceleration and velocity (see e.g., The Difference Calculus, Chapter 9 in Reference [[4]]):

(12.8)\[\mathbf{a}_n = \frac{ \mathbf{v}_{n+1/2} - \mathbf{v}_{n-1/2} } { t_{n+1/2} - t_{n-1/2}} ,\]

and

(12.9)\[\mathbf{v}_{n+1/2} = \frac{ \mathbf{d}_{n+1} - \mathbf{d}_{n} } { t_{n+1} - t_{n}} ,\]

where, as shown in Fig. 12.2, the time axis is discretized with notions of whole step configurations at times \(t_{n-1}, t_n, t_{n+1}\) and half-step configurations at times \(t_{n-1/2}, t_{n+1/2}, \dots\)

Fig. 12.2 Graphical construction of the central difference time integrator.

Rearranging, these difference formulas ((12.8) and (12.9)) can be converted into integration formulas:

(12.10)\[\begin{split}\begin{gathered} \mathbf{v}_{n+1/2} = \mathbf{v}_{n-1/2} + \frac{1}{2} \left( \Delta t_{n-1/2} + \Delta t_{n+1/2} \right) \mathbf{a}_n \\ \mathbf{d}_{n+1} = \mathbf{d}_n + \Delta t_{n+1/2} \mathbf{v}_{n+1/2} \end{gathered}\end{split}\]

Combining these integration formulas with the equilibrium equation evaluated at \(t_n\), we can express the algorithm as

(12.11)\[\begin{split}\begin{gathered} \mathbf{a}_n = \mathbf{M}^{-1} \left[ \mathbf{F}^{\mathrm{ext}}(\mathbf{d}_n) - \mathbf{F}^{\mathrm{int}}(\mathbf{d}_n)\right] \\ \mathbf{v}_{n+1/2} = \mathbf{v}_{n-1/2} + \frac{1}{2} \left( \Delta t_{n-1/2} + \Delta t_{n+1/2} \right) \mathbf{a}_n \\ \mathbf{d}_{n+1} = \mathbf{d}_n + \Delta t_{n+1/2} \mathbf{v}_{n+1/2} \end{gathered}\end{split}\]

The velocity and displacement updates emanate from the central difference approximations to the acceleration \(\mathbf{a}_n\) and velocity \(\mathbf{v}_{n+1/2}\), respectively, giving the algorithm its name. The velocity measures that are utilized by the algorithm are shifted by a half step (said to be centered at the half-step), whereas accelerations and displacements are centered at the whole step. Fig. 12.3 graphically reveals the simplicity of the explicit time integration scheme.

Fig. 12.3 Graphical representation of the central difference time integrator.

As already mentioned, explicit finite element schemes are only conditionally stable, meaning that they only remain stable when the time increment \(\Delta t\) is less than some critical limit. This limit, sometimes called the Courant stability limit (see Reference [[5]]), can be shown to be as follows

(12.12)\[\Delta t \leq \frac{2}{\omega},\]

where \(\omega\) is the highest natural frequency in the mesh. An important necessary step in the central difference explicit time integrator is the estimation of this highest natural frequency in the discretized problem. Explicit dynamics problems frequently involve large deformations with potentially significant geometric, material, and contact nonlinearities, all of which can cause significant changes in the critical time step. Therefore, estimation of the critical time step must be made repeatedly throughout the problem simulation. It is thus important that the this calculation be as accurate and efficient as possible to make the most of the explicit method.

12.3.1. Element-based Critical Time Step Estimate

Stable time step estimates for explicit finite element methods are traditionally based on the conservative estimate of the frequency:

(12.13)\[\omega = 2 \left( \frac{c}{h}\right)_{\mathrm{max}} ,\]

where \(c\) and \(h\) are the sound speed and characteristic mesh size, respectively, associated with the element in the mesh having the largest ratio of these two quantities. Combining (12.12) and (12.13) we find that

(12.14)\[\Delta t^{\mathrm{max}} = \left( \frac{h}{c}\right)_{\mathrm{min}} .\]

In other words, the time step may be no larger than the amount of time required for a sound wave to traverse the element in the mesh having the smallest transit time. Such an estimation of the critical time step is based solely on element level calculations and is, in fact, part of the element internal force calculation. This is due in large part to the estimate of the sound speed of the material, which is of a dilatational wave. The accuracy of directly applying this condition is limited in practice due to the arbitrary finite element geometries in a typical mesh because the definition of characteristic length is somewhat of an art for distorted elements. Alternatively, the stability limit as reported in Reference [[3]] is related to the maximum global eigenvalue, \(\lambda_{\text{max}}\):

(12.15)\[\Delta t^2 = \frac{4}{\lambda_{\text{max}}}.\]

Because the maximum element eigenvalue is an upper bound on the maximum global eigenvalue (Reference [[6]]), we can compute an element-based stable time step estimate using

(12.16)\[\Delta t^E = \frac{2}{\left.\sqrt{\lambda}\right|_{\text{max over }e }}.\]

Details of how this element-based time step is calculated for different elements are covered in the chapter on element formulations.

12.3.2. Nodal-based Critical Time Step Estimate

A method is now described in which the maximum element modal stiffness are used to estimate a maximum nodal stiffness which, when combined with the lumped nodal mass, gives a sharper upper bound on the maximum global eigenvalue.

Let \(\lambda_{\text{max}}\) denote the largest eigenvalue of the generalized problem

(12.17)\[\left( \mathbf{K} - \lambda \mathbf{M} \right) \mathbf{u} = 0\]

and \(\mathbf{u}_{\text{max}}\) the eigenvalue corresponding to \(\lambda_{\text{max}}\). In (12.17), \(\mathbf{K}\) is the stiffness matrix and \(\mathbf{M}\) the diagonal, lumped mass matrix. The Rayleigh quotient for the maximum eigenvalue is

(12.18)\[\lambda _{\mathrm{max}} = \frac {\mathbf{u}_{\mathrm{max}}^T \mathbf{K} \mathbf{u}_{\mathrm{max}} } { \mathbf{u}_{\mathrm{max}}^T \mathbf{M} \mathbf{u}_{\mathrm{max}} } .\]

Noting that the numerator of (12.18) is twice the strain energy \(S\) of the system when deformed into the mode shape \(\mathbf{u}_{\mathrm{max}}\), we can write

(12.19)\[2S = \mathbf{u}_{\mathrm{max}}^T \mathbf{K} \mathbf{u}_{\mathrm{max}} = \sum_{e=1}^{n_e} \left( \mathbf{u}_{\mathrm{max}}^e \right)^T \mathbf{K}^e \left( \mathbf{u}_{\mathrm{max}}^e \right) .\]

We observe that the eigenvalue problem for the element stiffness matrix \(\mathbf{K}^e\) may be stated as

(12.20)\[\mathbf{K}^e \phi ^e = k^e \phi ^e .\]

Consequently,

(12.21)\[(\mathbf{u}^e )^T \mathbf{K}^e \mathbf{u}^e \leq k_{\mathrm{max}}^e (\mathbf{u}^e )^T \mathbf{u}^e\]

for all \(\mathbf{u}^e\) where \(k_{\mathrm{max}}^e\) is the maximum eigenvalue (so called modal stiffness) of the element stiffness matrix. From this result, we define a global stiffness matrix \(\hat{\mathbf{K}}\) assembled from the element stiffness matrices \(\hat{\mathbf{K}}^e\) defined as

(12.22)\[\hat{\mathbf{K}}^e = k_{\mathrm{max}}^e \mathbf{I}^e\]

where \(\mathbf{I}^e\) is an ndofe by ndofe identity matrix (ndofe is the number of degrees of freedom in the element). Based on (12.18), (12.20) and (12.21),

(12.23)\[2S \leq \sum_{e=1}^{n_e} \left( \mathbf{u}_{\mathrm{max}}^e \right)^T \hat{\mathbf{K}}^e \left( \mathbf{u}_{\mathrm{max}}^e \right)\]

leading to

(12.24)\[\lambda _{\mathrm{max}} = \frac {\mathbf{u}_{\mathrm{max}}^T \mathbf{K} \mathbf{u}_{\mathrm{max}} } { \mathbf{u}_{\mathrm{max}}^T \mathbf{M} \mathbf{u}_{\mathrm{max}} } \leq \frac{\mathbf{u}_{\mathrm{max}}^T \hat{\mathbf{K}} \mathbf{u}_{\mathrm{max}} } { \mathbf{u}_{\mathrm{max}}^T \mathbf{M} \mathbf{u}_{\mathrm{max}} } = \hat{\lambda}_{\mathrm{max}} .\]

Given the mode shape \(\mathbf{u}_{\mathrm{max}}\), the expression for \(\hat{\lambda}_{\mathrm{max}}\) is easily evaluated since both \(\hat{\mathbf{K}}\) and \(\mathbf{M}\) are diagonal. Methods for predicting this mode shape have been developed for specific ‘template’ geometries (Reference [[7]], but for general finite element geometries this remains impractical.

Rather than directly calculating \(\hat{\lambda}_{\mathrm{max}}\), we seek an upper bound. To this end, we define the ratio for every node \(I\) as

(12.25)\[\hat{\lambda}^I = \frac {\hat{K}^I }{M^I} ,\]

where \(\hat{K}^I\) and \(M^I\) are the diagonal elements in the \(I^{\text{th}}\) row of \(\hat{\mathbf{K}}\) and \(\mathbf{M}\), respectively. Without loss of generality, the ratios are ordered such that \(\hat{\lambda}^m \geq \hat{\lambda}^{m-1} \geq \dots \geq \hat{\lambda}^1\), in which case (12.24) can be written as

(12.26)\[\hat{\lambda}_{\mathrm{max}} = \frac { \sum_I {\mathbf{u}_{\mathrm{max}^I}^T \hat{\mathbf{K}}^I } } { \sum_I {\mathbf{u}_{\mathrm{max}^I}^T \mathbf{M}^I }} = \hat{\lambda}^m \left\{ \frac{ 1+ \left[ \left( ( u_{\mathrm{max}}^{m-1} )^2 M^{m-1}\right) / \left( ( u_{\mathrm{max}}^{m} )^2 M^{m}\right)\right] \frac{\hat{\lambda}^{m-1}}{\hat{\lambda}^{m}} + \dots} { 1+ \left[ \left( ( u_{\mathrm{max}}^{m-1} )^2 M^{m-1}\right) / \left( ( u_{\mathrm{max}}^{m} )^2 M^{m}\right)\right] + \dots} \right\}\]

Since all the ratios \(\hat{\lambda}^{m-1} / \hat{\lambda}^{m}\) are less than or equal to one, it follows immediately that

(12.27)\[\hat{\lambda}_{\mathrm{max}} \leq \hat{\lambda}^m = \left. \frac{\hat{K}^I}{M^I} \right|_{\text{max over }I} ,\]

in which \(M^I\) is the lumped mass at node \(I\), and \(\hat{K}^I\) is the assembly of the maximum element modal stiffness at node \(I\), that is

(12.28)\[\hat{K}^I = \sum_{e \in e^I} {k_{\text{max}}^e} ,\]

where \(e^I\) is the set of elements that are connected to node \(I\).

(12.15), (12.24), and (12.27) lead to a nodal-based stable time step estimate:

(12.29)\[\Delta t^N = \frac{2} {\left. \sqrt{\frac{\hat{K}^I}{M^I} } \right|_{\text{max over }e } }.\]

Now we show that the nodal-based stable time step estimate is always greater than or equal to the element-based estimate. Following a similar procedure outlined in (12.26), we can write

(12.30)\[\hat{\lambda}^I = \frac{\hat{K}^I}{M^I} = \frac{\sum_{e \in e^I} k_{\text{max}}^e}{\sum_{e \in e^I} m^e} = \frac{ \frac{k_{\text{max}}^1}{m^1} + (m^2/m^1)\frac{k_{\text{max}}^2}{m^2} + \dots }{1+(m^2/m^1)+ \dots} \leq \frac{k_{\text{max}}^1}{m^1} ,\]

where the element eigenvalues \({k_{\text{max}}^e}/{m^e}\) are arranged in descending order, \({k_{\text{max}}^1}/{m^1} \geq {k_{\text{max}}^2}/{m^2} \geq \dots\). Thus the nodal-based estimate of the maximum eigenvalue at node \(I\) is bounded by the largest of all element eigenvalues connected to node \(I\). It follows from (12.30) that

(12.31)\[\hat{\lambda}^m = \left. \frac{\hat{K}^I}{M^I}\right|_{\text{max over }I} \leq \left. \frac{k_{\text{max}}^e}{m^e}\right|_{\text{max over }e \in e^I} \left. \frac{}{}\right|_{\text{max over }I} = \left. \frac{k_{\text{max}}^e}{m^e}\right|_{\text{max over }e} = \left. \lambda_{\text{max}}^e \right|_{\text{max over }e}\]

Since \(\hat{\lambda}^m \leq \left. \lambda_{\text{max}}^e \right|_{\text{max over }e}\), it follows directly from (12.15) that the nodal-based estimate is always greater than or equal to the element-based estimate.

The cost of the nodal-based estimate calculation includes the element eigenvalue analysis (which must be done in the case of the element based calculation) plus the cost of an assembly procedure every time step. (12.29) must be evaluated at each node as opposed to evaluating (12.16) for every element.

12.3.3. Lanczos-based Critical Time Step Estimate

The paper [[8]], which is reproduced here, demonstrates the cost-effective use of the Lanczos method for estimating the critical time step in an explicit, transient dynamics code. The Lanczos method can give a significantly larger estimate for the critical time-step than an element-based method (the typical scheme). However, the Lanczos method represents a more expensive method for calculating a critical time-step than element-based methods. Our paper shows how the additional cost of the Lanczos method can be amortized over a number of time steps and lead to an overall decrease in run-time for an explicit, transient dynamics code. We present an adaptive hybrid scheme that synthesizes the Lanczos-based and element-based estimates and allows us to run near the critical time-step estimate provided by the Lanczos method.

12.3.3.1. Introduction

Codes using explicit time integration techniques are important for simulating transient dynamics problems involving large deformation of solids with various nonlinear effects (contact, nonlinear materials, element death, etc.). The second order central difference operator used in explicit codes is stable if the time step is no larger than the critical time step. For most problems in solid mechanics, the critical time step is extremely small and the number of time steps required for a typical analysis is quite large. Therefore, the accurate, efficient, and reliable calculation of the critical time step is of fundamental importance.

The element-based method [[9]] is an efficient method for producing a critical time step estimate at every time step. However, it can produce a conservative estimate for the critical time step in many cases. The Lanczos [[10]] method is a reliable procedure for producing a time step that is the theoretical maximum value for a structure and is usually much better than the element-based estimate. The cost of obtaining a Lanczos based estimate will not offset the cost benefit of the increased value for the critical time step. Therefore, it is not feasible to call the Lanczos method at every explicit dynamics time step. In this paper we outline a cost-effective method for utilizing the Lanczos method (together with an element-based scheme) for the critical time step estimation.

Benson [[11]] investigates estimating the critical time step by using the power iteration. Parlett [[12]] presents analysis comparing the Lanczos method and power iteration. The Lanczos method provides a more rapid approximation, in terms of matrix-vector products, relative to the power iteration for approximating the largest eigenvalue as the relative separation of the largest eigenvalue decreases. Hence, we can expect the Lanczos method to require less matrix-vector products to approximate the critical time step to a specified tolerance. We also remark that in contrast to our paper, Benson [[11]] does not present a scheme that addresses two crucial issues when using the power iteration (or Lanczos method) for estimating the critical time step.

Two crucial issues must be addressed when using the Lanczos method to estimate the critical time step. First, the Lanczos-based time step estimate must be used for two to three times the number of explicit time integration steps required to recover the cost of the Lanczos method if we are to see a noticeable reduction in overall computation times for a problem. (We explore the cost of the Lanczos method in terms of internal force calculations in later sections.) Second, the Lanczos method provides an overestimate of the critical time step, and so we need an effective scheme to scale back the Lanczos-based critical time step estimate. We address both these issues and present an adaptive hybrid scheme that synthesizes the Lanczos-based and element-based estimates and allows us to run near the critical time-step estimate provided by the Lanczos method.

We also remark that in addition to the increased efficiency that can be achieved with the Lanczos-based time step, we also have the added benefit of increased accuracy. For explicit transient dynamic codes, using a time step as close as possible to the critical time step [[3]] gives the most accurate answer. Reducing the time step in an explicit transient dynamics code actually increases the error.

Our paper is organized as follows. Section Section 12.3.3.2 discusses the critical time step and motivates a Lanczos-based estimate. The Lanczos iteration and method are briefly introduced in section Section 12.3.3.3. A cost benefit analysis of the element-based and Lanczos-based approximations to the critical time is considered in section Section 12.3.3.4. A practical implementation within an explicit dynamics code is the subject of section Section 12.3.3.5. Several numerical examples are presented in section Section 12.3.3.6, and we provide our conclusions in section Section 12.3.3.7.

12.3.3.2. Critical time step

Let \(\mathbf{K}\) and \(\mathbf{M}\) be the stiffness and mass matrices arising in an explicit dynamics simulation so that \(\mathbf{M}\) is a diagonal matrix due to mass lumping. The critical time step for second order central time differencing is bounded from above by \(2\omega^{-2}_{\max}\) where \(\omega^2_{\max}\) is the largest eigenvalue of the generalized eigenvalue problem

(12.32)\[\mathbf{K} \mathbf{u} = \mathbf{M} \mathbf{u}\omega^2_{\max}, \quad \left( \mathbf{K}, \mathbf{M} \in \mathbb{R}^{n\times n} \right),\]

where we assume that \(\omega^2_{\max}\) is positive. An inexpensive [[13]] upper bound to \(\omega^2 ` is given by the maximum element eigenvalue :math:\)omega^2_{max,e}` over all the element eigenvalue problems

(12.33)\[\mathbf{K}^e \mathbf{u}^e = \mathbf{M}^e \mathbf{u}^e\omega^2_{e}, \quad \left( \mathbf{K}^e, \mathbf{M}^e \in \mathbb{R}^{n^e\times n^e} \right),\]

where \(n^e \ll n\). Therefore, \(\omega^{-2}_{\max,e} \leq \omega^{-2}_{\max}\) and we have a lower bound for the critical time step. The maximal element eigenvalue is typically computed analytically [[9]] for the finite elements that are typically used in transient dynamics.

The Lanczos method rapidly provides a lower bound \(\omega^2_{\max,L}\) to \(\omega^2_{\max}\) so that

(12.34)\[\omega^{-2}_{\max,e} \leq \omega^{-2}_{\max} \leq \omega^{-2}_{\max,L}.\]

In fact, the Lanczos iteration is sharp so that \(\omega^{-2}_{\max} \lessapprox \omega^{-2}_{\max,L}\) so that with care, an excellent approximation to the critical time step is computed for a modest cost. This approximation may be dramatically superior to the standard element based estimate. The details of a careful use of the Lanczos-based estimate is the subject of section Section 12.3.3.5.

12.3.3.3. Lanczos iteration

The Lanczos reduction rapidly provides approximations to the maximum and minimum eigenvalues of a symmetric \(\mathbf{A}\in \mathbb{R}^{n\times n}\), in particular the largest in magnitude eigenvalue. Suppose that

(12.35)\[\mathbf{A} \mathbf{Q}_j = \mathbf{Q}_j \mathbf{T}_j + \mathbf{f}_j \mathbf{e}^T_j,\]

is a Lanczos reduction of length \(j\) where \(\mathbf{f}_j \in\mathbb{R}^{n}\), and \(\mathbf{e}_j \in \mathbb{R}^{j}\) contains column \(j\) of the identity matrix \(\mathbf{I}_n \in \mathbb{R}^{n\times n}\). If we denote

\[\begin{split}\mathbf{T}_j = \left( \begin{array}{cccc} \alpha_1 & \beta_2 & \cdots & 0 \\ \beta_2 & \alpha_2 & \cdots & 0 \\ \vdots & & \ddots & \beta_j \\ 0 & \cdots & \beta_j & \alpha_j \\ \end{array} \right), \quad \alpha_i, \beta_i \in \mathbb{R}\end{split}\]

and

\[\begin{split}\mathbf{Q}_j = \left( \begin{array}{cccc} \mathbf{q}_1 & \mathbf{q}_2 & \cdots & \mathbf{q}_j \\ \end{array} \right), \quad \mathbf{q}_i \in \mathbb{R}^{n }\end{split}\]

then the familiar Lanczos three-term recurrence is recovered by equating column \(j\) of ((12.35)) to obtain

(12.36)\[\mathbf{f}_j = \mathbf{A} \mathbf{q}_j-\mathbf{q}_j\alpha_j-\mathbf{q}_{j-1}\beta_{j-1}^T.\]

Furthermore, because of the orthonormality of \(\mathbf{Q}_j\), we have

\[\begin{split}\alpha_j &=\mathbf{q}_j^T \mathbf{A}\mathbf{q}_j, \\ \mathbf{q}_{j+1}\beta_{j+1} &= \mathbf{f}_j, \\ \mathbf{q}_i^T \mathbf{f}_j & = 0, \quad i=1,\ldots,j\end{split}\]

and so \(\mathbf{q}_{j+1} = \mathbf{f}_j \beta_{j+1}^{-1}\), where we assume that \(\beta_{j+1}\) is non-zero. We define a Lanczos iteration to be that computing \(\mathbf{A} \mathbf{q}_j, \alpha_j,\beta_{j+1}\), and \(\mathbf{f}_j\). We define the Lanczos method that of computing \(m\) iterations and computing the largest in magnitude eigenvalue of \(\mathbf{T}_m\).

The largest eigenvalue of the symmetric tridiagonal matrix \(\mathbf{T}_j\) approximates the largest in magnitude eigenvalue of \(\mathbf{A}\). We can determine the quality of the approximation produced by an eigenpair of \(\mathbf{T}_j\). If we post multiply ((12.35)) by \(\mathbf{s}\) where \(\mathbf{T}_j\mathbf{s} = \mathbf{s} \theta\) (and \(\|\mathbf{s}\|=1\)), then

(12.37)\[\mathbf{A} (\mathbf{Q}_j \mathbf{s}) - (\mathbf{Q}_j \mathbf{s})\theta = \mathbf{f}_j (\mathbf{e}^T_j\mathbf{s}).\]

In words, the residual of the approximate eigenpair \((\mathbf{Q}_j\mathbf{s},\theta)\) is proportional to \(\mathbf{f}_j\) (note that \(\mathbf{e}^T_j\mathbf{s}\) is notation for the last component of \(\mathbf{s}\)). The implication is that we can easily monitor the quality of the approximation produced by the Lanczos method. If \(\theta\) is the largest in magnitude eigenvalue of \(\mathbf{T}_j\), then \(\theta \leq \omega_{\max}^2 \leq \|\mathbf{f}_j\|_2 \, |\mathbf{e}^T_j\mathbf{s}| + \theta\) (see [[12]] for a discussion). Hence,

(12.38)\[\frac{1}{\|\mathbf{f}_j\|_2 \, |\mathbf{e}^T_j\mathbf{s}| + \theta} \leq \omega_{\max}^{-2} \leq \frac{1}{\theta}.\]

We also remark that the norm of the residual is a non-increasing function of \(j\); again see [[12]].

The Lanczos iteration is adapted for computing the largest eigenvalues of ((12.32)) by replacing \(\mathbf{A}\) with \(\mathbf{M}^{-1}\mathbf{K}\) and computing an \(\mathbf{M}\)-orthonormal \(\mathbf{Q}_j\). This orthonormality is needed so that \(\mathbf{M}^{-1}\mathbf{K}\) is symmetric in the inner product induced by \(\mathbf{M}\). See [[14], [12]] for further discussion and implementations.

The cost of a careful implementation of a Lanczos iteration, :math:` j>1`, is one matrix-vector product with \(\mathbf{K}\) and \(\mathbf{M}^{-1}\), and two vector products and vector subtractions. Within an explicit dynamics code, the cost of computing a Lanczos vector is approximately the cost of an internal force calculation, represented by the matrix-vector product \(\mathbf{K}\mathbf{q}_j\). Therefore, we approximate the cost of computing the Lanczos-based time step estimate as

(12.39)\[m \tau\]

where \(m\) denotes the number of Lanczos iterations and \(\tau\) represents the CPU (central processor unit) time needed for an element-based explicit dynamics time integration step.

The Lanczos method only requires knowledge of \(\mathbf{K}\) via its application on a vector. If internal force calculations are used for the needed matrix-vector products, the Lanczos vectors \(\mathbf{q}_j\) are scaled so that they represent velocities associated with small strain. When these scaled vectors are sent to the internal force calculation, the internal force calculation becomes a matrix-vector product with a (constant) tangent stiffness matrix \(\mathbf{K}_T\).

12.3.3.4. Cost-Benefit Analysis

This section provides a simple model for assessing the cost of using the Lanczos method for computing an estimate of the critical time step. We assume that Lanczos-based time step is valid for \(n_{L}\) time integration steps. We address the important issue of the adapting the time step when we present the details for practical use of the Lanczos method in a subsequent section.

Denote by \(\Delta t_{L}\) and \(\Delta t_{e}\) the time steps estimate of the critical time step computed by the Lanczos and element-based methods, where the ratio \(\rho\) of \(\Delta t_{L}\) to \(\Delta t_{e}\) is at least as large as one because of ((12.34)). After \(n_{L}\) time steps, the dynamics simulation is advanced in time \(n_{L}\Delta t_{L}\). Let \(n_e\) be the number of element-based time steps so that \(n_e\Delta t_{e} \leq n_{L}\Delta t_{L} < (n_e+1)\Delta t_{e}.\) In terms of \(\rho\), we have the relationship

(12.40)\[n_e \leq \rho n_{L} < n_e+1,\]

so bounding the number of Lanczos-based explicit integration steps in terms of \(\rho\) and the number of element-based integration steps.

Let us examine the computational costs in terms of CPU time in performing the above \(n_{L}\) and \(n_{e}\) integration steps. Denote by \(\tau\) the CPU time for an element-based time integration step and assume that it is dominated by the cost of an internal force calculation. Using ((12.40)), the CPU time of \(n_{L}\) time integration steps is

(12.41)\[(n_{L} + m) \tau,\]

and the CPU time of \(n_{e}\) time integration steps is \(n_{e}\tau\). Equating these two CPU times, determines when the cost of both approaches is equivalent and results in the relationship

(12.42)\[\hat{n}_{e}= m + \hat{n}_L.\]

Using ((12.42)) within ((12.40)) gives

(12.43)\[\frac{m}{\rho-1} \leq \hat{n}_L < \frac{m+1}{\rho-1}\]

so bounding the minimum number of Lanczos-based time integration steps in terms of the number of Lanczos iterations and \(\rho\) so that the cost of the computing the Lanczos-based time step is amortized.

Our cost benefit analysis provides theb reak-even point at which the Lanczos method becomes cost-effective by overcoming the associated overhead. For example, let \(\rho = 1.25\) and \(m = 20\) so that \(\hat{n}_{L}\) is bounded from below by \(80\), and by ((12.41)) \(\hat{n}_{e} = 100\). Hence, the time integration with the Lanczos-based and element-based estimates of the critical time step give the same simulation time for the same CPU time. If we use the Lanczos-based time step \(\Delta t_{L}\) for more than \(80\) time integration steps, then the Lanczos-based approach is cost-effective.

A Lanczos-based critical time estimate is cost effective if \(m\) is small and \(\rho\) is not close to one. The size of \(m\) is dependent upon the ability of the Lanczos method to rapidly provide an accurate approximation to \(\omega^2_{\max}.\) If \(\rho\) approaches one, then the Lanczos-based critical time step approaches the element-based critical time step, implying that \(\hat{n}_L\) must increase to offset the cost of the \(m\) Lanczos iterations. Section Section 12.3.3.6 demonstrate that \(m\) is small and that \(\rho\) is not close to one for realistic problems.

Our section ends by considering the additional cost involved with contact. The addition of contact to an analysis can add computational costs to a time step that are as large as or larger than the internal force calculations. Therefore, for an analysis with contact, running at a larger time step than the element-based estimate can have an even greater impact on reducing CPU time for an analysis.

The above analysis is easily extended to the case where we have contact. If the CPU time of contact over a time step is some multiple \(\gamma\) of \(\tau\), then in analogy to ((12.42)) and ((12.43)), we have

\[(1 + \gamma) \hat{n}_{e} = m+ (1 + \gamma) \hat{n}_{L},\]

and

\[\frac{m}{(\rho-1)(1+\gamma)} \leq \hat{n}_L < \frac{m+1+\gamma}{(\rho-1)(1+\gamma)}\]

Again, for example, let \(\rho = 1.25\) and \(m = 20\) and assume the computational cost of contact calculations is the same as an internal force calculation so that \(\gamma=1\). Hence, the break-even point is \(\hat{n}_{L} = 40\) and \(\hat{n}_{e} = 50\). The additional cost of the contact calculations within the time integration reduces the break-even point over that with no contact (\(\gamma=0\)).

12.3.3.5. Using the Lanczos-based estimate

The previous section shows how the repeated use of a Lanczos-based time step estimate could be cost-effective within an explicit transient dynamics simulation. This section presents an adaptive scheme that combines the Lanczos-based estimate with an element-based estimate of the critical times-step over a number of explicit time integration steps.

Section 12.3.3.2 explained that the Lanczos method provides an approximation to the maximum eigenvalue of (12.32) from below so overestimating the critical time step. Therefore, we scale back the Lanczos-based time. The scheme to determine a scaled-back value employs the element-based time step estimate. Again, let \(\Delta t_{L}\) and \(\Delta t_{e}\) be the time steps computed by the Lanczos and element-based methods. The scaled back estimate for the critical time step, \(\Delta t_{s}\), is computed from the equation

\[\Delta t_{s} = \Delta t_{e} + f_{s} ( \Delta t_{L} - \Delta t_{e} ),\]

where \(f_{s}\) is a scale factor. (The value for \(f_{s}\) ranges from \(0.9\) to \(0.95\) for our problems—a rigorous estimate can be made by using ((12.38)).) This value of \(f_{s}\) results in \(\Delta t_{s}\) close to and slightly less than the critical time step. Once \(\Delta t_{s}\) is determined, the ratio

\[t_{r} = \frac{\Delta t_{s}}{\Delta t_{e}}\]

is computed. This ratio is then used to scale subsequent element-based estimates for the critical time step. If \(\Delta t_{e(n)}\) is the \(n^{th}\) element-based time step after the time step where the Lanczos method is computed, then the \(n^{th}\) time step computed is

\[\Delta t_{(n)} = t_{r} \Delta t_{e(n)}.\]

The ratio \(t_{r}\) is used until the next call to the Lanczos method. The next call to the Lanczos method is controlled by one of two mechanisms. First, the user can set the frequency with which the Lanczos method is called. The user can set a parameter so that the Lanczos method is called only once every \(n\) time steps. This number remains fixed throughout an analysis. Second, the user can control when the Lanczos method is called based on changes in the element-based time step. For this second method, the change in the element-based critical time step estimate is tracked. At the \(n^{th}\) step after the call to the Lanczos iteration, the element-based time step is \(\Delta t_{e(n)}\). If the value

\[\frac{|\Delta t_{e(n)} - \Delta t_{e} |}{\Delta t_{e}}\]

is greater than some limit set by the user, then the Lanczos method is called. If there is a small, monotonic change in the element-based time step over a large number of time integration steps, this second mechanism will result in the Lanczos method being computed. If there is a large, monotonic change in the element-based critical time step over a few time steps, the Lanczos method will also be called.

These two mechanisms for calling the Lanczos method may be combined resulting in an adaptive scheme for estimating the critical time step during an explicit transient dynamics simulation. For example, suppose the second mechanism, the mechanism based on a change in the element-based time step, results in a call to the Lanczos method. This resets the counter for the first mechanism, the mechanism using a set number of time steps between calls to the Lanczos iteration.

12.3.3.6. Numerical experiments

This method for reusing a Lanczos-based time step estimate has been implemented in Presto [[15]], and employed within a number of explicit dynamics simulations. We discuss several of these examples.

Example one: The Lanczos method has been used to obtain a critical time step estimate for a cubic block consisting solely of cubic elements—a \(10 \times 10 \times 10\) mesh of eight-node hexahedral elements. We know that, for a cubic eight-node hexahedral element, the element-based estimate is conservative by a factor of \(1/\sqrt{3}\). The Lanczos method yields a critical time estimate for this mesh that is \(\rho=\sqrt{3}\) (approximately 1.732) times larger than the element-based estimate. This is done by using 20 Lanczos vectors.

Example two: Critical time step estimates were made for two mechanical systems. The systems consisted of cylindrical metal cans containing a variety of components. Some of these components have relatively simple geometries, while other components have complex shapes. A number of the components with complex shapes are a foam material used to absorb impact loads. One component was modeled with approximately 250,000 degrees of freedom, and the other one was modeled with approximately 350,000 degrees of freedom. For both of these models, a good estimate for the maximum eigenvalue was obtained with the Lanczos method by computing only twenty Lanczos vectors. For the model with 250,000 degrees of freedom, an actual analysis was run. The value for \(\rho\) for this problem was 1.83. The break-even point for this case (\(n_{L} = 20\) and \(\rho = 1.83\)) is \(n_{e} = 45\). It was possible to use the same scale factor for 1700 time steps for this analysis, which is well above the break-even point. The extended use of the Lanczos based estimate reduced the computation cost by over 56%.

Example three: A study of a large-scale model involving 1.7 million nodes (5.1 million degrees of freedom) showed that only 45 Lanczos vectors were required to obtain a good estimate of the maximum eigenvalue. The value of \(\rho\) for this problems was 1.2. Use of this Lanczos based estimated for this problem would be extremely cost-effective.

12.3.3.7. Conclusions

The Lanczos method is cost-effective for estimating the critical time step in an explicit, transient dynamics code. The Lanczos method can give a significantly larger estimate for the critical time-step than an element-based method (the typical scheme). The adaptive hybrid scheme synthesizes the Lanczos-based and element-based estimates and allows us to run near the critical time-step estimate provided by the Lanczos method.

Not all problems will lend themselves reuse of one Lanczos-based estimate for thousands of time steps. However, if it is possible to use the Lanczos-based estimate for two to three times the number of time steps required for break-even, we begin to see a noticeable reduction in the total CPU time required for a problem.

In addition, to the increased efficiency we can achieve with the Lanczos iteration, we also have the added benefit of increased accuracy. For explicit transient dynamic codes, using a time step as close as possible to the critical time gives the most accurate answer. Reducing the time step in an explicit transient dynamics code actually increases the error.

12.4. Implicit Finite Element Methods

To introduce the concept of an implicit time finite element method, we examine the trapezoidal rule, which is simply the member of the Newmark family obtained by setting \(\alpha=0\), \(\beta=1/4\), and \(\gamma=1/2\). Substitution of these values into (12.6) yields

(12.44)\[\begin{split}\begin{gathered} \mathbf{M} \mathbf{a}_{n+1} + \mathbf{F}^{\mathrm{int}}(\mathbf{d}_{n+1}) = \mathbf{F}^{\mathrm{ext}}(\mathbf{d}_{n+1}) \\ \mathbf{d}_{n+1} = \mathbf{d}_n + \Delta t \mathbf{v}_n + \frac{\Delta t^2}{4} \left[ \mathbf{a}_n + \mathbf{a}_{n+1} \right] \\ \mathbf{v}_{n+1} = \mathbf{v}_n + \frac{\Delta t}{2} \left[ \mathbf{a}_n + \mathbf{a}_{n+1} \right] . \end{gathered}\end{split}\]

Insight into this method can be obtained by combining the first two equations in (12.44) and solving for \(\mathbf{d}_{n+1}\) to get

(12.45)\[\begin{split}\begin{gathered} \frac{4}{\Delta t^2} \mathbf{M} \mathbf{d}_{n+1} + \mathbf{F}^{\mathrm{int}}(\mathbf{d}_{n+1}) = \mathbf{F}^{\mathrm{ext}}(\mathbf{d}_{n+1}) + \mathbf{M} \left( \mathbf{a}_n + \Delta t \mathbf{v}_n + \frac{4}{\Delta t^2} \mathbf{d}_n \right) \\ \mathbf{a}_{n+1} = \frac{4}{\Delta t^2} \left( \mathbf{d}_{n+1} - \mathbf{d}_n \right) -\frac{4}{\Delta t} \mathbf{v}_{n} - \mathbf{a}_n \\ \mathbf{v}_{n+1} = \mathbf{v}_n + \frac{\Delta t}{2} \left[ \mathbf{a}_n + \mathbf{a}_{n+1} \right] \end{gathered}\end{split}\]

Solving the first equation is the most expensive procedure involved in updating the solution from \(t_n\) to \(t_{n+1}\). This equation is not only fully coupled, but also non-linear in general due to the internal force vector.

Note that we can write the first equation of (12.45) in terms of a dynamic incremental residual \(\mathbf{r}_{n+1}\) via

(12.46)\[\mathbf{r} ( \mathbf{d}_{n+1}) = \left[ \mathbf{F}^{\mathrm{ext}}(\mathbf{d}_{n+1}) + \mathbf{M} \left( \mathbf{a}_n + \Delta t \mathbf{v}_n + \frac{4}{\Delta t^2} \mathbf{d}_n \right) - \left( \frac{4}{\Delta t^2} \mathbf{M} \mathbf{d}_{n+1} + \mathbf{F}^{\mathrm{int}}(\mathbf{d}_{n+1}) \right) \right] = 0\]

This system has the same form as (11.10), which suggests that the same sort of nonlinear solution strategies are needed for implicit dynamic calculations as in quasistatics (Section 11). Equation solving is the topic of the next chapter, where we will discuss at some length the techniques used to solve (11.9) and (12.46) in Sierra/SolidMechanics, particularly for parallel computing.

---

[1]

H.M. Hilber, T.J.R. Hughes, and R.L. Taylor. Improved numerical dissipation for time integration algorithms in structural dynamics. Earthquake Engineering and Structural Dynamics, 5:283–292, 1977.

[2]

N.M. Newmark. A method of computation for structural dynamics. Journal of Engineering Mechanics, ASCE, 85(EM3):67–94, 1959.

[3] (1,2,3,4)

R.D. Krieg and S.W. Key. Transient shell response by numerical time integration. International Journal for Numerical Methods in Engineering, 7:273–286, 1973.

[4]

R.W. Hamming. Numerical Methods for Scientists and Engineers. Dover, 1973.

[5]

R. Courant, K. Friedrichs, and H. Lewy. Über die partiellen differenzengleichungen der mathematischen physik. Mathematische Annalen, 100:32–74, 1928.

[6]

D.P. Flanagan and T. Belytschko. Simultaneous relaxation in structural dynamics. Journal of the Engineering Mechanics Division, ASCE, 107:1039–1055, 1981.

[7]

G.L. Goudreau. Evaluation of discrete methods for the linear dynamic response of elastic and viscoelastic solids. Technical Report 69–15, U.C. Berkeley, Department of Civil Engineering, 1970.

[8]

J.R. Koteras and R.B. Lehoucq. Estimating the critical time-step in explicit dynamics using the Lanczos method. International Journal for Numerical Methods in Engineering, 69(13):2780–2788, 2006.

[9] (1,2)

D.P. Flanagan and T. Belytschko. Eigenvalues and stable time steps for the uniform strain hexahedron and quadrilateral. ASME Journal of Applied Mechanics, 51:35–40, 1984.

[10]

C. Lanczos. An iteration method for the solution of the eigenvalue problem of linear differential and integral operators. J. Research of the National Bureau of Standards, 45(4):255–282, October 1950. Research Paper 2133.

[11] (1,2)

D.J. Benson. Stable time step estimation for multi-material eulerian hydrocodes. Computer Methods in Applied Mechanics and Engineering, 167:191–205, 1998.

[12] (1,2,3,4)

B.N. Parlett. The symmetric eigenvalue problem. Number 20 in Classic in Applied Mathematics. SIAM, 1998. Reprint of “The symmetric eigenvalue problem”, Prentice-Hall, 1980.

[13]

B.M. Irons. Applications of a theorem on eigenvalues to finite element problems. Technical Report CR/132/70, University of Wales, Department of Civil Engineering, Swansea, U.K., 1970.

[14]

B. Nour-Omid. The Lanczos algorithm for solution of large generalized eigenproblems. In The Finite Element Method: Linear static and dynamic finite element analysis,Thomas J.R. Hughes, pages 582–600. Dover, 2000.

[15]

J.R. Koteras and A. S. Gullerud. Presto user’s guide version 1.05. Technical Report SAND2003-1089, Sandia National Laboratories, Albuquerque, NM, April 2003.