Explicit Dynamic Analysis

The explicit dynamics procedure performs a large number of small time increments efficiently. An explicit central-difference time integration rule is used; each increment is relatively inexpensive (compared to the direct-integration dynamic analysis procedure available in Abaqus/Standard) because there is no solution for a set of simultaneous equations. The explicit central-difference operator satisfies the dynamic equilibrium equations at the beginning of the increment, t; the accelerations calculated at time t are used to advance the velocity solution to time $t+\mathrm{\Delta }t/2$ and the displacement solution to time $t+\mathrm{\Delta }t$.

The explicit dynamic analysis procedure is based on the implementation of an explicit integration rule together with the use of diagonal (“lumped”) element mass matrices. A correction phase is used to implicitly treat connectors and constraints.

Predicted nodal accelerations, ${{\ddot{u}}_{\left(i\right)}^N}^{Pred}$ , at the beginning of increment number i of an explicit dynamic step are computed by:

where $M^{NJ}$ is the mass matrix, $P^J$ is the applied load vector, and $I^J$ is the explicit internal force vector. A lumped mass matrix is used because its inverse is simple to compute and because the vector multiplication of the mass inverse by the inertial force requires only n operations, where n is the number of degrees of freedom in the model. The internal force vector, $I^J$ , is assembled from contributions from the individual elements such that a global stiffness matrix need not be formed.

The equations of motion for the body are integrated using the explicit central-difference integration rule to compute predicted mid-increment nodal velocity (or angular velocity), ${{\dot{u}}_{\left(i+\frac{1}{2}\right)}^N}^{Pred}$ , and predicted end-of-increment displacement (or rotation), ${u_{\left(i+1\right)}^N}^{Pred}$ as follows:

The central-difference integration operator is explicit in the sense that the predicted kinematic state is advanced using known values of ${\dot{u}}_{\left(i-\frac{1}{2}\right)}^N$ and $u_{\left(i\right)}^N$ from the previous increment.

However, the predicted kinematic quantities do not account for implicitly treated connections and constraints. To account for these entities, implicit forces are computed considering the predicted configuration. Calculation of the implicit force vector, ${F_{Impl}}_{\left(i\right)}^J$ , often requires solution of nonlinear systems of equations each increment, which is typically efficient but can degrade performance for situations with overlapping constraints, as discussed in [add xref]. These implicit forces cause nodal acceleration corrections at increment i:

Acceleration corrections are applied to determine corrected kinematic quantities that satisfy connection and constraint behavior as follows:

Sometimes it is necessary to enforce the constraints implicitly, especially if there are multiple constraints applied at a node. The implicit enforcement of constraints can lead to a performance degradation (see Computational Cost for Constraint Solution in Abaqus/Explicit).

The explicit integration scheme in Abaqus/Explicit requires nodal mass or inertia to exist at all activated degrees of freedom (see Conventions) unless constraints are applied using boundary conditions. More precisely, a nonzero nodal mass must exist unless all activated translational degrees of freedom are constrained and nonzero rotary inertia must exist unless all activated rotational degrees of freedom are constrained. Nodes that are part of a rigid body do not require mass, but the entire rigid body must possess mass and inertia unless constraints are used. Nodes that belong to Eulerian elements also do not require mass, since the surrounding Eulerian elements might be void at some time during the simulation.

When degrees of freedom at a node are activated by elements with a nonzero mass density (for example, solid, shell, beam) or mass and inertia elements, a nonzero nodal mass or inertia occurs naturally from the assemblage of lumped mass contributions.

When degrees of freedom at a node are activated by elements with no mass (for example, spring, dashpot, or connector elements), care must be taken either to constrain the node or to add mass and inertia as appropriate.

The explicit procedure integrates through time by using many small time increments. The central-difference operator is conditionally stable, and the stability limit for the operator (with no damping) is given in terms of the highest frequency of the system as

With damping, the stable time increment is given by

