# Contact analysis within the bi-potential framework using cell-based smoothed finite element method

Qianwei Chen, Yan Li, Zhiqiang Feng, Huijian Chen

► To cite this version:

Qianwei Chen, Yan Li, Zhiqiang Feng, Huijian Chen. Contact analysis within the bi-potential framework using cell-based smoothed finite element method. 11th International Conference on Computational Methods (ICCM2020), Aug 2020, Online, France. hal-03254422

HAL Id: hal-03254422

https://hal.science/hal-03254422

Submitted on 8 Jun 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Contact analysis within the bi-potential framework using cell-based
smoothed finite element method

*Qianwei Chen¹, †Yan Li¹, Zhiqiang Feng¹,² and Huijian Chen¹

¹School of Mechanics and Engineering, Southwest Jiaotong University, Chengdu, China

² LMEE Univ-Evry, Université Paris-Saclay, Evry, France

*Presenting author: qwchen0110@163.com

†Corresponding author: liyan001de@gmail.com

**Abstract**

This paper presents a cell-based smoothed finite element method (CS-FEM) for solving two-
dimensional contact problems with the bi-potential formulation. The contact force and the
relative displacement are coupled with each other and solved by Uzawa algorithm. Three
contact states are investigated accurately. The CS-FEM is performed with six different kinds
of smoothing domains which are constructed by dividing the background element into
different regions. Three numerical examples are presented to verify the accuracy of the
method. The effect of the friction coefficient for the contact are also investigated. All
solutions agree well with reference values. The results produced by the CS-FEM are more
accurate than those of the traditional FEM. Besides, the CS-FEM can provide both upper
bound and lower bound solutions for the strain energy while using different smoothing
domains.

**Keywords:** Contact analysis, CS-FEM, Bi-potential formulation, Uzawa algorithm

**1. Introduction**

The contact problems exist everywhere in life. For example, the contact between the tire and
the road, the contact between the shaft and the bearing, the contact between the knife and the
flesh in virtual surgery, etc. They play important roles in many engineering applications [1][2].

There are two reasons contributed to the non-linearity of the contact problems [3][4]. One is
that both the contact surfaces and their positions change during the contact process. The other
one is the non-linearity of contact conditions, including the non-penetration of contact bodies
and the tangential friction conditions. The widely used contact algorithms in engineering
applications are the penalty function method [5]-[8], the Lagrange multiplier method [9]-[12]
and the linear complementarity techniques [13]-[16]. The bi-potential method proposed by De
Saxcé and Feng provides an effective tool for dissipative constitutive modelling [17][18].
Applying the augmented Lagrangian method to the law of contact, the implicit equation of the
projection on the Coulomb's cone is equivalent to the original contact inequality. There has
demonstrated that the Uzawa algorithm is more efficient when compared with the Newton
method [19]. There are many researches and applications in impact, hyperelasticity, wear and
other problems [20]-[23].

The finite element method (FEM) is one of the most effective numerical methods for solving
contact problems. Other widely used approaches include the boundary element methods
[24][25] and the meshless methods [26]-[28]. The smoothed finite element method (S-FEM)
proposed by Liu et al. is a weakened weak form method based on the G-space theory [29].
When using different kinds of smoothing domains which are based on cells, edges, nodes and
faces of the background elements, the cell-based smoothed finite element method (CS-FEM),
the edge-based smoothed finite element method, the node-based smoothed finite element method and the face-based smoothed finite element method are created [30]-[33]. So far, the S-FEM has been widely used in various fields, including the acoustics [34], materials science [35], vibration [36], fluid-structure coupling [37], and electromagnetics [38]. The S-FEM has also been applied to solve contact problems [16][26][39][40].

In this paper, the CS-FEM with quadrilateral elements is used to solve the two-dimensional contact problems within the bi-potential framework. The contact forces can be solved using the Uzawa algorithm. The rest of the paper is organized as follows: in Section 2, the governing equations and the contact criterion are introduced. Section 3 illustrates the smoothed finite element method, especially the CS-FEM. The smoothed Galerkin weak form is obtained by using the smoothing strain technique. In Section 4, the contact analysis is performed within the bi-potential framework. The Uzawa algorithm is applied to solve the contact force. Finally, three numerical examples are presented to examine the numerical accuracy of the proposed method.

## 2. Problem statement

### 2.1 Governing equations

Considering an elastic contact body with domain $\Omega_i$ and the force boundary condition $\Gamma_f^i$, the displacement boundary condition $\Gamma_u^i$, the contact boundary condition $\Gamma_c^i$ as shown in Fig. 1.

<figure><img src="image_2.png" /><figcaption>Figure 1. Contact body</figcaption></figure>

The governing equation of this contact body is:

$$ \nabla^T \sigma + \mathbf{f} = 0 \qquad (1) $$

with

