Application of the method of fundamental solutions to 2D and 3D Signorini problems

Engineering Analysis with Boundary Elements 58 (2015) 48–57

Contents lists available at ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

Application of the method of fundamental solutions to 2D and 3D Signorini problems Hongyan Zheng, Xiaolin Li n College of Mathematics Science, Chongqing Normal University, Chongqing 400047, PR China

art ic l e i nf o

a b s t r a c t

Article history: Received 12 January 2015 Received in revised form 24 February 2015 Accepted 14 March 2015 Available online 2 April 2015

This paper presents an application of the method of fundamental solutions (MFS) for the numerical solution of 2D and 3D Signorini problems. In our application, by using a projection technique to tackle the nonlinear Signorini boundary inequality conditions, the original Signorini problem is transformed into a sequence of linear elliptic boundary value problems and then solved by the MFS. Convergence and efficiency of the present MFS is proved theoretically and verified numerically. & 2015 Elsevier Ltd. All rights reserved.

Keywords: Method of fundamental solutions Signorini problem Projection iterative Convergence

1. Introduction Signorini problems come up in the modeling of many realistic science and engineering applications such as the shallow dam problem [1–4], the electropaint process [3–7], the unilateral contact problem [8,9], the free boundary problem [10,11], and so on. Besides, in the theory of variational inequalities [12,13], a broad class of problems arising in industry, social, economics, finance, pure and applied sciences give rise to Signorini problems. The numerical solution of Signorini problems has been frequently dominated by classical numerical methods such as the finite difference method and the finite element method (FEM) [1,6,8]. These methods require domain meshing which is often arduous, computationally expensive, and fraught with pitfalls. In Signorini problems, the boundary potential and its normal derivative alternate on the Signorini boundary in conjunction with certain nonlinear boundary inequality conditions. To obtain the solution in the domain, we need first to determine the number and position where the change from one type of boundary condition to the other occurs. Therefore, the primary focus in solving Signorini problems is on the Signorini boundary of the domain and thus, the boundary element method (BEM) is particularly suitable for the approximate solution of such problems. Some applications of the BEM for Signorini problems can be found in Refs. [2,3,9,11,14–17].

n

Corresponding author. E-mail address: [email protected] (X. Li).

http://dx.doi.org/10.1016/j.enganabound.2015.03.008 0955-7997/& 2015 Elsevier Ltd. All rights reserved.

The BEM reduces the computational dimensions of the original problem by one and gives a simple discretization of infinite domain problems [18], but it involves the generation of elements on the boundary surface and the computation of some complex singular integrals on boundary elements. In some cases, these processes can also be very difficult and computationally expensive. To alleviate the meshing-related issues, some boundary type meshless methods, such as the boundary node method (BNM) [19], the boundary point interpolation method (BPIM) [20], the hybrid BNM [21–25], the Galerkin BNM [26,27], the dual BNM [28,29] and the boundary element-free method (BEFM) [30–32], have been developed by introducing meshless shape functions into boundary integral equations. In recent years, the BNM [33], the BPIM [34] and the BEFM [35] have been extended to 2D Signorini problems. However, these boundary type meshless methods still involve the computation of boundary integrals. The method of fundamental solutions (MFS) [36,37] is a boundary type meshless method, in which the solution is approximated as a linear combination of fundamental solutions. The MFS eliminates the issue of computing integrals required in the BEM and the rest boundary type methods aforementioned. Being mathematically simple, integration-free, easy-to-program, truly meshless and the extensible to multidimensional problems make the MFS very attractive in solving boundary value problems. In 1998, Poullikkas et al. pioneered the application of the MFS to 2D Signorini problems [4]. In 2001, they also applied the MFS to a special 3D Signorini problem, known as the electropainting problem [5]. In their works, the original problem is reformulated as a nonlinear least-squares problem with nonlinear inequality constraints, and then solved by a special

H. Zheng, X. Li / Engineering Analysis with Boundary Elements 58 (2015) 48–57

least-squares minimization routine to accommodate constraints. However, the success of their application depends severely on the quality of the optimization algorithm used to handle the nonlinear constrained problem. In this paper, a new application of the MFS is developed for boundary-only analysis of Signorini problems. The numerical formulae are valid for 2D and 3D Signorini problems and also valid for both interior and exterior problems simultaneously. In our application, the nonlinear boundary inequality conditions are incorporated naturally into an iterative scheme by using a projection operator [12,13]. Then, the original Signorini problem is transformed into a sequence of well-posed linear boundary value problems. Finally, the MFS is used iteratively for solving these linear problems. The application of the MFS in this paper, differs from that in Refs. [4,5], involves only linear boundary value problems and linear system of algebraic equations, and thus avoids the design of a quite sophisticated and time-consuming optimization algorithm. As a result, the present application is expected to have higher computational speed and efficiency. Convergence and efficiency of the present MFS is also proved theoretically and verified numerically in detail. The following discussions begin with a detailed numerical implementation of the MFS for Signorini problems in Section 2. Then, Section 3 provides the associated iterative algorithm and convergence analysis. Finally, numerical examples and conclusions are given in Sections 4 and 5, respectively.

2. The MFS for Signorini problems Let Ω be a d-dimensional domain in R , where d ¼2 or d 3 denotes the spatial dimension. A generic
T point in R is denoted as x ¼ ðx1 ; x2 ; …; xd ÞT or y ¼ y1 ; y2 ; …; yd . Let Γ D , Γ N and Γ S a ∅ be three disjointed parts constituting the boundary Γ of Ω. Consider the following Signorini problem:

