7.4 Multiphasic Contact

Contact Integral

See Section Triphasic and Multiphasic Materials 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} \]

Gap Function

The gap function \(g\), representing the distance between the contact surfaces, is defined by

\[ \begin{equation} \mathbf{x}^{\left(2\right)}=\mathbf{x}^{\left(1\right)}+g\mathbf{n}^{\left(1\right)},\quad g=\left(\mathbf{x}^{\left(2\right)}-\mathbf{x}^{\left(1\right)}\right)\cdot\mathbf{n}^{\left(1\right)}\,.\label{eq656} \end{equation} \]

The linearization of variables associated with motion, pressure, and concentration, is given by

\[ \begin{equation} \begin{aligned}D\mathbf{x}^{\left(1\right)}= & \Delta\mathbf{u}^{\left(1\right)} & D\mathbf{x}^{\left(2\right)} & =\Delta\mathbf{u}^{\left(2\right)}+\mathbf{g}_{\alpha}^{\left(2\right)}\Delta\eta_{\left(2\right)}^{\alpha}\\ D\tilde{p}^{\left(1\right)}= & \Delta\tilde{p}^{\left(1\right)} & D\tilde{p}^{\left(2\right)} & =\Delta\tilde{p}^{\left(2\right)}+\frac{\partial\tilde{p}^{\left(2\right)}}{\partial\eta_{\left(2\right)}^{\alpha}}\Delta\eta_{\left(2\right)}^{\alpha}\\ D\tilde{c}^{\gamma\left(1\right)}= & \Delta\tilde{c}^{\gamma\left(1\right)} & D\tilde{c}^{\gamma\left(2\right)} & =\Delta\tilde{c}^{\gamma\left(2\right)}+\frac{\partial\tilde{c}^{\gamma\left(2\right)}}{\partial\eta_{\left(2\right)}^{\alpha}}\Delta\eta_{\left(2\right)}^{\alpha}\\ D\delta\mathbf{v}^{\left(1\right)} & =\mathbf{0} & D\mathbf{v}^{\left(2\right)} & =\frac{\partial\delta\mathbf{v}^{\left(2\right)}}{\partial\eta_{\left(2\right)}^{\alpha}}\Delta\eta_{\left(2\right)}^{\alpha}\\ D\delta\tilde{p}^{\left(1\right)} & =0 & D\delta\tilde{p}^{\left(2\right)} & =\frac{\partial\delta\tilde{p}^{\left(2\right)}}{\partial\eta_{\left(2\right)}^{\alpha}}\Delta\eta_{\left(2\right)}^{\alpha}\\ D\delta\tilde{c}^{\gamma\left(1\right)} & =0 & D\delta\tilde{c}^{\gamma\left(2\right)} & =\frac{\partial\delta\tilde{c}^{\gamma\left(2\right)}}{\partial\eta_{\left(2\right)}^{\alpha}}\Delta\eta_{\left(2\right)}^{\alpha} \end{aligned} \,,\label{eq657} \end{equation} \]

where

\[ \begin{equation} \Delta\eta_{\left(2\right)}^{\alpha}=\left(\Delta\mathbf{u}^{\left(1\right)}-\Delta\mathbf{u}^{\left(2\right)}\right)\cdot a^{\alpha\beta}\mathbf{g}_{\beta}^{\left(1\right)}-a^{\alpha\beta}g\mathbf{n}^{\left(1\right)}\cdot\frac{\partial\Delta\mathbf{u}^{\left(1\right)}}{\partial\eta_{\left(1\right)}^{\beta}}\,,\label{eq658} \end{equation} \]

with \(a^{\alpha\beta}=\left(A_{\alpha\beta}\right)^{-1}\) and \(A_{\alpha\beta}=\mathbf{g}_{\alpha}^{\left(1\right)}\cdot\mathbf{g}_{\beta}^{\left(2\right)}\).

Penalty Method

Let the normal component of the contact traction be described by the penalty function,

\[ \begin{equation} t_{n}=\begin{cases} \varepsilon_{n}g & g<0\\ 0 & g\geqslant0 \end{cases}\,,\label{eq659} \end{equation} \]

