7.10 Rigid-Deformable Coupling

In FEBio deformable bodies can be coupled with rigid bodies ¹. At these rigid-deformable interfaces, the coupling of nodal degrees of freedom of deformable elements that attach to rigid bodies requires a modification of the global stiffness matrix and residual vector. This section describes the coupling between rigid and deformable bodies.

The position of a node shared by any number of deformable finite elements is denoted by \(\mathbf{x}\) in the current configuration. If the node belongs to one or more deformable elements but is not connected to a rigid body, then \(\mathbf{x}\) is given in terms of the nodal displacement \(\mathbf{u}\) by \(\mathbf{x}=\mathbf{X}+\mathbf{u}\); the corresponding nodal virtual velocity is \(\delta\mathbf{v}\) and the linearization of \(\mathbf{x}\) along an incremental displacement is denoted by \(\Delta\mathbf{u}\).

The contribution to the virtual work of the nodal force \(\mathbf{f}^{a}\) at node \(a\) is given by

\[ \begin{equation} \delta G=\delta\mathbf{v}^{a}\cdot\mathbf{f}_{n+\alpha}^{a}\,,\label{eq:rdc-virtual-work} \end{equation} \]

whera \(\delta\mathbf{v}^{a}\) is the virtual velocity of node \(a\) and \(\mathbf{f}_{n+\alpha}^{a}\) is the global nodal force, evaluated at the intermediate time \(t_{n+\alpha}\) as \(\mathbf{f}_{n+\alpha}^{a}=\alpha\mathbf{f}_{n+1}^{a}+\left(1-\alpha\right)\mathbf{f}_{n}^{a}\). The linearization of this virtual work along the incremental displacement \(\Delta\mathbf{u}^{b}\) of node \(b\) is

\[ \begin{equation} D\delta G=-\alpha\delta\mathbf{v}^{a}\cdot\mathbf{K}^{ab}\cdot\Delta\mathbf{u}^{b},\label{eq:rdc-linearization-vw} \end{equation} \]

where \(\mathbf{K}^{ab}=\left(\partial\mathbf{f}^{a}/\partial\mathbf{x}^{b}\right)_{n+1}\) is the contribution to the global stiffness matrix from the interactions of the degrees of freedom of nodes \(a\) and \(b\).

Now we consider the cases when either node \(a\), or node \(b\), or both, are attached to a rigid body. Our objective is to determine how to modify the global residual vector and stiffness matrix to account for the coupling of deformable and rigid body degrees of freedom.

When node \(a\) is attached to rigid body \(a\), its position is given in terms of that rigid body's degrees of freedom by the general relation (6.3-18). The corresponding virtual velocity is given in (7.9-4), reproduced here as

\[ \begin{equation} \delta\mathbf{v}^{a}=\delta\mathbf{r}^{a}-\hat{\mathbf{z}}_{n+\alpha}^{a}\cdot\delta\boldsymbol{\theta}^{a}\,,\label{eq:rdc-rb-virtual-velocity} \end{equation} \]

where \(\mathbf{z}_{n+\alpha}^{a}=\boldsymbol{\Lambda}_{n+\alpha}^{a}\cdot\mathbf{Z}^{a}\) is the position of node \(a\) relative to the center of mass of rigid body \(a\), at the intermediate time \(t_{n+\alpha}\). Now, the contribution of the global nodal force \(\mathbf{f}_{n+\alpha}^{a}\) to \(\delta G\) must be modified from \eqref{eq:rdc-virtual-work} according to

\[ \begin{equation} \delta G=\delta\mathbf{v}^{a}\cdot\mathbf{f}_{n+\alpha}^{a}=\left[\begin{array}{cc} \delta\mathbf{r}^{a} & \delta\boldsymbol{\theta}^{a}\end{array}\right]\left[\begin{array}{c} \mathbf{f}_{n+\alpha}^{a}\\ \hat{\mathbf{z}}_{n+\alpha}^{a}\cdot\mathbf{f}_{n+\alpha}^{a} \end{array}\right]\,.\label{eq:rdc-rb-virtual-work} \end{equation} \]

In other words, the displacement degrees of freedom of node \(a\) should be eliminated from the global system of equations and replaced with the translation and rotation degrees of freedom of rigid body \(a\). The force vector \(\mathbf{f}_{n+\alpha}^{a}\) should be made to contribute to the translation degrees of freedom of the center of mass of rigid body \(a\), whereas the moment \(\mathbf{z}_{n+\alpha}^{a}\times\mathbf{f}_{n+\alpha}^{a}\) should contribute to the rotation degrees of freedom of the rigid body.

When node \(b\) is connected to rigid body \(b\), the incremental displacement \(\Delta\mathbf{u}^{b}\) should be replaced with the rigid body incremental motions,

\[ \begin{equation} \Delta\mathbf{u}^{b}=\Delta\mathbf{r}^{b}-\hat{\mathbf{z}}_{n+1}^{b}\cdot\Delta\boldsymbol{\theta}^{b}\,.\label{eq:rdc-rb-incremental-disp} \end{equation} \]

Now, the contribution \(\mathbf{K}^{ab}\) to the global stiffness matrix needs to be modified from \eqref{eq:rdc-linearization-vw} according to the three possible cases:

• Node \(a\) belongs to rigid body \(a\), node \(b\) belongs to flexible elements only,

  \[ \begin{equation} D\delta G=-\alpha\left[\begin{array}{cc} \delta\mathbf{r}^{a} & \delta\boldsymbol{\theta}^{a}\end{array}\right]\left[\begin{array}{c} \mathbf{K}^{ab}\\ \hat{\mathbf{z}}_{n+\alpha}^{a}\cdot\mathbf{K}^{ab} \end{array}\right]\left[\Delta\mathbf{u}^{b}\right]\,.\label{eq:rdc-arigid-bflex} \end{equation} \]

• Node \(b\) belongs to rigid body \(b\), node \(a\) belongs to flexible elements only,

  \[ \begin{equation} D\delta G=-\alpha\left[\delta\mathbf{v}^{a}\right]\left[\begin{array}{cc} \mathbf{K}^{ab} & -\mathbf{K}^{ab}\cdot\hat{\mathbf{z}}_{n+1}^{b}\end{array}\right]\left[\begin{array}{c} \Delta\mathbf{r}^{b}\\ \Delta\boldsymbol{\theta}^{b} \end{array}\right]\,.\label{eq:rdc-aflex-brigid} \end{equation} \]

• Node \(a\) belongs to rigid body \(a\), node \(b\) belongs to rigid body \(b\),

  \[ \begin{equation} D\delta G=-\alpha\left[\begin{array}{cc} \delta\mathbf{r}^{a} & \delta\boldsymbol{\theta}^{a}\end{array}\right]\left[\begin{array}{cc} \mathbf{K}^{ab} & -\mathbf{K}^{ab}\cdot\hat{\mathbf{z}}_{n+1}^{b}\\ \hat{\mathbf{z}}_{n+\alpha}^{a}\cdot\mathbf{K}^{ab} & -\hat{\mathbf{z}}_{n+\alpha}^{a}\cdot\mathbf{K}^{ab}\cdot\hat{\mathbf{z}}_{n+1}^{b} \end{array}\right]\left[\begin{array}{c} \Delta\mathbf{r}^{b}\\ \Delta\boldsymbol{\theta}^{b} \end{array}\right]\,.\label{eq:rdc-arigid-brigid} \end{equation} \]

---

1. Maker, B. N.. "Rigid bodies for metal forming analysis with NIKE3D." University of California, Lawrence Livermore Lab Rept, vol. UCRL-JC-119862, pp. 1-8 (1995). ↩