Implicit Dynamic Analysis Using Direct Integration
A direct-integration dynamic analysis in
Abaqus/Standard:
must be used when nonlinear dynamic response is being studied;
can be fully nonlinear (general dynamic analysis) or can be based on
the modes of the linear system (subspace projection method); and
can be used to study a variety of applications, including:
dynamic responses requiring transient fidelity and involving
minimal energy dissipation;
dynamic responses involving nonlinearity, contact, and moderate
energy dissipation; and
quasi-static responses in which considerable energy dissipation
provides stability and improved convergence behavior for determining an
essentially static solution.
General nonlinear dynamic analysis in
Abaqus/Standard
uses implicit time integration to calculate the transient dynamic or
quasi-static response of a system. The procedure can be applied to a broad
range of applications calling for varying numerical solution strategies, such
as the amount of numerical damping required to obtain convergence and the way
in which the automatic time incrementation algorithm proceeds through the
solution. Typical dynamic applications fall into three categories:
Transient fidelity applications, such as an analysis of satellite
systems, require minimal energy dissipation. In these applications small time
increments are taken to accurately resolve the vibrational response of the
structure, and numerical energy dissipation is kept at a minimum. These
stringent requirements tend to degrade convergence behavior for simulations
involving contact or nonlinearities.
Moderate dissipation applications encompass a more general range of
dynamic events in which a moderate amount of energy is dissipated by
plasticity, viscous damping, or other effects. Typical applications include
various insertion, impact, and forming analyses. The response of these
structures can be either monotonic or nonmonotonic. Accurate resolution of
high-frequency vibrations is usually not of interest in these applications.
Some numerical energy dissipation tends to reduce solution noise and improve
convergence behavior in these applications without significantly degrading
solution accuracy.
Quasi-static applications are primarily interested in determining a
final static response. These problems typically show monotonic behavior, and
inertia effects are introduced primarily to regularize unstable behavior. For
example, the statically unstable behavior may be due to temporarily
unconstrained rigid body modes or “snap-through” phenomena. Large time
increments are taken when possible to obtain the final solution at minimal
computational cost. Considerable numerical dissipation may be required to
obtain convergence during certain stages of the loading history.
Based on the classifications listed above, you should indicate the type of
application you are studying when performing a general dynamic analysis.
Abaqus/Standard
assigns numerical settings based on your classification of the application
type, and this classification can significantly affect a simulation. In some
cases accurate results can be obtained with more than one application-type
setting, in which case analysis efficiency should be considered. A general
trend is that—among the three classifications—the high-dissipation quasi-static
classification tends to result in the best convergence behavior and the
low-dissipation transient fidelity classification tends to have the highest
likelihood of convergence difficulty.
Diagnostics for Modeling Errors Associated with Mass Properties
Accurate representation of inertia properties is necessary for accurate
dynamic analyses. In some cases
Abaqus/Standard
provides diagnostic messages when it detects likely modeling errors associated
with the specification of inertia properties. The most common way of specifying
inertia properties is with material densities.
Abaqus/Standard
issues a warning message to the data (.dat) file if a
material density is omitted in a dynamic analysis (this warning is not issued
if the density is zero only for certain values of temperature or field
variables). Other methods of specifying inertia properties include:
point mass and rotary inertia definitions, and
constraining nodes without inertia themselves to nodes having inertia
properties defined.
In some circumstances
Abaqus/Standard
attempts to solve systems of equations involving effective inversion of the
global mass matrix to directly adjust velocities and accelerations during a
general dynamic analysis as described in
Initial Conditions
and
Intermittent Contact/Impact
below. These additional velocity and acceleration adjustments occur by default
only for transient fidelity application types as defined above. If the global
mass matrix is found to be singular,
Abaqus/Standard
issues an error message by default, because singular mass is an indication that
the mass properties are not realistic due to a modeling error.
Diagnostic feedback specific to the global mass matrix being singular is
typically not provided for quasi-static and moderate dissipation application
types, although warnings typically are issued regarding the lack of material
density. Singular mass is not necessarily detrimental to a quasi-static
analysis. For example, it would be reasonable to only define inertia properties
(such as density) in components or regions with temporary static instabilities
(such as initially unconstrained rigid body modes that become constrained once
contact occurs) in a quasi-static analysis.
You can control the course of action Abaqus/Standard takes on detecting a singular global mass matrix.
Numerical Details
The effect of the application-type classification on numerical aspects of
general dynamic analyses is described below. In most cases the settings
determined by the application type are sufficient to successfully perform an
analysis. However, detailed user controls are provided to override settings on
an individual basis.
Time Integration Methods
Abaqus/Standard
uses the Hilber-Hughes-Taylor time integration by default unless you specify
that the application type is quasi-static. The Hilber-Hughes-Taylor operator is
an extension of the Newmark -method.
Numerical parameters associated with the Hilber-Hughes-Taylor operator are
tuned differently for moderate dissipation and transient fidelity applications
(as discussed later in this section). The backward Euler operator is used by
default if the application classification is quasi-static.
These time integration operators are implicit, which means that the
operator matrix must be inverted and a set of simultaneous nonlinear dynamic
equilibrium equations must be solved at each time increment. This solution is
done iteratively using Newton's method. The principal advantage of these
operators is that they are unconditionally stable for linear systems; there is
no mathematical limit on the size of the time increment that can be used to
integrate a linear system. An unconditionally stable integration operator is of
great value when studying structural systems because a conditionally stable
integration operator (such as that used in the explicit method) can lead to
impractically small time steps and, therefore, a computationally expensive
analysis.
Marching through a simulation with a finite time increment size generally introduces some
degree of numerical damping. This damping differs from the material damping discussed in
Material Damping (and in many
cases these two forms of damping work well together). The amount of damping associated
with the time integration varies among the operator types (for example, the backward
Euler operator tends to be more dissipative than the Hilber-Hughes-Taylor operator) and
in many cases (such as with the Hilber-Hughes-Taylor operator) depends on settings of
numerical parameters associated with the operator. The ability of the operator to
effectively treat contact conditions is often of considerable importance with respect to
their usefulness. For example, some changes in contact conditions can result in
“negative damping” (nonphysical energy source) for many time integrators, which can be
very undesirable.
It is possible to override the time integrator implied by the
application-type classification; for example, you can perform a moderate
dissipation dynamic analysis using the backward Euler integrator. Changing the
default integrator is not generally recommended but may be useful in special
cases.
Additional Control over Integrator Parameters
Additional user controls enable modifications to settings of numerical
parameters associated with the Hilber-Hughes-Taylor operator (see
Hilber,
Hughes, and Taylor (1977) for descriptions of the numerical parameters).
The default parameter settings depend on the specified application type, as
indicated in
Table 1
(see
Czekanski,
El-Abbasi, and Meguid (2001) for the basis of these settings).
Table 1. Default parameters for the Hilber-Hughes-Taylor integrator.
Parameter
Application
Transient Fidelity
Moderate Dissipation
–0.05
–0.41421
0.275625
0.5
0.55
0.91421
These parameters can be adjusted or modified individually if the Hilber-Hughes-Taylor operator
is being used. If the default settings of these parameters correspond to the transient
fidelity settings shown in Table 1
and you explicitly modify the parameter alone, the other parameters are adjusted automatically to and . This relation provides control of the numerical damping associated
with the time integrator while preserving desirable characteristics of the integrator.
The numerical damping grows with the ratio of the time increment to the period of
vibration of a mode. Negative values of provide damping; whereas results in no damping (energy preserving) and is exactly the
trapezoidal rule (sometimes called the Newmark -method, with and ). The setting provides the maximum numerical damping. It gives a damping ratio of
about 6% when the time increment is 40% of the period of oscillation of the mode being
studied. Allowable values of , , and are: , , .
Default Incrementation Schemes
Automatic time incrementation is used by default for nonlinear dynamic
procedures. The main factors used to control adjustments to the time increment
size for an implicit dynamic procedure are the convergence behavior of the
Newton iterations and the accuracy of the time integration. The time increment
size may vary considerably during an analysis. Details of the time increment
control algorithm depend on the type of dynamic application you are studying.
The following factors are considered by default in the time increment
control algorithm if you specify a quasi-static–type application (the same
factors control the time increment size for purely static analyses):
The time increment size is reduced if an increment appears to be
diverging or if the convergence rate is slow.
The time increment size is fairly aggressively increased if rapid
convergence occurs in previous increments.
Analyses for moderate dissipation-type applications also use these same
factors, as well as a default upper bound on the time increment size equal to
one-tenth of the step duration.
The following factors are considered by default in the time increment
control algorithm if you specify a transient fidelity–type application:
The time increment size is reduced if an increment appears to be
diverging or if the convergence rate is slow.
The time increment size is reduced if changes in contact status are
detected during the first attempt of processing an increment. The new increment
size is set such that the end of the increment corresponds to the average time
of the contact status changes that were detected with the previous increment
size. (In such cases an additional very small time increment is used to enforce
compatibility of velocities and accelerations across active contact
interfaces.)
The time increment size is reduced if the half-increment residual
(out-of-balance force) halfway through a time increment exceeds the
half-increment residual tolerance, which is 10,000 times the time average force
for a contact analysis or 1000 times the time average force for an analysis
without contact.
The time increment is gradually increased if rapid convergence occurs
in previous increments.
The upper bound for the time increment size is equal to 1/100 of the
step duration.
Intermittent Contact/Impact
The second and third factors described in the preceding list often result
in very small time increment sizes for contact simulations that are performed
as a transient fidelity application (and the time increment size tends to
remain small due to the fourth factor). This problem can be avoided by
specifying a different application type or by using more detailed user
controls, as discussed below.
General Settings for the Time Increment Controls
A high-level user control over which factors are considered by the time increment control
algorithm can be used to override the defaults implied by the specified application type
for the analysis. Regardless of the application type you have specified, you can enforce
time increment controls associated with either quasi-static applications or transient
fidelity applications.
Controlling the Half-Increment Residual
Controls associated with the half-increment residual tolerance are
provided for tuning the time incrementation. These controls are intended for
advanced users and typically do not need to be modified.
Controlling Incrementation Involving Contact
By default, specifying a transient fidelity application typically results in reduced time
increment sizes on changes in contact status. An extra time increment with a very small
size is subsequently performed to enforce compatibility of velocities and accelerations
across active contact interfaces. Direct user control over these incrementation aspects
is available.
Direct Time Incrementation
You may directly specify the time increment size to be used. This approach is not generally
recommended but may be useful in special cases. The analysis ends if convergence
tolerances are not satisfied within the maximum number of iterations allowed.
It is possible to ignore convergence tolerances: the solution to an
increment is accepted after the specified maximum number of iterations allowed
even if convergence tolerances are not satisfied. Ignoring convergence
tolerances can result in highly nonphysical results and is not recommended
except by analysts with a thorough understanding of how to interpret results
obtained this way.
Default Amplitude for Loads
Loads such as applied forces or pressures are ramped on by default if you have selected the
quasi-static application classification; such ramping tends to enhance robustness
because the load increment size is proportional to the time increment size. For example,
if the Newton iterations are unable to converge for a particular time increment size,
the automatic time incrementation algorithm reduces the time increment size and restart
the Newton iterations with a smaller load increment considered.
For the other application classifications the dynamic procedure applies loads with a step
function by default such that the full load is applied in the first increment of the
step (regardless of the time increment size) and the load magnitude remains constant
over each step. Thus, if the first increment is unable to converge with the original
time increment size, reducing the time increment does not reduce the load increment by
default. In some cases the convergence behavior still improves on reducing the time
increment because the regularizing effect of inertia on the integration operators is
inversely proportional to the square of the time increment size. See Defining an Analysis for more information on default amplitude types for the
various procedures and how to override the default.
The “Subspace Projection” Method
The alternative approach provided in Abaqus/Standard for nonlinear dynamic problems is the “subspace projection” method. See Subspace dynamics for the theory
behind this method. In this method the modes of the linear system are extracted in an
eigenfrequency extraction step (Natural Frequency Extraction) prior to the
dynamic analysis and are used as a small set of global basis vectors to develop the
solution. These modes include eigenmodes and, if activated in the eigenfrequency extraction
step, residual modes. The method works well when the system exhibits mildly nonlinear
behavior, such as small regions of plastic yielding or rotations that are not small but not
too large.
This method can be very effective. As with the other direct integration methods, it is more
expensive in terms of computer time than the modal methods of purely linear dynamic
analysis, but it is often significantly less expensive than the direct integration of all of
the equations of motion of the model. However, since the subspace projection method is based
on the modes of the system, it is not accurate if there is extreme nonlinear response that
cannot be modeled well by the modes that form the basis of the solution.
Selecting the Modes on Which to Project
You can select the modes of the system on which the subspace projection is performed. The mode
numbers can be listed individually, or they can be generated automatically. If you choose
not to select the modes, all modes extracted in the prior frequency extraction step,
including residual modes if they were activated, are used in the subspace projection.
Numerical Implementation
The subspace projection method is implemented in
Abaqus/Standard
using the explicit (central difference) operator to integrate the equations of
motion written in terms of the modes of the linear system. This integration
method is particularly effective here because the modes are orthogonal with
respect to the mass matrix so that the projected system always has a diagonal
mass matrix.
A fixed time increment is used: this increment is the smaller of the time
increment that you specify or 80% of the stable time increment, which is
for the linear system, where
is the highest circular frequency of the modes that are used as the basis of
the solution. The 80% factor is intended as a safety factor so that any
increase in this highest frequency caused by nonlinear effects is less likely
to cause the integration to become unstable. The 80% is rather arbitrary; in
some cases it may be nonconservative. You must monitor the response—for
example, the energy balance—to ensure that the time increment is not causing
instability. Instability is a concern if the nonlinearities can stiffen the
system significantly, although in many practical cases such stiffening effects
are more prominent in increasing the lower frequencies of the system than in
affecting the highest frequencies that are likely to be retained to represent
the dynamic behavior accurately.
Accuracy of the Subspace Projection Method
The effectiveness of the subspace projection method depends on the value of
the modes of the linear system as a set of global interpolation functions for
the problem, which is a matter of judgment on your part—the same sort of
judgment as required when deciding if a particular mesh of finite elements is
sufficient. The method is valuable for mildly nonlinear systems and for cases
where it is easy to extract enough modes that you can be confident that they
describe the system adequately.
If nonlinear geometric effects are considered in the subspace dynamics step,
it is possible to perform a dynamic simulation for some time, reextract the
modes on the current stressed geometry by using another frequency extraction
step, and then continue the analysis with the new modes as the subspace basis
system. This procedure can improve the accuracy of the method in some cases.
Material Damping
You can introduce Rayleigh damping, as explained in Material Damping. This damping
acts in addition to numerical damping associated with the time integrator (discussed
previously).
Initial Conditions
Initial Conditions
describes all of the available initial conditions. Initial velocities must be
defined in global directions regardless of the use of nodal transformations
(see
Transformed Coordinate Systems).
If initial velocities are specified at nodes for which displacement boundary conditions are also
specified, the initial velocities are ignored at these nodes. However, if a displacement
boundary condition refers to an amplitude curve with an analytically defined time variation
(that is, excluding the piecewise linear tabular and equally spaced definitions), Abaqus/Standard computes the initial velocity for the nodes involved in the boundary condition as the
time derivative (evaluated at time zero) of the analytic variation.
When initial velocities are specified for dynamic analysis, they should be consistent with all of
the constraints on the model, especially time-dependent boundary conditions. Abaqus/Standard ensures that initial velocities are consistent with boundary conditions and with
multi-point and equation constraints but does not check for consistency with internal
constraints such as incompressibility of the material. In case of a conflict, boundary
conditions and multi-point constraints take precedence over initial conditions.
Specified initial velocities are used in a dynamic step only if it is the
first dynamic step in an analysis. If a dynamic step is not the first dynamic
step and there is an immediately preceding dynamic step, the velocities from
the end of the preceding step are used as the initial velocities for the
current step. If a dynamic step is not the first dynamic step and the
immediately preceding step is not a dynamic step, zero initial velocities are
assumed for the current step.
Controlling Calculation of Accelerations at the Beginning of a Dynamic Step
By default, Abaqus/Standard calculates accelerations at the beginning of the dynamic step for transient fidelity
applications. You can choose to bypass these acceleration calculations, in which case Abaqus/Standard assumes that initial accelerations for the current step are zero unless there is an
immediately preceding dynamic step. If the immediately preceding step is also a dynamic
step, bypassing the acceleration calculations causes Abaqus/Standard to use the accelerations from the end of the previous step to continue the new step. It
is appropriate to bypass the acceleration calculations if the loading has not changed
suddenly at the start of the dynamic step, but it is not correct if the loading at the
beginning of the first increment is significantly different from that at the end of the
previous step. In cases where large loads are applied suddenly, high-frequency noise due
to the bypass of the acceleration calculations may greatly increase the half-increment
residual.
Boundary Conditions
Boundary conditions can be applied to any of the displacement or rotation
degrees of freedom (1–6), to warping degree of freedom 7 in open-section beam
elements, to fluid pressure degree of freedom 8 for hydrostatic fluid elements,
or to acoustic pressure degree of freedom 8 for acoustic elements (Boundary Conditions).
Amplitude references can be used to prescribe time-varying boundary
conditions in a direct-integration dynamic step. Default amplitude variations
are described in
Defining an Analysis.
In direct time integration dynamic analysis, when a node with a prescribed motion is used in an
equation constraint or a multi-point constraint to control the motion of another node, the
equation or multi-point constraint is imposed correctly for the displacement and velocity of
the dependent node. However, the acceleration is not rigorously transmitted to the dependent
node, which may cause some high-frequency noise.
In the subspace projection method it is not currently possible to specify nonzero boundary
conditions directly. Instead, acceleration boundary conditions can be approximated by using
appropriate combinations of large point masses and concentrated loads. At the node where
such a boundary condition is desired, attach a large point mass that is approximately
105–106 times larger than the mass of the original model. In
addition, a concentrated load of magnitude equal to the product between the large point mass
and the desired acceleration must be specified in the direction of the approximated boundary
condition. Since the point mass is significantly larger than the mass of the model, the big
mass–concentrated load combination approximates the desired acceleration in the specified
direction accurately. Boundary conditions other than accelerations must be converted into
acceleration histories before they can be approximated.
Loads
The following loads can be prescribed in a dynamic analysis:
Concentrated nodal forces can be applied to the displacement degrees of
freedom (1–6); see
Concentrated Loads.
Distributed pressure forces or body forces can be applied; see
Distributed Loads.
The distributed load types available with particular elements are described in
Abaqus Elements Guide.
The following predefined fields can be specified in a dynamic analysis, as
described in
Predefined Fields:
Although temperature is not a degree of freedom in stress/displacement elements, nodal
temperatures can be specified as a predefined field. Any difference between the applied
and initial temperatures causes thermal strain if a thermal expansion coefficient is
given for the material (Thermal Expansion). The
specified temperature also affects temperature-dependent material properties, if any.
The values of user-defined field variables can be specified. These
values only affect field-variable-dependent material properties, if any.
Material Options
Most material models that describe mechanical behavior are available for use
in a dynamic analysis. The following material properties are not active during
a dynamic analysis: thermal properties (except for thermal expansion), mass
diffusion properties, electrical conductivity properties, and pore fluid flow
properties.
Rate-Dependent Yield and Friction in Quasi-Static Response
You can control whether to consider or ignore the strain rate–dependence of the yield
stress and the slip rate–dependence of the friction coefficient within the step in a
quasi-static response.
Elements
Other than generalized axisymmetric elements with twist, any of the
stress/displacement elements in
Abaqus/Standard
(including those with temperature, pressure, and electrical potential degrees
of freedom) can be used in a dynamic analysis. Inertia effects are ignored in
hydrostatic fluid elements, and the inertia of the fluid in pore pressure
elements is not taken into account.
Output
In addition to the usual output variables available in
Abaqus/Standard
(see
Abaqus/Standard Output Variable Identifiers),
the following variables are provided specifically for implicit dynamic
analysis:
Variables for a specified element set or for the entire model:
XC
Current coordinates of the center of mass.
XCn
Coordinate n of the center of mass
().
UC
Displacement of the center of mass.
UCn
Displacement component n of the center of mass
().
URCn
Rotation component n of the center of mass.
VC
Equivalent rigid body velocity components.
VCn
Component n of the equivalent rigid body velocity
().
VRCn
Component n of the equivalent rigid body
angular velocity ().
HC
Angular momentum about the center of mass.
HCn
Component n of the angular momentum about the
center of mass ().
HO
Angular momentum about the origin.
HOn
Component n of the angular momentum about the
origin ().
RI
Rotary inertia about the origin.
RIij
-component
of the rotary inertia about the origin ().
MASS
Mass.
VOL
Current volume.
Input File Template
HEADING
…
BOUNDARYData lines to specify zero-valued boundary conditionsINITIAL CONDITIONSData lines to specify initial conditionsAMPLITUDE, NAME=nameData lines to define amplitude variations
**
STEP (,NLGEOM)
Once NLGEOM is specified, it will be active in all subsequent steps.DYNAMICData line to control automatic time incrementationBOUNDARYData lines to describe zero-valued or nonzero boundary conditionsCLOAD and/or DLOAD and/or INCIDENT WAVEData lines to specify loadsTEMPERATURE and/or FIELDData lines to prescribe predefined fieldsCECHARGE and/or DECHARGE (if electrical potential degrees of
freedom are active)
Data lines to specify charges END STEP
References
Czekanski, A., N. El-Abbasi, and S. A. Meguid, “Optimal
Time Integration Parameters for Elastodynamic Contact
Problems,” Communications in Numerical
Methods in
Engineering, vol. 17, pp. 379–384, 2001.
Hilber,
H.M., T. J. R. Hughes, and R. L. Taylor, “Improved
Numerical Dissipation for Time Integration Algorithms in Structural
Dynamics,” Earthquake Engineering and
Structural
Dynamics, vol. 5, pp. 283–292, 1977.