where \(\varepsilon_{n}\) is a penalty factor associated with \(t_{n}\). Similarly, let

\[ \begin{equation} \begin{cases} w_{n}=\varepsilon_{p}\pi=\varepsilon_{p}\left(\tilde{p}^{\left(1\right)}-\tilde{p}^{\left(2\right)}\right) & t_{n}<0\\ \tilde{p}^{\left(i\right)}=\tilde{p}^{\ast} & t_{n}=0 \end{cases}\,,\label{eq660} \end{equation} \]

and

\[ \begin{equation} \begin{cases} j_{n}^{\alpha}=\varepsilon_{c}\chi^{\alpha}=\varepsilon_{c}\left(\tilde{c}^{\alpha\left(1\right)}-\tilde{c}^{\alpha\left(2\right)}\right) & t_{n}<0\\ \tilde{c}^{\alpha\left(i\right)}=\tilde{c}^{\alpha\ast} & t_{n}=0 \end{cases}\,,\label{eq661} \end{equation} \]

where \(\varepsilon_{p}\) and \(\varepsilon_{c}\) are penalty factors associated with \(w_{n}\equiv w_{n}^{\left(1\right)}\) and \(j_{n}^{\alpha}\equiv j_{n}^{\alpha\left(1\right)}\), respectively. It follows that

\[ \begin{equation} \begin{aligned}Dt_{n} & =\varepsilon_{n}\left(\Delta\mathbf{u}^{\left(2\right)}-\Delta\mathbf{u}^{\left(1\right)}+\mathbf{g}_{\alpha}^{\left(2\right)}\Delta\eta_{\left(2\right)}^{\alpha}\right)\cdot\mathbf{n}^{\left(1\right)}\\ Dw_{n} & =\varepsilon_{p}\left(\Delta\tilde{p}^{\left(1\right)}-\Delta\tilde{p}^{\left(2\right)}-\frac{\partial\tilde{p}^{\left(2\right)}}{\partial\eta_{\left(2\right)}^{\alpha}}\Delta\eta_{\left(2\right)}^{\alpha}\right)\\ Dj_{n}^{\gamma} & =\varepsilon_{c}\left(\Delta\tilde{c}^{\gamma\left(1\right)}-\Delta\tilde{c}^{\gamma\left(2\right)}-\frac{\partial\tilde{c}^{\gamma\left(2\right)}}{\partial\eta_{\left(2\right)}^{\alpha}}\Delta\eta_{\left(2\right)}^{\alpha}\right) \end{aligned} \,.\label{eq662} \end{equation} \]

Given these relations, it can be shown that the directional derivatives of the various terms appearing in the integrand of \(\delta G_{c}\) are

\[ \begin{equation} \begin{aligned} & -D\left(t_{n}\left(\delta\mathbf{v}^{\left(1\right)}-\delta\mathbf{v}^{\left(2\right)}\right)\cdot\mathbf{g}_{1}^{\left(1\right)}\times\mathbf{g}_{2}^{\left(1\right)}\right)=\\ & J_{\eta}^{\left(1\right)}\varepsilon_{n}\left(\delta\mathbf{v}^{\left(1\right)}-\delta\mathbf{v}^{\left(2\right)}\right)\cdot\left(\mathbf{n}^{\left(1\right)}\otimes\mathbf{n}^{\left(1\right)}\right)\cdot\left(\Delta\mathbf{u}^{\left(1\right)}-\Delta\mathbf{u}^{\left(2\right)}\right)\\ & -J_{\eta}^{\left(1\right)}t_{n}\frac{\partial\delta\mathbf{v}^{\left(2\right)}}{\partial\eta_{\left(2\right)}^{\alpha}}\cdot\left(\mathbf{n}^{\left(2\right)}\otimes\mathbf{g}_{\left(2\right)}^{\alpha}\right)\cdot\left(\Delta\mathbf{u}^{\left(1\right)}-\Delta\mathbf{u}^{\left(2\right)}\right)\\ & +t_{n}\left(\delta\mathbf{v}^{\left(1\right)}-\delta\mathbf{v}^{\left(2\right)}\right)\cdot\left(\mathbf{g}_{2}^{\left(1\right)}\times\frac{\partial\Delta\mathbf{u}^{\left(1\right)}}{\partial\eta_{\left(1\right)}^{1}}-\mathbf{g}_{1}^{\left(1\right)}\times\frac{\partial\Delta\mathbf{u}^{\left(1\right)}}{\partial\eta_{\left(1\right)}^{2}}\right) \end{aligned} \,,\label{eq665} \end{equation} \]

