Implicit dynamic analysis

ProductsAbaqus/Standard

Abaqus offers dynamic analysis options for both linear and nonlinear problems. In the case of purely linear systems, methods based on the eigenmodes of the system are almost always chosen because they can provide insight into the structure's behavior that is not otherwise available and because they are usually significantly more cost-effective than the direct integration methods that are usually used for nonlinear problems. The linear dynamic analysis methods provided in Abaqus/Standard are discussed in Modal dynamics. For mildly nonlinear dynamic analysis problems the “modal projection method” is provided. The basis of that method is to use the eigenmodes of the linear system (extracted from an eigenfrequency analysis) as a set of global Ritz functions—a set of global interpolation functions, in the terminology of the finite element method—whose amplitudes define the response. Abaqus/Standard provides direct time integration using the explicit, central difference operator for this option. For any more severely nonlinear case the dynamic response is obtained by direct time integration of all of the degrees of freedom of the finite element model.

The choice of operator used to integrate the equations of motion in a dynamic analysis is influenced by many factors. Abaqus/Standard is designed to analyze structural components, by which we mean that the overall dynamic response of a structure is sought, in contrast to wave propagation solutions associated with relatively local response in continua. Belytschko (1976) labels these “inertial problems” and classifies them by stating that “wave effects such as focusing, reflection, and diffraction are not important.” Structural problems are considered “inertial” because the response time sought is long compared to the time required for waves to traverse the structure.

Dynamic integration operators are broadly characterized as implicit or explicit. Explicit schemes, as used in Abaqus/Explicit, obtain values for dynamic quantities at $t + Δ t$ based entirely on available values at time t. The central difference operator, which is the most commonly used explicit operator for stress analysis applications, is only conditionally stable, the stability limit being approximately equal to the time for an elastic wave to cross the smallest element dimension in the model. Implicit schemes remove this upper bound on time step size by solving for dynamic quantities at time $t + Δ t$ based not only on values at t, but also on these same quantities at $t + Δ t$ . But because they are implicit, nonlinear equations must be solved. In structural problems implicit integration schemes usually give acceptable solutions with time steps typically one or two orders of magnitude larger than the stability limit of simple explicit schemes, but the response prediction will deteriorate as the time step size, $Δ t$ , increases relative to the period, T, of typical modes of response. See, for example, Hilber, Hughes, and Taylor (1977) and Hilber and Hughes (1978) for a discussion of such errors. Three factors should be considered when selecting the maximum allowable time step size: the rate of variation of the applied loading, the complexity of the nonlinear damping and stiffness properties, and the typical period of vibration of the structure. In general, a maximum increment versus period ratio $Δ t / T <$ 1/10 is a good rule of thumb for obtaining reliable results. Thus, the relative economy of the two techniques of integration depends on the stability limit of the explicit scheme, the ease with which the nonlinear equations can be solved for the implicit operator, the relative size of time increments that can provide acceptable accuracy with the implicit scheme compared to the stability limit of the explicit scheme, and the size of the model.

In Abaqus/Standard the time step for implicit integration can be chosen automatically on the basis of the “half-increment residual,” a concept introduced in Hibbitt and Karlsson (1979). By monitoring the values of equilibrium residuals at $t + Δ t / 2$ once the solution at $t + Δ t$ has been obtained, the accuracy of the solution can be assessed and the time step adjusted appropriately.

To discuss the dynamic procedures further, we write out the d'Alembert force in the overall equilibrium equation, Equation 1. The body force at a point, $f$ , can be written as an externally prescribed body force, $F$ , and a d'Alembert force:

f = F - ρ \ddot{u},

where $ρ$ is the current density of the material at this point and $u$ is the displacement of the point. The body force term in the virtual work equation is

\int_{V} f \cdot δ v d V = \int_{V} F \cdot δ v d V - \int_{V} ρ \ddot{u} \cdot δ v d V .

The d'Alembert term can be written more conveniently in terms of the reference volume and reference density, $ρ_{0}$ , as

