김경종
27 KiB
Procedures
\begin{array}{l} u ^ {N} = \sum_ {\alpha} \Phi_ {\alpha} ^ {N} q _ {\alpha}, \quad \dot {u} ^ {N} = \sum_ {\alpha} \Phi_ {\alpha} ^ {N} \dot {q} _ {\alpha}, \\ \ddot {u} ^ {N} = \sum_ {\alpha} \Phi_ {\alpha} ^ {N} \ddot {q} _ {\alpha}, \quad \delta u ^ {N} = \sum_ {\alpha} \Phi_ {\alpha} ^ {N} \delta q _ {\alpha}, \\ \end{array}
where q _ { \alpha } , \dot { q } _ { \alpha } , \ddot { q } _ { \alpha } , and \delta \boldsymbol { q } _ { o } represent generalized displacement, velocity, acceleration, and displacement variation, respectively. Substitution into the virtual work expression yields the formula for the acceleration associated with mode ® as
Equation 2.4.3-3
\ddot {q} _ {\alpha} = \frac {1}{m _ {\alpha}} \Phi_ {\alpha} ^ {N} \left(P ^ {N} - I ^ {N}\right) \qquad \mathrm{(nosumon} \alpha \mathrm{)},
where m _ { \alpha } is the generalized mass associated with mode ®:
m _ {\alpha} = \Phi_ {\alpha} ^ {N} M ^ {N M} \Phi_ {\alpha} ^ {M}.
From Equation 2.4.3-3 it is seen that the element residuals are projected onto the vector space spanned by the chosen number of eigenmodes. Having calculated the generalized acceleration for each mode, the generalized displacement and velocity are calculated with the central difference operator
Equation 2.4.3-4
\begin{array}{l} q _ {t + \Delta t} = q _ {t} + \ddot {q} _ {t} \Delta t ^ {2}, \\ \dot {q} _ {t + \Delta t} = \Delta q _ {t} / \Delta t + \ddot {q} _ {t} \Delta t / 2. \\ \end{array}
The nodal values for all kinematic variables are obtained using the formulæ in Equation 2.4.3-2.
When initial velocities are applied, either explicitly using the *INITIAL CONDITIONS option with TYPE=VELOCITY or implicitly by continuation of the previous dynamic step, the initial velocity vector \dot { u } _ { 0 } has to be projected into the eigenspace. This leads to an initial generalized velocity for the mode ® in the form
\dot {q} _ {\alpha} ^ {0} = \frac {1}{m _ {\alpha}} \Phi_ {\alpha} ^ {N} M ^ {N M} \dot {u} _ {0} ^ {M}.
As in standard explicit dynamic integration, the method is conditionally stable. The stability limit is determined by the highest eigenfrequency of the modes used in the analysis:
\Delta t _ {s t a b l e} = \frac {2}{\omega_ {\mathrm{max}}},
where \omega _ { \mathrm { m a x } } is the highest circular frequency of the eigenmodes that are used as the basis of the solution. Since we will generally use a relatively small number of modes, this stability limit is usually much less restrictive than the stability limit for standard explicit integration.
Procedures
Throughout the procedure a fixed time increment is used: the value is chosen as the smaller of the time increment specified by the user, or 80% of the stable time increment. The 80% factor is intended as a safety factor, so a small increase in the highest eigenfrequency caused by nonlinear effects will not cause the integration to become unstable.
The subspace dynamic procedure is activated by including the SUBSPACE parameter on the *DYNAMIC option. In addition, the VECTORS parameter can be used to specify the number of eigenvectors spanning the solution space. If this parameter is omitted, the number of vectors will be equal to the number of eigenvectors extracted in the *FREQUENCY procedure.
It is possible to perform a subspace dynamic simulation for some time and then reextract the modes based on the current, stressed geometry (by using another *FREQUENCY step), followed by continuation of the analysis with the new modes as the subspace basis system. This can improve the accuracy of the method in certain cases.
Note that the method is noniterative; hence, there are no tolerances required for the procedure.
2.4.4 Equivalent rigid body dynamic motion
It is often useful to obtain the equivalent rigid body motion of part of a model (or of the whole model): the position and translational velocity of the part's center of mass and its angular rotation and velocity about the same center of mass. ABAQUS/Standard provides such output, based on equivalent momentum. This section defines how these values are calculated.
Let V be the volume of a part for which the equivalent rigid body motion values are requested. The density of the part in its initial configuration is \rho _ { 0 } ( S _ { i } ) , where S _ { i , i } = 1 , 2 , 3 , are material coordinates in the part. The spatial position of a material particle in its initial configuration is \mathbf { X } ( S _ { i } ) and in the current configuration is \mathbf { x } ( S _ { i } ) , resulting in a displacement of \mathbf { u } ( S _ { i } ) . We wish to compute the current spatial position of the center of mass of the part, \mathbf { x } _ { \circ } . ; the translational velocity of an equivalent rigid body, \dot { \mathbf { x } } _ { \circ } ; the angular velocity of this equivalent rigid body, \omega _ { \circ } ; the translational displacement of an equivalent rigid body motion, \mathbf { u } _ { \circ } ; and the rotation of an equivalent rigid body motion around the center of mass, \phi . In these definitions an "equivalent rigid body" means a rigid body with the same mass distribution and the same translational and angular momentum as the actual deforming part in the current configuration.
For simplicity of notation we define some quantities. The mass of the part is
\mathcal {M} \stackrel {\mathrm{def}} {=} \int_ {V ^ {0}} \rho_ {0} d V ^ {0}.
Its first mass moment about the origin is
\mathcal {J} _ {\circ} \stackrel {\mathrm{def}} {=} \int_ {V ^ {0}} \rho_ {0} \mathbf {x} d V ^ {0}.
Its second mass moment about the origin is
Procedures
\pmb {\mathcal {I}} _ {\circ} \stackrel {\mathrm{def}} {=} \int_ {V ^ {0}} \rho_ {0} (\mathbf {x} \cdot \mathbf {x} \mathbf {I} - \mathbf {x x}) d V ^ {0},
where I is a unit matrix.
It is convenient to invoke the relation \mathbf { x } = N ^ { N } \mathbf { x } ^ { N } (the summation convention is assumed), where N ^ { N } are the finite element interpolation functions associated with each degree of freedom and \mathbf { x } ^ { N } is the vector of current nodal positions. We can now write
\boldsymbol {\mathcal {I}} _ {\circ} \stackrel {\mathrm{def}} {=} \int_ {V ^ {0}} \rho_ {0} \big (N ^ {N} N ^ {M} \mathbf {x} ^ {N} \cdot \mathbf {x} ^ {M} \mathbf {I} - N ^ {N} N ^ {M} \mathbf {x} ^ {N} \mathbf {x} ^ {M} \big) d V ^ {0}.
Recognizing that the primitive mass matrix is
M ^ {N M} \stackrel {\mathrm{def}} {=} \int_ {V ^ {0}} \rho_ {0} N ^ {N} N ^ {M} d V ^ {0},
we have
\pmb {\mathcal {I}} _ {\circ} \stackrel {\mathrm{def}} {=} M ^ {N M} (\mathbf {x} ^ {N} \cdot \mathbf {x} ^ {M} \mathbf {I} - \mathbf {x} ^ {N} \mathbf {x} ^ {M}).
We can immediately obtain
\mathbf {x} _ {\circ} = \mathcal {J} _ {\circ} / \mathcal {M}
and, by equating the translational momentum of the equivalent and the actual body,
\dot {\mathbf {x}} _ {\circ} = \int_ {V ^ {0}} \rho_ {0} \dot {\mathbf {x}} d V ^ {0} / \mathcal {M}.
The angular velocity of the part is defined by equating the angular momentum of the part and of the equivalent rigid body about the center of mass:
\int_ {V ^ {0}} \rho_ {0} (\mathbf {x} - \mathbf {x} _ {\circ}) \times \dot {\mathbf {x}} d V ^ {0} = \pmb {\mathscr {I}} \cdot \pmb {\omega} _ {\circ},
where
\boldsymbol {\mathcal {I}} = \boldsymbol {\mathcal {I}} _ {\circ} + \mathcal {M} \left(\mathbf {x} _ {\circ} \mathbf {x} _ {\circ} - \mathbf {x} _ {\circ} \cdot \mathbf {x} _ {\circ} \mathbf {I}\right)
is the second mass moment of the part about its center of mass.
ABAQUS uses the lumped mass formulation for low-order elements. As a consequence, the second mass moments of inertia can deviate from the theoretical values, especially for coarse meshes. Certain ABAQUS elements produce lumped or structural contributions to this second mass moment (rotary
Procedures
inertias) not shown in these equations.
This provides
\boldsymbol {\omega} _ {\circ} = \boldsymbol {\mathcal {I}} ^ {- 1} \cdot \left[ \boldsymbol {\mathcal {H}} _ {\circ} - \boldsymbol {\mathcal {J}} _ {\circ} \times \dot {\mathbf {x}} _ {\circ} \right],
where
\pmb {\mathcal {H}} _ {\circ} \stackrel {\mathrm{def}} {=} \int_ {V ^ {0}} \rho_ {0} \mathbf {x} \times \dot {\mathbf {x}} d V ^ {0} = M ^ {N M} \mathbf {x} ^ {N} \times \dot {\mathbf {x}} ^ {M}
is the angular momentum of the part about the origin.
The perceived translational motion of the center of mass in an equivalent rigid body motion is calculated as
\mathbf {u} _ {\circ} = \int_ {V ^ {0}} \rho_ {0} \mathbf {u} d V ^ {0} / \mathcal {M}.
The equivalent rigid body rotation of the part with respect to its center of mass requires some conceptual approximations as follows. Denote the relative positions of a material particle with respect to the center of mass in the undeformed configuration and in the deformed configuration X and x, respectively. Consider that the configurations are known and that the axis of rotation of the body is denoted by the unit vector p. A material particle sees such rotation relative to the center of mass as
\pmb {\phi} = \phi \mathbf {p} = \frac {\overline {{\mathbf {X}}} _ {p} \times \overline {{\mathbf {x}}} _ {p}}{| \overline {{\mathbf {X}}} _ {p} \times \overline {{\mathbf {x}}} _ {p} |} \tan^ {- 1} \left(\frac {| \overline {{\mathbf {X}}} _ {p} \times \overline {{\mathbf {x}}} _ {p} |}{\overline {{\mathbf {X}}} _ {p} \cdot \overline {{\mathbf {x}}} _ {p}}\right),
where the subscript p denotes the projection of a vector into a plane normal to p. We now generalize this concept by integration of the constituent parts. Define
\mathbf {a} = \int_ {V} \rho \overline {{\mathbf {X}}} _ {p} \times \overline {{\mathbf {x}}} _ {p} d V, \qquad b = \int_ {V} \rho \overline {{\mathbf {X}}} _ {p} \cdot \overline {{\mathbf {x}}} _ {p} d V.
The average Euler rotation then follows with the equation
\phi_ {a v} = \frac {\mathbf {a}}{| \mathbf {a} |} \tan^ {- 1} \left(\frac {| \mathbf {a} |}{b}\right).
These integrals are not easily calculated, but with the modifications below they can be expanded in such a way that the (initially unknown) current center of mass, x , only appears in products with the position of particles in the known current configuration.
A necessary condition for the validity of the intuitive generalization above is that if the part undergoes an arbitrary rigid body rotation, the formula returns the rotation. That can easily be proved as follows.
Procedures
Every material particle rotates exactly an angle \phi = \phi _ { a v } in such a way that
\mathbf {a} = \mathbf {p} \sin \phi \int_ {V} \rho | \overline {{\mathbf {X}}} _ {p} | | \overline {{\mathbf {x}}} _ {p} | d V, \qquad b = \cos \phi \int_ {V} \rho | \overline {{\mathbf {X}}} _ {p} | | \overline {{\mathbf {x}}} _ {p} | d V,
and, therefore,
{\frac {\mathbf {a}}{| \mathbf {a} |}} = \mathbf {p}, \qquad {\frac {| \mathbf {a} |}{b}} = \tan \phi .
In all of these equations the direction vector p is unknown. To determine p we consider the characteristics of the displacement field of a rigid body rotation. For such a field, 1) the displacement of a particle is orthogonal to the rotation vector, and 2) the displacement is orthogonal to the position vector at half the motion. In a deformable body context we try to determine p by forcing these two statements to be true in an average sense. Considering that we are looking strictly at a rotation with respect to the center of mass, its definition automatically ensures that the first statement is satisfied. The second condition can be written in the form
\int_ {V} \rho \left[ \mathbf {p} \times (\overline {{\mathbf {x}}} + \overline {{\mathbf {X}}}) \right] \times [ \overline {{\mathbf {x}}} - \overline {{\mathbf {X}}} ] d V = 0,
which is a homogeneous set of equations in the components of p with coefficients made out of integrals of known quantities, from which p can be solved. We can then calculate the projection of the old and new position onto the plane normal to p with
\overline {{\mathbf {X}}} _ {p} = \overline {{\mathbf {X}}} - \mathbf {p p} \cdot \overline {{\mathbf {X}}}, \quad \overline {{\mathbf {x}}} _ {p} = \overline {{\mathbf {x}}} - \mathbf {p p} \cdot \overline {{\mathbf {x}}},
and by simple substitution one then obtains
\mathbf {a} = \int_ {V} \rho \left[ \overline {{\mathbf {X}}} - \left(\overline {{\mathbf {X}}} \cdot \mathbf {p}\right) \mathbf {p} \right] \times \left[ \overline {{\mathbf {x}}} - \left(\overline {{\mathbf {x}}} \cdot \mathbf {p}\right) \mathbf {p} \right] d V = \mathbf {p} \mathbf {p} \cdot \int_ {V} \rho \overline {{\mathbf {X}}} \times \overline {{\mathbf {x}}} d V,
where the vector \int _ { V } \rho \mathbf { \overline { { X } } } \times \mathbf { \overline { { x } } } d V is easily calculated from available quantities. With p known, a becomes determined. The quantity b can be calculated using with the same expressions, which yields
b = \int_ {V} \rho [ \overline {{\mathbf {X}}} - (\overline {{\mathbf {X}}} \cdot \mathbf {p}) \mathbf {p} ] \cdot [ \overline {{\mathbf {x}}} - (\overline {{\mathbf {x}}} \cdot \mathbf {p}) \mathbf {p} ] d V = (\mathbf {I} - \mathbf {p p}): \int_ {V} \rho \overline {{\mathbf {X}}} \overline {{\mathbf {x}}} d V.
Once a and b are known, Á is readily determined.
The determination of the equivalent rigid body rotation is based on average particle translations.
Rotational degrees of freedom are ignored in the calculation of this variable; it is assumed that such rotations will produce motions of points that will measurably contribute to the calculation. However, it is possible to find pathological cases in which that would not be the case; for instance, if the part of the model considered consists of rotary inertia elements only, the calculated average rigid body rotation
will be calculated as zero, even if the elements have indeed rotated.
2.4.5 Explicit dynamic analysis
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 {\mathbf {u}} ^ {(i + \frac {1}{2})} = \dot {\mathbf {u}} ^ {(i - \frac {1}{2})} + \frac {\Delta t ^ {(i + 1)} + \Delta t ^ {(i)}}{2} \ddot {\mathbf {u}} ^ {(i)},
\mathbf {u} ^ {(i + 1)} = \mathbf {u} ^ {(i)} + \Delta t ^ {(i + 1)} \dot {\mathbf {u}} ^ {(i + \frac {1}{2})},
where u_ is velocity and uÄ is acceleration. The superscript (i) refers to the increment number and \begin{array} { r } { i - \frac { 1 } { 2 } } \end{array} and \begin{array} { r } { i + \frac { 1 } { 2 } } \end{array} refer to midincrement values. The central difference integration operator is explicit in that the kinematic state can be advanced using known values of \dot { \mathbf { u } } ^ { ( i - \frac { 1 } { 2 } ) } and \ddot { \mathbf { 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 {\mathbf {u}} ^ {(i)} = \mathbf {M} ^ {- 1} \cdot (\mathbf {F} ^ {(i)} - \mathbf {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 { \mathbf { u } } ^ { ( i + \frac { 1 } { 2 } ) } , \dot { \mathbf { 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 {\mathbf {u}} ^ {(i + 1)} = \dot {\mathbf {u}} ^ {(i + \frac {1}{2})} + \frac {1}{2} \Delta t ^ {(i + 1)} \ddot {\mathbf {u}} ^ {(i + 1)}.
The central difference operator is not self-starting because the value of the mean velocity \dot { \mathbf { 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 {\mathbf {u}} ^ {(+ \frac {1}{2})} = \dot {\mathbf {u}} ^ {(0)} + \frac {\Delta t ^ {(1)}}{2} \ddot {\mathbf {u}} ^ {(0)}.
Substituting this expression into the update expression for \dot { \mathbf { u } } ^ { ( i + \frac { 1 } { 2 } ) } yields the following definition of \dot { \mathbf { u } } ^ { ( - \frac { 1 } { 2 } ) } :
\dot {\mathbf {u}} ^ {(- \frac {1}{2})} = \dot {\mathbf {u}} ^ {(0)} - \frac {\Delta t ^ {(0)}}{2} \ddot {\mathbf {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
\Delta t \leq \frac {2}{\omega_ {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
\Delta t \leq \frac {2}{\omega_ {m a x}} (\sqrt {1 + \xi^ {2}} - \xi),
where \xi 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:
\Delta t = \frac {2}{\omega_ {m a x} ^ {e l e m e n t}},
\omega _ { m a x } ^ { e l e m e n t }
Procedures
increment is given by the minimum taken over all the elements. The above stability limit can be rewritten as
\Delta t = \min (\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 B _ { i I } is the element gradient operator. 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 { \lambda } and \hat { \mu } , are determined in the following manner. We define \Delta p as the increment in the equivalent pressure stress, \begin{array} { r } { p = - \frac { 1 } { 3 } \mathrm { t r a c e } ( \pmb { \sigma } ) } \end{array} , \Delta \mathbf { S } as the increment in the deviatoric stress, \Delta \epsilon _ { v o l } as the increment of volumetric strain, and \Delta \mathbf { e } as the deviatoric strain increment. We assume a hypoelastic stress-strain rule of the form
\Delta p = - \hat {K} \Delta \epsilon_ {v o l},
\Delta \mathbf {S} = 2 \hat {\mu} \Delta \mathbf {e},
where \hat { K } is the effective bulk modulus. The effective moduli can then be computed as
\hat {K} = \frac {- \Delta p}{\Delta \epsilon_ {v o l}},
\hat {\mu} = \frac {1}{2} \frac {\Delta \mathbf {S} : \Delta \mathbf {e}}{\Delta \mathbf {e} : \Delta \mathbf {e}},
\hat {\lambda} = \hat {K} - \frac {2}{3} \hat {\mu}.
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, \lambda _ { o } and 2 \mu _ { 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 {\mu} = \frac {1}{2} (3 (\lambda_ {o} + 2 \mu_ {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 {\lambda} + 2 \hat {\mu})}{\rho}}.
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} \rho c _ {d} L _ {e} \dot {\epsilon} _ {v o l},
where b _ { 1 } is a damping coefficient (default=.06), \rho is the current material density, c _ { d } is the current dilatational wave speed, L _ { e } is an element characteristic length, and \dot { \epsilon } _ { 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} = \rho (b _ {2} L _ {e} \dot {\epsilon} _ {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
\xi = b _ {1} - b _ {2} 2 \frac {L _ {e}}{c _ {d}} \mathrm{min} (0, \dot {\epsilon} _ {v o l}).
Linear bulk viscosity is always included in ABAQUS/Explicit. The parameters b _ { 1 } and b _ { 2 } can be redefined by using the *BULK VISCOSITY option. The default values are b _ { 1 } = 0 . 0 6 and b _ { 2 } = 1 . 2 . *BULK VISCOSITY defines history data and 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}}{1 2} \rho c _ {d} L \frac {\Delta \kappa}{\Delta t},
where b _ { 1 } is a damping coefficient (default = 0.06), h _ { 0 } is the original thickness, \rho 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 \Delta \kappa = \Delta \kappa _ { 1 1 } + \Delta \kappa _ { 2 2 } 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 {\lambda} + 2 \hat {\mu})}{\rho}}.
The resultant pressure moment mh, where h is the current thickness, is added to the direct components of the moment resultant.
2.5 Modal dynamics
2.5.1 Eigenvalue extraction
There are many important areas of structural analysis in which it is essential to be able to extract the eigenvalues of the system and, hence, obtain its natural frequencies of vibration or investigate possible bifurcations that may be associated with kinematic instabilities. For example, structural evaluation for seismic events is often based on linear analysis, using the structure's modes up to a limiting cut-off frequency, which is usually taken as 33 Hz (cycles/second). Once the modes are available, their orthogonality property allows the linear response of the structure to be constructed as the response of a number of single degree of freedom systems. This opens the way to several response evaluation methods that are computationally inexpensive and provide useful insight into the structure's dynamic