# 8.3. ECM Theory

Explicit transient dynamics is a well-established capability for modeling large deformations of structures. It is common practice in explicit transient dynamics to seek a balance between computational efficiency and accuracy. Mass scaling [[1]] has traditionally been used as an approach to increase the critical time step limit associated with the central difference time integrator. Unfortunately, this has an undesirable side effect of mass damping dynamic modal response over the entire frequency spectrum. To mitigate this effect, methods have been developed in which the damping is proportional to the frequency [[2]]. In Sierra/SM the Explicit Control Modes algorithm performs an efficient modal decomposition of the frequency spectrum, allowing mass damping only on the high frequency modes. Examples will be presented that demonstrate that this approach yields accurate low frequency response, while often using larger time steps due to the mass scaling the high frequency response.

## 8.3.1. Introduction

Finite element analysis of transient dynamic problems is a production capability in many application areas. In these analyses an important question to be addressed by the analyst is the choice of using an explicit or implicit time integrator. It is well understood that the central difference explicit time integrator is efficient per time step but is restricted to taking relatively small, critical (or stable) time steps [[3], [4]]. An implicit time integrator, specifically the Hilber–Hughes–Taylor (HHT) time integrator [[5]], with the proper choice of parameters has no such stability limit allowing larger time steps but produces a system of equations that need to be solved every time step. As noted in [[5]], a consequence of choosing a large time step for the implicit time integrator is that it produces numerical damping in all frequencies, but predominantly in the highest frequencies. How much damping and in what modes depends on the particular problem. Thus, the question of what time integrator to use is much more than one of efficiency. Certainly, one must know well the class of problems to be solved when making this choice.

Simulations requiring necessarily finer discretizations to accurately represent modal stiffness and resolve details in the stress field are amenable to Explicit Control Modes. For the explicit time integrator, this imposes a critical time step restriction that can be limiting. However for some—possibly many—analyses the structural response is in the lower frequency spectrum, i.e., the influence of the modal content inherent in fine mesh on the low-frequency dynamics is of interest and not necessarily the high-frequency dynamics themselves. More precisely, spatial resolution as opposed to temporal resolution is needed for many problems (this premise is one that we intend to support in the examples).

It seems appropriate, then, to consider an algorithmic approach that can improve the stability limit of the explicit time integrator. Most importantly, we recognize that this approach needs to be accurate for low-mode response and competitive with implicit dynamics.

## 8.3.2. Modal Decomposition Approach

The objective of this algorithm is to modally decompose the dynamics (in the context of an explicit transient dynamics time integrator) into low-frequency and high frequency response. Having this decomposition may provide options such as integrating the low-frequency modes with explicit time integration and the high-frequency modes with an implicit time integrator.

The decomposition algorithm is based on applying multigrid concepts within an explicit central difference time integrator. We will limit the algorithm to consider only one addition level of coarsening. Thus, in addition to the fine mesh or reference mesh, we introduce a set of coarse basis functions that will describe the low-mode response.

The vector of external nodal forces on the fine mesh is \(f^{ext}_{fm}\). Also the vector of internal nodal forces on the fine mesh \(f^{int}_{fm}\) is a obtained from the divergence of the stress. In this work we assume that there is no contact, in which case the nodal residual force is

\[r = f^{ext}_{fm} - f^{int}_{fm}.\]

Let \(M\) denote the diagonal, lumped [[6]], mass matrix for the fine mesh, and let \(x\) denote the vector of displacements on the fine mesh. Expressed in terms of the nodal displacements, the dynamic equilibrium equations is

(8.1)\[M \ddot{x} = r.\]

Let \(\Phi\) denote an interpolation (prolongation) matrix associated with a coarse space of functions. The number of rows in \(\Phi\) equals the number of rows in \(x\) or \(r\), while the number of columns in is typically smaller. The matrix can be obtained from either a coarse finite element mesh or by using an algebraic approach [[7][8]]. Given \(\Phi\), the acceleration on the fine mesh can be written as

(8.2)\[\ddot{x} = \Phi \ddot{q} + \ddot{x}_{hf}\]

where \(q\) is a vector of generalized displacements associated with the low frequency part of the response, and \(x_{hf}\) is a vector of displacements associated with the high frequency part of the response.

The task now is to derive the equation to accomplish this decomposition making use only of the residual vector, \(r\), and mass matrix, \(M\), on the fine mesh, recognizing that there are no properties on the coarse mesh in the usual finite element sense. As in the multigrid method, the modal stiffness of low-mode response and the corresponding mass matrix is obtained using a restriction operator of properties/quantities from the fine mesh.

The low and high frequencies are decoupled by imposing the \(M\) orthogonality,

\[\ddot{q}^T \Phi^T M \ddot{x}_{hf}=0,\]

of the high and low frequency displacements. The orthogonality condition holds for all \(\ddot{q}\) if and only if

(8.3)\[\Phi^T M \ddot{x}_{hf}=0.\]

(8.3) implies that the high frequency part of the residual is orthogonal to the coarse space spanned by the columns of \(\Phi\). The coarse mesh mass matrix is given by

\[M_c = \Phi^T M \Phi.\]

Substitution of (8.2) into (8.1), pre-multiplying by \(\Phi^T\) and making use of (8.3) leads to the low frequency equilibrium condition

(8.4)\[M_c \ddot{q} = \Phi^T r.\]

This way of obtaining a coarse system from \(\Phi\) is called Galerkin coarsening [[8]]. For reference, the coarse grid stiffness matrix \(K_c\) corresponding to the fine mesh tangent stiffness matrix \(K\) is given by

\[K_c = \Phi^T K \Phi.\]

Next the high frequency equilibrium equation is derived. Solving for \(\ddot{q}\) in (8.4) gives