where ${\xi }_{max}$ is the fraction of critical damping in the mode with the highest frequency. Contrary to our usual engineering intuition, introducing damping to the solution reduces the stable time increment. In Abaqus/Explicit a small amount of damping is introduced in the form of bulk viscosity to control high frequency oscillations. Physical forms of damping, such as dashpots or material damping, can also be introduced. Bulk viscosity and material damping are discussed below.

The current dilatational wave speed, $c_d$, is determined in Abaqus/Explicit by calculating the effective hypoelastic material moduli from the material's constitutive response. Effective Lamé's constants, $\widehat{\lambda }$ and $\widehat{G}=2\widehat{\mu }$, are determined in the following manner. Define $\mathrm{\Delta }p$ as the increment in the mean stress, $\mathrm{\Delta }\mathbf{S}$ as the increment in the deviatoric stress, $\mathrm{\Delta }{ϵ}_{vol}$ as the increment of volumetric strain, and $\mathrm{\Delta }\mathbf{e}$ as the deviatoric strain increment. We assume a hypoelastic stress-strain rule of the form

The effective moduli can then be computed as

For shell elements defined by a shell cross-section that requires numerical integration (see Using a Shell Section Integrated during the Analysis to Define the Section Behavior), the effective moduli for the section are computed by integrating the effective moduli at the section points through the thickness. These effective moduli represent the element stiffness and determine the current dilatational wave speed in the element as

where $\rho$ is the density of the material.

In an isotropic, elastic material the effective Lamé's constants can be defined in terms of Young's modulus, E, and Poisson's ratio, $\nu$, by

and

The time increment used in an analysis must be smaller than the stability limit of the central-difference operator. Failure to use a small enough time increment will result in an unstable solution. When the solution becomes unstable, the time history response of solution variables such as displacements will usually oscillate with increasing amplitudes. The total energy balance will also change significantly.

If the model contains only one material type, the initial time increment is directly proportional to the size of the smallest element in the mesh. If the mesh contains uniform size elements but contains multiple material descriptions, the element with the highest wave speed will determine the initial time increment.

In nonlinear problems—those with large deformations and/or nonlinear material response—the highest frequency of the model will continually change, which consequently changes the stability limit. Abaqus/Explicit has two strategies for time incrementation control: fully automatic time incrementation (where the code accounts for changes in the stability limit) and fixed time incrementation.

The use of small increments (dictated by the stability limit) is advantageous because it allows the solution to proceed without iterations and without requiring tangent stiffness matrices to be formed. It also simplifies the treatment of contact.

The explicit dynamics procedure is ideally suited for analyzing high-speed dynamic events, but many of the advantages of the explicit procedure also apply to the analysis of slower (quasi-static) processes. A good example is sheet metal forming, where contact dominates the solution and local instabilities might form due to wrinkling of the sheet.

The results in an explicit dynamics analysis are not automatically checked for accuracy as they are in Abaqus/Standard (Abaqus/Standard uses the half-increment residual). In most cases this is not of concern because the stability condition imposes a small time increment such that the solution changes only slightly in any one time increment, which simplifies the incremental calculations. While the analysis might take an extremely large number of increments, each increment is relatively inexpensive, often resulting in an economical solution. It is not uncommon for Abaqus/Explicit to take over 10⁵ increments for an analysis. The method is, therefore, computationally attractive for problems where the total dynamic response time that must be modeled is only a few orders of magnitude longer than the stability limit; for example, wave propagation studies or some “event and response” applications.

The computer time involved in running a simulation using explicit time integration with a given mesh is proportional to the time period of the event. The time increment based on the element-by-element stability estimate can be rewritten (ignoring damping) in the form

where the minimum is taken over all elements in the mesh, $L_e$ is a characteristic length associated with an element (see Explicit dynamic analysis), $\rho$ is the density of the material in the element, and $\widehat{\lambda }$ and $\widehat{\mu }$ are the effective Lamé's constants for the material in the element (defined above).