ΔuðxÞ ¼ 0;

xAΩ

ð1Þ

uðxÞ ¼ ϕðxÞ;

x A ΓD

ð2Þ

qðxÞ ¼ φðxÞ;

x A ΓN

ð3Þ

and uðxÞ ¼ ϕðxÞ;

qðxÞ o φðxÞ;

x A ΓS

ð4Þ

qðxÞ ¼ φðxÞ;

x A ΓS

ð5Þ

or uðxÞ o ϕðxÞ;

where u A H 1 ðΩÞ is the unknown function, q ¼ ∂u=∂n is the normal derivative of u, n ¼ ðn1 ; n2 ; …; nd ÞT is the unit outward normal to Γ , and ϕ A H 1=2 ðΓ D [ Γ S Þ and φ A H  1=2 ðΓ N [ Γ S Þ are prescribed boundary functions. To obtain the solution u in Ω, we need first to determine on which parts of Γ S boundary conditions (4) apply, and on the remaining parts boundary conditions (5) apply. Then, the solution u and its derivatives in Ω can be solved using a standard numerical method such as the FEM, the BEM and meshless methods. As a result, the primary focus in solving Signorini problem (1)–(5) is on the Signorini boundary of the domain and thus, the inequality Signorini boundary conditions (4) and (5) need to be tackled efficiently. For doing this, it can be verified firstly that conditions (4) and (5) are equivalent to the following boundary conditions: uðxÞ r ϕðxÞ;

qðxÞ r φðxÞ;

ðuðxÞ  ϕðxÞÞðqðxÞ  φðxÞÞ ¼ 0;

x A ΓS ð6Þ

Then, let [12,13] Pa≔minða; 0Þ;

aAR

where P: R-R  [ f0g is a projection operator. And let qðxÞ  φðxÞ ¼ P½ðqðxÞ  φðxÞÞ  αðuðxÞ  ϕðxÞÞ;

ð7Þ

x A ΓS

ð8Þ

where α is an arbitrary positive constant. In light of Eq. (7), from Eq. (8) we have q r φ. If q ¼ φ, then ðq  φÞ  αðu  ϕÞ Z0, and thus u r ϕ. Otherwise, if q o φ, then ðq  φÞ  αðu  ϕÞ o0, and thus recalling again Eq. (8) leads to u ¼ ϕ. As a result, from Eq. (8) we can deduce the Signorini boundary conditions (4) and (5), and thus deduce Eq. (6). On the other hand, from Eq. (6) we can deduce Eq. (8) immediately. Summarizing, we have shown that Eq. (6) is equivalent to Eq. (8). In view of Eq. (8), an implicit iterative scheme can be defined for numerical computation as qðk þ 1Þ ðxÞ ¼ φðxÞ þ P

h

  i qðkÞ ðxÞ  φðxÞ  α uðk þ 1Þ ðxÞ  ϕðxÞ ;

k ¼ 0; 1; 2; … x A Γ S

ð9Þ where the superscript ðkÞ denotes the value at the kth iteration. Since the projection operator P is also nonlinear, Eq. (9) cannot be implemented directly. In this study, this operator is checked point-wise. Let fxi gN i ¼ 1  Γ be a set of N boundary nodes which are necessary for the implementation of the MFS. After the kth iteration, both uðkÞ ðxi Þ and qðkÞ ðxi Þ are known for all xi A Γ S . Thus, if on Γ SN  Γ S inequality     ð10Þ qðkÞ ðxi Þ  φðxi Þ  α uðkÞ ðxi Þ  ϕðxi Þ 4 0 is true for all xi A Γ SN , and for xi A Γ SR ≔Γ S \Γ SN this inequality is false, then from Eq. (9) we use qðk þ 1Þ ðxi Þ ¼ φðxi Þ;

d

49

xi A Γ SN

ð11Þ

and

  qðk þ 1Þ ðxi Þ ¼ qðkÞ ðxi Þ  α uðk þ 1Þ ðxi Þ  ϕðxÞi ;

xi A Γ SR

ð12Þ

for the ðk þ 1Þth iteration. It should be stressed that the MFS is incidental to check and update the inequality given by Eq. (10). According to Eqs. (11) and (12), the original Signorini problem (1)–(5) is reduced to the following linear problem:

Δuðk þ 1Þ ðxÞ ¼ 0;

xAΩ

ð13Þ

uðk þ 1Þ ðxÞ ¼ ϕðxÞ;

x A ΓD

ð14Þ

qðk þ 1Þ ðxÞ ¼ φðxÞ;

x A ΓN

ð15Þ

qðk þ 1Þ ðxi Þ ¼ φðxi Þ;

xi A Γ SN

αuðk þ 1Þ ðxi Þ þ qðk þ 1Þ ðxi Þ ¼ qðkÞ ðxi Þ þ αϕðxi Þ;

ð16Þ xi A Γ SR

ð17Þ

where k ¼ 0; 1; 2; … The boundary value problem (13)–(17) is now solved using the MFS. With the help of the MFS, the solution of this problem can be approximated as þ 1Þ uðk þ 1Þ ðxÞ  uðk ðxÞ ¼ N

N X j¼1

  þ 1Þ aðk U x; yj ; j

x A Ω ¼ Ω [ Γ ; k ¼ 0; 1; 2; … þ 1Þ aðk j

ð18Þ

is the jth unknown coefficient at the ðk þ 1Þth where iteration, yj is the source point placed on a fictitious boundary Γ 0 , N is the number of source points, and   8  ln x  yj ; d ¼ 2 > <   1   1  U x; yj ¼ d¼3 2ðd  1Þπ > : x  y  ; j is the fundamental solution of Laplace's equation.

50

H. Zheng, X. Li / Engineering Analysis with Boundary Elements 58 (2015) 48–57

From Eq. (18), the normal derivative qðk þ 1Þ ðxÞ≔∂uðk þ 1Þ ðxÞ=∂nx can be approximated as þ 1Þ ðxÞ ¼ qðk þ 1Þ ðxÞ  qðk N

N X j¼1

  þ 1Þ aðk Q x; yj ; j

xAΓ

ð19Þ

    where Q x; yj ≔∂U x; yj =∂nx . þ 1Þ To obtain the unknown coefficient aðk , using Eq. (18) and j collocating the Dirichlet boundary condition (14) for all boundary nodes xi A Γ D yield N X j¼1

  þ 1Þ aðk U x i ; y j ¼ ϕðx i Þ j

for all xi A Γ D

ð20Þ

Then, using Eq. (19) and collocating the Neumann boundary condition (15) and (16) for all boundary nodes xi A Γ N [ Γ SN yield N X j¼1





þ 1Þ aðk Q xi ; yj ¼ φðxi Þ j

for all xi A Γ N [ Γ SN

ð21Þ

Finally, using Eqs. (18) and (19) and collocating the Robin boundary condition (17) for all boundary nodes xi A Γ SR yield N X j¼1

     þ 1Þ aðk αU xi ; yj þ Q xi ; yj ¼ qk ðxi Þ þ αϕðxi Þ j

for all xi A Γ SR ð22Þ

3. Check the boundary conditions on Γ S using Eq. (10) and determine Γ SN and Γ SR . 4. Form the linear algebra system given by Eqs. (20)–(22), and þ 1Þ solve the system to obtain aðk . j þ 1Þ ðk þ 1Þ 5. Obtain uðk ð x Þ and q ð x i i Þ for all nodes x i A Γ S from Eqs. N N (23) and (24), respectively. 6. Compute the relative error, sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 X  2 X  ðk þ 1Þ þ 1Þ eðuÞ ¼ uN ðxi Þ  uðkÞ = uðk ðx i Þ N ðx i Þ N xi A Γ S

xi A Γ S

7. Stop the algorithm if eðuÞ r τ. Otherwise, update k to k þ 1, and go to step 3. In Algorithm 3.1, the iterative process continues until the stop condition eðuÞ r τ holds. In the subsequent numerical examples, unless otherwise specified, the allowable tolerance is taken as τ ¼ 10  6 . Obviously, the number of iteration is determined by the stop condition.
 Lemma 3.1. Let uðkÞ ; qðkÞ be the analytical solution of the boundary value problem (13)–(17) and let ðu; qÞ be the analytical solution of the Signorini problem (1)–(5), then for any α 4 0, there holds   lim uðkÞ ; qðkÞ ¼ ðu; qÞ k-1

Eqs. (20)–(22) contain a set of coupled N linear algebraic equaþ 1Þ tions, which are solved together for N unknowns aðk . Then, the j ðk þ 1Þ potential u ðxÞ at an internal point x A Ω can be computed approximately using Eq. (18), and the MFS approximate values of uðk þ 1Þ ðxÞ and qðk þ 1Þ ðxÞ at any boundary point x A Γ can be computed respectively using Eqs. (18) and (19). Especially, uðk þ 1Þ ðxi Þ and qðk þ 1Þ ðxi Þ at boundary nodes xi A Γ S can be approximated respectively as uðk þ 1Þ ðxi Þ  uðNk þ 1Þ ðxi Þ ¼

N X j¼1



þ 1Þ aðk U xi ; yj j



for all xi A Γ S

ð23Þ

for all xi A Γ S

ð24Þ

and qðk þ 1Þ ðxi Þ  qðNk þ 1Þ ðxi Þ ¼

N X j¼1

  þ 1Þ aðk Q xi ; yj j

If boundary singularity exists, accuracy of the MFS solution decreases near the singular point. In Ref. [38], a modified MFS that incorporates the singularity has been developed and applied to some singular problems. This modified MFS gives more accurate results than the standard MFS and predicts the nature of the leading singular term.

3. Iterative algorithm and convergence analysis In this section, iterative algorithm of the MFS for Signorini problems (1)–(5) is concluded firstly. Then, convergence analysis of this algorithm is presented in detail. Algorithm 3.1. Choose a tolerance τ 4 0 and set k ¼0. 1. Select N collocation nodes fxi gN i ¼ 1 on the true boundary Γ and 0 place N source nodes fyj gN on a fictitious boundary Γ . j¼1 2. Assume initially that the boundary condition on Γ S is qð0Þ ¼ φ

ð25Þ

then solve the boundary value problem 8 ð0Þ > < Δu ¼ 0 in Ω uð0Þ ¼ ϕ on Γ D > : qð0Þ ¼ φ on Γ [ Γ N

S

ð0Þ

by the MFS to obtain u

on Γ S .



 ϕ 4 0 on Γ S , then Proof. In Eq. (9), if qðkÞ  φ
 α uðk þ 1Þ   φ ¼ qðk þ 1Þ . Thus, qðk þ 1Þ  qðkÞ þ α uðk þ 1Þ  ϕ o 0, and using q r φ yields qðk þ 1Þ Zq. Consequently,   D E qðk þ 1Þ  qðkÞ þ α uðk þ 1Þ  ϕ ; qðk þ 1Þ  q r0 ð26Þ ΓS

where the duality 〈; 〉Γ S is defined for functions w and v as R 〈w; v〉Γ S ¼ Γ S wv dΓ ,
and the 
norm J J Γ S is defined as pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi J v J Γ S ¼ 〈v; v〉Γ S . If qðkÞ  φ  α uðk þ 1Þ  ϕ r 0 on Γ S , then Eq.
 (9) is simplified as qðk þ 1Þ ¼ qðkÞ  α uðk þ 1Þ  ϕ . In this case, Eq. (26) is also true. Using Eq. (26) yields    ðk þ 1Þ Jq  q J 2Γ S ¼ J qðk þ 1Þ  qðkÞ  α uðk þ 1Þ  ϕ    þ qðkÞ  α uðk þ 1Þ  ϕ  q J 2Γ S     ¼ J qðkÞ  q  α uðk þ 1Þ  ϕ J 2Γ S  J qðk þ 1Þ  qðkÞ þ α uðk þ 1Þ  ϕ J 2Γ S   D E þ 2 qðk þ 1Þ  qðkÞ þ α uðk þ 1Þ  ϕ ; qðk þ 1Þ  q ΓS     ðkÞ ðk þ 1Þ 2 ðk þ 1Þ ðkÞ r Jq qα u  ϕ J ΓS  J q  q þ α uðk þ 1Þ  ϕ J 2Γ S D E ¼ J qðkÞ  q J 2Γ S  J qðk þ 1Þ  qðkÞ J 2Γ S  2α uðk þ 1Þ  ϕ; qðk þ 1Þ q

ΓS

ð27Þ From Eqs. (6) and (9) it follows that 〈u  ϕ; q  φ〉Γ S ¼ 0, u r ϕ and qðk þ 1Þ r φ. Then, applying Green's formula, in view of Eqs. (1) and (13), we have D E uðk þ 1Þ  ϕ; qðk þ 1Þ  q

D E ¼ uðk þ 1Þ  u; qðk þ 1Þ  q ΓS ΓS D E  ðk þ 1Þ þ u  ϕ; q φ  u  ϕ; q  φ Γ ΓS

  ¼ J ∇ uðk þ 1Þ  u J 2L2 Ω ð Þ D E þ u  ϕ; qðk þ 1Þ  φ

ΓS

Z0

S

H. Zheng, X. Li / Engineering Analysis with Boundary Elements 58 (2015) 48–57

51

4.1. Problem in an annular domain Thus, from Eq. (27) one gets J qðk þ 1Þ  q J 2Γ S r J qðkÞ q J 2Γ S  J qðk þ 1Þ  qðkÞ J 2Γ S r J qðkÞ  q J 2Γ S and then invoking Eq. (25) provides 1 X k¼0

J qðk þ 1Þ  qðkÞ J 2Γ S r

1  X k¼0

J qðkÞ  q J 2Γ S  J qðk þ 1Þ  q J 2Γ S



r J qð0Þ  q J 2Γ S ¼ J φ  q J 2Γ S Besides, according to Eqs. (15), and (25), it can

be shown

(17)   inductively that the sequence qðkÞ is bounded. Thus, qðkÞ is a Cauchy sequence and then we
can verify  that ðu; qÞ is the unique cluster point of the sequence uðkÞ ; qðkÞ . This ends the proof. □   ðkÞ of the In Algorithm 3.1, the approximate solution uðkÞ N ; qN boundary value problem (13)–(17) is computed using the MFS. As the MFS has been successfully applied to a large number of boundary value problems, many works on the theoretical analysis including convergent results have been developed [39–44]. These theoretical results indicate that the MFS achieves high accuracy such as exponential convergence under some conditions.   ðkÞ Lemma 3.2. Let uðkÞ be the MFS solution generated by N ;
qN Algorithm 3.1 and let uðkÞ ; qðkÞ be the analytical solution of the boundary value problem (13)–(17), then     ðkÞ ¼ uðkÞ ; qðkÞ lim uðkÞ N ; qN N-1

Finally, the convergence result of the current meshless MFS algorithm for the Signorini problem (1)–(5) follows straightforward from Lemmas 3.1 and 3.2.   ðkÞ be the MFS solution generated by Theorem 3.1. Let uðkÞ N ; qN Algorithm 3.1 and let ðu; qÞ be the analytical solution of the Signorini problem (1)–(5), then   ðkÞ lim ¼ ðu; qÞ uðkÞ N ; qN k-1;N-1

4. Numerical experiments Some numerical examples are presented in this section to demonstrate the efficiency and accuracy of the current meshless MFS algorithm and to confirm the convergence analysis. In this paper, all numerical codes have been implemented in MATLAB 7.0.1, and the computations have been tested on a desktop PC with an Intel 2.93 GHz CPU and 2 GB RAM.

The first example is a Signorini problem in an annular domain. The exact solution is [14] pffiffiffiffiffiffiffiffiffiffiffi 3 pffiffiffiffiffiffiffiffiffiffiffi uðxÞ ¼ Im A þ B sign x1 þ i A  B sign x2 ;

 x ¼ ðx1 ; x2 ÞT A Ω ¼ a o r o b qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  
2 2
2 2 x2  x2 2 and x1  x22 þ 4a14 r 4  a4 , B ¼ 14r2 2 ar 2 þ ar2 where A ¼ 2r12 qffiffiffiffiffiffiffiffiffiffiffiffiffiffi r ¼ x21 þ x22 . It can be verified that the solution can be simplified as pffiffiffiffiffiffiffiffiffiffiffi uðxÞ ¼ 2ðA þ 2BÞ A  B sign x2 And its derivatives can be expressed as

∂u ∂A ∂B pffiffiffiffiffiffiffiffiffiffiffi ¼2 þ2 A  B þ ðA þ 2BÞ ∂x1 ∂x1 ∂x1

pffiffiffiffiffiffiffiffiffiffiffi ∂A ∂B  =ð2 A  BÞ sign x2 ∂x1 ∂x1 and



∂u ∂A ∂B pffiffiffiffiffiffiffiffiffiffiffi ¼2 þ2 A  B þ ðA þ 2BÞ ∂x2 ∂x2 ∂x2

pffiffiffiffiffiffiffiffiffiffiffi ∂A ∂B  =ð2 A  BÞ sign x2 ∂x2 ∂x2

In this study, we use a ¼0.1 and b¼ 0.25 in computation. Besides, Dirichlet boundary conditions corresponding to the analytical solution are imposed on the outer circle Γ D ¼ fx A R2 : r ¼ bg, and Signorini boundary conditions are imposed on the inner circle as uðxÞ 4 ϕðxÞ; qðxÞ ¼ φðxÞ

 x A Γ S ¼ x A R2 : r ¼ a

or

uðxÞ ¼ ϕðxÞ;

qðxÞ 4 φðxÞ;

where ϕðxÞ ¼ minð0; uðxÞÞ and 8  6=a; x2 4 0 > <
 φðxÞ ¼  6 x21  x22 2 =a5 ; x2 r 0; x2 Z  jx1 j > : 0; x2 o  jx1 j Fig. 1 gives a comparison between the MFS solutions with the exact solutions for u and q on Γ S when 64 boundary nodes are used on each circle. Here, the results are shown as plots over the arc angle 2πθ. We can find from this figure that the numerical solutions are in excellent agreement with the exact solutions. For error estimation and convergence studies, the problem is solved with different numbers of boundary nodes. An error norm

Fig. 1. Results of (a) u and (b) q on the Signorini boundary Γ S .

52

H. Zheng, X. Li / Engineering Analysis with Boundary Elements 58 (2015) 48–57

is defined as ffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X


2 X

2 eðuÞ ¼ u xj uN xj = u xj xj A Γ S

xj A Γ S

where u and uN are exact and numerical solutions, respectively. Fig. 2 gives the log–log plot of error eðuÞ with respect to the boundary nodes number N. In this analysis, six kinds of nodal arrangement of 32, 64, 128, 256, 512 and 1024 boundary nodes are used. The errors of the Galerkin boundary element method (GBEM) [14], the boundary element-linear complementary method (BE-LCM) [15] and the boundary element-projection iterative method (BE-PIM) [16] are also given in this figure for comparison. We can observe that both error and convergence of the MFS are the best. In addition, the experimental convergence rate obtained using the GBEM, the BE-LCM, the BE-PIM and the MFS is 1.11, 1.01, 1.84 and 2.34, respectively. Tables 1 and 2 tabulate the corresponding number of iteration and CPU time, respectively. From Table 1 we can find that the number of iteration of the MFS is little influenced by the number of boundary nodes, and is much less than that of the GBEM and the BE-LCM. From Table 2 we can find that the CPU time required in the BE-LCM is about 30 times that in the MFS, and the CPU time required in the BE-PIM is about tenfold that in the MFS. Therefore, CPU time of the MFS is much less than that of the BE-LCM and the BE-PIM. In sum, the MFS outperforms the GBEM, the BE-LCM and

the BE-PIM in terms of accuracy, convergence, iterative number and CPU time for this experiment. Approximate values of u and its derivatives ∂u=∂x1 and ∂u=∂x2 at some inner points are also computed in this example. The MFS results for u along the inner circle r ¼ 0.2 are shown in Fig. 3 together with the analytical solutions. The results for ∂u=∂x1 and ∂u=∂x2 are shown in Fig. 4. These figures indicate that the numerical results at inner points match the analytical solutions very well for both u and its derivatives.

4.2. Shallow dam problem The problem considered here is the well-known shallow dam problem related to percolation in gently sloping beaches [1–4]. Fig. 5 depicts the schematic diagram. The mathematical model of

Fig. 3. Results of u on the inner circle r¼ 0.2.

Fig. 2. Convergence of the MFS and three BEMs for the Signorini problem in an annular domain.

Table 1 Number of iteration for the Signorini problem in an annular domain. N

32

64

128

256

GBEM BE-LCM BE-PIM MFS

24 56 4 6

28 91 4 6

36 141 5 6

48 203 6 7

512

1024 Fig. 4. Results of ∂u=∂x1 and ∂u=∂x2 on the inner circle r ¼0.2.

263 7 7

319 8 7

Table 2 CPU time (in seconds) for the Signorini problem in an annular domain. N

32

64

128

256

512

1024

BE-LCM BE-PIM MFS

0.047 0.078 0.016

0.125 0.141 0.031

0.453 0.438 0.047

2.407 1.782 0.203

16.953 8.141 0.656

125.907 40.188 3.875 Fig. 5. Schematic diagram for the shallow dam problem.

H. Zheng, X. Li / Engineering Analysis with Boundary Elements 58 (2015) 48–57

53

Fig. 6. Approximate solution of the shallow dam problem with (a) G1 ðx1 Þ, (b) G2 ðx1 Þ and (c) G3 ðx1 Þ.

Table 4

Table 3 Number of iteration for the shallow dam problem when N ¼ 64.

Number of iteration for the shallow dam problem when τ ¼ 10  6 .

τ

10  2

10  3

10  4

10  5

10  6

N

28

48

64

84

128

256

512

1024

The BPIM The MFS in Ref. [4] Our MFS

5 77 2

10 84 3

14 270 3

19 345 3

24 893 4

The BPIM The MFS in Ref. [4] Our MFS

25 72 4

25 119 4

24 893 4

24 946 5

25

25

25

24

7

5

5

5

this problem can be formulated as 8 ΔuðxÞ ¼ 0 in Ω ¼ f0r x1 r 1; 0r x2 r 1g > > > > > > < uðxÞ ¼ 0 on Γ D ¼ fx1 ¼ 0; 0 r x2 r 1g qðxÞ ¼ 0 on Γ N ¼ f0 r x1 r 1; x2 ¼ 0g > > > uðxÞ ¼ Gð1Þ on Γ D ¼ fx1 ¼ 1; 0r x2 r 1g > > > : uðxÞ o Gðx1 Þ; qðxÞ ¼ 0 or uðxÞ ¼ Gðx1 Þ; qðxÞ o 0;

Table 5 CPU time (in seconds) for the shallow dam problem when τ ¼ 10  6 .

on

Γ S ¼ f0 r x1 r 1; x2 ¼ 1g

where the given function Gðx1 Þ depicts the surface profile. The following three cases of surface profiles [1–4] are considered:

1 G1 ð x 1 Þ ¼ x1 ð1  x1 Þ  x1 2 G2 ð x 1 Þ ¼

N

28

48

64

84

128

256

512

1024

The BPIM Our MFS

0.078 0.016

0.125 0.016

0.141 0.031

0.157 0.063

0.250 0.094

0.672 0.204

2.469 0.750

10.734 3.625





1 5 5  x1 1  x 1 ð1  x 1 Þ  x 1 2 2 3

G3 ðx1 Þ ¼ sin ð12x1 Þ  2 Fig. 6 plots a comparison between our MFS solutions with the BPIM solutions [34] for u on Γ S when 32 boundary nodes are used on per side of the domain. Although no analytical solutions exist for this realistic physical problem, we can find that our MFS solutions are in good agreement with the BPIM solutions. In order to investigate the convergence behavior of our method, the problem is solved choosing different numbers of boundary nodes (N) and various values of the tolerance (τ). The numerical results for G1 ðx1 Þ are compared to those obtained using the MFS in Ref. [4] and the BPIM in Ref. [34]. When N ¼ 64, the number of iteration for various τ is tabulated in Table 3. When τ ¼ 10  6 , the number of iteration and CPU time for various N are tabulated in Tables 4 and 5, respectively. Tables 3 and 4 indicate that the number of iteration of our MFS is little influenced by the number of boundary nodes, and is much less than that of the MFS in Ref. [4] and the BPIM in Ref. [34]. From Table 5 we can find that the CPU time required in the BPIM is about threefold that in our MFS. Therefore, the convergence of our MFS is faster than that of the MFS in Ref. [4] and the BPIM in Ref. [34]. 4.3. 2D electropainting problem The well-known electropainting problem describes the coating process of a metal surface with paint [3,6]. This problem plays a key role in the application of corrosion-protection paints to vehicle

Fig. 7. Schematic diagram for the electropainting problem.

bodies in the motor manufacturing industry. The workpiece to be painted is immersed in a electrobath containing an electrolyte paint liquor with charged ions. A potential difference is established between the workpiece and the edge of the electrobath. Consequently, the paint particles become charged and are deposited on the surface of the workpiece. This realistic industrial problem can be modeled as a Signorini problem on a square domain ð  0:5; 0:5Þ  ð0; 1Þ for the electric potential u in the paint liquor. Owing to the symmetry of the problem, only the right-hand side of the square domain needs to be simulated as shown in Fig. 7. Then, the electric potential u satisfies

ΔuðxÞ ¼ 0; uðxÞ ¼ 1;

x A Ω ¼ f0 rx1 r 0:5; 0 r x2 r1g x A Γ D ¼ f0 r x1 r 0:5; x2 ¼ 0g

54

H. Zheng, X. Li / Engineering Analysis with Boundary Elements 58 (2015) 48–57

Fig. 8. Numerical results of the 2D electropainting problem for (a) ε ¼ 0:4, (b) ε ¼ 0:5, (c) ε ¼ 0:55 and (d) ε ¼ 0:7.

Fig. 9. Paint distributions of the 2D electropainting problem for (a) ε ¼ 0:4, (b) ε ¼ 0:5, (c) ε ¼ 0:55 and (d) ε ¼ 0:7.

H. Zheng, X. Li / Engineering Analysis with Boundary Elements 58 (2015) 48–57

55

Fig. 10. Paint distributions of the 3D electropainting problem for (a) ε ¼ 1, (b) ε ¼ 1:5, (c) ε ¼ 3 and (d) ε ¼ 6.

qðxÞ ¼ 0;

x A Γ N ¼ fx1 ¼ 0; 0 r x2 r 1g

uðxÞ 4 0;

qðxÞ ¼  ε

or

uðxÞ ¼ 0;

qðxÞ 4  ε;

x A Γ S ¼ Γ \ðΓ D [ Γ N Þ where the known constant ε is the critical cut-off current and is a property of the paint, and the Signorini boundary Γ S makes up the workpiece to be painted. This problem has a unique solution [6]. In this problem, the primary aim is to decide which parts of Γ S should be painted (u 40; q þ ε ¼ 0Þ and the corresponding thickness distribution. For various values of ε, Fig. 8 plots numerical results of the MFS and the BEM [16] for u on Γ S when 240 boundary nodes are used. Here, the results are shown as plots over the arc-length s, and s ¼0, s¼ 0.5, s ¼1.5 and s ¼2 correspond to the vertexes ð0; 0Þ, ð0:5; 0Þ, ð0:5; 1Þ and ð0; 1Þ, respectively. We can find from this figure that the results of the MFS are in excellent agreement with those of the BEM. Paint distributions over the domain for various ε are shown in Fig. 9. The results are obtained using 240 boundary nodes. Figs. 8(a) and 9(a) indicate that the workpiece is completely painted with the corner receiving the least amount of paint for a small ε. If ε is

increased, Figs. 8(b) and (c) and 9(b) and (c) indicate that the paint film near the corner becomes thinner and eventually unpainted. With a higher ε, Figs. 8(d) and 9(d) indicate that both the corner and the top boundary are unpainted. As a result, although no analytical solutions exist for this realistic industrial problem, the reliability of the MFS can be assessed based on the reasonable physical characteristics obtained. Additionally, we can find that the results of our MFS match those of the FEM [6], the MFS [4], and the BEM [3,15].

4.4. 3D electropainting problem To verify that the present MFS has good capability to tackle 3D Signorini problems, the last example embraces a 3D electropainting problem [5,7]. The mathematical model of this problem can be formulated as 8

 Δu ¼ 0 in Ω ¼  a rx1 r a;  b r x2 r b;  c r x3 rc > > >

 > > Γ ¼  a r x ra; x ¼ b;  c rx rc q ¼ 0 on > N 1 2 3 <

 u 4 0 q ¼  ε or u ¼ 0; q 4  ε on Γ S ¼  a rx1 ra; x2 ¼  b;  c r x3 rc > > > > u ¼ 1 on Γ D ¼ Γ \ðΓ N [ Γ S Þ > > :

56

H. Zheng, X. Li / Engineering Analysis with Boundary Elements 58 (2015) 48–57

where u is the electric potential and ε is the critical current necessary for paint deposition. In this study, a, b and c are set to be unity in computation. In this problem, Γ S can be considered as the car roof to be painted. On Γ S , the painted subregions are unknown in advance. Therefore, we need to decide the painted subregions and the corresponding paint thickness. Fig. 10 plots the paint distribution over Γ S , obtained for different values of ε. In this analysis, we used 150 nodes on the whole boundary. We can observe from Fig. 10 that the paint thickness is symmetric with respect to the axes of coordinates. Besides, the minimum paint thickness occurs near the center of the car roof, and the paint thickness increases as one moves from the center. For the smaller values of ε, the car roof is completely painted with the center receiving the least amount of paint, as shown in Fig. 10(a). If the value of ε is increased, the car roof center becomes thinner (Fig. 10(b)) and eventually unpainted (Fig. 10(c)). At an even higher value of ε, the unpainted area becomes larger (Fig. 10(d)). Therefore, although there are no analytical solutions available to verify the current MFS for such a complex industrial problem, the reliability can be judged from the reasonable physical characteristics.

5. Conclusion In this study, we present an application of the MFS for boundary-only analysis of 2D and 3D Signorini problems. In this numerical method, the nonlinear Signorini boundary inequality conditions are incorporated naturally into a linear iterative scheme and thus can be implemented easily. Then, the nonlinear Signorini problem is transformed into a sequence of linear elliptic boundary value problems and solved numerically by the MFS. The MFS is ideally suitable for Signorini problems since, in such problems, it is the Signorini boundary which is of prime interest. Convergence of the MFS for Signorini problems is proven theoretically. The theoretical result indicates that the MFS generates a sequence of strongly convergent approximations. Some numerical examples are given to demonstrate the accuracy, convergence and efficiency of the MFS for 2D and 3D Signorini problems. The numerical results are accurate and are in excellent agreement with analytical solutions and previously reported numerical results. The errors and CPU times are compared with those obtained by the BEMs. These comparisons indicate that the current MFS has higher computational precision and convergence rate, but has lower CPU time. Therefore, compared with the BEMs, the computational efficiency is improved markedly by the current MFS.

Acknowledgments This work was supported by the National Natural Science Foundation of China (No. 11471063), the Natural Science Foundation Project of CQ CSTC (No. cstc2014jcyjA00005) and the Educational Commission Foundation of Chongqing of China (No. KJ130604). References [1] Aitchison JM, Elliott CM, Ockendon JR. Percolation in gently sloping beaches. IMA J Appl Math 1983;30:269–87. [2] Karageorghis A. Numerical solution of a shallow dam problem by a boundary element method. Comput Methods Appl Mech Eng 1987;61:265–76.

[3] Aitchison JM, Poole MW. A numerical algorithm for the solution of Signorini problems. J Comput Appl Math 1998;94:55–67. [4] Poullikkas A, Karageorghis A, Georgiou G. The method of fundamental solutions for Signorini problems. IMA J Appl Math 1998;18:273–85. [5] Poullikkas A, Karageorghis A, Georgiou G. The numerical solution of threedimensional Signorini problems with the method of fundamental solutions. Eng Anal Bound Elem 2001;25:221–7. [6] Aitchison JM, Lacey AA, Shillor M. A model for an electropaint process. IMA J Appl Math 1984;33:17–31. [7] Poole PW, Aitchison JM. Numerical model of an electropaint process with applications to the automotive industry. IMA J Math Appl Bus Ind 1997;8:347–60. [8] Coorevits P, Hild P, Lhalouani K, Sassi T. Mixed finite element methods for unilateral problems: convergence analysis and numerical studies. Math Comput 2002;71:1–25. [9] 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. [10] Howison D, Morgan JD, Ockendon JR. A class of codimension two free boundary problems. SIAM Rev. 1997;39:221–53. [11] Chen JT, Hsiao CC, Lee YT. Study of free-surface seepage problem using hypersingular equations. Commun Numer Methods Eng 2007;23:755–69. [12] Noor MA. Projection iterative methods for extended general variational inequalities. J Appl Math Comput 2010;32:83–95. [13] Han JY, Xiu NH, Qi HD. Nonlinear complementarity theory and algorithm. Shanghai: Science and Technology Press; 2006. [14] Spann W. On the boundary element method for the Signorini problem of the Laplacian. Numer Math 1993;65:337–56. [15] Zhang S, Zhu J. The boundary element-linear complementary method for the Signorini problem. Eng Anal Bound Elem 2012;36:112–7. [16] Zhang S, Zhu J. A projection iterative algorithm boundary element method for the Signorini problem. Eng Anal Bound Elem 2013;37:176–81. [17] Zhang S. A projection iterative algorithm for the Signorini problem using the boundary element method. Eng Anal Bound Elem 2015;50:313–9. [18] Cheng Y, Ji X, He P. Infinite similar boundary element method for dynamic fracture mechanics. Acta Mech Sin 2004;36(1):43–8. [19] Mukherjee YX, Mukherjee S. The boundary node method for potential problems. Int J Numer Methods Eng 1997;40:797–815. [20] Liu GR, Gu YT. Boundary meshfree methods based on the boundary point interpolation methods. Eng Anal Bound Elem 2004;28:475–87. [21] Zhang J, Yao Z, Li H. A hybrid boundary node method. Int J Numer Methods Eng 2002;53:751–63. [22] Miao Y, He T, Yang Q, Zheng J. Multi-domain hybrid boundary node method for evaluating top-down crack in asphalt pavements. Eng Anal Bound Elem 2010;34:755–60. [23] Wang Q, Miao Y, Zhu H. A fast multipole hybrid boundary node method for composite materials. Comput Mech 2013;51(6):885–97. [24] Miao Y, Wang Q, Zhu H. A new fast multipole hybrid boundary node method for thermal analysis of 3D composites. Comput Mech 2014;53:77–90. [25] Miao Y, Sun T, Zhu H, Wang Q. A new model for the analysis of reinforced concrete members with a coupled HdBNM/FEM. Appl Math Model 2014;38:5582–91. [26] Li X, Zhu JA. Galerkin boundary node method and its convergence analysis. J Comput Appl Math 2009;230:314–28. [27] Li X. A meshless interpolating Galerkin boundary node method for Stokes flows. Eng Anal Bound Elem 2015;51:112–22. [28] Li X. Implementation of boundary conditions in BIEs-based meshless methods: a dual boundary node method. Eng Anal Bound Elem 2014;41:139–51. [29] Li X, Li S. Meshless boundary node methods for Stokes problems. Appl Math Model 2015;39:1769–83. [30] Cheng Y, Peng M. Boundary element-free method for elastodynamics. Sci Chin Ser G Phys Mech Astron 2005;48(6):641–57. [31] Peng M, Cheng Y. A boundary element-free method (BEFM) for twodimensional potential problems. Eng Anal Bound Elem 2009;33:77–82. [32] Li X. An interpolating boundary element-free method for three-dimensional potential problems. Appl Math Model 2014 http://dx.doi.org/10.1016/j.apm. 2014.10.071. [33] Wang Y, Li X. A meshless algorithm with moving least squares approximations for elliptic Signorini problems. Chin Phys B 2014;23(9):090202. [34] Ren Y, Li X. A meshfree method for Signorini problems using boundary integral equations. Math Probl Eng 2014:490127. [35] Li F, Li X. The interpolating boundary element-free method for unilateral problems arising in variational inequalities. Math Probl Eng 2014:518727. [36] Fairweather G, Karageorghis A. The method of fundamental solutions for elliptic boundary value problems. Adv Comput Math 1998;9:69–95. [37] Chen CS, Karageorghis A, Smyrlis YS. The method of fundamental solutions—a meshless method. Atlanta: Dynamic Publishers; 2008. [38] Karageorghis A. Modified methods of fundamental solutions for harmonic and biharmonic problems with boundary singularities. Numer Methods Part Differ Equ 1992;8:1–19. [39] Katsurada M, Okamoto H. The collocation points of the fundamental solution method for the potential problem. Comput Math Appl 1996;31:123–37.

H. Zheng, X. Li / Engineering Analysis with Boundary Elements 58 (2015) 48–57

[40] Mitic P, Rashed YF. Convergence and stability of the method of meshless fundamental solutions using an array of randomly distributed sources. Eng Anal Bound Elem 2004;28:143–53. [41] Li X. Convergence of the method of fundamental solutions for Poisson's equation on the unit sphere. Adv Comput Math 2008;28:169–282. [42] Li X, Zhu J. The method of fundamental solutions for nonlinear elliptic problems. Eng Anal Bound Elem 2009;33:322–9.

57

[43] Li X, Li S. A meshless method for nonhomogeneous polyharmonic problems using method of fundamental solution coupled with quasi-interpolation technique. Appl Math Model 2011;35:3698–709. [44] Li ZC. The method of fundamental solutions for annular shaped domains. J Comput Appl Math 2009;228:355–72.