Explicit dynamic analysis

Abaqus/Explicit provides explicit direct integration for performing dynamic analysis of problems in which inertia effects are considered.

This page discusses:

Stability
Bulk viscosity
Rotational bulk viscosity for shell elements

See Also

In Other Guides

Explicit Dynamic Analysis

ProductsAbaqus/Explicit

The explicit dynamics analysis procedure in Abaqus/Explicit is based upon the implementation of an explicit integration rule together with the use of diagonal or “lumped” element mass matrices. The equations of motion for the body are integrated using the explicit central difference integration rule

{\dot{u}}^{(i + \frac{1}{2})} = {\dot{u}}^{(i - \frac{1}{2})} + \frac{Δ t^{(i + 1)} + Δ t^{(i)}}{2} {\ddot{u}}^{(i)},

u^{(i + 1)} = u^{(i)} + Δ t^{(i + 1)} {\dot{u}}^{(i + \frac{1}{2})},

where $\dot{u}$ is velocity and $\ddot{u}$ is acceleration. The superscript $(i)$ refers to the increment number and $i - \frac{1}{2}$ and $i + \frac{1}{2}$ refer to midincrement values. The central difference integration operator is explicit in that the kinematic state can be advanced using known values of ${\dot{u}}^{(i - \frac{1}{2})}$ and ${\ddot{u}}^{(i)}$ from the previous increment. The explicit integration rule is quite simple but by itself does not provide the computational efficiency associated with the explicit dynamics procedure. The key to the computational efficiency of the explicit procedure is the use of diagonal element mass matrices because the inversion of the mass matrix that is used in the computation for the accelerations at the beginning of the increment is triaxial:

{\ddot{u}}^{(i)} = M^{- 1} \cdot (F^{(i)} - I^{(i)}),

where $M$ is the diagonal lumped mass matrix, $F$ is the applied load vector, and $I$ is the internal force vector. The explicit procedure requires no iterations and no tangent stiffness matrix.

Special treatment of the mean velocities ${\dot{u}}^{(i + \frac{1}{2})}$ , ${\dot{u}}^{(i - \frac{1}{2})}$ etc. is required for initial conditions, certain constraints, and presentation of results. For presentation of results, the state velocities are stored as a linear interpolation of the mean velocities:

{\dot{u}}^{(i + 1)} = {\dot{u}}^{(i + \frac{1}{2})} + \frac{1}{2} Δ t^{(i + 1)} {\ddot{u}}^{(i + 1)} .

The central difference operator is not self-starting because the value of the mean velocity ${\dot{u}}^{(- \frac{1}{2})}$ needs to be defined. The initial values (at time $t = 0$ ) of velocity and acceleration are set to zero unless they are specified by the user. We assert the following condition:

{\dot{u}}^{(+ \frac{1}{2})} = {\dot{u}}^{(0)} + \frac{Δ t^{(1)}}{2} {\ddot{u}}^{(0)} .

Substituting this expression into the update expression for ${\dot{u}}^{(i + \frac{1}{2})}$ yields the following definition of ${\dot{u}}^{(- \frac{1}{2})}$ :

{\dot{u}}^{(- \frac{1}{2})} = {\dot{u}}^{(0)} - \frac{Δ t^{(0)}}{2} {\ddot{u}}^{(0)} .

Stability

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 eigenvalue in the system as

Δ t \leq \frac{2}{ω_{m a x}} .

In Abaqus/Explicit a small amount of damping is introduced to control high frequency oscillations. With damping the stable time increment is given by

Δ t \leq \frac{2}{ω_{m a x}} (\sqrt{1 + ξ^{2}} - ξ),

where $ξ$ is the fraction of critical damping in the highest mode. Contrary to our usual engineering intuition, introducing damping to the solution reduces the stable time increment.