$$ \left\{ \begin{array}{ll} \mathbf{L}^T \sigma = \bar{\mathbf{t}} & \text{on } \Gamma_f \\ \mathbf{L}^T \sigma = \bar{\mathbf{r}} & \text{on } \Gamma_c \\ \mathbf{u} = \bar{\mathbf{u}} & \text{on } \Gamma_u \end{array} \right. \qquad (2) $$

where $\nabla$ is a differential operator which can be written for 2D problems as follows:

$$ \nabla = \begin{bmatrix} \partial/\partial x & 0 \\ 0 & \partial/\partial y \\ \partial/\partial y & \partial/\partial x \end{bmatrix} \qquad (3) $$
The stress $\sigma = [\sigma_{xx} \ \sigma_{yy} \ \tau_{xy}]^T$. The body force $\mathbf{f} = [f_x \ f_y]^T$ and displacements $\mathbf{u} = [u_x \ u_y]^T$. $\bar{\mathbf{t}}$, $\bar{\mathbf{u}}$ and $\bar{\mathbf{r}}$ are the loads, displacements and contact forces on the boundary, respectively. $\mathbf{L}$ is the matrix of the unit normal vectors defined by:

$$ \mathbf{L} = \begin{bmatrix} \partial x / \partial n & 0 \\ 0 & \partial y / \partial n \\ \partial y / \partial n & \partial x / \partial n \end{bmatrix}. \qquad (4) $$

## 2.2 Contact criterion

The Signorini condition mainly describes a normal contact relationship. It has the following three characteristics at each contact point, they are the geometric condition of non-penetration, the static condition of non-adhesion and the mechanical complementary condition. Therefore, the following formula expresses the relationship between the normal distance $x_n^\alpha$ and normal contact force $r_n^\alpha$ for any contact point $\alpha$. The Signorini condition is defined as:

$$ \text{Signor}(x_n^\alpha, r_n^\alpha) \Leftrightarrow \begin{cases} x_n^\alpha \ge 0 \\ r_n^\alpha \ge 0 \\ x_n^\alpha r_n^\alpha = 0 \end{cases} \qquad (5) $$

In order to satisfy the tangential criterion of contact, the Coulomb friction rule is adopted and:

$$ \text{Coul}(\mathbf{x}_t^\alpha, \mathbf{r}_t^\alpha) \Leftrightarrow \begin{cases} \|\mathbf{x}_t^\alpha\| = 0 & \|\mathbf{r}_t^\alpha\| \le \mu r_n^\alpha \\ \|\mathbf{x}_t^\alpha\| \ne 0 & \mathbf{r}_t^\alpha = -\mu r_n^\alpha \frac{\mathbf{x}_t^\alpha}{\|\mathbf{x}_t^\alpha\|} \end{cases} \qquad (6) $$

where $\mathbf{x}_t^\alpha$ and $\mathbf{r}_t^\alpha$ are the tangential components of the spacing vector and the contact force vector for contact point $\alpha$, respectively. $\mu$ is the friction coefficient.

Let $K_\mu$ represents the Coulomb's cone which is expressed as:

$$ K_{\mu} = \{ \mathbf{r} \in \mathbb{R}^3 \text{ such that } \|\mathbf{r}_{t}^{\alpha}\| - \mu r_{n}^{\alpha} \le 0 \} \qquad (7) $$

The contact satisfies a complex non-smooth dissipative law including the separating, sticking and sliding. By combining Eqs. (5) and (7), the contact criterion can be written as:

$$ \begin{align} \text{Separating: } & x_n^\alpha > 0 \quad \text{and} \quad r_n^\alpha = 0 \notag \\ \text{Sticking: } & \|\mathbf{x}_t^\alpha\| = 0 \quad \text{and} \quad \mathbf{r}_t^\alpha \in \text{int}(K_\mu) \tag{8} \\ \text{Sliding: } & \|\mathbf{x}_t^\alpha\| \ne 0 \quad \text{and} \quad \mathbf{r}_t^\alpha \in \text{bd}(K_\mu) \quad \text{with} \quad \mathbf{r}_t^\alpha = -\mu r_n^\alpha \frac{\mathbf{x}_t^\alpha}{\|\mathbf{x}_t^\alpha\|} \notag \end{align} $$

where $\text{int}(K_\mu)$ and $\text{bd}(K_\mu)$ represent the interior and the boundary of Coulomb's cone. Eq. (8) is called the Signorini-Coulomb conditions.

# 3. Cell-based smoothed finite element method

## 3.1 Galerkin weak form

For problem domain $\Omega$, by introducing a weight function $\varphi$, Eq. (1) can be transformed into an equivalent integral form:

$$ \int_{\Omega} \varphi^{\mathrm{T}} \nabla^{\mathrm{T}} \sigma d\Omega + \int_{\Omega} \varphi^{\mathrm{T}} \mathbf{f} d\Omega = 0 \qquad (9) $$
Using Green's divergence theorem, we can convert Eq. (9) to:

$$ \int_{\Omega} (\nabla \boldsymbol{\phi})^{\mathrm{T}} \mathbf{C} \varepsilon d\Omega - \int_{\Omega} \boldsymbol{\phi}^{\mathrm{T}} \mathbf{f} d\Omega - \int_{\Gamma_f} \boldsymbol{\phi}^{\mathrm{T}} \bar{\mathbf{t}} d\Gamma - \int_{\Gamma_c} \boldsymbol{\phi}^{\mathrm{T}} \bar{\mathbf{r}} d\Gamma = 0 \quad (10) $$

where **C** is the material matrix. The physical equation is **σ**=$\mathbf{C}\varepsilon$, where $\varepsilon$ is the vector of strain. The weight function takes the variation of displacement as **φ** = $\delta\mathbf{u}$ and the Galerkin weak form can be got:

$$ \int_{\Omega} (\delta \boldsymbol{\varepsilon})^{\mathrm{T}} \mathbf{C} \varepsilon d\Omega - \int_{\Omega} (\delta \mathbf{u})^{\mathrm{T}} \mathbf{f} d\Omega - \int_{\Gamma_f} (\delta \mathbf{u})^{\mathrm{T}} \bar{\mathbf{t}} d\Gamma - \int_{\Gamma_c} (\delta \mathbf{u})^{\mathrm{T}} \bar{\mathbf{r}} d\Gamma = 0 \quad (11) $$

where $\delta\boldsymbol{\varepsilon} = \nabla\delta\mathbf{u}$.

## 3.2 Smoothed Galerkin weak form

In S-FEM, the smoothed strain $\bar{\boldsymbol{\varepsilon}}$ is calculated by the strain smoothing operation as:

$$ \bar{\boldsymbol{\varepsilon}}(\bar{\mathbf{x}}) = \int_{\Omega_s} \boldsymbol{\varepsilon}(\mathbf{x}) W(\bar{\mathbf{x}} - \mathbf{x}) d\Omega \quad (12) $$

where $\Omega_s$ is the smoothing domain which can be created based on nodes, edges, cells of the background elements. In this work the CS-FEM model is adopted with six different kinds of smoothing regions as shown in Fig. 2. $W$ is a smoothing function that should be positive over the local support domain and satisfies the unity property $\int_{\Omega_s} W(\bar{\mathbf{x}} - \mathbf{x}) d\Omega = 1$.

<figure><img src="image_3.png" /><figcaption>Figure 2. Six different smoothing domains based on a quadrilateral element: 1SD, 2SDs, 3SDs, 4SDs, 8SDs, 16SDs</figcaption></figure>

Here the Heaviside-type smoothing function is adopted:

$$ W(\bar{\mathbf{x}} - \mathbf{x}) = \begin{cases} 1 / \bar{A} & \mathbf{x} \in \Omega_s \\ 0 & \mathbf{x} \notin \Omega_s \end{cases} \quad (13) $$

where $\bar{A}$ is the area of the smoothing domain. Substituting Eq. (13) into Eq. (12) and using the Green's divergence theorem, the smoothed strain can be obtained:

$$ \begin{aligned} \bar{\boldsymbol{\varepsilon}} &= \int_{\Omega_s} \nabla \mathbf{u} W(\bar{\mathbf{x}} - \mathbf{x}) d\Omega \\ &= -\int_{\Omega_s} \nabla W(\bar{\mathbf{x}} - \mathbf{x}) \mathbf{u} d\Omega + \frac{1}{\bar{A}} \int_{\Gamma_s} \mathbf{L} \mathbf{u} d\Gamma \quad (14) \\ &= \frac{1}{\bar{A}} \int_{\Gamma_s} \mathbf{L} \mathbf{u} d\Gamma \end{aligned} $$

where $\Gamma_s$ is the boundary of the smoothing domain. Note that only boundary integration instead of domain integral is required in this calculation.

By substituting Eq. (14) into Eq. (11), the following smoothed Galerkin weak form can be obtained:
$$ \int_{\Omega} (\delta \bar{\boldsymbol{\epsilon}})^{\mathrm{T}} \mathbf{C} \bar{\boldsymbol{\epsilon}} d\Omega - \int_{\Omega} (\delta \mathbf{u})^{\mathrm{T}} \mathbf{f} d\Omega - \int_{\Gamma_f} (\delta \mathbf{u})^{\mathrm{T}} \bar{\mathbf{t}} d\Gamma - \int_{\Gamma_c} (\delta \mathbf{u})^{\mathrm{T}} \bar{\mathbf{r}} d\Gamma = 0 \quad (15) $$

### 3.3 Discretized system of equations

Assume the displacement function as:

$$ \mathbf{u}(\mathbf{x}) = \mathbf{N}(\mathbf{x})\mathbf{d} \quad \mathbf{x} \in \Omega_s \quad (16) $$

where vector $\mathbf{N}(\mathbf{x})$ and $\mathbf{d}$ represent the shape function and nodal displacements respectively.

Substituting Eq. (16) into Eq. (14), the smoothing strain $\bar{\boldsymbol{\epsilon}}$ can be written as:

$$ \bar{\boldsymbol{\epsilon}} = \frac{1}{A} \int_{\Gamma_s} \mathbf{L} \mathbf{u} d\Gamma = \bar{\mathbf{B}} \mathbf{d} \quad (17) $$

where $\bar{\mathbf{B}}$ denotes the smoothed strain matrix which can be calculated as:

$$ \bar{\mathbf{B}} = \frac{1}{A} \int_{\Gamma_s} \mathbf{L} \mathbf{N}(\mathbf{x}) d\Gamma = \begin{bmatrix} \bar{b}_x & 0 \\ 0 & \bar{b}_y \\ \bar{b}_y & \bar{b}_x \end{bmatrix} \quad \mathbf{x} \in \Gamma_s \quad (18) $$

with

$$ \bar{b}_h = \frac{1}{A} \sum_{p=1}^{N_p} n_{h,p} \mathbf{N}(\mathbf{x}_p^G) l_p \quad h = x, y \quad (19) $$

where $N_p$ represents the number of segments of the boundary $\Gamma_s$, $n_{h,p}$ and $\mathbf{x}_p^G$ are respectively the outward unit normal and Gauss points on each segment, $l_p$ represents the length of the p-th segment.

Substituting Eq. (16) and Eq. (17) into Eq. (15), we can get the discretized algebraic system of equations:

$$ \sum_{e=1}^{N_e} \sum_{s=1}^{N_s} \bar{A}_s \bar{\mathbf{B}}^{\mathrm{T}} \mathbf{C} \bar{\mathbf{B}} \mathbf{d} - \sum_{e=1}^{N_e} \int_{\Omega_e} \mathbf{N}^{\mathrm{T}} \mathbf{f} d\Omega - \sum_{e=1}^{N_e} \int_{\Gamma_{ef}} \mathbf{N}^{\mathrm{T}} \bar{\mathbf{t}} d\Gamma - \sum_{e=1}^{N_e} \int_{\Gamma_{ec}} \mathbf{N}^{\mathrm{T}} \bar{\mathbf{r}} d\Gamma = 0 \quad (20) $$

where $N_e$ and $N_s$ represent the number of background elements and smoothing domains. $\Omega_e$ represents background element. $\Gamma_{ef}$ and $\Gamma_{ec}$ represent the boundary of forces and contact forces per element. For simplicity of description, we can use the following symbols:

$$ \mathbf{F}_{\mathrm{ext}} = \sum_{e=1}^{N_e} \int_{\Omega_e} \mathbf{N}^{\mathrm{T}} \mathbf{f} d\Omega + \sum_{e=1}^{N_e} \int_{\Gamma_{ef}} \mathbf{N}^{\mathrm{T}} \bar{\mathbf{t}} d\Gamma \quad (21) $$

$$ \mathbf{R} = \sum_{e=1}^{N_e} \int_{\Gamma_{ec}} \mathbf{N}^{\mathrm{T}} \bar{\mathbf{r}} d\Gamma \quad (22) $$

$$ \mathbf{F}_{\mathrm{int}} = \sum_{e=1}^{N_e} \sum_{s=1}^{N_s} \bar{A}_s \bar{\mathbf{B}}^{\mathrm{T}} \mathbf{C} \bar{\mathbf{B}} \mathbf{d} \quad (23) $$

where $\mathbf{F}_{\mathrm{ext}}$ is the external forces. $\mathbf{R}$ is the contact forces. $\mathbf{F}_{\mathrm{int}}$ is the internal forces.

The Eq. (20) can be rewritten as:

$$ \bar{\mathbf{K}} \mathbf{d} = \mathbf{F}_{\mathrm{ext}} + \mathbf{R} \quad (24) $$

where $\bar{\mathbf{K}}$ is stiffness matrix obtained by:

$$ \bar{\mathbf{K}} = \sum_{e=1}^{N_e} \sum_{s=1}^{N_s} \bar{A}_s \bar{\mathbf{B}}^{\mathrm{T}} \mathbf{C} \bar{\mathbf{B}} \quad (25) $$
Here Eq. (24) cannot be solved directly because both **d** and **R** are unknown quantities. In our work we compute the contact force **R** on the contact surface first. Then the displacements **d** can be solved by taking **R** as an external force.

## 4. The contact within the bi-potential framework

### 4.1 Contact kinematics

Considering two elastic bodies $\Omega_1$ and $\Omega_2$ coming into contact, as shown in Fig. 3. $P_1$ is the contact points on the boundary $\Gamma_1$. $P_2$ is the normal projection point of the contact point $P_1$ on the boundary $\Gamma_2$. For 2D problems, we can set up a local coordinate system based on the normal and tangential vector at point $P_1$. The initial gap **g** between $P_1$ and $P_2$ is determined by a contact collision detector as shown in Fig. 3:

<figure><img src="image_4.png" /><figcaption>Figure 3. Contact kinematics</figcaption></figure>

Assuming:

$$d(P_1) = \mathbf{\Phi}_1 d_1 \quad d(P_2) = \mathbf{\Phi}_2 d_2 \qquad (26)$$

where $d_1$ and $d_2$ denotes respectively the nodal displacements in $\Omega_1$ and $\Omega_2$. $\mathbf{\Phi}_1$ and $\mathbf{\Phi}_2$ are integration matrix in corresponding regions. If we consider the case with $N_c$ contact points $P_1^\alpha$ ($\alpha = 1 \cdots N_c$) on $\Gamma_1$, and $P_2^\alpha$ is the normal projection points of the contact points $P_1^\alpha$ on $\Gamma_2$. The relative position between $P_1^\alpha$ and $P_2^\alpha$ is:

$$d^\alpha = d(P_2^\alpha) - d(P_1^\alpha) \qquad (27)$$

In local coordinates it can be expressed as:

$$d^\alpha = x_i^\alpha \mathbf{t} + x_n^\alpha \mathbf{n} \qquad (28)$$

where

$$\mathbf{x}^\alpha = \begin{Bmatrix} x_i^\alpha \\ x_n^\alpha \end{Bmatrix} = \begin{Bmatrix} \mathbf{t}^\mathrm{T} \mathbf{d}^\alpha \\ \mathbf{n}^\mathrm{T} \mathbf{d}^\alpha \end{Bmatrix}. \qquad (29)$$

Substituting Eq. (26) and Eq. (27) to Eq. (29), we get:

$$\begin{Bmatrix} x_1^\alpha \\ x_n^\alpha \end{Bmatrix} = \begin{bmatrix} -\mathbf{t}^\mathrm{T} \mathbf{\Phi}_1 & 0 & \mathbf{t}^\mathrm{T} \mathbf{\Phi}_2 & 0 \\ 0 & -\mathbf{n}^\mathrm{T} \mathbf{\Phi}_1 & 0 & \mathbf{n}^\mathrm{T} \mathbf{\Phi}_2 \end{bmatrix} \begin{Bmatrix} d_{1x} \\ d_{1y} \\ d_{2x} \\ d_{2y} \end{Bmatrix} \qquad (30)$$

By defining the transformation matrix **T** as:

$$\mathbf{T}^\alpha = [-\mathbf{T}_1^\alpha | \mathbf{T}_2^\alpha] \quad \mathbf{T}_k^\alpha = \begin{bmatrix} \mathbf{t}^\mathrm{T} \mathbf{\Phi}_k & 0 \\ 0 & \mathbf{n}^\mathrm{T} \mathbf{\Phi}_k \end{bmatrix} \qquad (31)$$
Eq. (30) can be written as:

$$ \mathbf{x}^{\alpha} = \mathbf{T}^{\alpha} \mathbf{d} \qquad (32) $$

whose incremental form is:

$$ \mathbf{x}^{\alpha(i+1)} = \mathbf{T}^{\alpha} \mathbf{d}^{(i)} + \mathbf{g}^{\alpha} \qquad (33) $$

where $\mathbf{g}^{\alpha} = (0 \ g^{\alpha})$.

In the local reference coordinate system, the contact force $\mathbf{r}^{\alpha}$ can be defined as:

$$ \mathbf{r}^{\alpha} = r_{t}^{\alpha} \mathbf{t} + r_{n}^{\alpha} \mathbf{n} \qquad (34) $$

Applying the contact virtual work principle:

$$ \mathbf{r}^{\alpha\mathrm{T}} \delta\mathbf{x}^{\alpha} = \mathbf{R}^{\alpha\mathrm{T}} \delta\mathbf{d}^{\alpha} \qquad (35) $$

We have:

$$ \mathbf{R}^{\alpha} = \mathbf{T}^{\alpha\mathrm{T}} \mathbf{r}^{\alpha} \qquad (36) $$

The whole system of equations of $N_c$ contact points can be written as:

$$ \left\{ \begin{array}{l} \mathbf{R} = \mathbf{T}^{\mathrm{T}} \mathbf{r} \\ \mathbf{x} = \mathbf{T} \mathbf{d} + \mathbf{g} \end{array} \right. \qquad (37) $$

with $\mathbf{T} = [\mathbf{T}^1 \ \cdots \ \mathbf{T}^\alpha]^\mathrm{T}$, $\mathbf{x} = [\mathbf{x}^1 \ \cdots \ \mathbf{x}^{N_c}]^\mathrm{T}$, $\mathbf{r} = [\mathbf{r}^1 \ \cdots \ \mathbf{r}^{N_c}]^\mathrm{T}$, $\mathbf{g} = [\mathbf{g}^1 \ \cdots \ \mathbf{g}^{N_c}]^\mathrm{T}$.

## 4.2 The contact solution method

For a contact point, the bi-potential function of the contact law can be written as:

$$ b_c(-\mathbf{x}^{\alpha}, \mathbf{r}^{\alpha}) = \prod_{\mathbf{R}_-} (-x_n^{\alpha}) + \prod_{K_{\mu}} (\mathbf{r}^{\alpha}) + \mu r_n^{\alpha} \|\mathbf{x}_t^{\alpha}\| \qquad (38) $$

where $\mathbf{R}_- = (-\infty, 0]$, and $\prod_{K_{\mu}} (\mathbf{r}^{\alpha})$ is the indicator function and it can be described as: if $\mathbf{r}^{\alpha} \in K_{\mu}$, $\prod_{K_{\mu}} (\mathbf{r}^{\alpha}) = 0$; otherwise $\prod_{K_{\mu}} (\mathbf{r}^{\alpha}) = +\infty$.

The Eq. (38) can be further expressed as an implicit subnormal form:

$$ \begin{gathered} -\mathbf{x}^{\alpha} \in \frac{\partial b_c(-\mathbf{x}^{\alpha}, \mathbf{r}^{\alpha})}{\partial r^{\alpha}} \\ \mathbf{r}^{\alpha} \in \frac{\partial b_c(-\mathbf{x}^{\alpha}, \mathbf{r}^{\alpha})}{\partial (-x^{\alpha})} \end{gathered} \qquad (39) $$

So, using the augmented Lagrangian method [41], the Eq. (38) and Eq. (39) can be equivalent to the following projection form on the Coulomb's cone:

$$ \mathbf{r}^{\alpha} = \mathrm{Proj}_{K_{\mu}}(\mathbf{r}^{\alpha*}) \qquad (40) $$

where $\mathbf{r}^{\alpha*}$ represents the vector of augmented contact force, and it is given by:

$$ \mathbf{r}^{\alpha*} = \mathbf{r}^{\alpha} - \rho \mathbf{x}^{\alpha*} \quad \text{with} \quad \mathbf{x}^{\alpha*} = \mathbf{x}^{\alpha} + \mu \|\mathbf{x}_t^{\alpha}\| \mathbf{n}^{\alpha} \qquad (41) $$

where $\rho$ is the bi-potential coefficient. In this work, the reciprocal of the maximum value of the diagonal elements of the matrix $\mathbf{D}$ as given in Eq. (45) is selected.

By summarizing the above equations and Eq. (8), there will be three contact situations within the bi-potential framework. They are $\mathbf{r}^{\alpha*} \in K_{\mu}$ (sticking), $\mathbf{r}^{\alpha*} \in K_{\mu}^*$ (separating) and $\mathbf{r}^{\alpha*} \in \mathbb{R}^2 - (K_{\mu} \cup K_{\mu}^*)$ (sliding), respectively, as shown in Fig. 4.
<figure><img src="image_5.png" /><figcaption>Figure 4. Coulomb cone and contact projection operators</figcaption></figure>

Consequently, we can define the projection operation by:

$$
\begin{align}
\text{Separating : } \operatorname{Proj}_{K_{\mu}}(\mathbf{r}^{\alpha*}) &= 0 \quad \text{if } \mu \|\mathbf{r}_t^{\alpha*}\| < -r_n^{\alpha*} \nonumber \\
\text{Sticking : } \operatorname{Proj}_{K_{\mu}}(\mathbf{r}^{\alpha*}) &= \mathbf{r}^{\alpha*} \quad \text{if } \|\mathbf{r}_t^{\alpha*}\| < \mu r_n^{\alpha*} \tag{42} \\
\text{Sliding : } \operatorname{Proj}_{K_{\mu}}(\mathbf{r}^{\alpha*}) &= \mathbf{r}^{\alpha*} - \left( \frac{\|\mathbf{r}_t^{\alpha*}\| - \mu r_n^{\alpha*}}{1 + \mu^2} \right) \left( \frac{\mathbf{r}_t^{\alpha*}}{\|\mathbf{r}_t^{\alpha*}\|} - \mu \mathbf{n}^{\alpha} \right) \quad \text{else} \nonumber
\end{align}
$$

4.3 Equilibrium equations of contact points

The system of equations related to each contact point is:

$$
\left\{
\begin{array}{l}
\bar{\mathbf{K}}\mathbf{d} = \mathbf{F}_{\text{ext}} + \mathbf{T}^{\text{T}}\mathbf{r} \\
\mathbf{x} = \mathbf{T}\mathbf{d} + \mathbf{g} \\
\mathbf{r}^{\alpha} = \operatorname{Proj}_{K_{\mu}}(\mathbf{r}^{\alpha*})
\end{array}
\right.
\qquad (43)
$$

eliminating **d** leads to the following system of equations:

$$
\left\{
\begin{array}{l}
\mathbf{x} = \mathbf{T}\bar{\mathbf{K}}^{-1}\mathbf{T}^{\mathrm{T}}\mathbf{r} + \mathbf{T}\bar{\mathbf{K}}^{-1}\mathbf{F}_{\mathrm{ext}} + \mathbf{g} \\
\mathbf{r}^{\alpha} = \mathrm{Proj}_{K_{\mu}}(\mathbf{r}^{\alpha*})
\end{array}
\right.
\qquad (44)
$$

Here, we can define some variables to make the equation more intuitive.

$$
\left\{
\begin{array}{l}
\mathbf{D} = \mathbf{T}\bar{\mathbf{K}}^{-1}\mathbf{T}^{\mathrm{T}} \\
\hat{\mathbf{x}} = \mathbf{T}\bar{\mathbf{K}}^{-1}\mathbf{F}_{\mathrm{ext}} + \mathbf{g}
\end{array}
\right.
\qquad (45)
$$

where **D** is named as the global Delassus operator [42].

For a contact system containing $N_c$ contact points, the contact points are coupled with each other. So, Eq. (44) can be transformed into an implicit system of equations.

$$
\mathbf{\Psi}^{\alpha}(\mathbf{r}, \mathbf{x}) = \left\{ \begin{array}{l} \mathbf{x}^{\alpha} - \mathbf{D}^{\alpha\alpha} \mathbf{r}^{\alpha} - \sum_{\beta=1, \beta \neq \alpha}^{N_c} \mathbf{D}^{\alpha\beta} \mathbf{r}^{\beta} - \hat{\mathbf{x}}^{\alpha} \\ \\ \mathbf{r}^{\alpha} - \operatorname{Proj}_{K_{\mu}}(\mathbf{r}^{\alpha*}) \end{array} \right\} = 0 \qquad (46)
$$

Where $\mathbf{D}^{\alpha\beta} = \mathbf{T}^{\alpha}\bar{\mathbf{K}}^{-1}\mathbf{T}^{\mathrm{T}\beta}$, $\Psi(\mathbf{r}, \mathbf{x}) = [\Psi^1 \ \dots \ \Psi^{N_c}]^\mathrm{T} = 0$.

4.4 Uzawa algorithm

Uzawa algorithm is a local iterative algorithm to solve implicitly Eq. (46). The calculation procedure mainly includes predictor and corrector for contact forces.
$$ \textit{Predictor : } \mathbf{r}^{\alpha*(i+1)} = \mathbf{r}^{\alpha(i)} - \rho^{(i)}(\mathbf{x}^{\alpha(i)} + \mu \|\mathbf{x}_t^{\alpha(i)}\|\mathbf{n}^{\alpha}) \quad (47) $$

$$ \textit{Correntor : } \mathbf{r}^{\alpha(i+1)} = \text{Proj}_{K_{\mu}}(\mathbf{r}^{\alpha*(i+1)}) $$

where i and i+1 are iteration numbers.

**Algorithm** Uzawa algorithm to solve contact problem

assume: $\mathbf{r}^{(0)} = 0$, $\mathbf{g}^{(0)} = (0 \ g)$, $\mathbf{x}^{\alpha(0)} = \mathbf{T}\bar{\mathbf{K}}^{-1}\mathbf{F}_{\text{ext}} + \mathbf{g}^{(0)}$.

for (i = 1 to $\Theta$)

for ($\alpha$ = 1 to $N_c$)

$$ \mathbf{r}^{\alpha*(i+1)} = \mathbf{r}^{\alpha(i)} - \rho^{(i)}(\mathbf{x}^{\alpha(i)} + \mu \|\mathbf{x}_t^{\alpha(i)}\|\mathbf{n}^{\alpha}); $$

$$ \mathbf{r}^{\alpha(i+1)} = \text{Proj}_{K_{\mu}}(\mathbf{r}^{\alpha*(i+1)}); $$

$$ \mathbf{x}^{\alpha(i+1)} = \mathbf{D}^{\alpha\alpha}\mathbf{r}^{\alpha(i+1)} + \sum_{\beta=1,\beta \neq \alpha}^{N_c} \mathbf{D}^{\alpha\beta}\mathbf{r}^{\beta(i)} + \hat{\mathbf{x}}^{\alpha(i)}; $$

end for

$$ \text{if } \left( \frac{\|\mathbf{r}^{(i+1)} - \mathbf{r}^{(i)}\|}{\|\mathbf{r}^{(i+1)}\|} \le \varepsilon_g \right) $$

break;

end if

end for

where $\Theta$ is the total number of contact iteration. $\varepsilon_g$ is a user-defined convergence coefficient.

## 5. Numerical examples

In this section, we will use three numerical examples to examine the efficiency of the contact analysis while combining the CS-FEM and the bi-potential formulation. For Q4 element, the characteristic length can be defined [29]:

$$ h = \sqrt{A_{\Omega} / N_e} \quad (48) $$

where $A_{\Omega}$ is the area of whole body. The strain energy error $e_{\text{error}}$ and the convergence ratio k are defined as:

$$ e_{\text{error}} = \frac{|E_n - E_r|}{E_r} \quad (49) $$

$$ k = \frac{\sum_{i=1}^{N_m} \log_{10}(e_{\text{error}})_i \left[ \log_{10}(h)_i - \frac{1}{N_m} \sum_{i=1}^{N_m} \log_{10}(h)_i \right]}{\sum_{i=1}^{N_m} \left[ \log_{10}(h)_i - \frac{1}{N_m} \sum_{i=1}^{N_m} \log_{10}(h)_i \right]^2} \quad (50) $$

where $E_n$ and $E_r$ are the numerical solution and reference solution of the strain energy, respectively. $N_m$ is the number of mesh models. In this work $N_m = 6$.

### 5.1 Frictionless contact of the flat and cylindrical bodies

Consider the contact of a crossed cylindrical-flat arrangement as shown in Fig. 5 (a). The diameter of the cylindrical sample is 12 mm. The bottom is a block, and its thickness is 10 mm [41]. This problem can be simplified to a plane strain model as shown in Fig. 5 (b).
The displacements are uniformly distributed on top of cylindrical body $\Omega_1$ while a full constraint is applied at the bottom of the flat body $\Omega_2$.

<figure><img src="image_6.png" /><figcaption>Figure 5. (a) Model of the crossed cylindrical-flat; (b) Plane strain model</figcaption></figure>

In this example, the Young's modulus and Poisson's ratio of the two bodies are $E_1=2 \times 10^7$ MPa, $\nu_1=0.3$ and $E_2=2 \times 10^8$ MPa, $\nu_2=0.3$ respectively. When using different number of nodes and elements are listed in Table 1, the problem domain meshed with quadrangles are given in Fig. 6.

<table><caption>Table 1. The number of nodes and elements of two elastic bodies.</caption><thead><tr><th>No.</th><th>M1</th><th>M2</th><th>M3</th><th>M4</th><th>M5</th><th>M6</th><th>Ref.</th></tr></thead><tbody><tr><td>Nodes</td><td>345</td><td>767</td><td>1237</td><td>1977</td><td>3050</td><td>4272</td><td>36175</td></tr><tr><td>Elements</td><td>299</td><td>699</td><td>1151</td><td>1871</td><td>2920</td><td>4118</td><td>35755</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>Figure 6. Six different mesh models</figcaption></figure>

5.1.1 Convergence of strain energy solution
While choosing $U_y = -0.01$ mm. Fig. 7 (a) compares the solutions of CS-FEM with different number of smoothing domains. From which we observed that all the solutions converge to the reference solution with the number of degrees of freedom increases. The reference value of the strain energy is obtained by the FEM-Q4 with 35755 quadrilateral elements, as shown in Table 1. And Fig. 7 (b) shows the convergence of the strain energy with respect to the characteristic length. Note that although the accuracy of CS-FEM-1SD is higher, it is unstable when the constraint is insufficient.

<figure><img src="image_8.png" /><figcaption>Figure 7. Numerical solutions (a) and the convergence of strain energy (b)</figcaption></figure>

According to the above figures, several phenomena can be drawn: 1) As the number of degrees of freedom increases, the strain energy solution of CS-FEM-1SD converges from the lower bound and CS-FEM-2SD to 16SD converges from the upper bound. 2) When using the displacement boundary condition, the CS-FEM-1SD provides a lower bound solution while the other models provide upper bound solutions. All the results produced by CS-FEM are more accurate than those of the FEM-Q4. 3) These calculation methods are approximately linear convergent, CS-FEM-1SD has the largest convergence ratio 1.87, the others convergence ratios are from 1.52 to 1.56. 4) The accuracy of CS-FEM-1SD is almost 5 times higher than the FEM-Q4 when using the same quadrilateral elements. CS-FEM-4SD is more stable, and its accuracy is more than 40% higher than the FEM-Q4.

### 5.1.2 Hertz contact verification

In order to investigate the accuracy of the CS-FEM within the bi-potential framework to solve the contact problem, here we adopt the Hertz theoretical solution for comparison [43]. The contact pressure distribution of this model is:

$$p(x) = p_0 \left(1 - \frac{x^2}{a^2}\right)^{1/2} \qquad (51)$$

where $x$ is the distribution of nodes on the contact surface. $a$ and $p_0$ are the half-width of the contact area and the maximum contact pressure, respectively, which can be expressed as:

$$a = \left( \frac{4Pr}{\pi E^*} \right)^{1/2} \qquad (52)$$

$$p_0 = \left( \frac{PE^*}{\pi r} \right)^{1/2} \qquad (53)$$

where $P$ is the contact reaction force, $E^*$ is the composite modulus and it can be written as:
$$E^* = \left( \frac{1 - v_1^2}{E_1} + \frac{1 - v_2^2}{E_2} \right)^{-1} \qquad (54)$$

$r$ is the relative curvature:

$$r = \left( \frac{1}{r_1} + \frac{1}{r_2} \right)^{-1} \qquad (55)$$

The flat body $\Omega_2$ is an infinite plane and $r_2 \to \infty$. thus, we have $r = r_1$.

For comparison, we use three boundary conditions: $U_{y1}=-0.002 \text{ mm}$, $U_{y2}=-0.006 \text{ mm}$ and $U_{y3}=-0.01 \text{ mm}$. CS-FEM-4SD and the mesh M6 are adopted for discretization. Then, the contact reaction forces corresponding to different boundary conditions are $P_1=14983.471 \text{ N}$, $P_2=52873.824 \text{ N}$ and $P_3=96133.574 \text{ N}$ respectively. Fig. 8 (a) shows the normal contact stress and the analytical solution of Hertz contact. We can find that the contact reaction force and the normal contact stress will increase when using larger displacements on the boundary. All the solutions agree well with the theoretical solution. Fig. 8 (b) shows the tangential relative slip on the contact surface with respect to the contact length.

<figure><img src="image_9.png" /><figcaption>Figure 8. Normal contact stress (a) and tangential relative slip on the contact surface (b)</figcaption></figure>

### 5.1.3 Analysis of penetration

Fig. 9 shows the penetration of different methods for solving the contact problem within the bi-potential framework. Since the largest penetration produced by CS-FEM-1D is less than $3 \times 10^{-20} \text{ mm}$, it can be neglected in the calculation.

<figure><img src="image_10.png" /><figcaption>Figure 9. Penetration by using different methods</figcaption></figure>
## 5.2 Contact between two elastic bodies

Considering the two elastic bodies $\Omega_1$ and $\Omega_2$ as shown in Fig. 10. The displacements are uniformly distributed on top of cylindrical body $\Omega_1$ while a full constraint is applied at the bottom of the flat body $\Omega_2$.

<figure><img src="image_11.png" /><figcaption>Figure 10. Contact between two elastic bodies</figcaption></figure>

In this example, a plane strain model is used. When using different number of nodes and elements are listed in Table 2, the problem domain meshed with quadrangles are given in Fig. 11.

<table><caption>Table 2. The number of nodes and elements of elastic body $\Omega_1$</caption><thead><tr><th>No.</th><th>M1</th><th>M2</th><th>M3</th><th>M4</th><th>M5</th><th>M6</th><th>Ref.</th></tr></thead><tbody><tr><td>Nodes</td><td>66</td><td>231</td><td>496</td><td>861</td><td>1326</td><td>1891</td><td>51681</td></tr><tr><td>Elements</td><td>50</td><td>200</td><td>450</td><td>800</td><td>1250</td><td>1800</td><td>51200</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>Figure 11. Six different mesh models for two elastic bodies</figcaption></figure>

### 5.2.1 Convergence of strain energy solution

While choosing $U_y = -0.04 m$ and the friction coefficient $\mu=0.8$. Fig. 7 (a) compares the solutions of CS-FEM with different number of smoothing domains. From which we observed that all the solutions converge to the reference solution with the number of degrees of freedom increases. The reference value of the strain energy is obtained by the FEM-Q4 with 35755 quadrilateral elements, as shown in Table 1. Fig. 7 (b) shows the convergence of the strain
energy with respect to the characteristic length. Note that although the accuracy of CS-FEM-1SD is higher, it is unstable when the constraint is insufficient.

<figure><img src="image_13.png" /><figcaption>Figure 12. Numerical solutions (a) and the convergence of strain energy (b)</figcaption></figure>

According to the above description, this example can also get a similar phenomenon to Fig. 7: As the number of degrees of freedom increases, the strain energy solution of CS-FEM converges from the lower bound and the upper bound. When using the displacement boundary condition, the CS-FEM-1SD provides a lower bound solution while the other models provide upper bound solutions. All the results produced by CS-FEM are more accurate than those of the FEM-Q4. The difference is (1) The CS-FEM-1SD has the largest convergence ratio 1.69, the others convergence ratios are from 1.43 to 1.49. (2) When using the same quadrilateral elements, the accuracy of the CS-FEM-1SD is almost 6 times higher than the FEM-Q4. The CS-FEM-4SD is more stable, and its accuracy is more than 30% higher than the FEM-Q4.

### 5.2.2 Effect of friction coefficient

In the following contact analysis, we consider the effect of friction coefficients. The CS-FEM-4SD and M6 model are selected. The friction coefficients are selected as 0.1, 0.15, 0.2, 0.25.

<figure><img src="image_14.png" /><figcaption>Figure 13. Normal contact stress (a) and tangential contact stress (b)</figcaption></figure>

Fig. 13 (a) and (b) shows the normal and tangential contact stress of the contact point using different friction coefficients. We can find that the normal contact stress will become higher at the right end, and the tangential contact force will rise first, and then will become higher at the right end when the friction coefficient increases. Fig. 14 (a) and (b) shows the tangential
relative slip and the ratio of contact forces which we can observe the contact states (sticking or sliding) of different contact points.

<figure><img src="image_15.png" /><figcaption>Figure 14. Tangential relative slip (a) and the evolution of contact states (b)</figcaption></figure>

### 5.2.3 Analysis of penetration

Fig. 15 shows the penetration of different cases. Since the largest penetration produced by CS-FEM-1D is less than $8 \times 10^{-16} m$, it can be neglected in the calculation.

<figure><img src="image_16.png" /><figcaption>Figure 15. The comparison of penetration when using different methods</figcaption></figure>

## 5.3 Contact between two elastic bodies with large slips

Finally, we consider the two elastic bodies $\Omega_1$ and $\Omega_2$ with large slip as shown in Fig. 16. The displacements are uniformly distributed on top of $\Omega_1$ while a full constraint is applied at the bottom of $\Omega_2$.

In this example, a plane stress model is used. When using different number of nodes and elements as listed in Table 3, the problem domain meshed with quadrangles are given in Fig. 17.
<figure><img src="image_17.png" /><figcaption>Figure 16. Contact between two elastic bodies with large slip</figcaption></figure>

<table><caption>Table 3. The number of nodes and elements of two elastic bodies Ω1 and Ω2.</caption><thead><tr><th>No.</th><th>M1</th><th>M2</th><th>M3</th><th>M4</th><th>M5</th><th>M6</th><th>Ref.</th></tr></thead><tbody><tr><td>Nodes</td><td>65</td><td>214</td><td>449</td><td>770</td><td>1177</td><td>1670</td><td>27377</td></tr><tr><td>Elements</td><td>43</td><td>172</td><td>387</td><td>688</td><td>1075</td><td>1548</td><td>26875</td></tr></tbody></table>

<figure><img src="image_18.png" /><figcaption>Figure 17. Six different mesh models for two elastic bodies</figcaption></figure>

### 5.3.1 Convergence of strain energy solution

While choosing $U_x=100$ mm, $U_y=-15$ mm and the friction coefficient $\mu = 0.2$. In the following analysis, the multi-step loading is used with 100 load steps. Fig. 18 (a) compares the solutions of CS-FEM with different number of smoothing domains. From which we observed that all the solutions converge to the reference solution with the number of degrees of freedom increases. The reference value of the strain energy is obtained by the FEM-Q4 with 26875 quadrilateral elements, as shown in Table 3. And Fig. 18 (b) shows the convergence of the strain energy with respect to the characteristic length. Note that although the accuracy of CS-FEM-1SD is much higher than the other models, it could be unstable when the constraint is insufficient.
<figure><img src="image_19.png" /><figcaption>Figure 18. Numerical solutions (a) and the convergence of the strain energy (b)</figcaption></figure>

According to these figures, we can get the similar phenomenon with the previous two examples: 1) The CS-FEM-1SD has the least convergence ratio 1.0, the others convergence ratios are from 1.43 to 1.48. 2) When using the same quadrilateral elements, the accuracy of the CS-FEM-1SD is almost 5 times higher than the FEM-Q4, and the accuracy of the CS-FEM-4SD is more than 30% higher than the FEM-Q4. Note that the CS-FEM-4SD is more stable.