\[ \begin{equation} \begin{aligned} & -D\left(w_{n}\left(\delta\tilde{p}^{\left(1\right)}-\delta\tilde{p}^{\left(2\right)}\right)\left|\mathbf{g}_{1}^{\left(1\right)}\times\mathbf{g}_{2}^{\left(1\right)}\right|\right)=\\ & -J_{\eta}^{\left(1\right)}\varepsilon_{p}\left(\delta\tilde{p}^{\left(1\right)}-\delta\tilde{p}^{\left(2\right)}\right)\left(\Delta\tilde{p}^{\left(1\right)}-\Delta\tilde{p}^{\left(2\right)}\right)\\ & +J_{\eta}^{\left(1\right)}\left[\varepsilon_{p}\left(\delta\tilde{p}^{\left(1\right)}-\delta\tilde{p}^{\left(2\right)}\right)\frac{\partial\tilde{p}^{\left(1\right)}}{\partial\eta_{\left(1\right)}^{\alpha}}\mathbf{g}_{\left(1\right)}^{\alpha}+w_{n}\frac{\partial\delta\tilde{p}^{\left(2\right)}}{\partial\eta_{\left(2\right)}^{\alpha}}\mathbf{g}_{\left(2\right)}^{\alpha}\right]\cdot\left(\Delta\mathbf{u}^{\left(1\right)}-\Delta\mathbf{u}^{\left(2\right)}\right)\\ & +w_{n}\left(\delta\tilde{p}^{\left(1\right)}-\delta\tilde{p}^{\left(2\right)}\right)\mathbf{n}^{\left(1\right)}\cdot\left(\mathbf{g}_{2}^{\left(1\right)}\times\frac{\partial\Delta\mathbf{u}^{\left(1\right)}}{\partial\eta_{\left(1\right)}^{1}}-\mathbf{g}_{1}^{\left(1\right)}\times\frac{\partial\Delta\mathbf{u}^{\left(1\right)}}{\partial\eta_{\left(1\right)}^{2}}\right) \end{aligned} \,,\label{eq666} \end{equation} \]

\[ \begin{equation} \begin{aligned} & -D\left(j_{n}^{\gamma}\left(\delta\tilde{c}^{\gamma\left(1\right)}-\delta\tilde{c}^{\gamma\left(2\right)}\right)\left|\mathbf{g}_{1}^{\left(1\right)}\times\mathbf{g}_{2}^{\left(1\right)}\right|\right)=\\ & -J_{\eta}^{\left(1\right)}\varepsilon_{c}\left(\delta\tilde{c}^{\gamma\left(1\right)}-\delta\tilde{c}^{\gamma\left(2\right)}\right)\left(\Delta\tilde{c}^{\gamma\left(1\right)}-\Delta\tilde{c}^{\gamma\left(2\right)}\right)\\ & +J_{\eta}^{\left(1\right)}\varepsilon_{c}\left(\delta\tilde{c}^{\gamma\left(1\right)}-\delta\tilde{c}^{\gamma\left(2\right)}\right)\frac{\partial\tilde{c}^{\gamma\left(1\right)}}{\partial\eta_{\left(1\right)}^{\alpha}}\mathbf{g}_{\alpha}^{\left(1\right)}\cdot\left(\Delta\mathbf{u}^{\left(1\right)}-\Delta\mathbf{u}^{\left(2\right)}\right)\\ & +J_{\eta}^{\left(1\right)}j_{n}^{\gamma}\frac{\partial\delta\tilde{c}^{\gamma\left(1\right)}}{\partial\eta_{\left(1\right)}^{\alpha}}\mathbf{g}_{\alpha}^{\left(1\right)}\cdot\left(\Delta\mathbf{u}^{\left(1\right)}-\Delta\mathbf{u}^{\left(2\right)}\right)\\ & +j_{n}^{\gamma}\left(\delta\tilde{c}^{\gamma\left(1\right)}-\delta\tilde{c}^{\gamma\left(2\right)}\right)\mathbf{n}^{\left(1\right)}\cdot\left(\mathbf{g}_{2}^{\left(1\right)}\times\frac{\partial\Delta\mathbf{u}^{\left(1\right)}}{\partial\eta_{\left(1\right)}^{1}}-\mathbf{g}_{1}^{\left(1\right)}\times\frac{\partial\Delta\mathbf{u}^{\left(1\right)}}{\partial\eta_{\left(1\right)}^{2}}\right) \end{aligned} \label{eq667} \end{equation} \]

