Dynamics

For a body \(\mathcal{B}\), the conservation of linear and angular momentum respectively read

\[\frac{d}{dt}\int_{\mathcal{B}}\rho\dot{\mathbf{x}}dv=\int_{\partial\mathcal{B}}\mathbf{t}da+\int_{\mathcal{B}}\rho\mathbf{b}dv\]
\[\frac{d}{dt}\int_{\mathcal{B}}\left(\mathbf{x}-\bar{\mathbf{x}}\right)\times\rho\dot{\mathbf{x}}dv= \int_{\partial\mathcal{B}}\left(\mathbf{x}-\bar{\mathbf{x}}\right)\times\mathbf{t}da+ \int_{\mathcal{B}}\left(\mathbf{x}-\bar{\mathbf{x}}\right)\times\rho\mathbf{b}dv\]

where \(t\) is time, \(\rho\) is the mass density, \(\dot{\mathbf{x}}\) is the point velocity, \(\mathbf{t}\) is the surface traction, \(\mathbf{b}\) is the body force, and \(\bar{\mathbf{x}}\) is a selected point. All of the mentioned quantities are spatial and so is the integration domain \(\mathcal{B}\), being the current configuration of the body. Sections below rephrase these balance principles, as required, for each of the kinematic models available in Solfec.

Rigid dynamics

Linear momentum balance reads

\[m\ddot{\bar{\mathbf{x}}}=\int_{\partial\mathcal{B}}\mathbf{t}da+\int_{\mathcal{B}}\rho\mathbf{b}dv\]

The spatial angular momentum balance reads

\[\frac{d}{dt}\left(\mathbf{j}\mathbf{\omega}\right)=\mathbf{j}\dot{\mathbf{\omega}}+\mathbf{\omega}\times\mathbf{j}\mathbf{\omega}=\mathbf{m}\]

The referential angular momentum balance reads

\[\mathbf{J}\dot{\mathbf{\Omega}}+\mathbf{\Omega}\times\mathbf{J}\mathbf{\Omega}=\mathbf{\Lambda}^{T}\mathbf{m}\]

where \(\mathbf{m}\) is the torque

\[\mathbf{m}=\int_{\partial\mathcal{B}}\left(\mathbf{x}-\bar{\mathbf{x}}\right)\times\mathbf{t}da+\int_{\mathcal{B}}\left(\mathbf{x}-\bar{\mathbf{x}}\right)\times\rho\mathbf{b}dv\]

and \(\mathbf{j}\) and \(\mathbf{J}\) are

\[\mathbf{j}=\mathbf{\Lambda}\mathbf{J}\mathbf{\Lambda}^{T},\,\,\,\mathbf{J}=\mbox{tr}\left(\mathbf{E}_{0}\right)\mathbf{I}-\mathbf{E}_{0}\]

the spatial and the referential inertia tensors, respectively. \(\mathbf{I}\) is the \(3\times3\) identity matrix. \(\mathbf{E}_{0}\) is the referential Euler tensor

(22)\[\mathbf{E}_{0}=\int_{\mathcal{B}_{0}}\left(\mathbf{X}-\bar{\mathbf{X}}\right)\otimes\left(\mathbf{X}-\bar{\mathbf{X}}\right)\rho_{0}dV\]

where \(\mathcal{B}_{0}\) is the reference configuration, \(\rho_{0}\) is the referential mass density, and the tensor product operator \(\otimes\) makes a \(3\times3\) matrix out of two 3-vectors, e.g. \(\mathbf{a}=\mathbf{b}\otimes\mathbf{c}\) means \(a_{ij}=b_{i}c_{j}\).

In matrix notation these balance principles read

\[\mathbf{M}\dot{\mathbf{u}}=\mathbf{f}\]

where

(23)\[\begin{split}\mathbf{M}=\left[\begin{array}{cc} \mathbf{J}\\ & m\mathbf{I} \end{array}\right]\end{split}\]

and

(24)\[\begin{split}\mathbf{f}=\left[\begin{array}{c} \mathbf{\Lambda}^{T}\int_{\partial\mathcal{B}}\left(\mathbf{x}-\bar{\mathbf{x}}\right)\times\mathbf{t}da+ \mathbf{\Lambda}^{T}\int_{\mathcal{B}}\left(\mathbf{x}-\bar{\mathbf{x}}\right)\times\rho\mathbf{b}dv-\mathbf{\Omega}\times\mathbf{J}\mathbf{\Omega}\\ \int_{\partial\mathcal{B}}\mathbf{t}da+\int_{\mathcal{B}}\rho\mathbf{b}dv \end{array}\right]\end{split}\]

Pseudo–rigid dynamics

Linear momentum balance reads

\[m\ddot{\bar{\mathbf{x}}}=\int_{\partial\mathcal{B}}\mathbf{t}da+\int_{\mathcal{B}}\rho\mathbf{b}dv\]

The referential angular momentum balance reads

\[\ddot{\mathbf{F}}\mathbf{E}_{0}+V\bar{\mathbf{P}}=\int_{\partial\mathcal{B}_{0}}\mathbf{t}\otimes \left(\mathbf{X}-\bar{\mathbf{X}}\right)dA+\int_{\mathcal{B}_{0}}\rho_{o}\mathbf{b}\otimes\left( \mathbf{X}-\bar{\mathbf{X}}\right)dV\]

where \(\mathbf{E}_{0}\) is defined in (22), \(V\) is the volume of the domain, and \(\bar{\mathbf{P}}\) is the average referential first Piola stress, defined as