### 5.3.2 Effect of friction coefficient and load step

In the following contact analysis, the CS-FEM-4SD and M6 model are selected. The friction coefficients are 0.05, 0.1, 0.15, 0.2, respectively. Fig. 19 (a) and (b) shows the normal and tangential contact stress of the contact point using different friction coefficients at load step 100. We can find that as the friction coefficient increases, the upward trend of normal contact stress on the contact surface will be more obvious. The left end of the contact surface decreases while the right end of the contact surface increases. The tangential contact stress will become higher at the right end.

<figure><img src="image_20.png" /><figcaption>Figure 19. Normal contact stress (a) and tangential contact stress (b)</figcaption></figure>

We define the friction coefficient is 0.2 to calculate the normal and tangential contact stress of the contact point at different load steps as shown in Fig. 20. We can find that both normal and tangential contact stress increase steadily in the process of loading.
<figure><img src="image_21.png" /><figcaption>Figure 20. Normal contact stress (a) and tangential contact stress (b)</figcaption></figure>

Fig. 21 shows the ratio of contact forces which we can observe the contact states (sticking or sliding) of different contact points A, B and C. We can find that during the loading process, they are all in a sliding state, and as the loading progresses, the sliding will be more stable.

<figure><img src="image_22.png" /><figcaption>Figure 21. Evolution of contact states</figcaption></figure>