where \(J_{\eta}^{\left(1\right)}=\left|\mathbf{g}_{1}^{\left(1\right)}\times\mathbf{g}_{2}^{\left(1\right)}\right|\).

Discretization

The contact integral may be discretized as

\[ \begin{equation} \delta G_{c}=\sum\limits_{e=1}^{n_{e}^{\left(1\right)}}\sum\limits_{k=1}^{n_{\mbox{int}}^{\left(e\right)}}W_{k}J_{\eta}^{\left(1\right)}\left[t_{n}\left(\delta\mathbf{v}^{\left(1\right)}-\delta\mathbf{v}^{\left(2\right)}\right)\cdot\mathbf{n}^{\left(1\right)}+w_{n}\left(\delta\tilde{p}^{\left(1\right)}-\delta\tilde{p}^{\left(2\right)}\right)+\sum\limits_{\gamma}j_{n}^{\gamma}\left(\delta\tilde{c}^{\gamma\left(1\right)}-\delta\tilde{c}^{\gamma\left(2\right)}\right)\right]\,.\label{eq668} \end{equation} \]

The variables may be interpolated over each element face according to

\[ \begin{equation} \begin{aligned}\delta\mathbf{v}^{\left(1\right)} & =\sum\limits_{a=1}^{m^{\left(1\right)}}N_{a}^{\left(1\right)}\delta\mathbf{v}_{a}^{\left(1\right)} & \delta\mathbf{v}^{\left(2\right)} & =\sum\limits_{b=1}^{m^{\left(2\right)}}N_{b}^{\left(2\right)}\delta\mathbf{v}_{b}^{\left(2\right)}\\ \Delta\mathbf{u}^{\left(1\right)} & =\sum\limits_{c=1}^{m^{\left(1\right)}}N_{c}^{\left(1\right)}\Delta\mathbf{u}_{c}^{\left(1\right)} & \Delta\mathbf{u}^{\left(2\right)} & =\sum\limits_{d=1}^{m^{\left(2\right)}}N_{d}^{\left(2\right)}\Delta\mathbf{u}_{d}^{\left(2\right)}\\ \delta\tilde{p}^{\left(1\right)} & =\sum\limits_{a=1}^{m^{\left(1\right)}}N_{a}^{\left(1\right)}\delta\tilde{p}_{a}^{\left(1\right)} & \delta\tilde{p}^{\left(2\right)} & =\sum\limits_{b=1}^{m^{\left(2\right)}}N_{b}^{\left(2\right)}\delta\tilde{p}_{b}^{\left(2\right)}\\ \Delta\tilde{p}^{\left(1\right)} & =\sum\limits_{c=1}^{m^{\left(1\right)}}N_{c}^{\left(1\right)}\Delta\tilde{p}_{c}^{\left(1\right)} & \Delta\tilde{p}^{\left(2\right)} & =\sum\limits_{d=1}^{m^{\left(2\right)}}N_{d}^{\left(2\right)}\Delta\tilde{p}_{d}^{\left(2\right)}\\ \delta\tilde{c}^{\gamma\left(1\right)} & =\sum\limits_{a=1}^{m^{\left(1\right)}}N_{a}^{\left(1\right)}\delta\tilde{c}_{a}^{\gamma\left(1\right)} & \delta\tilde{c}^{\gamma\left(2\right)} & =\sum\limits_{b=1}^{m^{\left(2\right)}}N_{b}^{\left(2\right)}\delta\tilde{c}_{b}^{\gamma\left(2\right)}\\ \Delta\tilde{c}^{\gamma\left(1\right)} & =\sum\limits_{c=1}^{m^{\left(1\right)}}N_{c}^{\left(1\right)}\Delta\tilde{c}_{c}^{\gamma\left(1\right)} & \Delta\tilde{c}^{\gamma\left(2\right)} & =\sum\limits_{d=1}^{m^{\left(2\right)}}N_{d}^{\left(2\right)}\Delta\tilde{c}_{d}^{\gamma\left(2\right)} \end{aligned} \,.\label{eq669} \end{equation} \]