\int_{V_{0}} ρ_{0} \ddot{u} \cdot δ v d V_{0},

where $\ddot{u}$ is the acceleration field. When implicit integration is used, the equilibrium equations are written at the end of a time step (at time $t + Δ t$ ), and $\ddot{u}$ is calculated from the time integration operator. The interpolator approximates the displacement at a point as

u = N^{N} u^{N},

so that

\ddot{u} = N^{N} {\ddot{u}}^{N},

provided that $N^{N}$ are not displacement dependent. This is the case for most of the elements in Abaqus; the form taken for the d'Alembert force terms in those instances where it is not true (the Hermite cubic beams, B23 and B33) is discussed in detail in Elements. With this interpolation assumption, the d'Alembert force term is

- (\int_{V_{0}} ρ_{0} N^{N} \cdot N^{M} d V_{0}) {\ddot{u}}^{M};

that is, the consistent mass matrix times the accelerations of the nodal variables. The finite element approximation to equilibrium, Equation 2, is

\begin{matrix} (1) & M^{N M} {\ddot{u}}^{M} + I^{N} - P^{N} = 0, \end{matrix}

where

M^{N M} = \int_{V_{0}} ρ_{0} N^{N} \cdot N^{M} d V_{0}

is the consistent mass matrix,

I^{N} = \int_{V_{0}} β^{N} : σ d V_{0}

is the internal force vector, and

P^{N} = \int_{S} N^{N} \cdot t d S + \int_{V} N^{N} \cdot F d V

is the external force vector. In this context the terms “matrix” and “vector” refer to matrices and vectors in the space of the nodal variables $u^{N}$ .

The definition of the mass matrix introduced above is the “consistent” mass: the mass matrix obtained by consistent use of the interpolation. The first-order elements in Abaqus all use “lumped” mass, where the mass matrix is a diagonal matrix. The lumped matrix is obtained by adding each row of the consistent matrix onto the diagonal. For these first-order elements the lumped mass matrix gives more accurate results in numerical experiments that calculate the natural frequencies of simple models.

Implicit operators available in Abaqus/Standard for time integration of the dynamic problem include the operator defined by Hilber, Hughes, and Taylor (1977) and the backward Euler operator. The Hilber-Hughes-Taylor operator is a generalization of the Newmark operator with controllable numerical damping—the damping being most valuable in the automatic time stepping scheme, because the slight high-frequency numerical noise inevitably introduced when the time step is changed is removed rapidly by a small amount of numerical damping. The operator replaces the actual equilibrium equation (Equation 1) with a balance of d'Alembert forces at the end of the time step and a weighted average of the static forces at the beginning and end of the time step:

\begin{matrix} (2) & {M^{N M} {\ddot{u}}^{M} |}_{t + Δ t} + (1 + α) ({I^{N} |}_{t + Δ t} - {P^{N} |}_{t + Δ t}) - α ({I^{N} |}_{t} - {P^{N} |}_{t}) + {L^{N} |}_{t + Δ t} = 0, \end{matrix}

where ${L^{N} |}_{t + Δ t}$ is the sum of all Lagrange multiplier forces associated with degree of freedom N. The operator definition is completed by the Newmark formulae for displacement and velocity integration:

\begin{matrix} (3) & {u |}_{t + Δ t} = {u |}_{t} + {Δ t \dot{u} |}_{t} + Δ t^{2} ({(\frac{1}{2} - β) \ddot{u} |}_{t} + {β \ddot{u} |}_{t + Δ t}) \end{matrix}

and

\begin{matrix} (4) & {\dot{u} |}_{t + Δ t} = {\dot{u} |}_{t} + Δ t ({(1 - γ) \ddot{u} |}_{t} + {γ \ddot{u} |}_{t + Δ t}), \end{matrix}

with

β = \frac{1}{4} {(1 - α)}^{2}, γ = \frac{1}{2} - α and - \frac{1}{2} \leq α \leq 0.