The time increment from the global stability estimate might be somewhat larger, but for this discussion we will assume that the above inequality always holds (when the inequality does not hold, the solution time will be somewhat faster).

For linear, nonisotropic elastic materials this stability limit is further scaled down by the square root of the ratio of the effective material stiffness to the maximum material stiffness in one particular direction. Since this effectively means that the time increment can be no larger than the time required to propagate a stress wave across an element, the computer time involved in running a quasi-static analysis can be very large: the cost of the simulation is directly proportional to the number of time increments required.

The number of increments, n, required is $n=T/\mathrm{\Delta }t$ if $\mathrm{\Delta }t$ remains constant, where T is the time period of the event being simulated. (Even the element-by-element approximation of $\mathrm{\Delta }t$ will not remain constant in general, since element distortion will change $L_e$ and nonlinear material response will change the effective Lamé constants. But the assumption is sufficiently accurate for the purposes of this discussion.) Thus,

In a two-dimensional analysis refining the mesh by a factor of two in each direction will increase the run time in the explicit procedure by a factor of eight—four times as many elements and half the original time increment size. Similarly, in a three-dimensional analysis refining the mesh by a factor of two in each direction will increase the run time by a factor of sixteen.

In a quasi-static analysis it is expedient to reduce the computational cost by either speeding up the simulation or by scaling the mass. In either case the kinetic energy should be monitored to ensure that the ratio of kinetic energy to internal energy does not get too large—typically less than 10%.

Bulk viscosity introduces damping associated with volumetric straining. Its purpose is to improve the modeling of high-speed dynamic events (see Stability above for a discussion of the effect of damping on the stable time increment). Abaqus/Explicit contains two forms of bulk viscosity: linear and quadratic. Linear bulk viscosity is included by default in an Abaqus/Explicit analysis.

The bulk viscosity parameters $b_1$ and $b_2$ defined below can be redefined and can be changed from step to step. If the default values are changed in a step, the new values will be used in subsequent steps until they are redefined. Bulk viscosities defined this way apply to the whole model. For an individual element set the linear and quadratic bulk viscosities can be scaled by a factor by defining section controls (see Section Controls).

Defining inelastic material behavior, dashpots, etc. will introduce energy dissipation into a model. In addition to these mechanisms, general (“Rayleigh”) material damping can be introduced (see Material Damping). Adding damping to a model, especially stiffness proportional damping, ${\beta }_R$ , might significantly reduce the stable time increment.

Abaqus/Explicit writes critical elements (elements with the smallest stable time increments) and their stable time increment values to the output database at each summary increment. By default, the number of critical elements written to the output database is 10.

The deformation speed in an element is defined as the largest absolute value of all the deformation rate components of an element times the element characteristic length, $L_e$. You can request diagnostic information about the deformation speed within a step definition, as described below. In a multistep analysis diagnostic requests remain in effect until they are explicitly redefined.

There are some analyses in which it is useful to monitor the value of a variable at every increment. For example, in a force-driven analysis such as hydro-forming, the simulation time that is sufficient to model the completion of the physical process might depend on the magnitude of the displacement of a node or a group of nodes in the model. Another example is a drop test simulation where the postfailure response is not of interest. Monitoring the values of critical variables and halting the analysis when those variables exceed a given criterion can reduce computational expense and turnaround time.

For such problems Abaqus/Explicit allows output variables to be monitored during an analysis to verify whether or not their values have exceeded or fallen below user-specified values in specified element or node sets. Comparisons of specified element integration point variables, element section variables, or nodal variables with user-specified values are performed at every increment. At the first occurrence of a variable exceeding the user-specified bounds, the variable name, the associated element or node number, and the increment number are written to the status (.sta) file. In addition, you can request that the analysis be stopped and/or the output state be written in the increment following the one in which the variable has exceeded the user-specified bound. At the end of each step in which variables are monitored, the maximum, minimum, or absolute maximum value that each variable attains during the course of the analysis, along with the number of the element or node where the extreme value occurred, will be written to the status file.