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¶

$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¶

$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¶

$\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.