\[\bar{\mathbf{P}}=\partial_{\mathbf{F}}\Psi\left(\mathbf{F}\right)\]
\[\Psi=\frac{1}{4}\left[\mathbf{F}^{T}\mathbf{F}-\mathbf{I}\right]:\mathbf{C}:\left[\mathbf{F}^{T}\mathbf{F}-\mathbf{I}\right]\]
\[C_{ijkl}=\lambda\delta_{ij}\delta_{kl}+\mu\left[\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}\right]\]

where the hyperelastic Saint Venant – Kirchhoff material is adopted. In the above \(\lambda\) and \(\mu\) are Lam’e constants, while \(\delta_{ij}\) is the Kronecker delta. The Lam’e constants can be expressed in terms of the Young modulus \(E\) and the Poisson ratio \(\nu\) as

\[\lambda=\frac{E\nu}{\left(1+\nu\right)\left(1-2\nu\right)}\]
\[\mu=\frac{E}{2+2\nu}\]

In matrix notation these balance principles read

\[\mathbf{M}\dot{\mathbf{u}}=\mathbf{f}\]

where

(25)\[\begin{split}\mathbf{M}=\left[\begin{array}{cccc} \mathbf{E}_{0}\\ & \mathbf{E}_{0}\\ & & \mathbf{E}_{0}\\ & & & m\mathbf{I} \end{array}\right]\end{split}\]

and

(26)\[\begin{split}\mathbf{f}=\left[\begin{array}{c} \int_{\partial\mathcal{B}_{0}}\mathbf{t}\otimes\left(\mathbf{X}-\bar{\mathbf{X}}\right)dA+\int_{\mathcal{B}_{0}}\rho_{o}\mathbf{b}\otimes\left(\mathbf{X}-\bar{\mathbf{X}}\right)dV-V\bar{\mathbf{P}}\\ \int_{\partial\mathcal{B}}\mathbf{t}da+\int_{\mathcal{B}}\rho\mathbf{b}dv \end{array}\right]\end{split}\]

It should be noted, that it is the row–wise composition of \(\dot{\mathbf{F}}\) in \(\mathbf{u}\) (cf. Kinematics), which allows us to use the computationally convenient block–diagonal form of \(\mathbf{M}\) for pseudo–rigid bodies.

Finite–element dynamics

The linear momentum balance reads

\[\mathbf{M}\dot{\mathbf{u}}+\mathbf{f}_{\text{int}}\left(\mathbf{q}\right)+ \mathbf{f}_{\text{damp}}\left(\mathbf{q},\mathbf{u}\right)=\mathbf{f}_{\text{ext}}\left(\mathbf{q}\right)\]

where \(\mathbf{M}\) is the mass matrix, \(\mathbf{f}_{\text{int}}\) is a vector of internal forces, \(\mathbf{f}_{\text{damp}}\) is a vector of damping forces, and \(\mathbf{f}_{\text{ext}}\) is a vector of external forces. These matrices and vectors are defined as follows

(27)\[\mathbf{M}=\int_{\mathcal{B}_{0}}\rho_{0}\mathbf{N}^{T}\mathbf{N}dV\]
(28)\[\mathbf{f}_{\text{int}}=\int_{\mathcal{B}_{0}}\left\{ \partial\mathbf{N}/\partial\mathbf{X}\right\} ^{T}\colon\mathbf{P}dV\]
(29)\[\mathbf{K}=\partial\mathbf{f}_{\text{int}}/\partial\mathbf{q}\]
\[\mathbf{f}_{\text{damp}}=\eta\mathbf{K}\]
(30)\[\mathbf{f}_{\text{ext}}=\int_{\mathcal{B}}\rho\mathbf{N}^{T}\mathbf{b}dv+\int_{\partial\mathcal{B}}\mathbf{N}^{T}\mathbf{t}da\]

where, except for \(\mathbf{f}_{\text{ext}}\), the so called total Lagrangian notation was used. The contraction of the strain matrix \(\mathbf{B}=\left\{ \partial\mathbf{N}/\partial\mathbf{X}\right\}^{T}\) with the first Piola stress tensor, \(\mathbf{B}:\mathbf{P}=\sum_{ij}B_{ij}P_{ij}\), creates the nodal components of the internal force vector. The derivative of the internal force with respect to displacements, \(\partial\mathbf{f}_{\text{int}}/\partial\mathbf{q}\), is customarily called the stiffness matrix, \(\mathbf{K}\). Stiffness proportional damping is used in Solfec, hence \(\mathbf{f}_{\text{damp}}=\eta\mathbf{K}\), where \(\eta\ge0\).

In matrix notation we simply have

\[\mathbf{M}\dot{\mathbf{u}}=\mathbf{f}\]

where

\[\mathbf{f}=\mathbf{f}_{\text{ext}}-\mathbf{f}_{\text{int}}-\mathbf{f}_{\text{damp}}\]

Implementation

Dynamics is implement in bod.c (rigid, pseudo–rigid) and fem.c (finite–element) files. Mass and stiffness matrices \(\mathbf{M}\) and \(\mathbf{K}\), and the damping factor \(\eta\), are declared in bod.h as follows:

struct general_body
{
  /* ... */

  MX *M;            /* inertia operator */

  MX *K;            /* stiffness operator */

  double damping;   /* stiffness proportional damping */

  /* ... */
}

Assembling of (23) is in bod.c:rig_dynamic_inverse.
Evaluation of (24) is in bod.c:rig_static_force.
Assembling of (25) is in bod.c:prb_dynamic_explicit_inverse.
Evaluation of (26) is in bod.c:prb_dynamic_force.
Assembling of diagonalized (27) is in fem.c:diagonal_inertia.
Evaluation of (28) is in fem.c:internal_force.
Assembling of (29) is in fem.c:tangent_stiffness.
Evaluation of (30) is in fem.c:external_force.