Then,

\[ \begin{equation} \begin{aligned}\delta G_{c} & =\sum\limits_{e=1}^{n_{e}^{\left(1\right)}}\sum\limits_{k=1}^{n_{\mbox{int}}^{\left(e\right)}}W_{k}J_{\eta}^{\left(1\right)}\left(\sum\limits_{a=1}^{m^{\left(1\right)}}\left[\begin{array}{cccc} \delta\mathbf{v}_{a}^{\left(1\right)} & \delta\tilde{p}_{a}^{\left(1\right)} & \delta\tilde{c}_{a}^{\alpha\left(1\right)} & \delta\tilde{c}_{a}^{\beta\left(1\right)}\end{array}\right]\cdot\left[\begin{array}{c} \mathbf{f}_{a}^{\left(1\right)}\\ w_{a}^{\left(1\right)}\\ j_{a}^{\alpha\left(1\right)}\\ j_{a}^{\beta\left(1\right)} \end{array}\right]\right.\\ & +\left.\sum\limits_{b=1}^{m^{\left(2\right)}}\left[\begin{array}{cccc} \delta\mathbf{v}_{b}^{\left(2\right)} & \delta\tilde{p}_{b}^{\left(2\right)} & \delta\tilde{c}_{b}^{\alpha\left(2\right)} & \delta\tilde{c}_{b}^{\beta\left(2\right)}\end{array}\right]\cdot\left[\begin{array}{c} \mathbf{f}_{b,k}^{\left(2\right)}\\ w_{b,k}^{\left(2\right)}\\ j_{b,k}^{\alpha\left(2\right)}\\ j_{b,k}^{\beta\left(2\right)} \end{array}\right]\right)\,, \end{aligned} \label{eq670} \end{equation} \]

where

\[ \begin{equation} \begin{aligned}\mathbf{f}_{a}^{\left(1\right)} & =t_{n}N_{a}^{\left(1\right)}\mathbf{n}^{\left(1\right)} & \mathbf{f}_{b,k}^{\left(2\right)} & =-t_{n}N_{b}^{\left(2\right)}\mathbf{n}^{\left(1\right)}\\ w_{a}^{\left(1\right)} & =w_{n}N_{a}^{\left(1\right)} & w_{b,k}^{\left(2\right)} & =-w_{n}N_{b}^{\left(2\right)}\\ j_{a}^{\gamma\left(1\right)} & =j_{n}^{\gamma}N_{a}^{\left(1\right)} & j_{b,k}^{\gamma\left(2\right)} & =-j_{n}^{\gamma}N_{b}^{\left(2\right)} \end{aligned} \,.\label{eq671} \end{equation} \]

Similarly,

