Direct integration dynamic analysis

In a direct integration dynamic analysis, activated by the *DYNAMIC keyword, the equation of motion is integrated in time using the $\alpha$-method developed by Hilber, Hughes and Taylor [65] (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 [24]. The parameter $\alpha$ lies in the interval [-1/3,0] and controls the high frequency dissipation: $\alpha$=0 corresponds to the classical Newmark method inducing no dissipation at all, while $\alpha$=-1/3 corresponds to maximum dissipation (default is $\alpha$=-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 [24]. 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 [24] in which the $[K]$ 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 $\alpha$-method, which in CalculiX is also used for explicit calculations (unless the massless contact method is selected), the maximum time step $\Delta t_{cr}$ is dictated by:

$$\Delta t_{cr} \approx \frac{\Omega _{cr}}{\omega_{max}},$$ (516)

where $\Omega _{cr}$ is given by Equation (2.477) in [24] and $\omega_{max}$, which is the highest natural frequency of an element, satisfies for volumetric elements:

$$\omega_{max} \approx 2 \frac{c}{h_{min}},$$ (517)

where $c$ is the velocity of sound for the material at stake and $h_{min}$ is the minimum height of the element. For an isotropic material $c$ satisfies [68]:

$$c= \sqrt{\frac{E(1-\nu )}{(1+\nu)(1-2 \nu) \rho } }= \sqrt{\frac{\lambda + 2 \mu }{\rho } } .$$ (518)

It corresponds to the wave speed of longitudinal waves, which are faster than transversal waves, for which the speed amounts to $\sqrt{\mu/\rho }$. For the special case of single crystal materials, which are anistropic materials characterized by three independent elastic contants, the derivation is more intricate [68]. 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

$$f=\cos(k x - \omega t),$$ (519)

in which $k$ is the wave number and $\omega$ the angular frequency. For constant values of $k$ and $\omega$ the wave has a constant amplitude for $k x -\omega t$ constant or a wave speed satisfying:

$$c=\frac{dx}{dt} = \frac{\omega }{k}.$$ (520)

Now, adding to this wave a wave with slightly different wave number and angular frequency one obtains:

$$g := \cos(k x - \omega t) + \cos[(k+\delta k) x - (\omega+ \delta \omega) t]$$

$$= \cos \left [ \left (A + \frac{\delta A}{2} \right ) - \frac{\delt... ...os \left [ \left (A + \frac{\delta A}{2} \right ) + \frac{\delta A}{2} \right ]$$

$$= 2 \cos (A+ \frac{\delta A}{2}) \cos \frac{\delta A}{2},$$ (521)

where

$$A := k x - \omega t,$$ (522)

and

$$\delta A := \delta k x - \delta \omega t.$$ (523)

The result in the last line of Equation (521) is the original wave ( $\cos (A+ \delta A/2) \approx \cos A$) with the so-called phase velocity

$$c_{ph} = \frac{\omega + \delta \omega/2}{k + \delta k/2} \approx \frac{\omega }{k}$$ (524)

modulated by a wave whose velocity amounts to the so-called group velocity:

$$c_g = \frac{\delta \omega }{\delta k} \approx \frac{\partial \omega }{\partial k}.$$ (525)

The function $\omega$ versus $k$ 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:

$$u_j = A p_j e^{i(k_m x_m - \omega t)},$$ (526)

where $\boldsymbol{u}$ is the displacement vector, A is the amplitude, $\boldsymbol{p}$ is the polarization vector and $\boldsymbol{k}$ is the wave number vector, is substituted in the general homogeneous equation of motion:

$$\sigma_{ij,i} = \rho \ddot{u}_j,$$ (527)

where $\boldsymbol{\sigma }$ 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

$$\sigma_{ij} = \Sigma_{ijkl} \epsilon_{kl},$$ (528)

the definition of linear strain

$$\epsilon_{kl}= (u_{k,l}+u_{l,k})/2,$$ (529)

and the symmetry properties of the elasticity tensor

$$\Sigma_{ijkl} = \Sigma_{jikl} = \Sigma_{ijlk} = \Sigma_{klij},$$ (530)

one obtains the following form for the equation of motion:

$$\Sigma_{ijkl} u_{k,li} = \rho \ddot{u}_j.$$ (531)

Sustitution of the wave equation now yields:

$$- \Sigma_{ijkl}k_l k_i p_k = - \rho \omega^2 p_j.$$ (532)

Setting $k_l := k n_l$ and $k_i := k n_i$ one obtains the eigenvalue problem

$$[\Sigma_{ijkl} n_l n_i] \{p_k \} = \rho c^2 \{p_j\},$$ (533)

where $c=\omega/k$ is the phase velocity we are looking for. Since (denoting the matrix on the left hand side by $C$):

$$C_{jk} = \Sigma_{ijkl} n_l n_i = \Sigma_{ljki} n_l n_i = \Sigma_{kilj} n_l n_i = \Sigma_{ikjl} n_l n_i = C_{kj}$$ (534)

the left hand matrix is symmetric. Furthermore, for an arbitrary vector $\boldsymbol{x}$ one obtains:

$$\Sigma_{ijkl} n_i x_j x_k n_l = \Sigma_{ijkl} x_i n_j x_k n_l \ge 0,$$ (535)

since $\Sigma_{ijkl} \epsilon_{ij} \epsilon_{kl} = \sigma_{kl} \epsilon_{kl} \ge 0$ (linear elastic energy) and $x_i n_j$ can be considered as the components of a fictitious strain tensor $\boldsymbol{x} \otimes \boldsymbol{n}$. This now means that $C$ is positive definite. Therefore, the eigenvalues are real and positive.

For an isotropic linear elastic material the elasticity tensor satisfies ([24]):

$$\Sigma_{ijkl}=\lambda \delta_{ij}\delta_{kl} + \mu (\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}),$$ (536)

and component $jk$ of the matrix $C$ amounts to:

$$C_{jk}= \Sigma_{ijkl}n_l n_i = (\lambda+\mu)n_j n_k + \mu \delta_{jk}.$$ (537)

Therefore, $C$ can be thought of consisting of the following row vectors (or column vectors, since it is a symmetric matrix):

$$C= \begin{bmatrix}(\lambda + \mu)n_1 \boldsymbol{n} + \mu \boldsy... ...e}_2 \\ (\lambda + \mu)n_3 \boldsymbol{n} + \mu \boldsymbol{e}_3 \end{bmatrix}.$$ (538)

Setting the determinant of the eigenvalue matrix to zero amounts to:

$$\vert C - \xi I \vert = [ (\lambda + \mu)n_1 \boldsymbol{n} + (\m... ...{e}_2] \times [(\lambda + \mu)n_3 \boldsymbol{n} + \mu \boldsymbol{e}_3] \} =0.$$ (539)

Multiplying out the different terms leads to:

$$(\lambda+\mu)n_1 (\mu-\xi)^2 \boldsymbol{n} \cdot(\boldsymbol{e}_2 \times \boldsymbol{e}_3) + (\lambda+\mu)n_2 (\mu-\xi)^2 \boldsymbol{e}_1 \cdot(\boldsymbol{n} \times \boldsymbol{e}_3)$$

$$+ \;\; (\lambda+\mu)n_3 (\mu-\xi)^2 \boldsymbol{e}_1 \cdot(\boldsymbol{e}_2 \times \boldsymbol{n}) + (\mu-\xi)^3 \boldsymbol{e}_1 \cdot(\boldsymbol{e}_2 \times \boldsymbol{e}_3)=0.$$ (540)

Using identities such as $\boldsymbol{e}_1 \cdot(\boldsymbol{e}_2 \times \boldsymbol{n})=\boldsymbol{n} \cdot(\boldsymbol{e}_1 \times \boldsymbol{e}_2)$ finally leads to:

$$(\lambda+\mu)(\mu-\xi)^2+(\mu-\xi)^3=0,$$ (541)

which leads to a double root $\xi=\mu$ and a single root $\xi=\lambda+2 \mu$. Since $\xi=\rho c^2$ the corresponding phase velocities are $c=\sqrt{\mu/\rho }$ and $c=\sqrt{(\lambda+2 \mu)/\rho }$. The polarization vectors are the eigenvectors of the system and are obtained by substituting the eigenvalues in the eigenvalue system $C = \xi I$. This leads to the polarization vectors $(n_2,-n_1,0)$ and $(n_3,0,-n_1)$ for the double root and $(n_1,n_2,n_3)$ for the single root. The former eigenvectors are orthogonal to the propagation direction of the wave $\boldsymbol{n}$ 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 $\omega$ w.r.t. the wave number vector $\boldsymbol{k}$ is needed. Multiplying the eigenvalue Equation [533] with the normalized polarization vector $p_j$ yields:

$$\Sigma_{ijkl} k_l k_i p_k p_j = \rho \omega^2.$$ (542)

Taking the derivative w.r.t. $k_q$ leads to:

$$\Sigma_{ijkl} \frac{\partial (k_l k_i)}{\partial k_q} p_k p_j = 2 \rho \omega \frac{\partial \omega }{\partial k_q},$$ (543)

or

$$\Sigma_{ijkl} [k_l \delta_{iq} + \delta_{lq} k_i] p_k p_j = 2 \rho \omega \frac{\partial \omega }{\partial k_q}.$$ (544)

This amounts to:

$$\Sigma_{qjkl} k_l p_k p_j + \Sigma_{ijkq} k_i p_k p_j = 2 \rho \omega \frac{\partial \omega }{\partial k_q}.$$ (545)

Since $\Sigma_{ijkq}= \Sigma_{kqij}= \Sigma_{qkij}= \Sigma_{qkji}$ the term $\Sigma_{ijkq} k_i p_k p_j$ equals $\Sigma_{qjkl} k_l p_j p_k$, i.e. it is identical with the first term in the equation. Therefore:

$$\Sigma_{qjkl} k_l p_k p_j = \rho \omega \frac{\partial \omega }{\partial k_q},$$ (546)

or

$$\Sigma_{ijkl} k_l p_k p_j = \rho \omega \frac{\partial \omega }{\partial k_i},$$ (547)

which amounts to:

$$(c_g)_i = \frac{\partial \omega }{\partial k_i} = \frac{\Sigma_{ijkl} p_j p_k n_l }{\rho c}.$$ (548)

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 [536] leads to:

$$\frac{\partial \omega }{\partial k_i} = \frac{1}{\rho c} [\lambda p_i p_k n_k + \mu p_i p_l n_l + \mu n_i ].$$ (549)

For the longitudinal wave, knowing that $c_L:=c=\sqrt{\lambda + 2 \mu /\rho }$, $\boldsymbol{n} \Vert \boldsymbol{p}$, $\Vert\boldsymbol{n}\Vert=1$ and $\Vert \boldsymbol{p} \Vert=1$ one gets

$$\frac{\partial \omega }{\partial k_i} = \frac{(\lambda+2 \mu )}{\rho c_L} n_i = c_L n_i.$$ (550)

For the transversal waves $c_T:=c=\sqrt{\mu/\rho }$ and $\boldsymbol{n} \perp \boldsymbol{p}$. Therefore:

$$\frac{\partial \omega }{\partial k_i} = \frac{\mu }{\rho c_T} n_i = c_T n_i.$$ (551)

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 $\boldsymbol{n}$. So for a given wave vector Equation [533] has to be solved to yield the phase velocities $c$ and the corresponding polarization vectors $\boldsymbol{p}$. The latter ones have to be substituted in Equation [548] 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 $\sigma_{ijkl}$. This velocity is then used to calculate $\omega_{max}$ and $\Delta t_{cr}$.

For the contact spring elements, the idealization of a spring with spring constant $k$ connecting two nodes with nodal masses $M_1$ and $M_2$ is used. For such a system one obtains:

$$\omega_{max}=\sqrt{\frac{k(m_1+m_2)}{m_1 m_2} }$$ (552)

where $m_1=M_1/2$ and $m_2=M_2/2$. 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 selective mass scaling can be used [73]. It was introduced in CalculiX in the course of a Master Thesis [23]. Applying selective mass scaling the time increment can be increased, which reduces the overall computational time. The selective mass scaling procedure starts from the lumped mass matrix, which for linear elements and one global coordinate direction looks like (examplary for a C3D8-element):

$$M=\frac{m^e}{8} \boldsymbol{I},$$ (553)

where $m^e$ is the total mass of the element. Now, this mass matrix is augmented by a multiple of itself after removal of the rigid body translational modes:

$$\Delta M = \frac{m^e}{8} (\boldsymbol{I} - \boldsymbol{n} \otimes \boldsymbol{n}),$$ (554)

where (writing $\boldsymbol{n}$ as a column vector)

$$\boldsymbol{n}^T = \{ 1\; 1\; 1\; 1\; 1\; 1\; 1\; 1 \}/\sqrt{8}.$$ (555)

This leads to:

$$\Delta M = \frac{m}{8} \begin{bmatrix}7/8 & -1/8 & \cdots & -1/8 \\ -1/8 & 7/8 & & \\ \vdots & & \ddots & \\ -1/8 & & & 7/8 \end{bmatrix}.$$ (556)

Now, the division by 8 in the denominator of the matrix is approximated by a division by 7 to yield 1 on the diagonal and a scaling factor $\beta$ is introduced:

$$\Delta M = \frac{\beta m}{7 \cdot 8} \begin{bmatrix}7 & -1 & \cdots & -1 \\ -1 & 7 & & \\ \vdots & & \ddots & \\ -1 & & & 7 \end{bmatrix}.$$ (557)

In conclusion, the diagonal elements of matrix $M$ are now modified by a factor $1+\beta$. Since the natural frequency $\omega \sim 1/\sqrt{m}$ and the time increment $\Delta t \sim 1/\omega$, one can conclude that

$$\frac{\Delta t_{\text{target}}}{\Delta t_{\text{stable}}} = \sqrt{\frac{m_{\text{target}}}{m} } = \sqrt{1+\beta },$$ (558)

from which $\beta=(\Delta t_{\text{target}}/\Delta t_{\text{stable}})^2-1$. For other element types “7” should be replaced by $n-1$ and “8” by $n$, where n is the number of nodes in the element. The resulting mass matrix is now not diagonal any more and a LU-decomposition has to be performed. Since the mass does not change during the calculation, this decomposition has to be performed only once at the start of the calculation.

In CalculiX, selective mass scaling is triggered by specifying the minimum time increment $\Delta t_{\text{target}}$ which the user is willing to tolerate underneath the *DYNAMIC keyword (third parameter). If for any volumetric element the increment size $\Delta t_{\text{stable}}$ calculated by CalculiX (based on the wave speed) is less than the minimum, the mass of this element is automatically selectively scaled. In that case 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. If the mass of too many elements is redistributed, the lowest eigenfrequencies might be affected, leading to wrong results. Therefore, the maximum allowed time increment should be limited by a maximum allowed change of the lowest eigenfrequency by, e.g. 1 %. It is the responsability of the user to check this by performing a *FREQUENCY calculation in which the desired time step is provided as the fourth entry underneath the *FREQUENCY card. In addition, the value of $\alpha$ can also be specified.

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 selective mass scaling nor spring stiffness reduction is applied.

The following damping options are available:

• Rayleigh damping by means of the *DAMPING keyword card underneath a *MATERIAL card. It assumes the damping matrix to be a linear combination of the problem's stiffness matrix and mass matrix. Although possibly defined for only one material, the coefficients of the linear combination apply to the whole model. For explicit calculations the damping matrix is allowed to be mass matrix proportional only.
• Dashpot damping by defining dashpot elements (cf. Section 6.2.39; for implicit dynamic calculations only).
• Contact damping by defining a contact damping constant and, optionally, a tangent fraction using the *CONTACT DAMPING keyword card (implicit dynamic calculations only).

For explicit dynamic calculations an additional hard contact formulation has been coded within a procedure characterized by:

• the formulation of the contact condition as a set-valued force law [102].
• a reduction of the dynamic equations to static equations (discarding the inertia) for all contact degrees of freedom (master and slave).
• a Verlet time integration scheme.

From now on the method will be called the massless explicit dynamics method. Its implementation closely follows the frictional flow diagram in [69]. 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:

• no nonzero boundary conditions are allowed.
• mass scaling cannot be used

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.

In the massless explicit dynamics method selective mass scaling is also used. It is triggered in exactly the same way as for the $\alpha$-method, i.e. by specifying a minimum time increment exceeding the stable time increment. In contrast to the $\alpha$-method the change of mass matrix is added to the unlumped mass. Selective mass scaling is particularly advantageous in the presence of few very small elements which lead to an extremely small stable time increment.

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.