<caption>Figure 22. Penetration</caption>

### 5.3.3 Analysis of penetration

In the penetration analysis, the CS-FEM-4SD model is selected. The friction coefficient is 0.2. Fig. 22 shows the penetration of the three points A, B and C during the loading process within the bi-potential framework. We can find that the contact force is different, and the corresponding penetration is also different. This is because the accuracy of the contact criterion is controlled by contact forces rather than displacements. We also find that the penetration can reach the order of $1.5 \times 10^{-17}$ mm, which shows a very high precision of our approach.

## 6. Conclusion

In this work, the CS-FEM is applied with six different smoothing domains based on quadrilateral background elements for solving contact problems. The smoothed Galerkin weak form and the discrete form are derived at first. Then multi-points coupling effect is considered within the bi-potential framework with the exact Signorini-Coulomb condition. The Uzawa algorithm is used to solve the contact force. This method can effectively simulate the three states of separating, sliding and sticking.
Through the study of numerical examples, we can get the following conclusions:

(1) As the number of degrees of freedom increases, the strain energy solution obtained by CS-FEM with different smoothing domains can converge to the reference solution with different smoothing domains. CS-FEM-1SD converges from the lower bound and CS-FEM-2SD to 16SD converges from the upper bound.

(2) When using the displacement boundary condition, the CS-FEM-1SD provides a lower bound solution while the other models provide upper bound solutions. All the results produced by CS-FEM are more accurate than those of the FEM-Q4.