\[ \begin{equation} \begin{aligned}-D\delta G_{c} & =\sum\limits_{e=1}^{n_{e}^{\left(1\right)}}\sum\limits_{k=1}^{n_{\mbox{int}}^{\left(e\right)}}W_{k}J_{\eta}^{\left(1\right)}\times\left(\sum\limits_{a=1}^{m^{\left(1\right)}}\left[\begin{array}{cccc} \delta\mathbf{v}_{a}^{\left(1\right)} & \delta\tilde{p}_{a}^{\left(1\right)} & \delta\tilde{c}_{a}^{\alpha\left(1\right)} & \delta\tilde{c}_{a}^{\beta\left(1\right)}\end{array}\right]\cdot\right.\\ & \left(\sum\limits_{c=1}^{m^{\left(1\right)}}\left[\begin{array}{cccc} \mathbf{K}_{ac}^{\left(1,1\right)} & 0 & 0 & 0\\ \mathbf{g}_{ac}^{\left(1,1\right)} & g_{ac}^{\left(1,1\right)} & 0 & 0\\ \mathbf{h}_{ac}^{\alpha\left(1,1\right)} & 0 & h_{ac}^{\alpha\alpha\left(1,1\right)} & 0\\ \mathbf{h}_{ac}^{\beta\left(1,1\right)} & 0 & 0 & h_{ac}^{\beta\beta\left(1,1\right)} \end{array}\right]\cdot\left[\begin{array}{c} \Delta\mathbf{u}_{c}^{\left(1\right)}\\ \Delta\tilde{p}_{c}^{\left(1\right)}\\ \Delta\tilde{c}_{c}^{\alpha\left(1\right)}\\ \Delta\tilde{c}_{c}^{\beta\left(1\right)} \end{array}\right]\right.\\ & +\left.\sum\limits_{d=1}^{m_{k}^{\left(2\right)}}\left[\begin{array}{cccc} \mathbf{K}_{ad,k}^{\left(1,2\right)} & 0 & 0 & 0\\ \mathbf{g}_{ad,k}^{\left(1,2\right)} & g_{ad,k}^{\left(1,2\right)} & 0 & 0\\ \mathbf{h}_{ad,k}^{\alpha\left(1,2\right)} & 0 & h_{ad,k}^{\alpha\alpha\left(1,2\right)} & 0\\ \mathbf{h}_{ad,k}^{\beta\left(1,2\right)} & 0 & 0 & h_{ad,k}^{\beta\beta\left(1,2\right)} \end{array}\right]\cdot\left[\begin{array}{c} \Delta\mathbf{u}_{d}^{\left(2\right)}\\ \Delta\tilde{p}_{d}^{\left(2\right)}\\ \Delta\tilde{c}_{d}^{\alpha\left(2\right)}\\ \Delta\tilde{c}_{d}^{\beta\left(2\right)} \end{array}\right]\right)\\ & +\sum\limits_{b=1}^{m_{k}^{\left(2\right)}}\left[\begin{array}{cccc} \delta\mathbf{v}_{b,k}^{\left(2\right)} & \delta\tilde{p}_{b,k}^{\left(2\right)} & \delta\tilde{c}_{b,k}^{\alpha\left(2\right)} & \delta\tilde{c}_{b,k}^{\beta\left(2\right)}\end{array}\right]\cdot\\ & \left(\sum\limits_{c=1}^{m^{\left(1\right)}}\left[\begin{array}{cccc} \mathbf{K}_{bc,k}^{\left(2,1\right)} & 0 & 0 & 0\\ \mathbf{g}_{bc,k}^{\left(2,1\right)} & g_{bc,k}^{\left(2,1\right)} & 0 & 0\\ \mathbf{h}_{bc,k}^{\alpha\left(2,1\right)} & 0 & h_{bc,k}^{\alpha\alpha\left(2,1\right)} & 0\\ \mathbf{h}_{bc,k}^{\beta\left(2,1\right)} & 0 & 0 & h_{bc,k}^{\beta\beta\left(2,1\right)} \end{array}\right]\cdot\left[\begin{array}{c} \Delta\mathbf{u}_{c}^{\left(1\right)}\\ \Delta\tilde{p}_{c}^{\left(1\right)}\\ \Delta\tilde{c}_{c}^{\alpha\left(1\right)}\\ \Delta\tilde{c}_{c}^{\beta\left(1\right)} \end{array}\right]\right.\\ & +\left.\left.\sum\limits_{d=1}^{m_{k}^{\left(2\right)}}\left[\begin{array}{cccc} \mathbf{K}_{bd,k}^{\left(2,2\right)} & 0 & 0 & 0\\ \mathbf{g}_{bd,k}^{\left(2,2\right)} & g_{bd,k}^{\left(2,2\right)} & 0 & 0\\ \mathbf{h}_{bd,k}^{\alpha\left(2,2\right)} & 0 & h_{bd,k}^{\alpha\alpha\left(2,2\right)} & 0\\ \mathbf{h}_{bd,k}^{\beta\left(2,2\right)} & 0 & 0 & h_{bd,k}^{\beta\beta\left(2,2\right)} \end{array}\right]\cdot\left[\begin{array}{c} \Delta\mathbf{u}_{d}^{\left(2\right)}\\ \Delta\tilde{p}_{d}^{\left(2\right)}\\ \Delta\tilde{c}_{d}^{\alpha\left(2\right)}\\ \Delta\tilde{c}_{d}^{\beta\left(2\right)} \end{array}\right]\right)\right)\,, \end{aligned} \label{eq672} \end{equation} \]