(8.5)\[\ddot{q} =M^{-1}_c \Phi^T r.\]

From (8.1), (8.2), and (8.5) determine the high frequency equilibrium condition

(8.6)\[M \ddot{x}_{hf} = r - M \Phi M^{-1}_c \Phi^T r.\]

At this point no approximations have been made. To sum up, substituting (8.5) and (8.6) into (8.2) leads to

(8.7)\[\ddot{x} = \underbrace{ \Phi M^{-1}_c \Phi^T r}_\text{low frequency}+ \underbrace{ M^{-1} \left( r - M \Phi M^{-1}_c \Phi^T r\right)}_\text{high frequency}\]

The lumped mass matrix is required to obtain the most accurate approximation properties for the explicit central difference time integrator [[9], [4]]. Thus, given that we are integrating the low-frequency response with central difference, a lumped representation is needed. It is unclear that the argument for finite elements and the fine mesh extend to the Galerkin coarse problem. The lumping is done by applying the restriction operator to the diagonal lumped fine mesh mass matrix,

\[M_c = \Phi^T M \Phi.\]

## 8.3.3. Explicit-Explicit Partitioning

First we consider explicit time integration for the low-frequency modes. The critical time step for integrating these modes is constructed, again via projection of nodal quantities on the fine mesh. If \(\Delta t_{cr}^{cm}\) denotes the critical time step for the coarse mesh, then a node-based estimate [[3]] is given by

> \[\Delta t_{cr}^{cm} = \min_{\text{coarse}~ \text{nodes}} 2 \sqrt{ \frac{ \Phi^T M }{ \Phi^T K^{max} } }\]

where \(K^{max}\) is a vector that contains the maximum modal stiffness for each node of the fine mesh. Details of the calculation of the maximum modal stiffness can be found in [[3]].

Next, we wish to make use of the assumption that the high-frequency dynamics are negligible. The accelerations represented by the second term in (8.7) correspond to those high-frequency modes. The idea is to replace \(M^{-1}\) in the second term of (8.7) by \(\tilde{M}^{-1}\), where

(8.8)\[\tilde{M}= M \alpha,\]

for a diagonal matrix \(\alpha\) that contains a scale factor for each node of the fine mesh. These scale factors are greater than \(1\) wherever the nodal based time step at a fine mesh node is smaller than the critical time step on the coarse mesh.

(8.9)\[\begin{split}\alpha_i = \begin{cases} \frac {K^{max}_i}{4M_i} \left( \Delta t_{cr}^{cm} \right)^2 & \text{if } 2 \sqrt{\frac{ M_i}{K^{max}_i}} > \Delta t_{cr}^{cm} \\ 1 & \text{otherwise} \end{cases}\end{split}\]

Consequently, the mass scaling produced by (8.9) is applied only to the high-frequency modes that could not otherwise be integrated stably with the central difference time integrator at the critical time step on the coarse mesh. The net result for the acceleration, \(\ddot{x}\), on the reference mesh is of the form

(8.10)\[\ddot{x} = \underbrace{ \Phi M^{-1}_c \Phi^T r}_\text{low frequency}+ \underbrace{ \tilde{M}^{-1} \left( r - M \Phi M^{-1}_c \Phi^T r\right)}_\text{mass-damped high frequency}\]

## 8.3.4. Energy Ratio: A Measure of Approximation

Kinetic energy calculations can be performed for the low frequency and high frequency contributions separately. Using time integrated acceleration components in (16), the kinetic energy in the low frequencies is,

(8.11)\[KE_{lf} = \frac{1}{2} M \| v_{lf} + \Delta t \Phi M^{-1}_c \Phi^T r \|^2\]

Likewise, the kinetic energy in the high frequencies is,

(8.12)\[KE_{hf} = \frac{1}{2} M \| v_{hf} + \Delta t \tilde{M}^{-1} \left( r - M \Phi M^{-1}_c \Phi^T r\right)\| ^2\]

With the kinetic energy quantities, an energy ratio is computed as follows,

\[{\rm Energy}~{\rm Ratio} = \frac{KE_{lf}}{KE_{lf} + KE_{hf}}\]

Obviously, the time integrated estimates of the kinetic energies require additional memory yet they provide a useful measure for the approximations being made with the explicit-explicit modal filtering. When there is little or no approximation made using a mass-damped high frequency response the energy ratio is asymptotically approaching unity. In contrast, when the approximation error is significant, the energy ratio is well below one.

[1]

T. Belytschko, W.K. Liu, and B. Moran. Nonlinear Finite Elements for Continua and Structures. John Wiley & Sons, 2000.

[2]

L. Olovsson, K. Simonsson, and M. Unosson. Selective mass scaling for explicit finite element analyses. International Journal for Numerical Methods in Engineering, 63(10):1436–1445, 2005.

[3] (1,2,3)

M. Heinstein, F. Mello, and C. Dohrmann. A nodal-based stable time step predictor for transient dynamcis with explicit time integration. Am Soc Mech Eng Press Vessels, Piping Division, 343:225–229, 1996.

[4] (1,2)

T.J.R. Hughes. The Finite Element Method: Linear static and dynamic finite element analysis. Dover, 2000. Reprint of “The Finite Element Method”, Prentice-Hall, 1987.

[5] (1,2)

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.

[6]

E. Hinton, T. Rock, and O.C. Zienkiewicz. A note on mass lumped and related processes in the finite element method. Earthquake Engineering and Structural Dynamics, 4(3):245–249, 1976.

[7]

A. Toselli and O. Widlund. Domain Decomposition Methods – Algorithms and Theory. Springer Series in Computational Mathematics. Springer, 2005.

[8] (1,2)

A. Brandt. General highly accurate algebraic coarsening. Elect. Trans. Numer. Anal., 10:1–20, 2000.

[9]

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.