(3) The solution of CS-FEM-1SD is the most accurate. When using the same background elements, its accuracy more than 5 times higher than those of the FEM-Q4. But this model could be unstable when the constraint is insufficient. CS-FEM-4SD has a higher stability in solving contact problems, and its accuracy can be more than 30-40% relative to the FEM-Q4 solution when using the same background elements.

## Acknowledgement

We gratefully acknowledge the financial support of the National Natural Science Foundation of China (Grant No. 1177020241), the National Key Research and Development Project of China (Grant No. 2017YFB0703202) and the Fundamental Research Funds for the Central Universities (Grant No.2682019CX42).

## References

[1] Wriggers, P. (2002) Computational Contact Mechanics, John Wiley & Sons.

[2] Johnson, K. L. (1985) Contact mechanics, *Journal of Tribology* **108**, 464.

[3] Zienkiewicz, O. C. and Taylor, R. L. (2005) *The Finite Element Method*.

[4] Bathe, K. J. and Chaudhary, A. (1985) A solution method for planar and axisymmetric contact problems. *International Journal for Numerical Methods in Engineering* **21**, 65-88.

[5] Weyler, R., Oliver, J., Sain, T., and Cante, J. C. (2012) On the contact domain method: a comparison of penalty and lagrange multiplier implementations, *Computer Methods in Applied Mechanics & Engineering* **205**, 68-82.

