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

HAL
open science

# A projection iterative algorithm boundary element method for the Signorini problem

Shougui Zhang, Jialin Zhu

► To cite this version:

Shougui Zhang, Jialin Zhu. A projection iterative algorithm boundary element method for the Signorini problem. Engineering Analysis with Boundary Elements, 2013, 37 (1), pp. 176-181. 10.1016/j.enganabound.2012.08.010. hal-01363584

HAL Id: hal-01363584

https://hal.science/hal-01363584

Submitted on 10 Sep 2016

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A projection iterative algorithm boundary element method for the Signorini problem

Shougui Zhang

College of Mathematics and Statistics, Chongqing University, Chongqing 401331, PR China

College of Mathematics Science, Chongqing Normal University, Chongqing 401331, PR China

Jialin Zhu

College of Mathematics and Statistics, Chongqing University, Chongqing 401331, PR China

We propose a projection iterative algorithm based on a fixed point equation for solving a certain class of Signorini problem. The satisfaction of the Signorini boundary conditions is verified in a projection iterative manner, and at each iterative step, an elliptic mixed boundary value problem is solved by a boundary element method which is suitable for any domain. We prove the convergence of the algorithm by the property of projection. The advantage of this algorithm is that it is easy to be implemented and converge quickly. Some numerical results show the accuracy and effectiveness of the algorithm.

# 1. Introduction

Domain discretization methods such as the finite element method (FEM) and the finite differential method (FDM) have been extensively applied for the numerical solution of Signorini problems [1–5]. However, these methods require the discretization of entire domain which results in high computational cost. As the Signorini conditions are given on the boundary of the domain, boundary discretization methods are suitable for the solution of such problems. Some applications of method of

fundamental solution (MFS) and boundary element method (BEM) schemes for Signorini problems can be found in [6–20].

It is worth noting that some simple and practical iterative algorithms linked with the BEM have been applied to Signorini problems [18–20]. For example, the switching algorithm with the BEM for Signorini problems of elliptic differential equations is provided in [18]. Signorini problems for the Laplace equation are solved by a decomposition–coordination method (DCM) in [19], which is based on boundary variational inequality formulations. The boundary element-linear complementary method (BE-LCM) has been employed to solve Signorini problems, and a projected successive over-relaxation iterative method is applied to the problem effectively in [20].

In recent years, the projection iterative algorithm has become a popular method for finding numerical solutions of Signorini
problems [4,5,21-26]. The semi-smooth Newton method for Signorini problems which bases on the projection algorithm was proposed in [4,5]. Some projection iterative methods were proposed for variational inequalities in [21-26].

In this paper, we propose a projection iterative algorithm based on a fixed point equation. The equivalence between the Signorini boundary condition considered and the fixed point equation on the boundary is discussed in Section 2, where we propose a new iterative algorithm for the fixed point problem and the proof of convergence of the algorithm is given in Appendix. The implementation of BEM is introduced in Section 3. Some numerical results are presented in Section 4, and we conclude the paper in Section 5.

## 2. An iterative algorithm for the Signorini problem

For the sake of simplicity, let us consider, as a model problem, a Signorini problem in two dimensions. Let $\Omega$ be a bounded domain in $\mathbb{R}^2$ with the smooth boundary $\Gamma = \Gamma_D \cup \Gamma_N \cup \Gamma_S$, and $\Gamma_S \neq \emptyset$. We want to find the solution $u$ satisfying

$$ egin{cases} 
\Delta u = f \text{ in } \Omega \\ 
u = g \text{ on } \Gamma_D \\ 
\frac{\partial u}{\partial n} = q \text{ on } \Gamma_N 
\end{cases} 
\quad (1) $$

and the following Signorini boundary conditions

$$ u \geq h, \quad \frac{\partial u}{\partial n} \geq 0, \quad (u-h) \frac{\partial u}{\partial n} = 0, \text{ on } \Gamma_S \quad (2) $$

where $f \in L^2(\Omega)$, $g \in H^{1/2}(\Gamma_D)$, $h \in H^{1/2}(\Gamma_S)$ and $q \in H^{-1/2}(\Gamma_N)$ are all given functions. It is well known that this problem has a unique solution.

In order to obtain our algorithm, we transfer the problem of satisfaction for Signorini boundary conditions considered to a fixed point problem. For $a \in \mathbb{R}$ we define the projection operator

$$ [a]_+ = \max(a, 0). $$

Then the Signorini boundary conditions (2) are equivalent to the following fixed point problem

$$ (u-h)-\left[(u-h)-c\frac{\partial u}{\partial n}\right]_+ = 0, \text{ on } \Gamma_S \quad (3) $$

where $c$ is any positive number [22,26]. By (3) we have $\min_{u \in H^1(\Omega)} u(u/\partial n) = 0$, which is equivalent to $u - h \ge 0$, $\partial u/\partial n \ge 0$ and $(u - h)\partial u/\partial n = 0$. Clearly the fixed point Eq. (3) is equivalent to the Signorini boundary conditions (2) [26].

According to the fixed point Eq (3), we propose an implicit scheme on $\Gamma_S$ as follows:

$$ u^{(k+1)} = h + \left[ (u^{(k)} - h) - c \frac{\partial u^{(k+1)}}{\partial n} \right]_+ , \quad k = 0, 1, 2, \dots \quad (4) $$

In the numerical iterative algorithm, the projection operation $[ ]_+$ needs be verified point-wise only on $\Gamma_S$. For any initial value $u^{(0)} \ge h$ and any positive constant $c > 0$, we can prove that the sequence $\{u_k\}$ generated by the algorithm converges to the unique solution of the problem, see Appendix.

