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

See Section 2.7↑ for a review of multiphasic materials. The contact interface is defined between surfaces $\gamma^{\left(1\right)}$ and $\gamma^{\left(2\right)}$ . Due to continuity requirements on the traction and fluxes, the external virtual work resulting from contact tractions $\mathbf{t}^{\left(i\right)}$ , solvent fluxes $w_{n}^{\left(i\right)}$ and solute fluxes $j_{n}^{\alpha\left(i\right)}$ for solute $\alpha$ ($i=1,2)$ , may be combined into the contact integral $$\begin{equation} \begin{aligned}\delta G_{c} & =\int_{\gamma^{\left(1\right)}}\left(\delta\mathbf{v}^{\left(1\right)}-\delta\mathbf{v}^{\left(2\right)}\right)\cdot\mathbf{t}^{\left(1\right)}\,da^{\left(1\right)}\\ & +\int_{\gamma^{\left(1\right)}}\left(\delta\tilde{p}^{\left(1\right)}-\delta\tilde{p}^{\left(2\right)}\right)w_{n}^{\left(1\right)}\,da^{\left(1\right)}\\ & +\sum\limits _{\alpha}\int_{\gamma^{\left(1\right)}}\left(\delta\tilde{c}^{\alpha\left(1\right)}-\delta\tilde{c}^{\alpha\left(2\right)}\right)j_{n}^{\alpha\left(1\right)}\,da^{\left(1\right)} \end{aligned} \,.\label{eq651} \end{equation}$$ In the current implementation, only frictionless contact is taken into consideration, so that the contact traction has only a normal component, $\mathbf{t}^{\left(i\right)}=t_{n}\mathbf{n}^{\left(i\right)}$ . To evaluate and linearize $\delta G_{c}$ , define the covariant basis vectors on each surface as $$\begin{equation} \mathbf{g}_{\alpha}^{\left(i\right)}=\frac{\partial\mathbf{x}^{\left(i\right)}}{\partial\eta_{\left(i\right)}^{\alpha}},\quad\alpha=1,2,\label{eq652} \end{equation}$$ where $\mathbf{x}^{\left(i\right)}$ represents the spatial position of points on $\gamma^{\left(i\right)}$ , and $\eta_{\left(i\right)}^{\alpha}$ represent the parametric coordinates of that point. The unit outward normal on each surface is then given by $$\begin{equation} \mathbf{n}^{\left(i\right)}=\frac{\mathbf{g}_{1}^{\left(i\right)}\times\mathbf{g}_{2}^{\left(i\right)}}{\left|\mathbf{g}_{1}^{\left(i\right)}\times\mathbf{g}_{2}^{\left(i\right)}\right|}.\label{eq653} \end{equation}$$ Now the contact integral may be rewritten as $$\begin{equation} \begin{aligned}\delta G_{c} & =\int_{\gamma^{\left(1\right)}}\left(\delta\mathbf{v}^{\left(1\right)}-\delta\mathbf{v}^{\left(2\right)}\right)t_{n}\mathbf{g}_{1}^{\left(1\right)}\times\mathbf{g}_{2}^{\left(1\right)}\,d\eta_{\left(1\right)}^{1}d\eta_{\left(1\right)}^{2}\\ & +\int_{\gamma^{\left(1\right)}}\left(\delta\tilde{p}^{\left(1\right)}-\delta\tilde{p}^{\left(2\right)}\right)w_{n}^{\left(1\right)}\left|\mathbf{g}_{1}^{\left(1\right)}\times\mathbf{g}_{2}^{\left(1\right)}\right|\,d\eta_{\left(1\right)}^{1}d\eta_{\left(1\right)}^{2}\\ & +\sum\limits _{\alpha}\int_{\gamma^{\left(1\right)}}\left(\delta\tilde{c}^{\alpha\left(1\right)}-\delta\tilde{c}^{\alpha\left(2\right)}\right)j_{n}^{\alpha\left(1\right)}\left|\mathbf{g}_{1}^{\left(1\right)}\times\mathbf{g}_{2}^{\left(1\right)}\right|\,d\eta_{\left(1\right)}^{1}d\eta_{\left(1\right)}^{2} \end{aligned} \,,\label{eq654} \end{equation}$$ and the linearization $D\delta G_{c}$ of $\delta G_{c}$ has the form $$\begin{equation} D\delta G_{c}=\sum\limits _{i=1}^{2}D\delta G_{c}\left[\Delta\mathbf{u}^{\left(i\right)}\right]+D\delta G_{c}\left[\Delta\tilde{p}^{\left(i\right)}\right]+\sum\limits _{\alpha}D\delta G_{c}\left[\Delta\tilde{c}^{\alpha\left(i\right)}\right]\,.\label{eq655} \end{equation}$$