In a direct integration dynamic analysis, activated by the
*DYNAMIC keyword, the equation of motion is integrated in
time using the -method developed by Hilber, Hughes and Taylor
[62] (unless the massless contact method is selected; this method is treated
at the end of this section). The method is implemented exactly as described in
[23]. The parameter
lies in the interval [-1/3,0] and
controls the high frequency dissipation:
=0 corresponds to the
classical Newmark method inducing no dissipation at all, while
=-1/3
corresponds to maximum dissipation (default is
=-0.05). The user can choose between an implicit and explicit version of the algorithm. The implicit version (default) is unconditionally stable.
In the explicit version, triggered by the parameter EXPLICIT on the *DYNAMIC
keyword card, the mass matrix is lumped, and a forward integration scheme is
used so that the solution can be calculated without solving a system of
equations. This corresponds to section 2.11.5 in [23]. The mass
matrix only has to be set up once at the start of each step and no stiffness
matrix is needed. Indeed, the terms in equation (2.475) in
[23] in which the matrix is used correspond to the internal
forces. They can be calculated directly from the stresses without need to set
up the stiffness matrix. Therefore, each increment is much faster than with the implicit
scheme. Furthermore, in the explicit method no iterations are
performed, so each increment consists of exactly one iteration. However, the
explicit scheme is only conditionally stable: in the
-method, which in
CalculiX is also used for explicit calculations (unless the massless contact
method is selected), the maximum
time step
is dictated by:
![]() |
(523) |
where
is given by Equation (2.477) in [23] and
, which is the highest natural frequency of an element,
satisfies for volumetric elements:
![]() |
(524) |
where is the velocity of sound for the material at stake and
is the minimum height of the
element. For an isotropic material
satisfies [65]:
![]() |
(525) |
It corresponds to the wave speed of longitudinal waves, which are faster than
transversal waves, for which the speed amounts to
. For the special case of single crystal materials, which are anistropic
materials characterized by three independent elastic contants, the derivation
is more intricate [65]. Here, the derivation will be given for
general anisotropic materials.
First, the difference between the phase velocity and group velocity of a wave is explained. A one-dimensional wave is described by
![]() |
(526) |
in which is the wave number and
the angular frequency. For
constant values of
and
the wave has a constant amplitude for
constant or a wave speed satisfying:
![]() |
(527) |
Now, adding to this wave a wave with slightly different wave number and angular frequency one obtains:
where
![]() |
(529) |
and
![]() |
(530) |
The result in the last line of Equation (528) is the original wave
(
) with the so-called phase velocity
![]() |
(531) |
modulated by a wave whose velocity amounts to the so-called group velocity:
![]() |
(532) |
The function versus
is called the dispersion relationship. If
this relationship is linear the phase and group velocity coincide. If not,
they differ.
To obtain the phase velocity in three dimensions for an arbitrary material the expression for a three-dimensional planar wave:
![]() |
(533) |
where
is the displacement vector, A is the amplitude,
is the polarization vector and
is the wave
number vector, is substituted in the general homogeneous equation of motion:
![]() |
(534) |
where
is the stress tensor, a comma denotes
differentiation w.r.t. space and a dot denotes differentiation
w.r.t. time. Indeed, using the elastic relationship
![]() |
(535) |
the definition of linear strain
![]() |
(536) |
and the symmetry properties of the elasticity tensor
![]() |
(537) |
one obtains the following form for the equation of motion:
![]() |
(538) |
Sustitution of the wave equation now yields:
![]() |
(539) |
Setting
and
one obtains the eigenvalue problem
where
is the phase velocity we are looking for. Since (denoting
the matrix on the left hand side by
):
![]() |
(541) |
the left hand matrix is symmetric. Furthermore, for an arbitrary vector
one obtains:
![]() |
(542) |
since
(linear elastic energy) and
can be considered as the
components of a fictitious strain tensor
. This now means that
is positive definite. Therefore, the
eigenvalues are real and positive.
For an isotropic linear elastic material the elasticity tensor satisfies ([23]):
and component of the matrix
amounts to:
![]() |
(544) |
Therefore, can be thought of consisting of the following row vectors (or
column vectors, since it is a symmetric matrix):
![]() |
(545) |
Setting the determinant of the eigenvalue matrix to zero amounts to:
![]() |
(546) |
Multiplying out the different terms leads to:
![]() |
![]() |
![]() |
|
![]() |
![]() |
![]() |
(547) |
Using identities such as
finally leads to:
![]() |
(548) |
which leads to a double root and a single root
. Since
the corresponding phase velocities are
and
. The polarization
vectors are the eigenvectors of the system and are obtained by substituting
the eigenvalues in the eigenvalue system
. This leads to the
polarization vectors
and
for the double root and
for the single root. The former eigenvectors are orthogonal to
the propagation direction of the wave
and therefore these are
transversal waves, whereas the latter eigenvector is parallel to the wave
vector and corresponds to a longitudinal wave.
In order to obtain an expression for the group velocity of the wave the
derivative of the angular frequency w.r.t. the wave number vector
is needed. Multiplying the eigenvalue Equation [540] with the normalized polarization vector
yields:
![]() |
(549) |
Taking the derivative w.r.t. leads to:
![]() |
(550) |
or
![]() |
(551) |
This amounts to:
![]() |
(552) |
Since
the term
equals
, i.e. it is
identical with the first term in the equation. Therefore:
![]() |
(553) |
or
![]() |
(554) |
which amounts to:
The latter equation expresses the group velocity as a function of the polarization vector, the wave vector and the phase velocity.
Substituting the isotropic elasticity tensor from Equation [543] leads to:
![]() |
(556) |
For the longitudinal wave, knowing that
,
,
and
one gets
![]() |
(557) |
For the transversal waves
and
. Therefore:
![]() |
(558) |
Consequently, for an isotropic elastic material the phase and group velocity coincide. For an anisotropic material, such as a single crystal, this is not necessarily the case.
Now, coming back to the original question of a stable time step in an explicit
dynamic procedure it is clear that the maximal group speed is to be taken. For
an isotropic material this is the longitudinal wave speed. For an anisotropic
material the group speed may depend on the wave vector
. So for a given wave vector Equation
[540] has to be solved to yield the phase velocities
and the corresponding polarization vectors
. The latter ones
have to be substituted in Equation [555] to find the group
velocities. Now, the wave vector direction has to be varied so as to find the
maximum group velocity feasible for a material characterized by the specific
elastic tensor
. This velocity is then used to calculate
and
.
For the contact spring elements, the idealization of a spring with spring
constant connecting two nodes with nodal masses
and
is
used. For such a system one obtains:
![]() |
(559) |
where and
. For node-to-face penalty contact and
face-to-face penalty contact the nodal mass on the facial sides can be
obtained by using the shape functions at the spring location.
To accelerate explicit dynamic calculations mass scaling can be used [70]. It was introduced in CalculiX in the course of a Master Thesis [22]. Applying mass scaling the time increment can be increased, which reduces the overall computational time. In the mass scaling procedure the lumped mass matrix of each element is augmented by a fully populated matrix. The resulting computational increase is, however, neglegible, since the LU-decomposition has only to be performed once at the start of the calculation. Mass scaling is triggered by specifying the minimum time increment which the user wants to allow underneath the *DYNAMIC keyword (third parameter). If for any volumetric element the increment size calculated by CalculiX (based on the wave speed) is less than the minimum, the mass of this element is automatically scaled and redistributed such that the total mass of the element does not change. This is obtained by moving mass from the off-diagonal positions of the element mass matrix onto the diagonal. If any mass scaling takes place, a message is printed and the elements for which the mass was redistributed are stored in file “jobname_WarnElementMassScaled.nam”. This file can be read in any active cgx-session by typing “read jobname_WarnElementMassScaled.nam inp” and the elements can be appropriately visualized. For spring elements the mimimum time increment specified by the user is obtained by reducing the spring stiffness. CalculiX stores the maximum spring stiffness reduction to standard output.
Without a minimum time increment no mass scaling nor spring stiffness reduction is applied.
The following damping options are available:
For explicit dynamic calculations an additional hard contact formulation has been coded within a procedure characterized by:
From now on the method will be called the massless contact method. Its implementation closely follows the frictional flow diagram in [66]. From this publication it is clear that the contactless stiffness matrix is needed. The submatrix related to the contact degrees of freedom (master and slave) is used to set up and solve an inclusion problem in an implicit way. The other submatrices are used to calculate the right hand side of the dynamic equations. The left hand side of these equations is made up of a combination of the mass and damping matrix. It is assumed that these latter matrices do not change during the step, therefore, they can be factorized once at the beginning of the step.
If the NLGEOM parameter is not used on the *STEP card and if all materials are linear the stiffness matrix is calculated only once at the beginning of the step, else it is calculated in each increment which substantially increases the computational time.
Limitations right now include:
The method is triggered by the option TYPE=MASSLESS on all *CONTACT PAIR cards in a *DYNAMIC, EXPLICIT procedure only. In a *STATIC or *DYNAMIC implicit step the contact definition is replaced by TYPE=NODE TO SURFACE. Note that within one input deck only one type of *CONTACT PAIR is allowed. Since the contact is hard, the only parameter is the friction coefficient.
For all dynamic calculations (implicit dynamics, explicit dynamics with penalty contact or explicit dynamics with massless contact) a energy balance can be requested. For implicit dynamics this is done by default, for explicit dynamics the balance is calculated if the user has requested the output variable ENER underneath a *EL PRINT, *EL FILE or *ELEMENT OUTPUT keyword. This increases the computational time.