where

\[ \begin{equation} \begin{aligned}\mathbf{K}_{ac}^{\left(1,1\right)} & =N_{a}^{\left(1\right)}\left(\varepsilon_{n}N_{c}^{\left(1\right)}\mathbf{N}^{\left(1\right)}+t_{n}\mathbf{A}_{c}^{\left(1\right)}\right)\\ \mathbf{K}_{ad,k}^{\left(1,2\right)} & =-\varepsilon_{n}N_{a}^{\left(1\right)}N_{d}^{\left(2\right)}\mathbf{N}^{\left(1\right)}\\ \mathbf{K}_{bc,k}^{\left(2,1\right)} & =-N_{c}^{\left(1\right)}\left(\varepsilon_{n}N_{b}^{\left(2\right)}\mathbf{N}^{\left(1\right)}+t_{n}\mathbf{M}_{b}^{\left(2\right)}\right)-t_{n}N_{b}^{\left(2\right)}\mathbf{A}_{c}^{\left(1\right)}\\ \mathbf{K}_{bd,k}^{\left(2,2\right)} & =N_{d}^{\left(2\right)}\left(\varepsilon_{n}N_{b}^{\left(2\right)}\mathbf{N}^{\left(1\right)}+t_{n}\mathbf{M}_{b}^{\left(2\right)}\right) \end{aligned} \,,\label{eq673} \end{equation} \]

\[ \begin{equation} \begin{aligned}\mathbf{g}_{ac}^{\left(1,1\right)} & =N_{a}^{\left(1\right)}\left(\varepsilon_{p}N_{c}^{\left(1\right)}\mathbf{p}^{\left(1\right)}-w_{n}\mathbf{A}_{c}^{\left(1\right)}\cdot\mathbf{n}^{\left(1\right)}\right)\\ \mathbf{g}_{ad,k}^{\left(1,2\right)} & =-\varepsilon_{p}N_{a}^{\left(1\right)}N_{d}^{\left(2\right)}\mathbf{p}^{\left(1\right)}\\ \mathbf{g}_{bc,k}^{\left(2,1\right)} & =N_{c}^{\left(1\right)}\left(-\varepsilon_{p}N_{b}^{\left(2\right)}\mathbf{p}^{\left(1\right)}+w_{n}\mathbf{m}_{b}^{\left(2\right)}\right)+w_{n}N_{b}^{\left(2\right)}\mathbf{A}_{c}^{\left(1\right)}\cdot\mathbf{n}^{\left(1\right)}\\ \mathbf{g}_{bd,k}^{\left(2,2\right)} & =N_{d}^{\left(2\right)}\left(\varepsilon_{p}N_{b}^{\left(2\right)}\mathbf{p}^{\left(1\right)}-w_{n}\mathbf{m}_{b}^{\left(2\right)}\right) \end{aligned} \,,\label{eq674} \end{equation} \]

\[ \begin{equation} \begin{aligned}g_{ac}^{\left(1,1\right)} & =-\varepsilon_{p}N_{a}^{\left(1\right)}N_{c}^{\left(1\right)} & g_{ad,k}^{\left(1,2\right)} & =\varepsilon_{p}N_{a}^{\left(1\right)}N_{d}^{\left(2\right)}\\ g_{bc,k}^{\left(2,1\right)} & =\varepsilon_{p}N_{b}^{\left(2\right)}N_{c}^{\left(1\right)} & g_{bd,k}^{\left(2,2\right)} & =-\varepsilon_{p}N_{b}^{\left(2\right)}N_{d}^{\left(2\right)} \end{aligned} \,,\label{eq675} \end{equation} \]