The time incrementation scheme in Abaqus/Explicit is fully automatic and requires no user intervention. Abaqus/Explicit uses an adaptive algorithm to determine conservative bounds for the highest element frequency. An estimate of the highest eigenvalue in the system can be obtained by determining the maximum element dilatational mode of the mesh. The stability limit based upon this highest element frequency is conservative in that it will give a smaller stable time increment than the true stability limit that is based upon the maximum frequency of the entire model. In general, constraints such as boundary conditions and contact have the effect of compressing the eigenvalue spectrum, which the element by element estimates do not take into account. Abaqus/Explicit contains a global estimation algorithm, which determines the maximum frequency of the entire model. This algorithm continuously updates the estimate for the maximum frequency. Abaqus/Explicit initially uses the element by element estimates. As the step proceeds, the stability limit will be determined from the global estimator once the algorithm determines that the accuracy of the global estimation is acceptable. The global estimation algorithm is not used when any of the following capabilities are included in the model: fluid elements, JWL equation of state, infinite elements, material damping, dashpots, thick shells (thickness to characteristic length ratio larger than 0.92) or thick beams (thickness to length ratio larger than 1.0), and nonisotropic elastic materials with temperature and field variable dependency.

A trial stable time increment is computed for each element in the mesh using the following expression:

Δ t = \frac{2}{ω_{m a x}^{e l e m e n t}},

where $ω_{m a x}^{e l e m e n t}$ is the element maximum eigenvalue. A conservative estimate of the stable time increment is given by the minimum taken over all the elements. The above stability limit can be rewritten as

Δ t = m i n (\frac{L_{e}}{c_{d}}),

where $L_{e}$ is the characteristic element dimension and $c_{d}$ is the current effective, dilatational wave speed of the material. The characteristic element dimension is derived from an analytic upper bound expression for the maximum element eigenvalue. Considering the 4-node uniform strain quadrilateral (CPE4R), the characteristic element dimension is

L_{e} = \frac{A}{\sqrt{B_{i I} B_{i I}}},

where $A$ is the element area and $B_{i I}$ is the element gradient operator (see Solid isoparametric quadrilaterals and hexahedra). Similar characteristic element dimensions are derived for all element types found in Abaqus/Explicit.

The current dilatational wave speed is determined in Abaqus/Explicit by calculating the effective hypoelastic material moduli from the material's constitutive response. Effective Lamé's constants, $\hat{λ}$ and $\hat{μ}$ , are determined in the following manner. We define $Δ p$ as the increment in the equivalent pressure stress, $p = - \frac{1}{3} trace (σ)$ , $Δ S$ as the increment in the deviatoric stress, $Δ ϵ_{v o l}$ as the increment of volumetric strain, and $Δ e$ as the deviatoric strain increment. We assume a hypoelastic stress-strain rule of the form

Δ p = - \hat{K} Δ ϵ_{v o l},

Δ S = 2 \hat{μ} Δ e,

where $\hat{K}$ is the effective bulk modulus. The effective moduli can then be computed as

\hat{K} = \frac{- Δ p}{Δ ϵ_{v o l}},

\hat{μ} = \frac{1}{2} \frac{Δ S : Δ e}{Δ e : Δ e},

\hat{λ} = \hat{K} - \frac{2}{3} \hat{μ} .

If the strain increments are insignificant, these relationships will not yield numerically meaningful results. In this circumstance Abaqus/Explicit sets the effective Lamé's constants to the initial values for the material, $λ_{o}$ and $2 μ_{o}$ . In the case where the volumetric strain increment is significant but the deviatoric stress increment is not, the effective shear modulus is estimated to be

2 \hat{μ} = \frac{1}{2} (3 (λ_{o} + 2 μ_{o}) - 3 \hat{K}) .

These effective moduli represent the element stiffness and determine the current dilatational wave speed in the element as

c_{d} = \sqrt{\frac{(\hat{λ} + 2 \hat{μ})}{ρ}} .

Bulk viscosity

