Kinematics

A kinematic model defines an allowed type of deformation for an individual body. Solfec includes three kinematic models:

• rigid body
• pseudo–rigid body
• finite–element body

A kinematic model is selected by specifying the kind paramter of the BODY command.

Rigid body

The motion of a rigid body reads

(1)\[\mathbf{x}\left(\mathbf{X},t\right)=\mathbf{\Lambda}\left(t\right)\left(\mathbf{X}-\bar{\mathbf{X}}\right)+\bar{\mathbf{x}}\left(t\right)\]

where \(\mathbf{\Lambda}\left(t\right)\) is a \(3\times3\) rotation matrix, \(\bar{\mathbf{X}}\) is a selected referential point, and \(\bar{\mathbf{x}}\left(t\right)\) is a spatial point. Hence \(\bar{\mathbf{x}}\left(t\right)=\mathbf{x}\left(\bar{\mathbf{X}},t\right)\) is the motion of the selected point \(\bar{\mathbf{X}}\). The term \(\mathbf{\Lambda}\left(t\right)\left(\mathbf{X}-\bar{\mathbf{X}}\right)\) represents the rotation of \(\mathbf{X}\) about the point \(\bar{\mathbf{X}}\). Representing rotations, \(\mathbf{\Lambda}\) must be orthogonal: \(\mathbf{\Lambda}^{T}\mathbf{\Lambda}=\mathbf{I}\), where \(\mathbf{I}\) is the \(3\times3\) identity matrix. Thus, the rigidity condition follows \(\left\Vert \mathbf{x}-\bar{\mathbf{x}}\right\Vert =\left\Vert \mathbf{X}-\bar{\mathbf{X}}\right\Vert\), that is the length of a line segment before and after rotation is preserved.

The velocity of the spatial point reads

(2)\[\dot{\mathbf{x}}\left(\mathbf{X},t\right)=\mathbf{\Lambda}\left(t\right)\hat{\mathbf{\Omega}}\left(t\right) \left(\mathbf{X}-\bar{\mathbf{X}}\right)+\dot{\bar{\mathbf{x}}}\left(t\right)=\hat{\mathbf{\omega}}\left(t\right) \mathbf{\Lambda}\left(t\right)\left(\mathbf{X}-\bar{\mathbf{X}}\right)+\dot{\bar{\mathbf{x}}}\left(t\right)\]

where \(\mathbf{\omega}\) is the spatial angular velocity vector, \(\mathbf{\Omega}\) is the referential angular velocity vector, and \(\dot{\bar{\mathbf{x}}}\) is the velocity of the spatial image of \(\bar{\mathbf{X}}\). W note that

\[\mathbf{\omega}=\mathbf{\Lambda}\mathbf{\Omega}\]

The hat operator \(\hat{\mathbf{y}}\) creates an anti–symmetric \(3\times3\) matrix out of a 3–vector as follows

\[\begin{split}\hat{\mathbf{y}}=\left[\begin{array}{ccc} 0 & -y_{3} & y_{2}\\ y_{3} & 0 & -y_{1}\\ -y_{2} & y_{1} & 0 \end{array}\right]\end{split}\]

In vector notation the configuration \(\mathbf{q}\) and the velocity \(\mathbf{u}\) can be expressed as

(3)\[\begin{split}\mathbf{q}=\left[\begin{array}{c} \Lambda_{11}\\ \Lambda_{21}\\ ...\\ \bar{x}_{1}\\ \bar{x}_{2}\\ \bar{x}_{3} \end{array}\right],\,\,\,\mathbf{u}=\left[\begin{array}{c} \Omega_{1}\\ \Omega_{2}\\ \Omega_{3}\\ \dot{\bar{x}}_{1}\\ \dot{\bar{x}}_{2}\\ \dot{\bar{x}}_{3} \end{array}\right]\end{split}\]

The configuration \(\mathbf{q}\) has 12 components and the velocity \(\mathbf{u}\) has 6 components.

Pseudo–rigid body

The motion of a pseudo–rigid body reads

(4)\[\mathbf{x}\left(t\right)=\mathbf{F}\left(t\right)\left(\mathbf{X}-\bar{\mathbf{X}}\right)+\bar{\mathbf{x}}\left(t\right)\]