Hilber, Hughes, and Taylor (1977) present cogent arguments for the use of Equation 2–Equation 4 for integrating structural dynamics problems. The main appeal of the operator is its controllable numerical damping and the form this damping takes, slowly growing at low frequencies, with more rapid growth in damping at high frequencies. Control over the amount of numerical damping is provided by the parameter $α$ : with $α = 0$ , there is no damping and the operator is the trapezoidal rule (Newmark, $β = 1 / 4$ ); while with $α = - 1 / 3$ , significant damping is available. This operator is used primarily because the slight numerical damping it offers is needed in the automatic time stepping scheme. Each time step change introduces some slight noise or “ringing” into the solution; a little numerical damping ( $α = - 0.05$ seems a good choice) quickly removes this high-frequency noise without having any significant effect on the meaningful, lower frequency response. An energy content output is available and should be printed to monitor the overall energy balance. This has been done in dynamic examples in the Abaqus Example Problems Guide and shows that the numerical dissipation is always quite small (less than 1% of the total energy).

The integration of rotations during implicit dynamic calculations is done to preserve accuracy in cases where the rotary inertia is different in different directions in a body. For this purpose the accelerations are integrated in the body axis system, so that Newmark's formula gives the change in velocity as

{{\dot{ϕ}}^{α} |}_{t + Δ t} = {{\dot{ϕ}}^{α} |}_{t} + Δ t [{γ {\ddot{ϕ}}^{α} |}_{t + Δ t} + {(1 - γ) {\ddot{ϕ}}^{α} |}_{t}],

where ${\dot{ϕ}}^{α}$ is the angular velocity of the node and ${\ddot{ϕ}}^{α}$ is its angular acceleration, both taken in the current direction of the $α$ body axis at time t or at time $t + Δ t$ ; $Δ t$ is the time increment.

In the global system this is

{\dot{ϕ} |}_{t + Δ t} = {Δ t γ \ddot{ϕ} |}_{t + Δ t} + [{{e^{α} |}_{t + Δ t} e^{α} |}_{t}] \cdot [{\dot{ϕ} |}_{t} + {Δ t (1 - γ) \ddot{ϕ} |}_{t}],

where $e^{α} = e^{α} (ϕ)$ are orthonormal base vectors defining the body axis system. Since these are orthonormal vectors, this can be rewritten as

\begin{matrix} (5) & {\dot{ϕ} |}_{t + Δ t} = {Δ t γ \ddot{ϕ} |}_{t + Δ t} + Δ C \cdot [{\dot{ϕ} |}_{t} + {Δ t (1 - γ) \ddot{ϕ} |}_{t}], \end{matrix}

where $Δ C$ is the incremental rotation matrix. Let $Δ θ$ be the increment in rotation from time t to $t + Δ t$ . Then $Δ C$ is

Δ C = exp [Δ \hat{θ}],

where $Δ \hat{θ}$ is the skew-symmetric matrix with axial vector $Δ θ$ . See Rotation variables for a discussion of rotation variables, rotation matrices, and the exponential mapping of a skew-symmetric matrix.

Newmark's formula for the time integration of the rotation increment in the body axis system gives

Δ θ^{α} = {Δ t {\dot{ϕ}}^{α} |}_{t} + Δ t^{2} [{(\frac{1}{2} - β) {\ddot{ϕ}}^{α} |}_{t} + {β {\ddot{ϕ}}^{α} |}_{t + Δ t}] .

Since $Δ C \cdot Δ θ = exp [Δ \hat{θ}] \cdot Δ θ = Δ θ$ , the components $Δ θ^{α}$ are the same relative to the body axes at time t or $t + Δ t$ . In the global system this is

Δ θ = {Δ t^{2} β \ddot{ϕ} |}_{t + Δ t} + Δ C \cdot [{Δ t \dot{ϕ} |}_{t} + {Δ t^{2} (\frac{1}{2} - β) \ddot{ϕ} |}_{t}] .