The process of the iterative algorithm is the following [5]:

Step 1: Assume initially that the boundary condition on $\Gamma_S$ is

$$ u^{(0)} = h \quad \text{on } \Gamma_S \quad (5) $$

Then solve the mixed boundary value problem (1) with the boundary condition on $\Gamma_S$ (5) to obtain the normal derivative $\partial u^{(0)}/\partial n$ on $\Gamma_S$.

Choose a constant number $c > 0$, and a tolerance $\tau > 0$.

Step 2: Verify the following boundary conditions on $\Gamma_S$. If on some part of $\Gamma_S$, which we noted by $\Gamma_{SD}^{(k)}$

$$ (u^{(k)} - h) - c \frac{\partial u^{(k)}}{\partial n} < 0 \quad \text{on } \Gamma_{SD}^{(k)} $$

Then we pose the Dirichlet boundary condition on the set of $\Gamma_{SD}^{(k)}$

$$ u^{(k+1)} = h. \quad (6) $$

and the remaining part of $\Gamma_S$, we noted $\Gamma_{SR}^{(k)}$ where the following condition is imposed.

$$ u^{(k+1)} = u^{(k)} - c \frac{\partial u^{(k+1)}}{\partial n}. \quad (7) $$

Step 3: Solve the following elliptic boundary value problem

