## 2 Problem Setting
### Equations of Motion  
The contact conditions together with the relevant boundary conditions complement the local balance of momentum and angular momentum (\(\sigma=\sigma^{\intercal}\)). The inertia effects are neglected and a quasi-static process will be considered to deal with Coulomb friction. Hence, the equations of motion are written as follows:  
\[\left\{ \begin{array}{ll}\mathbf{\
abla\cdot P+f=0}&\text{in}\ \ \ \Omega\\ \mathbf{P\,N=F_{0}}&\text{on}\ \ \ \Gamma_{f}\\ \mathbf{u}=\mathbf{u}_{d}&\text{on}\ \ \Gamma_{u}\\ g_{n}\geq 0,\ \ \ \ \Pi_{n}\leq 0,\ \ \ \ g_{n}\Pi_{n}=0&\text{on}\ \ \Gamma_{c}\\ \|\mathbf{\Pi}_{t}\|\leq\mu|\Pi_{n}|,\ \ \ \ \ \|\dot{\mathbf{u}}_{t}\|\mathbf{\Pi}_{t}-\mu|\Pi_{n}|\dot{ \mathbf{u}}_{t}=0,\ \ \ \ \ \|\mathbf{s}\|\|\mathbf{\Pi}_{t}-\mu|\Pi_{n}|\mathbf{s}\|=0&\text{on}\ \ \Gamma_{c}\end{array} \right. \tag{19}\]  
\(\mathbf{f}\) defines the volume external forces, \(\mathbf{F}_{0}\) is a prescribed traction ([PERSON] boundary conditions) and \(\mathbf{u}_{d}\) is a prescribed displacement ([PERSON] boundary conditions). In order to obtain the weak formulation, which is the starting point of the finite element method, the first equation of (19) is multiplied by a test function \(\mathbf{v}\) and integration by parts is then performed. The weak formulation is stated as follows:  
\[\text{Find }\mathbf{u}\in\mathcal{U}=\left\{\mathbf{u}\in H^{\text{I}}(\Omega) |\ \mathbf{u}=\mathbf{u}_{d}\ \text{ on }\ \Gamma_{u}\right\}\]  
\[\int_{\Omega}\mathbf{P}:\mathbf{\
abla v\;d\Omega}-\int_{\Omega}\mathbf{f\cdot v \;d\Omega}-\int_{\partial\Omega}\mathbf{P\mathbf{N\cdot v\;d\Gamma}}=0\ \ \ \ \ \forall \mathbf{v}\in\mathcal{V} \tag{20}\]  
where \(\mathcal{V}=\left\{\mathbf{u}\in H^{\text{I}}(\Omega)|\ \mathbf{u}=\mathbf{0}\ \text{ on }\ \Gamma_{u}\right\}\). The integral defined on the boundary are separated as follows:  
\[\int_{\partial\Omega}\mathbf{P\mathbf{N\cdot v\;d\Gamma}}=\int_{\Gamma_{f}} \mathbf{F}_{0}\mathbf{\cdot v\;d\Gamma}+\int_{\Gamma_{c}}\mathbf{P\mathbf{N \cdot v\;d\Gamma}} \tag{21}\]  
The first term takes into account the Neumann boundary condition. The second one is related to the contact constraint and could be written as follows  
\[\int_{\Gamma_{c}}\mathbf{P\mathbf{N\cdot v\;d\Gamma}}=\int_{\Gamma_{c}} \left(\Pi_{n}\mathbf{n}+\mathbf{\Pi}_{t}\right)\mathbf{\cdot}\left(v_{n} \mathbf{n}+\mathbf{v}_{t}\right)\,\mathrm{d\Gamma}\]  
\[=\int_{\Gamma_{c}}\left(\Pi_{n}v_{n}+\mathbf{\Pi}_{t}\mathbf{\cdot v}_{t}\right) \,\mathrm{d\Gamma}\]  
The weak formulation could now be written as follows:  
\[\left\{ \begin{array}{ll}\text{Find }\mathbf{u}\in\mathcal{U}\\ \int_{\Omega}\mathbf{P}:\mathbf{\
abla v\;d\Omega}-\int_{\Omega}\mathbf{f\cdot v \;d\Omega}-\int_{\Gamma_{f}}\mathbf{F}_{0}\mathbf{\cdot v\;d\Gamma}=\int_{ \Gamma_{c}}\Pi_{n}v_{n}+\mathbf{\Pi}_{t}\mathbf{\cdot v_{t}\;d\Gamma}&\forall \mathbf{v}\in\mathcal{V}\\ g_{n}\geq 0,\ \ \ \ \ \Pi_{n}\leq 0,\ \ \ \ \ g_{n}\Pi_{n}=0&\text{on }\ \Gamma_{c}\\ \|\mathbf{\Pi}_{t}\|\leq\mu|\Pi_{n}|&\|\dot{\mathbf{u}}_{t}\|\mathbf{\Pi}_{t}- \mu|\Pi_{n}|\dot{\mathbf{u}}_{t}=0&\text{on }\ \Gamma_{c}\end{array} \right. \tag{23}\]  
Several formulations can be applied to incorporate the contact constraints into the variational formulation. Most standard finite element codes that are able to handle contact problems use either the penalty or the Lagrange multiplier method. Each of the methods has its own advantages and disadvantages. The methods are designed to fulfill the constraint equations in the normal direction in the contact interface. For the tangential part, a constitutive relation is needed when stick/slip motion occurs [19]. The penalty method is used to treat the normal contact condition. An active set strategy is employed to obtain the active contact zone. The normal contact condition contribution to the weak formulation is approximated as follows [20].