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

A similar procedure can now be used to calculate the discretized contact stiffness matrix. The linearization of the contact integral is repeated here: $$\begin{equation} \begin{aligned}\Delta G^{c}\left(\boldsymbol{\varphi},\mathbf{w}\right) & =\sum\limits _{e=1}^{N_{sel}}\int\limits _{\Gamma^{\left(1\right)e}}\Delta\left(t_{N}\delta g\right)\,d\Gamma\\ & =\sum\limits _{e=1}^{N_{sel}}\sum\limits _{i=1}^{N_{int}^{e}}w_{i}j\left(\boldsymbol{\xi}_{i}\right)\Delta\left(t_{N}\delta g\right)\left(\boldsymbol{\xi}_{i}\right) \end{aligned} \,.\label{eq584} \end{equation}$$ Using matrix notation we can rewrite equation (7.1.5-1) as, $$\begin{equation} \Delta W^{c}\left(\boldsymbol{\varphi},\mathbf{w}\right)=\sum\limits _{e=1}^{N_{sel}}\sum\limits _{i}^{N_{int}^{e}}w_{i}j\left(\mathbf{\xi}_{i}\right)\delta\boldsymbol{\Phi}\cdot\mathbf{k}^{c}\Delta\boldsymbol{\Phi}\,,\label{eq585} \end{equation}$$ where $\delta\boldsymbol{\Phi}$ is as above and $\Delta\boldsymbol{\Phi}$ similar to $\delta\boldsymbol{\Phi}$ with $\delta$ replaced with $\Delta$ , and $\mathbf{k}^{c}$ is $$\begin{equation} \begin{aligned}\mathbf{k}^{c} & =\varepsilon_{N}H\left(\lambda_{N}^{k}+\varepsilon_{N}g\right)\mathbf{NN}^{T}+t_{N}\left\{ g\left[m^{11}\mathbf{\bar{N}}_{1}\mathbf{\bar{N}}_{1}^{T}\right.\right.\\ & +m^{12}\left(\mathbf{\bar{N}}_{1}\mathbf{\bar{N}}_{2}^{T}+\mathbf{\bar{N}}_{2}\mathbf{\bar{N}}_{1}^{T}\right)+\left.m^{22}\mathbf{\bar{N}}_{2}\mathbf{\bar{N}}_{2}^{T}\right]-\mathbf{D}_{1}\mathbf{N}_{1}^{T}\\ & -\mathbf{D}_{2}\mathbf{N}_{2}^{T}-\mathbf{N}_{1}\mathbf{D}_{1}^{T}-\mathbf{N}_{2}\mathbf{D}_{2}^{T}+\left.\kappa_{12}\left(\mathbf{D}_{1}\mathbf{D}_{2}^{T}+\mathbf{D}_{2}\mathbf{D}_{1}^{T}\right)\right\} \end{aligned} \,,\label{eq586} \end{equation}$$ where, $$\begin{equation} \begin{aligned}\mathbf{N}= & \left[\begin{array}{c} \boldsymbol{\nu}\\ -N_{1}\left(\bar{\boldsymbol{\xi}}\right)\boldsymbol{\nu}\\ \vdots\\ -N_{4}\left(\bar{\boldsymbol{\xi}}\right)\boldsymbol{\nu} \end{array}\right]\,, & \mathbf{T}_{\alpha} & =\left[\begin{array}{c} \boldsymbol{\tau}_{\alpha}\\ -N_{1}\left(\bar{\boldsymbol{\xi}}\right)\boldsymbol{\tau}_{\alpha}\\ \vdots\\ -N_{4}\left(\bar{\boldsymbol{\xi}}\right)\boldsymbol{\tau}_{\alpha} \end{array}\right]\,, & \mathbf{N}_{\alpha} & =\left[\begin{array}{c} \mathbf{0}\\ -N_{1,\alpha}\left(\bar{\boldsymbol{\xi}}\right)\boldsymbol{\nu}\\ \vdots\\ -N_{4,\alpha}\left(\bar{\boldsymbol{\xi}}\right)\boldsymbol{\nu} \end{array}\right]\,.\end{aligned} \label{eq587} \end{equation}$$ The following vectors are also defined which depend on the vectors of (7.1.5-4): $$\begin{equation} \begin{aligned}\mathbf{D}_{1} & =\frac{1}{\det\mathbf{A}}\left[A_{22}\left(\mathbf{T}_{1}+g\mathbf{N}_{1}\right)-A_{12}\left(\mathbf{T}_{2}+g\mathbf{N}_{2}\right)\right]\\ \mathbf{D}_{2} & =\frac{1}{\det\mathbf{A}}\left[A_{11}\left(\mathbf{T}_{2}+g\mathbf{N}_{2}\right)-A_{12}\left(\mathbf{T}_{1}+g\mathbf{N}_{1}\right)\right]\\ \mathbf{\bar{N}}_{1} & =\mathbf{N}_{1}-\kappa_{12}\mathbf{D}_{2}\\ \mathbf{\bar{N}}_{2} & =\mathbf{N}_{2}-\kappa_{12}\mathbf{D}_{1} \end{aligned} \,,\label{eq588} \end{equation}$$ where the matrix $\mathbf{A}$ is defined as, $$\begin{equation} A_{ij}=m_{ij}+g\kappa_{ij}\,.\label{eq589} \end{equation}$$ Here, $m_{ij}=\boldsymbol{\tau}_{i}\cdot\boldsymbol{\tau}_{j}$ is the surface metric tensor and $\kappa_{ij}=\boldsymbol{\nu}\cdot\boldsymbol{\varphi}_{t,ij}^{\left(2\right)}\left(\mathbf{\bar{Y}}\right)$ denotes the components of the surface curvature at $\bar{\boldsymbol{\xi}}$ .