[6] Biotteau, E. and Ponthot, J. P. (2012) Modeling frictional contact conditions with the penalty method in the extended finite element framework, *Journal of Business Ethics*.

[7] Okabe, M. and Kikuchi, N. (2010) Penalty method to a two-body contact problem, *Engineering Mechanics*, 700-703.

[8] Belytschko, T. and Neal, M. O. (1991) Contact-impact by the pinball algorithm with penalty and lagrangian methods, *International Journal for Numerical Methods in Engineering* **31**, 547-572.

[9] Carpenter, N. J., Taylor, R. L., and Katona, M. G. (1991) Lagrange constraints for transient finite element surface contact, *International Journal for Numerical Methods in Engineering* **32**, 103-128.

[10] Pitkaranta, J. (1980) The finite element method with lagrange multipliers, *Mathematics of Computation* **37**, 13-30.

[11] Baillet, L., and Sassi, T. (2002) Finite element method with lagrange multipliers for contact problems with friction, *Comptes Rendus Mathematique* **334**, 917-922.

[12] H. Barbosa. (1991) The finite element method with lagrange multipliers on the boundary: circumventing the Babuka-Brezzi condition, *Computer Methods in Applied Mechanics & Engineering* **85**, 109-128.

[13] Farnood Gholami, Mostafa Nasri, József Kövecses, and Marek Teichmann. (2016) A linear complementarity formulation for contact problems with regularized friction. *Mechanism & Machine Theory* **105**, 568-582.