\[ \begin{equation} \begin{aligned}\mathbf{h}_{ac}^{\gamma\left(1,1\right)} & =N_{a}^{\left(1\right)}\left(N_{c}^{\left(1\right)}\varepsilon_{c}\mathbf{q}^{\gamma\left(1\right)}-\mathbf{A}_{c}^{\left(1\right)}\cdot j_{n}^{\gamma}\mathbf{n}^{\left(1\right)}\right)\\ \mathbf{h}_{ad,k}^{\gamma\left(1,2\right)} & =-N_{a}^{\left(1\right)}N_{d}^{\left(2\right)}\varepsilon_{c}\mathbf{q}^{\gamma\left(1\right)}\\ \mathbf{h}_{bc,k}^{\gamma\left(2,1\right)} & =N_{b}^{\left(2\right)}\left(-N_{c}^{\left(1\right)}\varepsilon_{c}\mathbf{q}^{\gamma\left(1\right)}+\mathbf{A}_{c}^{\left(1\right)}\cdot j_{n}^{\gamma}\mathbf{n}^{\left(1\right)}\right)+N_{c}^{\left(1\right)}j_{n}^{\gamma}\mathbf{m}_{b}^{\left(2\right)}\\ \mathbf{h}_{bd,k}^{\gamma\left(2,2\right)} & =N_{b}^{\left(2\right)}N_{d}^{\left(2\right)}\varepsilon_{c}\mathbf{q}^{\gamma\left(1\right)}-N_{d}^{\left(2\right)}j_{n}^{\gamma}\mathbf{m}_{b}^{\left(2\right)} \end{aligned} \,,\label{eq676} \end{equation} \]

\[ \begin{equation} \begin{aligned}h_{ac}^{\gamma\gamma\left(1,1\right)} & =-\varepsilon_{c}N_{a}^{\left(1\right)}N_{c}^{\left(1\right)}\\ h_{ad,k}^{\gamma\gamma\left(1,2\right)} & =\varepsilon_{c}N_{a}^{\left(1\right)}N_{d}^{\left(2\right)}\\ h_{bc,k}^{\gamma\gamma\left(2,1\right)} & =\varepsilon_{c}N_{b}^{\left(2\right)}N_{c}^{\left(1\right)}\\ h_{bd,k}^{\gamma\gamma\left(2,2\right)} & =-\varepsilon_{c}N_{b}^{\left(2\right)}N_{d}^{\left(2\right)} \end{aligned} \,,\label{eq677} \end{equation} \]

and

\[ \begin{equation} \begin{aligned}\mathbf{N}^{\left(1\right)} & =\mathbf{n}^{\left(1\right)}\otimes\mathbf{n}^{\left(1\right)} & \mathbf{A}_{c}^{\left(1\right)} & =\frac{1}{J_{\eta}^{\left(1\right)}}\boldsymbol{\mathcal{A}}\left\{ \frac{\partial N_{c}^{\left(1\right)}}{\partial\eta_{\left(1\right)}^{1}}\mathbf{g}_{2}^{\left(1\right)}-\frac{\partial N_{c}^{\left(1\right)}}{\partial\eta_{\left(1\right)}^{2}}\mathbf{g}_{1}^{\left(1\right)}\right\} \\ \mathbf{M}_{b}^{\left(2\right)} & =\mathbf{n}^{\left(2\right)}\otimes\mathbf{m}_{b}^{\left(2\right)} & \mathbf{m}_{b}^{\left(2\right)} & =\frac{\partial N_{b}^{\left(2\right)}}{\partial\eta_{\left(2\right)}^{\alpha}}\mathbf{g}_{\left(2\right)}^{\alpha}\\ \mathbf{p}^{\left(1\right)} & =\frac{\partial\tilde{p}^{\left(1\right)}}{\partial\eta_{\left(1\right)}^{\alpha}}\mathbf{g}_{\left(1\right)}^{\alpha} & \mathbf{q}^{\gamma\left(1\right)} & =\frac{\partial\tilde{c}^{\gamma\left(1\right)}}{\partial\eta_{\left(1\right)}^{\alpha}}\mathbf{g}_{\left(1\right)}^{\alpha} \end{aligned} \,.\label{eq678} \end{equation} \]