$\newcommand{\Dev}{\operatorname{Dev}}$ $\newcommand{\dev}{\operatorname{dev}}$ $\newcommand{\grad}{\operatorname{grad}}$ $\newcommand{\Grad}{\operatorname{Grad}}$ $\newcommand{\divg}{\operatorname{div}}$ $\newcommand{\Ei}{\operatorname{Ei}}$ $\newcommand{\tr}{\operatorname{tr}}$ $\newcommand{\Divg}{\operatorname{Div}}$ $\newcommand{\cay}{\operatorname{cay}}$ $\newcommand{\rot}{\operatorname{rot}}$

Consider a damper (e.g., a dashpot) inserted between points on two rigid bodies. The force generated by this dashpot is $$\begin{equation} \mathbf{f}^{\left(1\right)}=\varepsilon_{c}\dot{\mathbf{c}}\,.\label{eq:rbd-1} \end{equation}$$ where $\dot{\mathbf{c}}$ is the relative velocity between insertion points, $$\begin{equation} \dot{\mathbf{c}}=\mathbf{v}_{n+\alpha}^{\left(2\right)}-\mathbf{v}_{n+\alpha}^{\left(1\right)}\,.\label{eq:rbd-2} \end{equation}$$ This relative velocity is related to the rigid body degrees of freedom by $$\begin{equation} \begin{aligned}\mathbf{x}_{n+\alpha}^{\left(i\right)} & =\mathbf{r}_{n+\alpha}^{\left(i\right)}+\mathbf{z}_{n+\alpha}^{\left(i\right)}\\ \mathbf{v}_{n}^{\left(i\right)} & =\dot{\mathbf{r}}_{n}^{\left(i\right)}+\boldsymbol{\omega}_{n}\times\mathbf{z}_{n}^{\left(i\right)}\\ \mathbf{v}_{n+1}^{\left(i\right)} & =\dot{\mathbf{r}}_{n+1}^{\left(i\right)}+\boldsymbol{\omega}_{n+1}\times\mathbf{z}_{n+1}^{\left(i\right)}\\ \mathbf{v}_{n+\alpha}^{\left(i\right)} & =\alpha\mathbf{v}_{n+1}^{\left(i\right)}+\left(1-\alpha\right)\mathbf{v}_{n}^{\left(i\right)} \end{aligned} \,.\label{eq:rbd-3} \end{equation}$$ Since the dashpot may freely rotate about its insertion points, we set $\mathbf{m}^{\left(1\right)}=\mathbf{0}$ . It follows that $$\begin{equation} \boxed{\begin{aligned}\mathbf{f}^{\left(1\right)} & =\varepsilon_{c}\left(\mathbf{v}_{n+\alpha}^{\left(2\right)}-\mathbf{v}_{n+\alpha}^{\left(1\right)}\right)\\ \mathbf{m}^{\left(1\right)} & =\mathbf{0} \end{aligned} }\,,\label{eq:rbd-4} \end{equation}$$ and the virtual work resulting from the dashpot is $$\begin{equation} \delta G=\left[\begin{array}{cccc} \delta\mathbf{r}^{\left(1\right)} & \delta\boldsymbol{\omega}^{\left(1\right)} & \delta\mathbf{r}^{\left(2\right)} & \delta\boldsymbol{\omega}^{\left(2\right)}\end{array}\right]\left[\begin{array}{c} \mathbf{f}^{\left(1\right)}\\ \hat{\mathbf{z}}_{n+\alpha}^{\left(1\right)}\cdot\mathbf{f}^{\left(1\right)}\\ -\mathbf{f}^{\left(1\right)}\\ -\hat{\mathbf{z}}_{n+\alpha}^{\left(2\right)}\cdot\mathbf{f}^{\left(1\right)} \end{array}\right]\,.\label{eq:rbd-5} \end{equation}$$ The linearization of $\delta G$ may be written as $$\begin{equation} \begin{aligned}-D\delta G & =\left[\begin{array}{cccc} \delta\mathbf{r}^{\left(1\right)} & \delta\boldsymbol{\omega}^{\left(1\right)} & \delta\mathbf{r}^{\left(2\right)} & \delta\boldsymbol{\omega}^{\left(2\right)}\end{array}\right]\varepsilon_{c}\alpha\left[\begin{array}{cccc} \mathbf{A} & -\mathbf{B}^{\left(1\right)} & -\mathbf{A} & \mathbf{B}^{\left(2\right)}\\ \hat{\mathbf{z}}_{n+\alpha}^{\left(1\right)}\cdot\mathbf{A} & -\hat{\mathbf{z}}_{n+\alpha}^{\left(1\right)}\cdot\mathbf{B}^{\left(1\right)} & -\hat{\mathbf{z}}_{n+\alpha}^{\left(1\right)}\cdot\mathbf{A} & \hat{\mathbf{z}}_{n+\alpha}^{\left(1\right)}\cdot\mathbf{B}^{\left(2\right)}\\ -\mathbf{A} & \mathbf{B}^{\left(1\right)} & \mathbf{A} & -\mathbf{B}^{\left(2\right)}\\ -\hat{\mathbf{z}}_{n+\alpha}^{\left(2\right)}\cdot\mathbf{A} & \hat{\mathbf{z}}_{n+\alpha}^{\left(2\right)}\cdot\mathbf{B}^{\left(1\right)} & \hat{\mathbf{z}}_{n+\alpha}^{\left(2\right)}\cdot\mathbf{A} & -\hat{\mathbf{z}}_{n+\alpha}^{\left(2\right)}\cdot\mathbf{B}^{\left(2\right)} \end{array}\right]\left[\begin{array}{c} \Delta\mathbf{r}^{\left(1\right)}\\ \Delta\boldsymbol{\theta}^{\left(1\right)}\\ \Delta\mathbf{r}^{\left(2\right)}\\ \Delta\boldsymbol{\theta}^{\left(2\right)} \end{array}\right]\end{aligned} \,,\label{eq:rbd-6} \end{equation}$$ where $$\begin{equation} \left[\mathbf{K}\right]=\varepsilon_{c}\alpha\left[\begin{array}{cccc} \mathbf{A} & -\mathbf{B}^{\left(1\right)} & -\mathbf{A} & \mathbf{B}^{\left(2\right)}\\ \hat{\mathbf{z}}_{n+\alpha}^{\left(1\right)}\cdot\mathbf{A} & -\hat{\mathbf{z}}_{n+\alpha}^{\left(1\right)}\cdot\mathbf{B}^{\left(1\right)} & -\hat{\mathbf{z}}_{n+\alpha}^{\left(1\right)}\cdot\mathbf{A} & \hat{\mathbf{z}}_{n+\alpha}^{\left(1\right)}\cdot\mathbf{B}^{\left(2\right)}\\ -\mathbf{A} & \mathbf{B}^{\left(1\right)} & \mathbf{A} & -\mathbf{B}^{\left(2\right)}\\ -\hat{\mathbf{z}}_{n+\alpha}^{\left(2\right)}\cdot\mathbf{A} & \hat{\mathbf{z}}_{n+\alpha}^{\left(2\right)}\cdot\mathbf{B}^{\left(1\right)} & \hat{\mathbf{z}}_{n+\alpha}^{\left(2\right)}\cdot\mathbf{A} & -\hat{\mathbf{z}}_{n+\alpha}^{\left(2\right)}\cdot\mathbf{B}^{\left(2\right)} \end{array}\right]\,,\label{eq:rbd-7} \end{equation}$$ and $$\begin{equation} \begin{aligned}\mathbf{A} & =\frac{\gamma}{\beta\Delta t}\mathbf{I}\\ \mathbf{B}^{\left(i\right)} & =\frac{\gamma}{\beta\Delta t}\hat{\mathbf{z}}_{n+1}^{\left(i\right)}\cdot\mathbf{T}^{T}\left(\boldsymbol{\theta}^{\left(i\right)}\right)+\hat{\boldsymbol{\omega}}_{n+1}\cdot\hat{\mathbf{z}}_{n+1}^{\left(i\right)} \end{aligned} \,.\label{eq:rbd-8} \end{equation}$$