In some situations it is useful to connect two non-conforming meshes together. This can be done by defining a tied contact interface. In FEBio, the tied contact works very similar to the sliding contact interface. We need to define a slave surface and a master surface, where it is assumed that the slave surface nodes will be tied to the master surface faces.
Gap Function
Just as in sliding contact, we need to define a gap function that measures the distance between the slave and master surface. In order to do that, we first define the projection of a slave node to the master surface.
\[ \begin{equation} \mathbf{\bar{Y}}\left(\mathbf{X}\right)=\arg\min\limits_{\mathbf{Y}\in\Gamma^{\left(2\right)}}\left\Vert \mathbf{X}-\mathbf{Y}\right\Vert \,.\label{eq679} \end{equation} \]
This definition is similar to that of the sliding interface, except that now the projection is done in the material reference frame. This implies that the projection only needs to be calculated once, at the beginning of the analysis. We can now proceed to the definition of the gap function.
\[ \begin{equation} \mathbf{g}\left(\mathbf{X}\right)=\boldsymbol{\varphi}^{\left(1\right)}\left(\mathbf{X}\right)-\boldsymbol{\varphi}^{\left(2\right)}\left(\mathbf{\bar{Y}}\left(\mathbf{X}\right)\right)\,.\label{eq680} \end{equation} \]
An important observation is that the gap function is now a vector quantity since the gap needs to be closed in all direction, not just the normal direction as is the case in sliding contact.
With the definition of the gap function at hand (equation \eqref{eq680}), we can define the contribution to the virtual work equation from the tied contact reaction forces.
\[ \begin{equation} W_{t}=\int\limits_{\Gamma_{c}}\mathbf{T}\cdot\delta\mathbf{g}\,d\Gamma\,.\label{eq681} \end{equation} \]
Here, \(\mathbf{T}\) is the reaction force that enforces the constraint \(\mathbf{g}\left(\mathbf{X}\right)=0\). Since we anticipate the use of an augmented Lagrangian formalism, we can write this reaction force as follows,
\[ \begin{equation} \mathbf{T}=\boldsymbol{\lambda}+\varepsilon\mathbf{g}\,.\label{eq682} \end{equation} \]
The vector quantity \(\boldsymbol{\lambda}\) is the Lagrangian multiplier and \(\varepsilon\) is a penalty factor.
Since equation \eqref{eq681} is nonlinear we need to calculate the linearization. For tied contact, this is simply given by the following equation.
\[ \begin{equation} \Delta W_{t}=\int\limits_{\Gamma_{c}}\varepsilon\Delta\mathbf{g}\cdot\delta\mathbf{g}\,d\Gamma\,.\label{eq683} \end{equation} \]
Where
\[ \begin{equation} \delta\mathbf{g}=\mathbf{w}^{\left(1\right)}-\mathbf{w}^{\left(2\right)}\label{eq684} \end{equation} \]
and
\[ \begin{equation} \Delta\mathbf{g}=\Delta\boldsymbol{\varphi}^{\left(1\right)}\left(\mathbf{X}\right)-\Delta\boldsymbol{\varphi}^{\left(2\right)}\left(\mathbf{\bar{Y}}\left(\mathbf{X}\right)\right)\,.\label{eq685} \end{equation} \]
We also introduced the notation \(\mathbf{w}^{\left(i\right)}=\delta\boldsymbol{\varphi}^{\left(i\right)}\).
The discretization of \eqref{eq683} will lead to a contribution to the stiffness matrix. Notice that due to symmetry between \(\delta\mathbf{g}\) and\(\Delta\mathbf{g}\) this matrix will be symmetric.
The contact integral \eqref{eq681} can be discretized as follows. First, we split the integration over all the slave surface elements.
\[ \begin{equation} W_{t}=\sum\limits_{e=1}^{nel}\int\limits_{\Gamma_{c}^{\left(e\right)}}\mathbf{T}\cdot\delta\mathbf{g}\,d\Gamma^{\left(e\right)}\,.\label{eq686} \end{equation} \]
The integration can be approximated by a quadrature rule,
\[ \begin{equation} W_{t}=\sum\limits_{e=1}^{nel}\sum\limits_{i=1}^{N_{int}^{\left(e\right)}}w^{i}j\left(\xi_{i}\right)\mathbf{T}\left(\xi_{i}\right)\cdot\delta\mathbf{g}\left(\xi_{i}\right)\,.\label{eq687} \end{equation} \]
If we use a nodally integrated elements, we have
\[ \begin{equation} \begin{aligned}\mathbf{w}^{\left(1\right)}\left(\xi_{i}\right) & =\mathbf{c}_{i}^{\left(1\right)}\,,\\ \mathbf{w}^{\left(2\right)}\left(\xi_{i}\right) & =\sum\limits_{j}N_{j}\left(\bar{\xi}_{i}\right)\mathbf{c}_{j}^{\left(2\right)}\,, \end{aligned} \label{eq688} \end{equation} \]
so that,
\[ \begin{equation} \delta\mathbf{g}\left(\xi_{i}\right)=\mathbf{c}_{i}^{\left(1\right)}-\sum\limits_{j}N_{j}\left(\bar{\xi}_{i}\right)\mathbf{c}_{j}^{\left(2\right)}\,.\label{eq689} \end{equation} \]
We can now write the contact integral \eqref{eq686} in its final form,
\[ \begin{equation} W_{t}=\sum\limits_{e=1}^{nel}\sum\limits_{i=1}^{N_{int}^{\left(e\right)}}w^{i}j\left(\xi_{i}\right)\left(\mathbf{N}\left(\xi_{i}\right)\mathbf{T}\left(\xi_{i}\right)\right)\cdot\delta\boldsymbol{\Phi}\,,\label{eq690} \end{equation} \]
where
\[ \begin{equation} \delta\boldsymbol{\Phi}^{T}\left(\xi_{i}\right)=\left[\begin{array}{ccccc} \mathbf{c}_{i}^{\left(1\right)} & \mathbf{c}_{1}^{\left(2\right)} & \mathbf{c}_{2}^{\left(2\right)} & \cdots & \mathbf{c}_{n}^{\left(2\right)}\end{array}\right]\,,\label{eq691} \end{equation} \]
\[ \begin{equation} \mathbf{N}\left(\xi_{i}\right)=\left[\begin{array}{cccc} \mathbf{I} & -\mathbf{N}_{1} & \cdots & -\mathbf{N}_{n}\end{array}\right]\,,\label{eq692} \end{equation} \]
and
\[ \begin{equation} \mathbf{N}_{i}=\left[\begin{array}{ccc} N_{i} & 0 & 0\\ 0 & N_{i} & 0\\ 0 & 0 & N_{i} \end{array}\right]\,.\label{eq693} \end{equation} \]
For the linearized tied contact integral \eqref{eq683}, a similar discretization procedure leads to,
\[ \begin{equation} \Delta W_{t}=\sum\limits_{e=1}^{nel}\sum\limits_{i=1}^{N_{int}^{\left(e\right)}}w^{i}j\left(\xi_{i}\right)\Delta\boldsymbol{\Phi}\cdot\mathbf{K}_{c}\delta\boldsymbol{\Phi}\,,\label{eq694} \end{equation} \]
where
\[ \begin{equation} \mathbf{K}_{c}=\varepsilon\mathbf{N}^{T}\mathbf{N}\,.\label{eq695} \end{equation} \]