where \(\mathbf{x}\) is the current image of a referential point \(\mathbf{X}\), \(\mathbf{F}\) is a spatially homogeneous deformation gradient \(\left(\mathbf{F}=\partial\mathbf{x}/\partial\mathbf{X}\right), \bar{\mathbf{X}}\) is a selected referential point and \(\bar{\mathbf{x}}=\bar{\mathbf{x}}\left(t\right)\) is the current image of \(\bar{\mathbf{X}}\). Deformation gradient \(\mathbf{F}\) is an invertible and orientation preserving \(\left(\det\left(\mathbf{F}\right)>0\right)\) \(3\times3\) matrix.

The velocity of the spatial point reads

(5)\[\dot{\mathbf{x}}\left(t\right)=\dot{\mathbf{F}}\left(t\right)\left(\mathbf{X}-\bar{\mathbf{X}}\right)+\dot{\bar{\mathbf{x}}}\left(t\right)\]

where \(\dot{\mathbf{F}}\) is the velocity of the deformation gradient and \(\dot{\bar{\mathbf{x}}}\) is the velocity of the spatial image of \(\bar{\mathbf{X}}\).

In vector notation the configuration \(\mathbf{q}\) and the velocity \(\mathbf{u}\) can be expressed as

(6)\[\begin{split}\mathbf{q}=\left[\begin{array}{c} F_{11}\\ F_{12}\\ ...\\ \bar{x}_{1}\\ \bar{x}_{2}\\ \bar{x}_{3} \end{array}\right],\,\,\,\mathbf{u}=\left[\begin{array}{c} \dot{F}_{11}\\ \dot{F}_{12}\\ ...\\ \dot{\bar{x}}_{1}\\ \dot{\bar{x}}_{2}\\ \dot{\bar{x}}_{3} \end{array}\right]\end{split}\]

Both, the configuration \(\mathbf{q}\) and the velocity \(\mathbf{u}\) are 12–component vectors.

Finite–element body

The finite–element motion reads

(7)\[\mathbf{x}\left(\mathbf{X},t\right)=\mathbf{X}+\mathbf{N}\left(\mathbf{X}\right)\mathbf{q}\left(t\right)\]

where \(\mathbf{x}\) is the current image of the referential point \(\mathbf{X}\), \(\mathbf{N}\left(\mathbf{X}\right)\) is a matrix of shape functions, and \(\mathbf{q}\left(t\right)\) is a vector of nodal displacements.

The velocity of the spatial point reads

(8)\[\dot{\mathbf{x}}\left(\mathbf{X},t\right)=\mathbf{N}\left(\mathbf{X}\right)\dot{\mathbf{q}}\left(t\right)\]

where \(\dot{\mathbf{q}}\) is the vector of nodal \(x,y,z\) velocities.

In vector notation, the configuration \(\mathbf{q}\) and the velocity \(\mathbf{u}\) can be expressed as

(9)\[\begin{split}\mathbf{q}=\left[\begin{array}{c} q_{1x}\\ q_{1y}\\ q_{1z}\\ ...\\ q_{nx}\\ q_{ny}\\ q_{nz} \end{array}\right],\,\,\,\mathbf{u}=\dot{\mathbf{q}}\end{split}\]

Both, the configuration \(\mathbf{q}\) and the velocity \(\mathbf{u}\) have size \(3\times n\), where \(n\) is the number of nodes in a finite–element mesh.

The matrix \(\mathbf{N}\) in (7) and (8) has the following form

\[\begin{split}\mathbf{N}=\left[\begin{array}{ccccccc} N_{1} & & & ... & N_{n}\\ & N_{1} & & ... & & N_{n}\\ & & N_{1} & ... & & & N_{n} \end{array}\right]\end{split}\]

where \(N_{i}\) are nodal shape functions, juxtaposed from element shape functions meeting at coincident mesh nodes.

Shape functions

Fig. 4 Element types in Solfec: tetrahedron, pyramid, wedge, hexahedron.

Finite–element types available in Solfec are depicted in Fig. 4. Shape functions for those elements are summarised in Table 25.

Isoparametric mapping is used to go back and forth between the natural element coordinates \(\xi_{1}\), \(\xi_{2}\), \(\xi_{3}\), local to individual element domains, and the referential mesh coordinates \(\mathbf{Y}\)