[14] Li, J., Zhang, H., and Pan, S. (2010) A second-order cone linear complementarity approach for contact problems with orthotropic friction law, *Chinese Journal of Solid Mechanics* **31**, 109-118.

[15] Li, J., Pan, S., and Zhang, H. (2009) A second-order cone linear complementarity approach for three-dimensional frictional contact problems, *Lixue Xuebao/Chinese Journal of Theoretical and Applied Mechanics* **41**, 869-877.

[16] Yue, J. H., Liu, G. R., Li, M., and Niu, R. P. (2018) A cell-based smoothed finite element method for multi-body contact analysis using linear complementarity formulation, *International Journal of Solids and Structures*, 110-126.

[17] De Saxcé, G. and Feng, Z. Q. (1991) New inequality and functional for contact with friction: the implicit standard material approach, *Mechanics of Structures and Machines* **19**, 301-325.

[18] De Saxcé, G. and Feng, Z. Q. (1998) The bipotential method: A constructive approach to design the complete contact law with friction and improved numerical algorithms, *Mathematical and Computer Modelling* **28**, 225-245.

[19] Joli, P., and Feng, Z. Q. (2008) Uzawa and Newton algorithms to solve frictional contact problems within the bipotential framework, *International Journal for Numerical Methods in Engineering* **73**, 317-330.
[20] Feng, Z. Q., Magnain, B., and Cros, J. M. (2006) Solution of large deformation impact problems with friction between Blatz-Ko hyperelastic bodies, *International Journal of Engineering Science* **44**, 113-126.