Solving for the unknown velocity and acceleration at time $t + Δ t$ , the velocity update equation is

{\dot{ϕ} |}_{t + Δ t} = \frac{γ}{Δ t β} Δ θ + Δ C \cdot [{(1 - \frac{γ}{β}) \dot{ϕ} |}_{t} + {Δ t (1 - \frac{γ}{2 β}) \ddot{ϕ} |}_{t}],

and the acceleration update equation becomes

{\ddot{ϕ} |}_{t + Δ t} = {\frac{1}{Δ t γ} \dot{ϕ} |}_{t + Δ t} + Δ C \cdot [{(1 - \frac{1}{γ}) \ddot{ϕ} |}_{t} - {\frac{1}{Δ t γ} \dot{ϕ} |}_{t}] .

The backward Euler operator solves the equilibrium equation (Equation 1) at the end of the time step and updates the displacement and velocity using

{u |}_{t + Δ t} = {u |}_{t} + {Δ t \dot{u} |}_{t + Δ t}

and

{\dot{u} |}_{t + Δ t} = {\dot{u} |}_{t} + {Δ t \ddot{u} |}_{t + Δ t} .

The automatic time stepping for dynamic problems in Abaqus is based on the half-increment residual first proposed in Hibbitt and Karlsson (1979). The concept is quite appealing intuitively. Satisfaction of Equation 1 at the end of each time step (or actually of the Hilber-Hughes-Taylor form, Equation 2) ensures equilibrium in the discrete sense of the finite element model at these points in time but does not say anything about the quality of equilibrium at intermediate time points. The idea of the half-increment residual is to calculate the equilibrium residual error (the left-hand side of Equation 1) at some intermediate time point (chosen as $t + Δ t / 2$ ) and to assess the error in the dynamic response prediction by the magnitude of that error.