$$ 
\left\{
\begin{array}{l}
\Delta u^{(k+1)} = 0 \text{ in } \Omega \\
u^{(k+1)} = g \text{ on } \Gamma_D \\
\frac{\partial u^{(k+1)}}{\partial n} = q \text{ on } \Gamma_N \\
u^{(k+1)} = h \text{ on } \Gamma_{SD}^{(k)} \\
\frac{\partial u^{(k+1)}}{\partial n} = \frac{1}{c}(u^{(k)} - u^{(k+1)}) \text{ on } \Gamma_{SR}^{(k)}
\end{array}
\right.
\quad (8) $$

Then we obtain $u^{(k+1)}$ and its normal derivative $\partial u^{(k+1)} / \partial n = 1/c(u^{(k)} - u^{(k+1)})$ on the boundary $\Gamma_S$. Or we solve another mixed boundary value problem

$$ 
\left\{
\begin{array}{l}
\Delta u^{(k+1)} = 0 \text{ in } \Omega \\
u^{(k+1)} = g \text{ on } \Gamma_D \\
\frac{\partial u^{(k+1)}}{\partial n} = q \text{ on } \Gamma_N \\
u^{(k+1)} = h \text{ on } \Gamma_{SD}^{(k)} \\
u^{(k+1)} = u^{(k)} - c \frac{\partial u^{(k+1)}}{\partial n} \text{ on } \Gamma_{SR}^{(k)}
\end{array}
\right.
\quad (9) $$

then obtain the normal derivatives $\partial u^{(k+1)} / \partial n$ and $u^{(k+1)} = u^{(k)} - c(\partial u^{(k+1)} / \partial n)$ on the boundary $\Gamma_S$.

Step 4: If $\|u^{(k+1)} - u^{(k)}\|_{\Gamma_S} \le \tau \|u^{(k+1)}\|_{\Gamma_S}$, the iteration is terminated. Otherwise, we set $k = k+1$ and go back to Step 2.

## 3. Application of the algorithm with the boundary element method

Since the numerical method needed in the iteration process to solve the boundary value (8) or (9) should be convenient to update the data on the boundary and to check the inequalities on the boundary, so that the BEM is more appropriate for the algorithm [27,28]. We need only boundary meshing, and the checking of the inequalities is performed point-wise on the boundary nodes only in the implementation of BEM.

For simplicity, we may assume that $f=0$. If $f \ne 0$, domain integrations are involved. While a dual reciprocity method [29-31] or multiple reciprocity method [32] are used, the domain integrals are avoidable. Otherwise, if we can find a particular solution of the Poisson equation [33], the problem can be transformed to a problem of Laplace equation.

The boundary integral equation relating the potential and normal derivative on the boundary $\Gamma$ is the following

$$ 
\alpha(y)u(y) = \int_{\Gamma} \frac{\partial u(x)}{\partial n_x} u^*(x, y) ds_x \\ 
- \int_{\Gamma} u(x) \frac{\partial u^*(x, y)}{\partial n_x} ds_x, \quad y \in \Gamma \quad (10) $$

where $n_x$ is the outward unit normal to the boundary $\Gamma$ and $u^*(x, y) = -1/2\pi (\ln|x-y|)$ is the fundamental solution of Laplace
equation, $\alpha(y) = \theta(y)/2\pi$, and $\theta(y)$ is the angle at the boundary point $y \in \Gamma$, $\alpha = 1/2$ for any constant element.

When we use constant approximation for $u^{(k)}$ and $\partial u^{(k)}/\partial n$ on each boundary element, the checking of the inequalities point-wise on nodes element by element along $\Gamma_S$ in step 2 is easy to implement. Each element on $\Gamma_S$ clearly belongs to either $\Gamma_{SD}^{(k)}$ or $\Gamma_{SR}^{(k)}$.

At the $(k+1)$th iterative step, the boundary data $u^{(k+1)}$ and $\partial u^{(k+1)}/\partial n$ are linked by the discrete form of the boundary integral equation (10) while the nodes $y_i$ on $\Gamma_D$, $\Gamma_N$ and $\Gamma_{SD}^{(k)}$, which reads as

$$ 
\begin{aligned}
& \frac{1}{2} u^{(k+1)}(y_i) - \sum_{j=1}^{N} \frac{1}{2\pi} \left[ \int_{\Gamma_j} \frac{\partial \ln |x-y_i|}{\partial n_x} ds_x \right] u^{(k+1)}(x_j) \\
& = - \sum_{j=1}^{N} \frac{1}{2\pi} \left[ \int_{\Gamma_j} \ln|x-y_i| ds_x \right] \frac{\partial u^{(k+1)}(x_j)}{\partial n_x}, \quad y_i \in \Gamma_D \cup \Gamma_N \cup \Gamma_{SD}^{(k)}
\end{aligned}
\quad (11) $$

where $N$ denotes the total number of boundary nodes on $\Gamma$. If the Robin boundary conditions on $\Gamma_{SR}^{(k)}$ is applied

$$ \frac{\partial u^{(k+1)}}{\partial n} = \frac{u^{(k)} - u^{(k+1)}}{c} \quad (12) $$

the corresponding discrete form while $y_i \in \Gamma_{SR}^{(k)}$ reads

$$ 
\begin{aligned}
& \frac{1}{2} u^{(k+1)}(y_i) - \sum_{j=1}^{N} \frac{1}{2\pi} \left[ \int_{\Gamma_j} \frac{\partial \ln |x-y_i|}{\partial n_x} ds_x \right] u^{(k+1)}(x_j) \\
& \quad - \sum_{j=1}^{N} \frac{1}{2c\pi} \left[ \int_{\Gamma_j} \ln|x-y_i| ds_x \right] u^{(k+1)}(x_j) \\
& = - \sum_{j=1}^{N} \frac{1}{2c\pi} \left[ \int_{\Gamma_j} \ln|x-y_i| ds_x \right] u^{(k)}(x_j), \quad y_i \in \Gamma_{SR}^{(k)}
\end{aligned}
\quad (13) $$

If the Robin boundary conditions on $\Gamma_{SR}^{(k)}$ is given by another form

$$ u^{(k+1)} = u^{(k)} - c \frac{\partial u^{(k+1)}}{\partial n} \quad (14) $$

the corresponding discrete form while $y_i \in \Gamma_{SR}^{(k)}$ reads

$$ 
\begin{aligned}
& \frac{c}{2} \frac{\partial u^{(k+1)}(y_i)}{\partial n} - \sum_{j=1}^{N} \frac{c}{2\pi} \left[ \int_{\Gamma_j} \frac{\partial \ln |x-y_i|}{\partial n_x} ds_x \right] \frac{\partial u^{(k+1)}(x_j)}{\partial n} \\
& \quad - \sum_{j=1}^{N} \frac{1}{2\pi} \left[ \int_{\Gamma_j} \ln|x-y_i| ds_x \right] \frac{\partial u^{(k+1)}(x_j)}{\partial n} \\
& = \frac{u^{(k)}(y_i)}{2} - \sum_{j=1}^{N} \frac{1}{2\pi} \left[ \int_{\Gamma_j} \frac{\partial \ln |x-y_i|}{\partial n_x} ds_x \right] u^{(k)}(x_j), \quad y_i \in \Gamma_{SR}^{(k)}
\end{aligned}
\quad (15) $$

When the boundary conditions needed on boundary $\Gamma_D \cup \Gamma_N \cup \Gamma_{SD}^{(k)}$ for $(k+1)$th step and $u^{(k)}(x_j)$ on $\Gamma_{SR}^{(k)}$ being determined, we can solve the linear equations (11) and (13) or (11) and (15) at $(k+1)$th step, which make up well-posed problems to obtain the unknown $u^{(k+1)}$ and its normal derivative $\partial u^{(k+1)}/\partial n$ at the nodes on $\Gamma_{SR}^{(k)}$, so that we can reset the sets $\Gamma_{SD}^{(k+1)}$ and $\Gamma_{SR}^{(k+1)}$ through the checking by Step 2.

In the iterative process, we found in the numerical tests that the value of $c$ has a little effect on the number of iteration steps in general.

## 4. Numerical examples

In this Section, we choose $\tau = 10^{-6}$ and present two numerical experiments. In the first example an explicit solution is known in an annular domain. In the second example the analytical solution is not known in a square domain.

### 4.1. The Signorini problem for the Laplace equations (13,19)

We consider a test example in the annular domain $\Omega = \{x \in \mathbb{R}^2 : a < |x| < b\}$ with a Dirichlet boundary condition on the boundary $\Gamma_D = \{x \in \mathbb{R}^2 : |x| = b\}$ and a complementary condition on the Signorini boundary $\Gamma_S = \{x \in \mathbb{R}^2 : |x| = a\}$. The analytical solution in the domain $\Omega$ is given by the following complex function
$$ u(x_1, x_2) = \operatorname{Im}w(x_1 + i x_2)^3 $$

where

$$ 
\begin{aligned}
w(x_1 + i x_2) &= \sqrt{\frac{1}{2} \left( \frac{x_1^2 - x_2^2}{r^2} \right)^2 + \frac{1}{4} \left( \frac{r^2}{a^2} - \frac{a^2}{r^2} \right)^2 + \frac{1}{4} \frac{x_1^2 - x_2^2}{r^2} \left( \frac{r^2}{a^2} + \frac{a^2}{r^2} \right)} \quad \text{sign } x_1 \\
&\quad i \sqrt{\frac{1}{2} \left( \frac{x_1^2 - x_2^2}{r^2} \right)^2 + \frac{1}{4} \left( \frac{r^2}{a^2} - \frac{a^2}{r^2} \right)^2 - \frac{1}{4} \frac{x_1^2 - x_2^2}{r^2} \left( \frac{r^2}{a^2} + \frac{a^2}{r^2} \right)} \quad \text{sign } x_2
\end{aligned}
$$

and $r = \sqrt{x_1^2 + x_2^2}$, $|x| \ge a$. Then we can easily obtain all boundary conditions from the analytical solution.

The Signorini boundary conditions on $\Gamma_S$ are prescribed as

$$ u \ge h, \quad \frac{\partial u}{\partial n} \ge \varphi, \quad (u-h)\left(\frac{\partial u}{\partial n} - \varphi\right) = 0 $$

where

$$ h(x_1, x_2) = \min(0, u(x_1, x_2)) $$

$$ 
\varphi(x_1, x_2) = \begin{cases} -\frac{6}{a}(x_2 > 0, |x| = a) \\ -\frac{6}{a^2}(x_1^2 - x_2^2)^2 (x_2 \le 0, x_2 \ge -|x_1|, |x| = a) \\ 0 (x_2 < -|x_1|, |x| = a) \end{cases} $$

The analytical solution on the Signorini boundary $\Gamma_S$ is

$$ u(x_1, x_2) = -\sqrt{\max\left(0, \frac{x_2^2 - x_1^2}{a^2}\right)^3} \quad \text{sign } x_2, $$

$$ \frac{\partial u}{\partial n}(x_1, x_2) = -\frac{6}{a^3} \sqrt{\max\left(0, \frac{x_1^2 - x_2^2}{a^2}\right)} |x_1| x_2. $$

Using the BEM, this problem has been solved by the linear complementary method [13] and the DCM [19], respectively.

We take $a=0.1$ and $b=0.25$ and use the parameterizations $t \to (\arccos 2\pi t, -\arcsin 2\pi t)$ and $t \to (b\cos 2\pi t, \sin 2\pi t)$. We chose $c=10,000$ and apply our method to this problem on a uniform grid for $t$. The discretization includes 32 boundary elements on $\Gamma_S$ and 32 boundary elements on $\Gamma_D$, so $N=64$. The exact solution

<figure><img src="image_3.png" /><figcaption>Fig. 1. Analytical and approximate solution for the potential on $\Gamma_S$.</figcaption></figure>
and the numerical results for the potential on $\Gamma_S$ are shown in Fig. 1. The results for the normal derivative are shown in Fig. 2. It can be observed that our results are in a good agreement with the exact solution.

By choosing different numbers of boundary elements $N=16, 32, 64, 128$, and $256$ on $\Gamma$, we investigate the convergence behavior of our method. Our numerical results are shown in Tables 1 and 2. Table 1 indicates that numerical results converge well as the number of boundary elements $N$ increases. In Table 2, we choose different values of $c$ and compare the number of iterations between the DCM in Ref. [19] and our method. Because of using the BEM for both methods, each step of our method takes approximately the same computational cost as the DCM. Then, we can observe that the convergence of our method is faster than that of the DCM. In addition, the results show that the different values of $c$ do not affect the iteration steps of the algorithm.

<figure><img src="image_4.png" /><figcaption>Fig. 2. Analytical and approximate solution for the normal derivative on $\Gamma_S$.</figcaption></figure>

<table><caption>**Table 1**
Convergence with $N$ of the solution for the potential on $\Gamma_S$ for $c=10,000$.</caption><thead><tr><th rowspan="2">t</th><th colspan="5">Approximate solution</th><th rowspan="2">Exact solution</th></tr><tr><th>N=16</th><th>N=32</th><th>N=64</th><th>N=128</th><th>N=256</th></tr></thead><tbody><tr><td>0.125</td><td>-0.1829</td><td>-0.0795</td><td>-0.0334</td><td>-0.0137</td><td>-0.0056</td><td>0</td></tr><tr><td>0.25</td><td>0.9400</td><td>0.9778</td><td>0.9919</td><td>0.9971</td><td>0.9990</td><td>1</td></tr><tr><td>0.375</td><td>-0.1829</td><td>-0.0795</td><td>-0.0334</td><td>-0.0137</td><td>-0.0056</td><td>0</td></tr><tr><td>0.5</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>0.625</td><td>0.1829</td><td>0.0795</td><td>0.0334</td><td>0.0137</td><td>0.0056</td><td>0</td></tr><tr><td>0.75</td><td>-0.9400</td><td>-0.9778</td><td>-0.9919</td><td>-0.9971</td><td>-0.9990</td><td>-1</td></tr><tr><td>0.875</td><td>0.1829</td><td>0.0795</td><td>0.0334</td><td>0.0137</td><td>0.0056</td><td>0</td></tr><tr><td>1.0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr></tbody></table>

<table><caption>**Table 2**
Iteration with increasing $c$ and $N$.</caption><thead><tr><th rowspan="2">N</th><th colspan="6">The method in this paper</th><th rowspan="2">The method in Ref. [19]</th></tr><tr><th>c=1</th><th>c=10</th><th>c=100</th><th>c=1000</th><th>c=10,000</th><th>c=100,000</th></tr></thead><tbody><tr><td>16</td><td>3</td><td>3</td><td>3</td><td>3</td><td>3</td><td>3</td><td>21</td></tr><tr><td>32</td><td>4</td><td>4</td><td>4</td><td>4</td><td>4</td><td>4</td><td>24</td></tr><tr><td>64</td><td>4</td><td>4</td><td>4</td><td>4</td><td>4</td><td>4</td><td>28</td></tr><tr><td>128</td><td>5</td><td>5</td><td>5</td><td>5</td><td>5</td><td>5</td><td>36</td></tr><tr><td>256</td><td>6</td><td>6</td><td>6</td><td>6</td><td>6</td><td>6</td><td>48</td></tr></tbody></table>

## 4.2. The Signorini problem for the Poisson equation

In this example, we consider the following Signorini problem for Poisson's equation.

$$ \Delta u = -\cos\left(\frac{\pi}{2} + \pi x_1\right) - 1 \quad \text{in } \Omega = (0,1) \times (0,1) $$

$$ u = 0, \text{ on } \Gamma_D = \{(x_1, x_2) : x_1 = 0, 0 \le x_2 \le 1\} \\ \cup \{(x_1, x_2) : x_1 = 1, 0 \le x_2 \le 1\} $$

$$ \frac{\partial u}{\partial n} = -7x_1(1-x_1), \quad \text{on } \Gamma_N = \{(x_1, x_2) : 0 \le x_1 \le 1, x_2 = 1\} $$

with the following Signorini boundary conditions on $\Gamma_S = \{(x_1, x_2): 0 \le x_1 \le 1, x_2 = 0\}$

$$ u \le h, \quad \frac{\partial u}{\partial n} \le 0 \quad \text{and} \quad (u-h) \frac{\partial u}{\partial n} = 0 \qquad (16) $$

where $h(x_1) = 5x_1(1-x_1)(0.5-x_1)\max(x_1, 1-x_1)$. This problem was solved by the FDM and semi-Newton methods in [5]. Since $f(x)$ is simple in this example, a particular solution $u_P$ of the Poisson equation may be determined analytically as the following

$$ u_P = \frac{1}{\pi^2} \cos\left(\frac{\pi}{2} + \pi x_1\right) - \frac{x_2^2}{2} $$

then the original problem can be transformed into the following homogeneous equation through the substitution $u_H=u-u_P$; that is,

$$ \Delta u_H = 0 \quad \text{in } \Omega = (0,1) \times (0,1) $$

$$ u_H = -\frac{1}{\pi^2} \cos\left(\frac{\pi}{2} + \pi x_1\right) + \frac{x_2^2}{2}, \quad \text{on } \Gamma_D $$

$$ \frac{\partial u_H}{\partial n} = -7x_1(1-x_1) + x_2, \quad \text{on } \Gamma_N $$

with the following Signorini boundary conditions on $\Gamma_S$

$$ u_H \le 5x_1(1-x_1)(0.5-x_1) \max(x_1, 1-x_1) - \frac{1}{\pi^2} \cos\left(\frac{\pi}{2} + \pi x_1\right), \quad \frac{\partial u_H}{\partial n} \le -x_2 $$

and

$$ \left[ u_H - 5x_1(1-x_1)(0.5-x_1) \max(x_1, 1-x_1) + \frac{1}{\pi^2} \cos\left(\frac{\pi}{2} + \pi x_1\right) \right] \left[ \frac{\partial u_H}{\partial n} + x_2 \right] = 0. $$

Then we apply our method to solve the Signorini problem of homogeneous equation and obtain $u_H$. The final numerical solution of the original problem can be obtained by adding the homogeneous solution and the particular solution; that is, $u=u_H+u_P$.

<figure><img src="image_5.png" /><figcaption>Fig. 3. Approximate solution for the potential and $h(x_1)$ on $\Gamma_S$.</figcaption></figure>
The numerical results of the potential and its normal derivative for $c=10,000$ are plotted in Figs. 3 and 4 as mesh size $h=1/40$ ($N=160$). It can be seen that our results are in a good agreement with the Signorini boundary conditions (16).

In order to investigate the convergence behavior of our method, we solve the problem by choosing $c=10,000$ and different

numbers of boundary elements $N=40, 80, 160, 320$, and $640$ on $\Gamma$, and our numerical results are shown in Tables 3 and 4. We note that numerical results converge well as the number of boundary elements $N$ increases. Table 4 compares the performance of the method in [5] and our method by choosing different mesh sizes $h=1/10, 1/20, 1/40, 1/80$ and $1/160$. Here, Iter, CPU and DOF refer to the number of iterations, CPU times (sec) and the degree of freedom, respectively. The results show that, when degree of freedom of the problem is large, our method is more efficient than the method in [5].

Here, we investigate that the value of $c$ affects the convergence of the algorithm again. The results of different $c$ and $N$ are presented in Table 5 and show that the number of iterations decreases slowly as $c$ increase.

## 5. Conclusion

In this paper, we proposed an algorithm coupling BEM and a projection iteration for the solution of Signorini problems and proved its convergence. The advantage of our algorithm is that it converges rapidly and it can be implemented easily for Signorini problems defined in domain with arbitrary shape. We have put some Signorini problems for the Laplace and Poisson equations to numerical tests, and numerical results agree well with numerical solutions obtained by other methods.

## Acknowledgements

This work was supported by Educational Commission Foundation of Chongqing of China under the grant number KJ100621, the National Natural Science Foundation of China under the grant number 11101454 and the Natural Science Foundation Project of CQ CSTC under the grant number cstc2011jjA30003.

## Appendix. Convergence analysis of the algorithm

Let $\langle \cdot, \cdot \rangle$ denote the extension of the usual $L^2(\Gamma)$ scalar product to $H^{1/2}(\Gamma) \times H^{-1/2}(\Gamma)$. To investigate the convergence of the iterative algorithm in Section 2, we will need the following lemmas [26].

**Lemma 1.** Let $u^*$ be the solution of the Signorini problem, then

$$ \left\langle \frac{\partial u^{(k+1)}}{\partial n}, u^{(k)} - u^* \right\rangle_{\Gamma_s} \geq \left\langle \frac{\partial u^{(k+1)}}{\partial n}, u^{(k)} - u^{(k+1)} \right\rangle_{\Gamma_s} $$

**Proof.** Since $u^*$ is the solution of the Signorini problem, we consider the $(k+1)$th iteration and obtain

$$ \begin{aligned} \left\langle \frac{\partial u^*}{\partial n}, u^{(k+1)} - u^* \right\rangle_{\Gamma_s} &= \left\langle \frac{\partial u^*}{\partial n}, (u^{(k+1)} - h) + (h - u^*) \right\rangle_{\Gamma_s} \\ &= \left\langle \frac{\partial u^*}{\partial n}, u^{(k+1)} - h \right\rangle_{\Gamma_s} - \left\langle \frac{\partial u^*}{\partial n}, u^* - h \right\rangle_{\Gamma_s} = \left\langle \frac{\partial u^*}{\partial n}, u^{(k+1)} - h \right\rangle_{\Gamma_s} \geq 0. \end{aligned} $$

Applying Green's formula, we have

$$ \begin{aligned} \left| \nabla(u^{(k+1)} - u^*) \right|^2_{L^2(\Omega)} &= \left\langle \frac{\partial(u^{(k+1)} - u^*)}{\partial n}, u^{(k+1)} - u^* \right\rangle_{\Gamma_s} \\ &= \left\langle \frac{\partial u^{(k+1)}}{\partial n}, u^{(k+1)} - u^* \right\rangle_{\Gamma_s} - \left\langle \frac{\partial u^*}{\partial n}, u^{(k+1)} - u^* \right\rangle_{\Gamma_s} \geq 0. \end{aligned} $$

Consequently,

$$ \left\langle \frac{\partial u^{(k+1)}}{\partial n}, u^{(k+1)} - u^* \right\rangle_{\Gamma_s} \geq \left\langle \frac{\partial u^*}{\partial n}, u^{(k+1)} - u^* \right\rangle_{\Gamma_s} \geq 0 $$

<figure><img src="image_6.png" /><figcaption>Fig. 4. Approximate solution for the normal derivative on $\Gamma_s$.</figcaption></figure>

<table><caption>**Table 3**
Convergence with $N$ of the solution for the potential on $\Gamma_s$.</caption><thead><tr><th rowspan="2">(x<sub>1</sub>, x<sub>2</sub>)</th><th colspan="5">Approximate solution</th></tr><tr><th>N=40</th><th>N=80</th><th>N=160</th><th>N=320</th><th>N=640</th></tr></thead><tbody><tr><td>(0.1,0)</td><td>-0.0102</td><td>-0.0099</td><td>-0.0098</td><td>-0.0098</td><td>-0.0098</td></tr><tr><td>(0.2,0)</td><td>-0.0261</td><td>-0.0256</td><td>-0.0255</td><td>-0.0255</td><td>-0.0255</td></tr><tr><td>(0.3,0)</td><td>-0.0443</td><td>-0.0436</td><td>-0.0435</td><td>-0.0435</td><td>-0.0435</td></tr><tr><td>(0.4,0)</td><td>-0.0626</td><td>-0.0616</td><td>-0.0615</td><td>-0.0615</td><td>-0.0615</td></tr><tr><td>(0.5,0)</td><td>-0.0808</td><td>-0.0793</td><td>-0.0791</td><td>-0.0791</td><td>-0.0790</td></tr><tr><td>(0.6,0)</td><td>-0.1032</td><td>-0.1003</td><td>-0.1000</td><td>-0.0999</td><td>-0.0999</td></tr><tr><td>(0.7,0)</td><td>-0.1470</td><td>-0.1470</td><td>-0.1470</td><td>-0.1470</td><td>-0.1470</td></tr><tr><td>(0.8,0)</td><td>-0.1920</td><td>-0.1920</td><td>-0.1920</td><td>-0.1920</td><td>-0.1920</td></tr><tr><td>(0.9,0)</td><td>-0.1620</td><td>-0.1620</td><td>-0.1620</td><td>-0.1620</td><td>-0.1620</td></tr></tbody></table>

<table><caption>**Table 4**
Performance of the semi-smooth Newton method and the proposed method.</caption><thead><tr><th rowspan="2">h</th><th colspan="3">The method in Ref. [5]</th><th colspan="3">The method in this paper</th></tr><tr><th>Iter</th><th>CPU</th><th>DOF (N)</th><th>Iter</th><th>CPU</th><th>DOF (N)</th></tr></thead><tbody><tr><td>1/10</td><td>3</td><td>0.047</td><td>9 × 11</td><td>3</td><td>0.188</td><td>4 × 10</td></tr><tr><td>1/20</td><td>3</td><td>0.125</td><td>19 × 21</td><td>3</td><td>0.312</td><td>4 × 20</td></tr><tr><td>1/40</td><td>4</td><td>4.797</td><td>39 × 41</td><td>5</td><td>0.859</td><td>4 × 40</td></tr><tr><td>1/80</td><td>6</td><td>18.98</td><td>79 × 81</td><td>5</td><td>3.141</td><td>4 × 80</td></tr><tr><td>1/160</td><td>9</td><td>75.42</td><td>159 × 161</td><td>6</td><td>15.05</td><td>4 × 160</td></tr></tbody></table>

<table><caption>**Table 5**
Iteration with increasing $c$ and $N$.</caption><thead><tr><th>c/N</th><th>1</th><th>10</th><th>100</th><th>1000</th><th>10,000</th><th>100,000</th></tr></thead><tbody><tr><td>40</td><td>11</td><td>7</td><td>5</td><td>4</td><td>3</td><td>3</td></tr><tr><td>80</td><td>10</td><td>7</td><td>5</td><td>4</td><td>3</td><td>3</td></tr><tr><td>160</td><td>10</td><td>7</td><td>6</td><td>5</td><td>5</td><td>5</td></tr><tr><td>320</td><td>10</td><td>7</td><td>6</td><td>5</td><td>5</td><td>5</td></tr><tr><td>640</td><td>10</td><td>7</td><td>6</td><td>6</td><td>6</td><td>6</td></tr></tbody></table>
Then

$$
\left\langle \frac{\partial u^{(k+1)}}{\partial n}, u^{(k)} - u^* \right\rangle_{\Gamma_s} = \left\langle \frac{\partial u^{(k+1)}}{\partial n}, u^{(k)} - u^{(k+1)} + u^{(k+1)} - u^* \right\rangle_{\Gamma_s} \\
\geq \left\langle \frac{\partial u^{(k+1)}}{\partial n}, u^{(k)} - u^{(k+1)} \right\rangle_{\Gamma_s} .
$$

**Lemma 2.** Let $u^*$ be the solution of the Signorini problem, then

$$
\|u^{(k+1)} - u^*\|_{\Gamma_s}^2 \le \|u^{(k)} - c \frac{\partial u^{(k+1)}}{\partial n} - u^*\|_{\Gamma_s}^2 - \|u^{(k+1)} - u^{(k)} + c \frac{\partial u^{(k+1)}}{\partial n}\|_{\Gamma_s}^2
$$

**Proof.** According to the algorithm (4), we use the projection and have

$$
\left\langle u^{(k+1)} - (u^{(k)} - c \frac{\partial u^{(k+1)}}{\partial n}), u^{(k+1)} - u^* \right\rangle_{\Gamma_s} \le 0.
$$

Hence,

$$
\begin{align*}
\|u^{(k+1)} - u^*\|_{\Gamma_s}^2 &= \|u^{(k+1)} - \left(u^{(k)} - c \frac{\partial u^{(k+1)}}{\partial n}\right) + \left(u^{(k)} - c \frac{\partial u^{(k+1)}}{\partial n}\right) - u^*\|_{\Gamma_s}^2 \\
&= \left\| \left(u^{(k)} - c \frac{\partial u^{(k+1)}}{\partial n}\right) - u^* \right\|_{\Gamma_s}^2 + \left\| u^{(k+1)} - \left(u^{(k)} - c \frac{\partial u^{(k+1)}}{\partial n}\right) \right\|_{\Gamma_s}^2 \\
&\quad + 2 \left\langle u^{(k+1)} - \left(u^{(k)} - c \frac{\partial u^{(k+1)}}{\partial n}\right), \left(u^{(k)} - c \frac{\partial u^{(k+1)}}{\partial n}\right) - u^* \right\rangle_{\Gamma_s} \\
&= \left\| \left(u^{(k)} - c \frac{\partial u^{(k+1)}}{\partial n}\right) - u^* \right\|_{\Gamma_s}^2 - \left\| u^{(k+1)} - \left(u^{(k)} - c \frac{\partial u^{(k+1)}}{\partial n}\right) \right\|_{\Gamma_s}^2 \\
&\quad + 2 \left\langle u^{(k+1)} - \left(u^{(k)} - c \frac{\partial u^{(k+1)}}{\partial n}\right), u^{(k+1)} - u^* \right\rangle_{\Gamma_s} \\
&\le \left\| \left(u^{(k)} - c \frac{\partial u^{(k+1)}}{\partial n}\right) - u^* \right\|_{\Gamma_s}^2 - \left\| u^{(k+1)} - \left(u^{(k)} - c \frac{\partial u^{(k+1)}}{\partial n}\right) \right\|_{\Gamma_s}^2.
\end{align*}
$$

Using the above Lemmas, we can prove the following conver-
gence theorem for the iterative algorithm.

**Theorem 3.** Let {$u^{(k)}$} be the sequence generated by the algorithm of Section 2. Then $u^{(k)}$ converges to the unique solution $u^*$ of the Signorini problem as $k \to \infty$.

**Proof.** Since $u^*$ is the solution of the Signorini problem, by the Lemmas and the projection principle we have

$$
\begin{align*}
\|u^{(k+1)} - u^*\|_{\Gamma_s}^2 &\le \|u^{(k)} - u^*\|_{\Gamma_s}^2 - 2c \left\langle \frac{\partial u^{(k+1)}}{\partial n}, u^{(k)} - u^* \right\rangle_{\Gamma_s} \\
&\quad + 2c \left\langle \frac{\partial u^{(k+1)}}{\partial n}, u^{(k)} - u^{(k+1)} \right\rangle_{\Gamma_s} - \|u^{(k+1)} - u^{(k)}\|_{\Gamma_s}^2 \\
&\le \|u^{(k)} - u^*\|_{\Gamma_s}^2 - 2c \left\langle \frac{\partial u^{(k+1)}}{\partial n}, u^{(k)} - u^{(k+1)} \right\rangle_{\Gamma_s} \\
&\quad + 2c \left\langle \frac{\partial u^{(k+1)}}{\partial n}, u^{(k)} - u^{(k+1)} \right\rangle_{\Gamma_s} - \|u^{(k+1)} - u^{(k)}\|_{\Gamma_s}^2 \\
&= \|u^{(k)} - u^*\|_{\Gamma_s}^2 - \|u^{(k+1)} - u^{(k)}\|_{\Gamma_s}^2.
\end{align*}
$$

Clearly, the sequence {$u^{(k)}$} is bounded, and

$$
\sum_{k=0}^{\infty} \|u^{(k+1)} - u^{(k)}\|_{\Gamma_s}^2 \le \sum_{k=0}^{\infty} \|u^{(k)} - u^*\|_{\Gamma_s}^2 - \|u^{(k+1)} - u^*\|_{\Gamma_s}^2 < +\infty
$$

implying that {$u^{(k)}$} is a Cauchy sequence. Let $\lim_{k \to \infty} u^{(k)} = u^{(\infty)}$, then
$u^{(\infty)}$ satisfies the fixed point Eq. (3). So the sequence {$u^{(k)}$}
converges to the unique solution $u^*$ of the Signorini problem,
completing the proof.

**References**

[1] Aitchison JM, Lacey AA, Shillor M. A model for an electropaint process. IMA J Appl Math 1984;33:17-31.

[2] Aitchison JM, Elliott CM, Ockendon JR. Percolation in gently sloping beaches. IMA J Appl Math 1983;30:269-87.

[3] Poole MW, Aitchison JM. Numerical model of an electropaint process with applications to the automotive industry. IMA J Math Appl Bus Ind 1997;8:347-60.

[4] Ito K, Kunisch K. Semi-smooth Newton methods for variational inequalities of the first kind. Math Model Numer Anal 2003;37:41-62.

[5] Ito K, Raleigh, Kunisch K, Graz. Semi-smooth Newton methods for the Signorini problem. Appl Math 2008;53:455-6.

[6] Bogomolny A. Fundamental solutions method for elliptic boundary value problems. J Numer Anal 1985;22:644-69.

[7] Poullikkas A, Karageorghis A, Georgiou G. The method of fundamental solutions for Signorini problems. IMA J Numer Annual 1998;18:273-85.

[8] Poullikkas A, Karageorghis A, Georgiou G. The numerical solution of three-dimensional Signorini problems with the method of fundamental solution. Eng Anal Boundary Elem 2001;25:221-7.

[9] Karageorghis A. Numerical solution of a shallow dam problem by a boundary element method. Comput Methods Appl Mech Eng 1987;61:265-76.

[10] Panagiotopoulos PD. Multivalued boundary integral equations for inequality problems. The convex case Acta Mech 1987;70:145-67.

[11] Han HD. A direct boundary element method for Signorini problems. Math Comput 1990;55:115-28.

[12] Han HD, Guan Z, Yu CQ. The canonical boundary element analysis of a model for an electropaint process—The canonical boundary element approximation of a Signorini problem. J Appl Math of Chin Univ 1994;3:101-11.

[13] Leontiev A, Huacasi W, Herskovits J. An optimization technique for the solution of the Signorini problem using the boundary element method. Struct Multi Optim 2002;24:72-7.

[14] Piermatei Filho O, Leontiev A. An optimization approach for unconfined seepage problem with semipermeable conditions. Struct Multi Optim 2009;39:581-8.

[15] Maischak M, Stephan EP. Adaptive hp-versions of BEM for Signorini problems. Appl Numer Math 2005;54:425-49.

[16] Panzeca T, Salerno M, Terravecchia S, Zito L. The symmetric boundary element method for unilateral contact problems. Comput Methods Appl Mech Eng 2008;197:2667-79.

[17] Polizzotto C. Variational boundary-integral-equation approach to unilateral contact problems in elasticity. Comp Mech 1993;13:100-15.

[18] Aitchison JM, Poole MW. A numerical algorithm for the solution of Signorini problems. J Comput Appl Math 1998;94:55-67.

[19] Spann W. On the boundary element method for the Signorini problem of the Laplacian. Numer Math 1993;65:337-56.

[20] Zhang SG, Zhu JL. The boundary element-linear complementary method for the Signorini problem. Eng Anal Boundary Elem 2012;36:112-7.

[21] Noor MA. New approximation schemes for general variational inequalities. J Math Anal Appl 2000;251:217-29.

[22] Noor MA. Some developments in general variational inequalities. Appl Math Comput 2004;52:199-277.

[23] Noor MA. Extended general variational inequalities. Appl Math Lett 2009;22:182-6.

[24] Noor MA. Projection iterative methods for extended general variational inequalities. J Appl Math Comput 2010;32:83-95.

[25] Noor MA. Some iterative schemes for general mixed variational inequalities. J Appl Math Comput 2010;34:57-70.

[26] Han JY, Xiu NH, Qi HD. Nonlinear Complementarity Theory and Algorithm. Shanghai: Shanghai Science and Technology Press; 2006.

[27] Brebbia CA. The Boundary Element Method for Engineers. London: Pentech Press; 1978.

[28] Zhu JL, Yuan ZQ. Boundary Element Analysis. Beijing: Science Press; 2009.

[29] Farcas A, Elliott L, Ingham DB, Lesnic D. The dual reciprocity boundary element method for solving Cauchy problems associated to the Poisson equation. Eng Anal Boundary Elem 2003;27:955-62.

[30] Wrobel LC, Brebbia CA. The dual reciprocity boundary element formulation for nonlinear diffusion problems. Compt Mech Eng 1987;65:147-64.

[31] Jumarhon B, Amini S, Chen K. On the boundary element dual reciprocity method. J Eng Anal Boundary Elem 1997;20:205-11.

[32] Nowak AJ, Neves AC. The Multiple Reciprocity Boundary Element Method. Computational Mechanics Publication; 1994.

[33] Alves CJS, Chen CS. A new method of fundamental solutions applied to nonhomogeneous elliptic problems. Adv Comput Math 2005;23:125-42.