[21] Feng, Z. Q., Peyraut, F., and Labed, N. (2003) Solution of large deformation contact problems with friction between Blatz-Ko hyperelastic bodies, *International Journal of Engineering Science* **41**, 2213-2225.

[22] Ning, P., Feng, Z. Q., Quintero, J. A., Zhou, Y.J., and Peng, L. (2018) Uzawa algorithm to solve elastic and elastic-plastic fretting wear problems within the bipotential framework, *Computational Mechanics* **62**, 1327-1341.

[23] Peng, L., Feng, Z. Q., Joli, P., Renaud, C., and Xu, W. Y. (2019) Bi-potential and co-rotational formulations applied for real time simulation involving friction and large deformation. *Computational Mechanics* **64**, 611-623.

[24] Eck, C., Steinbach, O., and Wendland, W. L. (1999) A symmetric boundary element method for contact problems with friction, *Mathematics and Computers in Simulation* **50**, 43-61.

[25] Zhu, C. (1995) A new boundary element/mathematical programming method for contact problems with friction, *Communications in Numerical Methods in Engineering* **11**, 683-690.

[26] Li, Y., Liu, G. R., Dai, K. Y., Luan, M., Zhong, Z. H., Li, G. Y., and Han, X. (2007) Contact analysis for solids based on linearly conforming radial point interpolation method, *Computational Mechanics* **39**, 537-554.

[27] Rabczuk, T., and Belytschko, T. (2007) A three-dimensional large deformation meshfree method for arbitrary evolving cracks, *Computer Methods in Applied Mechanics and Engineering* **196**, 2777-2799.

[28] Rabczuk, T., and Ren, H. (2017) A peridynamics formulation for quasi-static fracture and contact in rock, *Engineering Geology*, 42-48.

[29] Liu, G. R., and Trung, N. T. (2010) *Smoothed Finite Element Methods*.

[30] Liu, G. R., Dai, K. Y., and Nguyen, T. T. (2007) A smoothed finite element method for mechanics problems, *Computational Mechanics* **39**, 859-877.

[31] Liu, G. R., Nguyen, T. T., Dai, K. Y., and Lam, K. Y. (2007) Theoretical aspects of the smoothed finite element method (SFEM), *International Journal for Numerical Methods in Engineering* **71**, 902-930.

[32] Thaihoang, C., Nguyenthanh, N., Nguyenxuan, H., Rabczuk, T., and Bordas, S. (2011) A cell — based smoothed finite element method for free vibration and buckling analysis of shells, *Ksce Journal of Civil Engineering* **15**, 347-361.

[33] Nguyenthoi, T., Rabczuk, T., Lamphat, T., Hohuu, V., and Phungvan, P. (2014) Free vibration analysis of cracked Mindlin plate using an extended cell-based smoothed discrete shear gap method (XCS-DSG3), *Theoretical and Applied Fracture Mechanics*, 150-163

[34] Li, E., He, Z. C., Xu, X., and Liu, G. R. (2015) Hybrid smoothed finite element method for acoustic problems, *Computer Methods in Applied Mechanics and Engineering*, 664-688.

[35] Nguyenxuan, H., Tran, L. V., Nguyenthoi, T., and Vudo, H. C. (2011) Analysis of functionally graded plates using an edge-based smoothed finite element method, *Composite Structures* **93**, 3019-3039.

[36] Dai, K. Y., and Liu, G. R. (2007) Free and forced vibration analysis using the smoothed finite element method (SFEM), *Journal of Sound and Vibration* **301**, 803-820.

[37] He, T. (2018) Towards straightforward use of cell-based smoothed finite element method in fluid-structure interaction, *Ocean Engineering*, 350-363.

[38] Li, S., Cui, X., and Li, G. (2017) Multi-physics analysis of electromagnetic forming process using an edge-based smoothed finite element method, *International Journal of Mechanical Sciences*, 244-252.

[39] Li, Y., Zhang, G., Liu, G. R., Huang, Y. N., and Zong, Z. (2013) A contact analysis approach based on linear complementarity formulation using smoothed finite element methods, *Engineering Analysis with Boundary Elements* **37**, 1244-1258.

[40] Li, Y., Liu, G. R., and Zhang, G. (2011) An adaptive NS/ES-FEM approach for 2D contact problems using triangular elements, *Finite Elements in Analysis and Design* **47**, 256-275.

[41] Klarbring, A. (1992) Mathematical programming and augmented Lagrangian methods for frictional contact problems, *Proceedings Contact Mechanics International Symposium*, 409-422.

[42] Francavilla, A. and Zienkiewicz, O. C. (1975) A note on numerical computation of elastic contact problems, *International Journal for Numerical Methods in Engineering* **9**, 913-924.

[43] Mccoll, I. R., Ding, J., and Leen, S. B. (2004) Finite element simulation and experimental validation of fretting wear, *Wear* **256**, 1114-1127.