The half-increment residual is based on the assumption that the accelerations vary linearly over the time interval (this is the basis of Newmark's formulae) so that, for any nodal displacement or rotation component, u:

{\ddot{u} |}_{τ} = {(1 - τ) \ddot{u} |}_{t} + {τ \ddot{u} |}_{t + Δ t}, 0 \leq τ \leq 1.

Having already solved for the state at $t + Δ t$ , this equation, together with Newmark's formulae now written for the time interval from t to $t + Δ t$ , requires that

\begin{matrix} {Δ u |}_{τ} & = {τ^{3} Δ u |}_{t + Δ t} + {τ (1 - τ^{2}) Δ t \dot{u} |}_{t} + {τ^{2} (1 - τ) \frac{Δ t^{2}}{2} \ddot{u} |}_{t}, \\ {\dot{u} |}_{τ} & = {\frac{γ}{β τ Δ t} Δ u |}_{τ} + {(1 - \frac{γ}{β}) \dot{u} |}_{t} + {(1 - \frac{γ}{2 β}) τ Δ t \ddot{u} |}_{t}, and \\ {\ddot{u} |}_{τ} & = {\frac{1}{β τ^{2} Δ t^{2}} Δ u |}_{τ} - {\frac{1}{β τ Δ t} \dot{u} |}_{t} + {(1 - \frac{1}{2 β}) \ddot{u} |}_{t}, \end{matrix}

where

{Δ u |}_{t + Δ t} = {u |}_{t + Δ t} - {u |}_{t}

is the increment in displacement obtained for the time step, $Δ t$ .

With these equations it is possible to evaluate the equilibrium residual at any time within the step. Presumably, if the solution is accurate, this residual will be small compared to significant forces in the problem. The residual at the end of the time step is

{R^{N} |}_{t + Δ t} = {M^{N M} {\ddot{u}}^{M} |}_{t + Δ t} + (1 + α) ({I^{N} |}_{t + Δ t} - {P^{N} |}_{t + Δ t}) - α ({I^{N} |}_{t} - {P^{N} |}_{t}) + {L^{N} |}_{t + Δ t} .

The residual at the start of the step is

{R^{N} |}_{t} = {M^{N M} {\ddot{u}}^{M} |}_{t} + (1 + α) ({I^{N} |}_{t} - {P^{N} |}_{t}) - α ({I^{N} |}_{t^{-}} - {P^{N} |}_{t^{-}}) + {L^{N} |}_{t},

where $t^{-}$ is the time at the start of the previous time step during normal time stepping analysis or $t^{-} = t$ if this is the first increment after an initial acceleration or impact calculation. Then the residual at $t + Δ t / 2$ is defined as

\begin{matrix} {R^{N} |}_{t + Δ t / 2} \overset{def}{=} & {M^{N M} {\ddot{u}}^{M} |}_{t + Δ t / 2} + (1 + α) ({I^{N} |}_{t + Δ t / 2} - {P^{N} |}_{t + Δ t / 2}) \\ - \frac{1}{2} α ({I^{N} |}_{t} - {P^{N} |}_{t} + {I^{N} |}_{t^{-}} - {P^{N} |}_{t^{-}}) + {L^{N} |}_{t + Δ t / 2}, \end{matrix}

where the ${I^{N} |}_{t + Δ t / 2},$ etc. are computed for conditions at time $(t + Δ t / 2)$ and

{L^{N} |}_{t + Δ t / 2} \overset{def}{=} \frac{1}{2} ({L^{N} |}_{t + Δ t} + {L^{N} |}_{t}) .

The “half-increment residual,” ${R |}_{t + Δ t / 2}$ , is defined as the magnitude of the largest entry in ${R^{N} |}_{t + Δ t / 2}$ and provides a measure of accuracy of the time-stepping solution.

The motivation behind the calculation of the half-increment residual is to provide a measure of the accuracy of the solution for a given time step. Numerical tests show that it provides a sensitive accuracy check on dynamic solutions and suggest that, if P is a typical magnitude of real forces in an undamped elastic system (for which the high-frequency response must be modeled reasonably accurately), then

if ${R |}_{t + Δ t / 2} \approx 0.1 P$ consistently, the time stepping solution has high accuracy;
if ${R |}_{t + Δ t / 2} \approx P$ consistently, the time stepping solution has moderately good accuracy;
if ${R |}_{t + Δ t / 2} \approx 10 P$ consistently, the time stepping solution is rather coarse.

Problems with large amounts of natural dissipation of energy, such as elastic-plastic systems, are usually less sensitive to time step choice than purely elastic problems, because the energy that appears in higher frequency modes is quickly dissipated. In such cases maximum half-increment residuals in the range of 1–10 times typical forces indicate quite acceptable accuracy for most studies, and even values of 10–100 times typical forces can give useful results for primary effects, such as overall deformation. Thus, the method can offer relatively cost-effective solutions for highly dissipative systems for which we require only moderately accurate prediction of the overall response.

The half-increment residual is the basis of the adaptive time incrementation scheme. If the half-increment residual is small, it indicates that the accuracy of the solution is high and the time step can be increased safely; conversely, if the half-increment calculation shows the solution is coarse, the time step used in the solution should be reduced. The algorithm is described in detail in Time integration accuracy in transient problems. The algorithm is purely empirical, but experience shows that it works quite well (Hibbitt and Karlsson, 1979), most especially in initially excited problems with high dissipation, such as impulsively loaded problems (or problems with short duration forcing) with extensive plasticity. In these cases the scheme is economic because the time step naturally increases as the solution progresses and the high-frequency response is dissipated.

The above observations are based on using the Hilber-Hughes-Taylor operator with $α = - 0.05.$ The slight numerical damping that the operator introduces removes the noise that inevitably enters the solution when the time step is changed. Even in quite lengthy problems the overall energy totals, computed on the basis of physical mechanisms in the model, balance well. The energy balance calculation is useful in assessing a solution—for example, the extent to which energy is dissipated by plasticity can be measured—and it is recommended that the user request occasional printout of the energy balance calculation when doing any analysis with Abaqus.