김경종
30 KiB
Procedures
techniques of integration depends on the stability limit of the explicit scheme, on the ease with which the nonlinear equations can be solved for the implicit operator, on 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 on the size of the model.
In ABAQUS/Standard the time step for implicit integration can be chosen automatically on the basis of the "half-step residual," a concept introduced in Hibbitt and Karlsson (1979). By monitoring the values of equilibrium residuals at t + \Delta t / 2 once the solution at t + \Delta 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 2.1.1-1. The body force at a point, f , can be written as an externally prescribed body force, F, and a d'Alembert force:
\mathbf {f} = \mathbf {F} - \rho \ddot {\mathbf {u}},
where \rho 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} \mathbf {f} \cdot \delta \mathbf {v} d V = \int_ {V} \mathbf {F} \cdot \delta \mathbf {v} d V - \int_ {V} \rho \ddot {\mathbf {u}} \cdot \delta \mathbf {v} d V.
The d'Alembert term can be written more conveniently in terms of the reference volume and reference density, \rho _ { 0 } , as
\int_ {V _ {0}} \rho_ {0} \ddot {\mathbf {u}} \cdot \delta \mathbf {v} d V _ {0},
where 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 + \Delta t ) , and uÄ is calculated from the time integration operator. The interpolator approximates the displacement at a point as
\mathbf {u} = \mathbf {N} ^ {N} u ^ {N},
so that
\ddot {\mathbf {u}} = \mathbf {N} ^ {N} \ddot {u} ^ {N},
provided that { \bf N } ^ { N } are not displacement dependent. This is the case for most of the elements in ABAQUS--in those instances where it is not true (the Hermite cubic beams, B23 and B33), the form taken for the d'Alembert force terms is discussed in detail in Chapter 3, "Elements." With this interpolation assumption, the d'Alembert force term is
- \left(\int_ {V _ {0}} \rho_ {0} \mathbf {N} ^ {N} \cdot \mathbf {N} ^ {M} d V _ {0}\right) \ddot {u} ^ {M};
Procedures
that is, the consistent mass matrix times the accelerations of the nodal variables. The finite element approximation to equilibrium, Equation 2.1.1-2, is
Equation 2.4.1-1
M ^ {N M} \ddot {u} ^ {M} + I ^ {N} - P ^ {N} = 0,
where
M ^ {N M} = \int_ {V _ {0}} \rho_ {0} \mathbf {N} ^ {N} \cdot \mathbf {N} ^ {M} d V _ {0}
is the consistent mass matrix,
I ^ {N} = \int_ {V _ {0}} \pmb {\beta} ^ {N}: \pmb {\sigma} d V _ {0}
is the internal force vector, and
P ^ {N} = \int_ {S} \mathbf {N} ^ {N} \cdot \mathbf {t} d S + \int_ {V} \mathbf {N} ^ {N} \cdot \mathbf {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 uN .
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.
The implicit operator used for time integration of the dynamic problem is the operator defined by Hilber, Hughes, and Taylor (1978). This operator is a single parameter 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 2.4.1-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:
Equation 2.4.1-2
M ^ {N M} \ddot {u} ^ {M} | _ {t + \Delta t} + (1 + \alpha) (I ^ {N} | _ {t + \Delta t} - P ^ {N} | _ {t + \Delta t}) - \alpha (I ^ {N} | _ {t} - P ^ {N} | _ {t}) + L ^ {N} | _ {t + \Delta t} = 0,
where L ^ { N } | _ { t + \Delta 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:
Equation 2.4.1-3
Procedures
u | _ {t + \Delta t} = u | _ {t} + \Delta t \dot {u} | _ {t} + \Delta t ^ {2} \left(\left(\frac {1}{2} - \beta\right) \ddot {u} | _ {t} + \beta \ddot {u} | _ {t + \Delta t}\right)
and
\dot {u} | _ {t + \Delta t} = \dot {u} | _ {t} + \Delta t \big ((1 - \gamma) \ddot {u} | _ {t} + \gamma \ddot {u} | _ {t + \Delta t} \big),
Equation 2.4.1-4
with
\beta = \frac {1}{4} (1 - \alpha) ^ {2}, \quad \gamma = \frac {1}{2} - \alpha \quad \text { and } \quad - \frac {1}{3} \leq \alpha \leq 0.
Hilber, Hughes, and Taylor (1978) present cogent arguments for the use of Equation 2.4.1-2-Equation 2.4.1-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, \beta = 1 / 4 ) , while with \alpha = - 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, since each time step change introduces some slight noise or "ringing" into the solution, and a little numerical damping ( \alpha = - 0 . 0 5 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 all of the dynamic examples in the ABAQUS Example Problems Manual 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 {\phi} ^ {\alpha} | _ {t + \Delta t} = \dot {\phi} ^ {\alpha} | _ {t} + \Delta t \left[ \gamma \ddot {\phi} ^ {\alpha} | _ {t + \Delta t} + (1 - \gamma) \ddot {\phi} ^ {\alpha} | _ {t} \right],
where \dot { \phi } ^ { \alpha } is the angular velocity of the node and \ddot { \phi } ^ { \alpha } is its angular acceleration, both taken in the current direction of the ® body axis at time t or at time t + \Delta t ; \Delta t is the time increment.
In the global system this is
\dot {\pmb {\phi}} | _ {t + \Delta t} = \Delta t \gamma \ddot {\pmb {\phi}} | _ {t + \Delta t} + [ \mathbf {e} ^ {\alpha} | _ {t + \Delta t} \mathbf {e} ^ {\alpha} | _ {t} ] \cdot \left[ \dot {\pmb {\phi}} | _ {t} + \Delta t (1 - \gamma) \ddot {\pmb {\phi}} | _ {t} \right],
where { \bf e } ^ { \alpha } = { \bf e } ^ { \alpha } ( \phi ) are orthonormal base vectors defining the body axis system. Since these are orthonormal vectors, this can be rewritten as
Equation 2.4.1-5
Procedures
\dot {\pmb {\phi}} | _ {t + \Delta t} = \Delta t \gamma \ddot {\pmb {\phi}} | _ {t + \Delta t} + \Delta \mathbf {C} \cdot \left[ \dot {\pmb {\phi}} | _ {t} + \Delta t (1 - \gamma) \ddot {\pmb {\phi}} | _ {t} \right],
where \Delta \mathbf { C } is the incremental rotation matrix. Let \Delta \pmb { \theta } be the increment in rotation from time t to t + \Delta t . Then \Delta \mathbf { C } is
\Delta \mathbf {C} = \exp [ \Delta \hat {\pmb {\theta}} ],
where \Delta \hat { \pmb { \theta } } is the skew-symmetric matrix with axial vector \Delta \pmb { \theta } . . See ``Rotation variables,'' Section 1.3.1, 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
\Delta \theta^ {\alpha} = \Delta t \dot {\phi} ^ {\alpha} | _ {t} + \Delta t ^ {2} \left[ \left(\frac {1}{2} - \beta\right) \ddot {\phi} ^ {\alpha} | _ {t} + \beta \ddot {\phi} ^ {\alpha} | _ {t + \Delta t} \right].
Since \Delta \mathbf { C } \cdot \Delta \pmb { \theta } = \exp [ \Delta \hat { \pmb { \theta } } ] \cdot \Delta \pmb { \theta } = \Delta \pmb { \theta } , the components \Delta \theta ^ { \alpha } are the same relative to the body axes at time t \ \mathrm { o r } \ t + \Delta t . In the global system this is
\Delta \pmb {\theta} = \Delta t ^ {2} \beta \ddot {\pmb {\phi}} | _ {t + \Delta t} + \Delta \mathbf {C} \cdot \left[ \Delta t \dot {\pmb {\phi}} | _ {t} + \Delta t ^ {2} \left(\frac {1}{2} - \beta\right) \ddot {\pmb {\phi}} | _ {t} \right].
Solving for the unknown velocity and acceleration at time t + \Delta t , , the velocity update equation is
\dot {\pmb {\phi}} | _ {t + \Delta t} = \frac {\gamma}{\Delta t \beta} \Delta \pmb {\theta} + \Delta \mathbf {C} \cdot \left[ \left(1 - \frac {\gamma}{\beta}\right) \dot {\pmb {\phi}} | _ {t} + \Delta t \left(1 - \frac {\gamma}{2 \beta}\right) \ddot {\pmb {\phi}} | _ {t} \right],
and the acceleration update equation becomes
\ddot {\pmb {\phi}} | _ {t + \Delta t} = \frac {1}{\Delta t \gamma} \dot {\pmb {\phi}} | _ {t + \Delta t} + \Delta \mathbf {C} \cdot \left[ \left(1 - \frac {1}{\gamma}\right) \ddot {\pmb {\phi}} | _ {t} - \frac {1}{\Delta t \gamma} \dot {\pmb {\phi}} | _ {t} \right].
The automatic time stepping for dynamic problems in ABAQUS is based on the half-step residual first proposed in Hibbitt and Karlsson (1979). The concept is quite appealing intuitively. Satisfaction of Equation 2.4.1-1 at the end of each time step (or actually of the Hilber-Hughes-Taylor form, Equation 2.4.1-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-step residual is to calculate the equilibrium residual error (the left-hand side of Equation 2.4.1-1) at some intermediate time point (chosen as t + \Delta t / 2 ) and to assess the error in the dynamic response prediction by the magnitude of that error.
The half-step 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:
Procedures
\ddot {u} | _ {\tau} = (1 - \tau) \ddot {u} | _ {t} + \tau \ddot {u} | _ {t + \Delta t}, \quad 0 \leq \tau \leq 1.
Having already solved for the state at t + \Delta t , this equation, together with Newmark's formulae now written for the time interval from t \ \mathrm { t o } \ t + \Delta t , requires that
\begin{array}{l} \Delta u | _ {\tau} = \tau^ {3} \Delta u | _ {t + \Delta t} + \tau (1 - \tau^ {2}) \Delta t \dot {u} | _ {t} + \tau^ {2} (1 - \tau) \frac {\Delta t ^ {2}}{2} \ddot {u} | _ {t}, \\ \dot {u} | _ {\tau} = \frac {\gamma}{\beta \tau \Delta t} \Delta u | _ {\tau} + \left(1 - \frac {\gamma}{\beta}\right) \dot {u} | _ {t} + \left(1 - \frac {\gamma}{2 \beta}\right) \tau \Delta t \ddot {u} | _ {t}, \quad \mathrm{and} \\ \ddot {u} | _ {\tau} = \frac {1}{\beta \tau^ {2} \Delta t ^ {2}} \Delta u | _ {\tau} - \frac {1}{\beta \tau \Delta t} \dot {u} | _ {t} + \left(1 - \frac {1}{2 \beta}\right) \ddot {u} | _ {t}, \\ \end{array}
where
\Delta u | _ {t + \Delta t} = u | _ {t + \Delta t} - u | _ {t}
is the increment in displacement obtained for the time step, \Delta 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 + \Delta t} = M ^ {N M} \ddot {u} ^ {M} | _ {t + \Delta t} + (1 + \alpha) (I ^ {N} | _ {t + \Delta t} - P ^ {N} | _ {t + \Delta t}) - \alpha (I ^ {N} | _ {t} - P ^ {N} | _ {t}) + L ^ {N} | _ {t + \Delta t}.
The residual at the start of the step is
R ^ {N} | _ {t} = M ^ {N M} \ddot {u} ^ {M} | _ {t} + (1 + \alpha) (I ^ {N} | _ {t} - P ^ {N} | _ {t}) - \alpha (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 + \Delta t / 2 is defined as
\begin{array}{l} R ^ {N} | _ {t + \Delta t / 2} \stackrel {\mathrm{def}} {=} M ^ {N M} \ddot {u} ^ {M} | _ {t + \Delta t / 2} + (1 + \alpha) (I ^ {N} | _ {t + \Delta t / 2} - P ^ {N} | _ {t + \Delta t / 2}) \\ - \frac {1}{2} \alpha (I ^ {N} | _ {t} - P ^ {N} | _ {t} + I ^ {N} | _ {t ^ {-}} - P ^ {N} | _ {t ^ {-}}) + L ^ {N} | _ {t + \Delta t / 2}, \\ \end{array}
where the I ^ { N } \big | _ { t + \Delta t / 2 } ; etc. are computed for conditions at time ( t + \Delta t / 2 ) and
L ^ {N} | _ {t + \Delta t / 2} \stackrel {\mathrm{def}} {=} \frac {1}{2} (L ^ {N} | _ {t + \Delta t} + L ^ {N} | _ {t}).
The "half-step residual," R \rvert _ { t + \Delta t / 2 } , is defined as the magnitude of the largest entry in R ^ { N } \vert _ { t + \Delta t / 2 } and provides a measure of accuracy of the time-stepping solution.
Procedures
The motivation behind the calculation of the half-step 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 + \Delta t / 2 } \approx 0 . 1 P consistently, the time stepping solution has high accuracy;
if R \rvert _ { t + \Delta t / 2 } ¼ P consistently, the time stepping solution has moderately good accuracy;
if R | _ { t + \Delta t / 2 } \approx 1 0 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-step 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-step residual is the basis of the adaptive time incrementation scheme. If the half-step 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-step calculation shows the solution is coarse, the time step used in solution should be reduced. The algorithm is described in detail in ``Time integration accuracy in transient problems,'' Section 8.3.4 of the ABAQUS/Standard User's Manual. 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 \alpha = - 0 . 0 5 : 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.
2.4.2 Intermittent contact/impact
Many dynamic problems involve intermittent contact, often with severe impact occurring when the structures hit. The contact algorithms in ABAQUS/Standard are designed to handle such cases. An example is the pipe whip problem, in which the fluid escaping from a ruptured pipe causes pipe motion and possible impact with restraints along the pipe. Depending on the geometry of the pipe and the position of the postulated break, severe impact can occur since the pipe may acquire high velocity before hitting a motionless restraint.
Impact conditions with hard contact are modeled in ABAQUS/Standard with the assumption that, at
Procedures
the time of impact, the two impacting surfaces instantaneously acquire the same velocity in the direction of the impact. This "fully plastic" impact concept is essential to any discrete model, the accuracy of the physical representation of local effects being dependent on the spatial model adopted. It should be emphasized that the plastic impact assumption is local: it is assumed that, at impact, energy is dissipated by some mechanism that is not modeled and whose spatial and temporal scale is infinitesimal compared to the discrete model. In the case of a beam model of a pipe hitting some other structure, this would presumably be the local plastic deformation of the pipe wall where it impacts the other structure. The points that acquire the same velocity sometimes separate just after the impact: part of the dynamic contact algorithm is designed to handle this case.
The local plastic impact concept is quite different from the simple "assumed coefficient of restitution" methods sometimes used in scoping analyses in place of detailed analysis of the history of contact forces. The method discussed here is specifically designed to allow as much--or as little--of this local detail to be obtained, as required. The usual time integration procedures in ABAQUS maintain the energy balance in terms of the energy mechanisms of the discrete model, so the instantaneous jumps that occur in the velocity and acceleration at impact imply that some other system of equations must govern the solution during impact. Thus, we view the impact events as separate from the usual time stepping and develop a set of impulse equations that allow the propagation of the solution over these instants of time and, hence, provide initial conditions from which the normal time stepping solution can continue to the next such event.
To develop the governing equations for impact, assume that, at some time t _ { 0 } , a part I of the surface of two bodies, A and B , comes into contact. Denote the velocities and accelerations of corresponding parts just before impact as
\dot {u} _ {A I} ^ {-}, \dot {u} _ {B I} ^ {-}, \ddot {u} _ {A I} ^ {-}, \ddot {u} _ {B I} ^ {-}
and
\dot {u} _ {A I} ^ {+}, \dot {u} _ {B I} ^ {+}, \ddot {u} _ {A I} ^ {+}, \ddot {u} _ {B I} ^ {+}
just after impact.
The fully plastic impact assumption requires that corresponding points acquire the same velocity and acceleration in the direction of impact immediately after impact so that, at time t _ { 0 } ^ { + } :
Equation 2.4.2-1
\begin{array}{l} \mathbf {n} \cdot \dot {\mathbf {u}} _ {A I} = \mathbf {n} \cdot \dot {\mathbf {u}} _ {B I}, \\ \mathbf {n} \cdot \ddot {\mathbf {u}} _ {A I} = \mathbf {n} \cdot \ddot {\mathbf {u}} _ {B I}. \\ \end{array}
Here n is the current normal to the interface surface I.
Writing \Delta \dot { \mathbf { u } } = \dot { \mathbf { u } } ^ { + } - \dot { \mathbf { u } } ^ { - } to describe the velocity "jump" at t0 at a point, we have
Equation 2.4.2-2
Procedures
\mathbf {n} \cdot \left(\dot {\mathbf {u}} _ {A I} ^ {-} + \Delta \dot {\mathbf {u}} _ {A I}\right) = \mathbf {n} \cdot \left(\dot {\mathbf {u}} _ {B I} ^ {-} + \Delta \dot {\mathbf {u}} _ {B I}\right).
The component of the force per unit area between the bodies across I in the normal direction n, N ^ { I } , must satisfy
Equation 2.4.2-3
N _ {A} ^ {I} = - N _ {B} ^ {I} \quad \mathrm{for} \quad t \geq t _ {0}
N _ {A} ^ {I} = N _ {B} ^ {I} = 0 \quad \mathrm{for} \quad t < t _ {0}.
Since finite velocity jumps are occurring at the time of impact in infinitesimal time (compared to the time scale of the simulation), during the infinitesimal interval t _ { 0 } ^ { - } to t _ { 0 } ^ { + } , N ^ { I } dominates all other forces in the system except for the d'Alembert forces. This simplifies the virtual work equation during t _ { 0 } ^ { - } to t _ { 0 } ^ { + } to
\sum_ {A, B} \Bigl [ \int_ {V} \rho \ddot {\mathbf {u}} \cdot \delta \mathbf {u} d V \Bigr ] + \int_ {I} N _ {A} \mathbf {n} \cdot \delta \mathbf {u} _ {A} d S + \int_ {I} N _ {B} \mathbf {n} \cdot \delta \mathbf {u} _ {B} d S = 0.
Integrating from t _ { 0 } ^ { - } to t _ { 0 } ^ { + } ,
Equation 2.4.2-4
\sum_ {A, B} \left[ \int_ {V} \int_ {t _ {0} ^ {-}} ^ {t _ {0} ^ {+}} \rho \ddot {\mathbf {u}} \cdot \delta \mathbf {u} d V \right] + \int_ {I} \int_ {t _ {0} ^ {-}} ^ {t _ {0} ^ {+}} \left(N _ {A} \mathbf {n} \cdot \delta \mathbf {u} _ {A} + N _ {B} \mathbf {n} \cdot \delta \mathbf {u} _ {B}\right) d S = 0.
But n \mathbf { \delta } \cdot \delta \mathbf { u } _ { A } = \mathbf { n } \cdot \delta \mathbf { u } _ { B } at t _ { 0 } , and N _ { B } = - N _ { A } , so the second term in Equation 2.4.2-4 is zero. However, the constraint (Equation 2.4.2-2) must be satisfied by augmenting Equation 2.4.2-4 with a Lagrange multiplier term with H as the multiplier:
\sum_ {A, B} \left[ \int_ {V} \rho \Delta \dot {\mathbf {u}} \cdot \delta \mathbf {u} d V \right] + \int_ {I} \delta \left[ H \mathbf {n} \cdot \left(\dot {\mathbf {u}} _ {A I} ^ {-} + \Delta \dot {\mathbf {u}} _ {A I} - \dot {\mathbf {u}} _ {B I} ^ {-} - \Delta \dot {\mathbf {u}} _ {B I}\right) \right] d S.
The first term has been integrated over t _ { 0 } ^ { - } to t _ { 0 } ^ { + } to give the velocity jump term. Taking the variation in the second term and noting that n cannot rotate in the infinitesimal time interval because there is no discontinuity in displacement, we obtain
Equation 2.4.2-5
\begin{array}{l} \sum_ {A, B} \left[ \int_ {V} \rho \Delta \dot {\mathbf {u}} \cdot \delta \mathbf {u} d V \right] + \int_ {I} H \mathbf {n} \cdot \left(\delta \mathbf {u} _ {A I} - \delta \mathbf {u} _ {B I}\right) d S \\ + \int_ {I} \left(\Delta \dot {\mathbf {u}} _ {A I} - \Delta \dot {\mathbf {u}} _ {B I}\right) \cdot \mathbf {n} \delta H d S \\ = - \int_ {I} \left(\dot {\mathbf {u}} _ {A I} ^ {-} - \dot {\mathbf {u}} _ {B I} ^ {-}\right) \cdot \mathbf {n} \delta H d S. \\ \end{array}
This equation is the impulse condition, which can be solved for the velocity jump \Delta \dot { \mathbf { u } } at all nodes. The solution also provides the impulse per unit area, H: the time integral of the pressure between the surfaces over the infinitesimal time of impact. The equilibrium equation written at time t _ { 0 } ^ { + } and including the constraint that \mathbf { n } \cdot \left( \ddot { \mathbf { u } } _ { A I } - \ddot { \mathbf { u } } _ { B I } \right) = 0 can then be used to obtain the initial accelerations immediately after impact.
The two sets of equations (for the velocity jumps and the initial accelerations at t _ { 0 } ^ { + } ) require solution of the mass matrix, augmented by the constraint, with different right-hand sides for the initial acceleration and the velocity jump equations. When the elements attached to the nodes that impact have consistent mass matrices, the equations will give velocity and acceleration jumps throughout that part of the model and not just at the impacting nodes themselves. From these initial conditions at t _ { 0 } ^ { + } the usual time stepping equations can continue, augmented by the constraint that \begin{array} { r } { \int _ { I } \mathbf { n } \cdot \left( \dot { \mathbf { u } } _ { A } - \dot { \mathbf { u } } _ { A } \right) d S = 0 } \end{array} , which is imposed by a Lagrange multiplier. Since this multiplier represents the interface pressure, its value is monitored for possible separation: if a negative value is seen (that is, if a state of tension exists between A and B across I), the constraint is removed: this requires another solution for the initial accelerations just after removal of the constraint, although there will be no velocity jump at that time. Separation can occur immediately after impact. If so, the equilibrium equations must be solved again without the constraint to find the corresponding accelerations.
In any real problem, impact and separation will occur at some intermediate point in a time step. To accommodate this, ABAQUS/Standard first solves the time step by ignoring impact, then estimates (by linear interpolation) the average time of impact or separation of all points that change in the increment, t _ { I } , and again solves the increment to that time. All surface contact changes are assumed to take place at that interpolated time point, with the surface clearances adjusted as necessary. For reasonable time steps this geometric adjustment is slight.
Typically, the time step used in the solution following a severe impact event is one or two orders of magnitude smaller than that preceding the event. As the high frequency noise generated by impact dissipates (through plasticity and the artificial damping introduced by the parameter ® in the time integration operator), the time step is expanded by the automatic time stepping algorithm and the solution proceeds accordingly.
Soft contact conditions, modeled with the *SURFACE BEHAVIOR, PRESSURE-OVERCLOSURE option, treat impact as an "elastic" event that does not destroy any kinetic energy. The change in velocity is determined by the amount of interpenetration. Thus, with soft contact the standard implicit integration procedure is used and no impulse equations need be solved. Under truly high velocity impact conditions, this may result in nodes bouncing back immediately after impact. This may lead to excessive contact chattering, resulting in convergence problems and small time increments. For this reason, softened contact is not recommended in "true" impact calculations. However, in certain dynamic calculations where impact effects are not critical--such as sheet forming, drop forging, or rolling operations--soft contact can work well because cutbacks due to impact calculations are avoided. The soft contact constraint is enforced via Lagrange multipliers. If the soft contact constraint compatibility is not satisfied within the given tolerance, a severe discontinuity iteration is performed.
2.4.3 Subspace dynamics
An alternative approach provided in ABAQUS/Standard for analysis of nonlinear dynamic problems is the "subspace projection" method. This method uses direct, explicit integration of the dynamic equations of equilibrium written in terms of a vector space spanned by a number of eigenvectors. The eigenmodes of the system, extracted in a *FREQUENCY step prior to the dynamic analysis, are used as the global basis vectors. This method can be very effective for systems with mild nonlinearities that do not substantially change the mode shapes. The method cannot be used for contact problems. 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. This method is implemented in ABAQUS/Standard using the explicit (central difference) operator to integrate the equations of motion projected onto the eigenmodes. In this procedure the explicit integration method is particularly effective because the eigenmodes are orthogonal with respect to the mass matrix so that the projected system always has a diagonal mass matrix and because the stability limit is determined by the highest eigenfrequency associated with the modes used in the analysis and not by the highest eigenfrequency of the structure.
Let us write the finite element approximation to the virtual work equations as
Equation 2.4.3-1
\delta u ^ {N} \bigg (M ^ {N M} \ddot {u} ^ {M} + I ^ {N} - P ^ {N} \bigg) = 0,
where
M ^ {N M} = \int_ {V _ {0}} \rho_ {0} \mathbf {N} ^ {N} \cdot \mathbf {N} ^ {M} d V _ {0}
is the consistent mass matrix,
I ^ {N} = \int_ {V _ {0}} \pmb {\beta} ^ {N}: \pmb {\sigma} d V _ {0}
is the internal force vector, and
P ^ {N} = \int_ {S} \mathbf {N} ^ {N} \cdot \mathbf {t} d S + \int_ {V} \mathbf {N} ^ {N} \cdot \mathbf {F} d V
is the external force vector.
The nodal displacements, velocities, accelerations, and the variations in displacement are expressed in terms of eigenmodes:
Equation 2.4.3-2