\[\begin{split}\mathbf{Y}\left(\xi_{1},\xi_{2},\xi_{3}\right)=\mathbf{N}\left(\xi_{1},\xi_{2},\xi_{3}\right)\left[\begin{array}{c} \mathbf{Y}_{1}\\ ...\\ \mathbf{Y}_{n} \end{array}\right]\end{split}\]

where \(\mathbf{Y}_{i}\) are referential coordinates of mesh nodes. The function \(\mathbf{N}\left(\mathbf{X}\right)\) in (7) and (8) has the following form

\[\mathbf{N}\left(\mathbf{X}\right)=\mathbf{N}\left(\mathbf{Y}^{-1}\left(\mathbf{X}\right)\right)\]

Newton iterations are required to solve \(\mathbf{X}=\mathbf{Y}\left(\xi_{1},\xi_{2},\xi_{3}\right)\) and find \(\xi_{1}\), \(\xi_{2}\), \(\xi_{3}\) for a given referential point \(\mathbf{X}\).

Table 25 Finite–element shape functions.

Tetrahedron

\(N_{1}=1-\left(\xi_{1}+\xi_{2}+\xi_{3}\right)\) \(N_{2}=\xi_{1}\) \(N_{3}=\xi_{2}\) \(N_{4}=\xi_{3}\)

Pyramid

\(r= \xi_{1}\xi_{2}\xi_{3}/\left(1-\xi_{3}\right)\text{ if }\xi_{3}\ne1\text{ or }0\text{ otherwise}\) \(N_{1}=0.25\left(\left(1+\xi_{1}\right)\left(1+\xi_{2}\right)-\xi_{3}+r\right)\) \(N_{2}=0.25\left(\left(1-\xi_{1}\right)\left(1+\xi_{2}\right)-\xi_{3}-r\right)\) \(N_{3}=0.25\left(\left(1-\xi_{1}\right)\left(1-\xi_{2}\right)-\xi_{3}+r\right)\) \(N_{4}=0.25\left(\left(1+\xi_{1}\right)\left(1-\xi_{2}\right)-\xi_{3}-r\right)\) \(N_{5}=\xi_{3}\)

Wedge

\(N_{1}=0.5\left(1-\xi_{1}-\xi_{2}\right)\left(1-\xi_{3}\right)\) \(N_{2}=0.5\xi_{1}\left(1-\xi_{3}\right)\) \(N_{3}=0.5\xi_{2}\left(1-\xi_{3}\right)\) \(N_{4}=0.5\left(1-\xi_{1}-\xi_{2}\right)\left(1+\xi_{3}\right)\) \(N_{5}=0.5\xi_{1}\left(1+\xi_{3}\right)\) \(N_{6}=0.5\xi_{2}\left(1+\xi_{3}\right)\)

Hexahedron

\(N_{1}=0.125\left(1-\xi_{1}\right)\left(1-\xi_{2}\right)\left(1-\xi_{3}\right)\) \(N_{2}=0.125\left(1+\xi_{1}\right)\left(1-\xi_{2}\right)\left(1-\xi_{3}\right)\) \(N_{3}=0.125\left(1+\xi_{1}\right)\left(1+\xi_{2}\right)\left(1-\xi_{3}\right)\) \(N_{4}=0.125\left(1-\xi_{1}\right)\left(1+\xi_{2}\right)\left(1-\xi_{3}\right)\) \(N_{5}=0.125\left(1-\xi_{1}\right)\left(1-\xi_{2}\right)\left(1+\xi_{3}\right)\) \(N_{6}=0.125\left(1+\xi_{1}\right)\left(1-\xi_{2}\right)\left(1+\xi_{3}\right)\) \(N_{7}=0.125\left(1+\xi_{1}\right)\left(1+\xi_{2}\right)\left(1+\xi_{3}\right)\) \(N_{8}=0.125\left(1-\xi_{1}\right)\left(1+\xi_{2}\right)\left(1+\xi_{3}\right)\)

Implementation

Kinematic models are implement in bod.c (rigid, pseudo–rigid) and fem.c (finite–element) files. Configuration and velocity vectors are declared in bod.h as follows:

struct general_body
{
  enum {OBS, RIG, PRB, FEM} kind; /* obstacle, rigid, pseudo-rigid, finite element */

  /* ... */

  double *conf,    /* configuration */
         *velo;    /* velocity */

  /* ... */
}