Bulk viscosity introduces damping associated with the volumetric straining. Its purpose is to improve the modeling of high-speed dynamic events.

There are two forms of bulk viscosity in Abaqus/Explicit. The first is found in all elements and is introduced to damp the “ringing” in the highest element frequency. This damping is sometimes referred to as truncation frequency damping. It generates a bulk viscosity pressure, which is linear in the volumetric strain:

p_{b v 1} = b_{1} ρ c_{d} L_{e} {\dot{ϵ}}_{v o l},

where $b_{1}$ is a damping coefficient (default=.06), $ρ$ is the current material density, $c_{d}$ is the current dilatational wave speed, $L_{e}$ is an element characteristic length, and ${\dot{ϵ}}_{v o l}$ is the volumetric strain rate.

The second form of bulk viscosity pressure is found only in solid continuum elements (except CPS4R). This form is quadratic in the volumetric strain rate:

p_{b v 2} = ρ {(b_{2} L_{e} {\dot{ϵ}}_{v o l})}^{2},

where $b_{2}$ is a damping coefficient (default=1.2) and all other quantities are as defined for the linear bulk viscosity. The quadratic bulk viscosity is applied only if the volumetric strain rate is compressive.

The quadratic bulk viscosity pressure will smear a shock front across several elements and is introduced to prevent elements from collapsing under extremely high velocity gradients. Consider a simple one element problem in which the nodes on one side of the element are fixed and the nodes on the other side have an initial velocity in the direction of the fixed nodes. If the initial velocity is equal to the dilatational wave speed of the material, the element—without the quadratic bulk viscosity—would collapse to zero volume in one time increment (because the stable time increment size is precisely the transit time of a dilatational wave across the element). The quadratic bulk viscosity pressure will introduce a resisting pressure that will prevent the element from collapsing.

The bulk viscosity pressure is not included in the material point stresses because it is intended as a numerical effect only—it is not considered to be part of the material's constitutive response. The bulk viscosity pressures are based upon the dilatational mode of each element. The fraction of critical damping in the dilatational mode of each element is given by

ξ = b_{1} - b_{2}^{2} \frac{L_{e}}{c_{d}} \min (0, {\dot{ϵ}}_{v o l}) .

Linear bulk viscosity is always included in Abaqus/Explicit. The parameters $b_{1}$ and $b_{2}$ can be redefined by the user. The default values are $b_{1} = 0.06$ and $b_{2} = 1.2$ . The bulk viscosity parameters can be changed from step to step. If the default values are changed in a step, the new values will be used in any subsequent steps unless they are redefined.

Rotational bulk viscosity for shell elements

For the displacement degrees of freedom, bulk viscosity introduces damping associated with volumetric straining. Linear bulk viscosity or truncation frequency damping is used to damp the high frequency ringing that leads to unwanted noise in the solution or spurious overshoot in the response amplitude. For the same reason, in shells the high frequency ringing in the rotational degrees of freedom is damped with linear bulk viscosity acting on the mean curvature strain rate. This damping generates a bulk viscosity “pressure moment,” m, which is linear in the mean curvature strain rate:

m = b_{1} \frac{h_{0}^{2}}{12} ρ c_{d} L \frac{Δ κ}{Δ t},

where $b_{1}$ is a damping coefficient (default = 0.06), $h_{0}$ is the original thickness, $ρ$ is the mass density, $c_{d}$ is the current dilatational wave speed, L is the characteristic length used for rotary inertia and transverse shear stiffness scaling, and $Δ κ = Δ κ_{11} + Δ κ_{22}$ is twice the increment in mean curvature. The dilatational wave speed is given in terms of the effective Lamé constants as

c_{d} = \sqrt{\frac{(\hat{λ} + 2 \hat{μ})}{ρ}} .

The resultant pressure moment $m h$ , where h is the current thickness, is added to the direct components of the moment resultant.