Numerical solution of two-dimensional mixed problems with variable coefficients by the boundary-domain integral and integro-differential equation methods

Engineering Analysis with Boundary Elements 35 (2011) 1279–1287

Contents lists available at ScienceDirect

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

Numerical solution of two-dimensional mixed problems with variable coefficients by the boundary-domain integral and integro-differential equation methods M.A. AL-Jawary , L.C. Wrobel School of Engineering and Design, Brunel University, Uxbridge UB8 3PH, UK

a r t i c l e i n f o

abstract

Article history: Received 25 February 2011 Accepted 16 June 2011

This paper presents a formulation of the boundary-domain integral equation (BDIE) and the boundarydomain integro-differential equation (BDIDE) methods for the numerical solution of two-dimensional mixed boundary-value problems (BVP) for a second-order linear elliptic partial differential equation (PDE) with variable coefficients. The methods use a specially constructed parametrix (Levi function) to reduce the BVP to a BDIE or BDIDE. The numerical formulation of the BDIDE employs an approximation for the boundary fluxes in terms of the potential function within the domain cells; therefore, the solution is fully described in terms of the variation of the potential function along the boundary and domain. Linear basis functions localised on triangular elements and standard quadrature rules are used for the calculation of boundary integrals. For the domain integrals, we have implemented Gaussian quadrature rules for two dimensions with Duffy transformation, by mapping the triangles into squares and eliminating the weak singularity in the process. Numerical examples are presented for several simple problems with square and circular domains, for which exact solutions are available. It is shown that the present method produces accurate results even with coarse meshes. The numerical results also show that high rates of convergence are obtained with mesh refinement. & 2011 Elsevier Ltd. All rights reserved.

Keywords: Mixed boundary-value problem Boundary-domain integral equation method Boundary-domain integro-differential equation method Variable coefficients Parametrix

1. Introduction The boundary element method (BEM) has become a powerful method for the numerical solution of BVPs, due to its ability (at least for problems with constant coefficients) of reducing a BVP for a linear PDE defined in a domain to an integral equation defined on the domain boundary, leading to a simplified discretisation process with boundary elements only. The main requirement for the reduction of the PDE to a boundary integral equation (BIE) is that a fundamental solution to the PDE must be available. Such fundamental solutions are well known for many PDEs with constant coefficients, see [1–4], but are not generally available when the coefficients of the original PDE are variable, except for some special cases [5,6]. The coefficients in the mathematical model of a physical problem typically correspond to the material parameters of the problem. In heterogeneous media the material parameters may vary with position and/or time. The governing equation of a physical problem in heterogeneous media is therefore likely to involve variable

 Corresponding author.

E-mail addresses: [email protected] (M.A. AL-Jawary), [email protected] (L.C. Wrobel). 0955-7997/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2011.06.004

coefficients. For this reason, there is a demand on the development of accurate and efficient numerical methods able to deal with the spatial variations of material coefficients [7,8]. Several BEM techniques have been proposed to treat problems with variable coefficients, e.g. by cell discretisation [7] or the dual reciprocity method (DRM) [8]. In this work, a parametrix (Levi function) is adopted, which is usually available. This allows a reduction of the mathematical problem to a boundary-domain integral or integro-differential equation (BDIE or BDIDE) [9,10]. A BDIE and a BDIDE formulations to solve problems with variable coefficients are presented in [9] using specially constructed localised parametrices to reduce a BVP with variable coefficients to a localised boundary-domain integral or integro-differential equation (LBDIE or LBDIDE). The use of specially constructed localised parametrices leads to sparsely populated systems of linear algebraic equations. An implementation of the LBDIE method for the numerical solution of a second-order linear elliptic PDE with variable coefficients is presented in [10], although the formulation is restricted to Neumann boundary-value problems. Boundary-domain integral equation methods have also been developed by Skerget et al. [11–13] for the solution of non-linear fluid mechanics problems described by the Navier–Stokes equations. Skerget’s formulation treats all the non-linear terms as body forces, which are included in the boundary integral equations as

1280

M.A. AL-Jawary, L.C. Wrobel / Engineering Analysis with Boundary Elements 35 (2011) 1279–1287

a domain integral evaluated by discretising the body into cells. Another related formulation developed by Popov and Power [14], named the dual reciprocity-multi-domain (DRM-MD) approach, combines the DRM with domain decomposition, leading to substantial improvements in the accuracy and convergence of the DRM formulation 5 for complex problems. The DRM-MD formulation has been applied to the solution of the Navier–Stokes equation [15] and to flow and solute transport in fractured porous media [16]. A further boundary-domain integral equation technique is the Analog Equation Method (AEM) of Katsikadelis [17], which has been applied to the solution of several elasticity problems, mostly related to plate bending. Katsikadelis and Nerantzaki [18] extended the AEM to a boundary-only method which decomposes the solution into a homogeneous part and a particular solution of the non-homogeneous one, and then obtained the particular solution via a radial basis function expansion of the domain term. In the present paper, the BDIE and BDIDE formulations proposed by Mikhailov [9,10] are extended to the treatment of mixed BVPs. The numerical algorithms developed here do not use the concept of localisation as in [9,10], but rather use global mesh-based discretisations. The paper also discusses different techniques to deal with the discontinuity of the flux at corners, by testing different positions of the collocation points. Numerical solutions of several test examples are included to validate the methods.

2. Reduction of the mixed BVP to a BDIE/BDIDE Let us consider the following stationary heat transfer BVP in an isotropic inhomogeneous medium for a two-dimensional body O, with prescribed temperature uðxÞ on part @D O of the boundary @O and prescribed heat flux tðxÞ on the remaining @N O part of @O, i.e. we consider the second-order linear elliptic PDE [9,10]:   2 X @ @uðxÞ ðLuÞðxÞ :¼ aðxÞ ð1Þ ¼ f ðxÞ, x A O @xi @xi i¼1 with the mixed boundary conditions uðxÞ ¼ uðxÞ, TuðxÞ ¼ tðxÞ,

x A @D O

ð2Þ

x A @N O

ð3Þ

where O is a bounded domain, u(x) the temperature, a(x) a known variable thermoconductivity coefficient, f(x) a known heat source, T a surface flux operator, ½TuðxÞ :¼ aðxÞð@u=@nÞðxÞ, n(x) the external normal vector to the boundary @O, and uðxÞ and tðxÞ are known functions. The BVP (1–3) appears particularly when modelling stationary heat transfer, elastostatics, electrostatics, and diffusion problems for functionally graded materials, as well as in flow in porous media. The Green formula for the differential operator L has the form Z Z ½uLvvLu dO ¼ ½uTvvTu dG ð4Þ O

@O

where u and v are arbitrary functions. Let L be a linear operator and F(x,y) be its fundamental solution, i.e. Lx Fðx,yÞ ¼ dðxyÞ where d is the Dirac delta function. Then one could take vðxÞ ¼ Fðx,yÞ, identify u(x) with a solution of Eq. (1), and thus arrive at the third Green identity Z Z cðyÞuðyÞ ½uðxÞTx Fðx,yÞFðx,yÞTuðxÞ dGðxÞ ¼ Fðx,yÞf ðxÞ dOðxÞ @O

O

ð5Þ

where 8 1 > > > < 0 cðyÞ ¼ > aðyÞ > > : 2p

if y A O if y 2 =O

ð6Þ

if y A @O and O  R2

where aðyÞ is the interior angle at a point y of the boundary @O, particularly, c(y)¼1/2 if y is a smooth point on the boundary. Substituting the boundary condition in the Green identity equation (5) and applying it for yA @O, we arrive at a direct boundary integral equation [1–4]. For partial differential operators with variable coefficients, like L in Eq. (1), a fundamental solution is generally not available in explicit form. However, a parametrix is often available, which is a function Pðx,yÞ satisfying the equation [9,10] Lx Pðx,yÞ ¼ dðxyÞ þRðx,yÞ

ð7Þ

The fundamental solution of the operator with frozen coefficients aðxÞ ¼ aðyÞ corresponding to the operator L defined in (1) can be used as a parametrix, in the two-dimensional case [9,10], 1 ln jxyj ð8Þ 2paðyÞ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where jxyj ¼ ðx1 y1 Þ2 þ ðx2 y2 Þ2 . Substituting Eq. (8) in Eq. (7), we obtain    2 X @ @ 1 ln jxyj ¼ dðxyÞ þ Rðx,yÞ aðxÞ @xi @xi 2paðyÞ i¼1

Pðx,yÞ ¼

By applying the product rule for differentiation, we get "    # 1 @aðxÞ @ 1 aðxÞ @2 1    ln jxyj þ lnjxyj aðyÞ @xi @xi 2p aðyÞ @x2i 2p i¼1 2 X

¼ dðxyÞ þRðx,yÞ P2 2 2 Now, since aðxÞ ¼ aðyÞ and i ¼ 1 ð@ =@xi Þ½ð1=2pÞln jxyj ¼ dðxyÞ, we have   2   aðxÞ X @2 1  dðxyÞ ¼ ln jxyj 2 aðyÞ i ¼ 1 @xi 2p and Rðx,yÞ ¼

  2 X 1 @aðxÞ @ 1   ln jxyj aðyÞ @xi @xi 2p i¼1

  @ 1 1 @r x y lnjxyj ¼ ¼ i 2i , @xi 2p 2pr @xi 2pr

r ¼ jxyj

Hence, the remainder Rðx,yÞ will then be Rðx,yÞ ¼

2 X

xi yi @aðxÞ , 2 @x 2 p aðyÞjxyj i i¼1

x,y A R2

ð9Þ

which has only a weak singularity at x ¼y. Substituting Pðx,yÞ for v(x) in Eq. (4) and taking u(x) as a solution to Eq. (1), we obtain the integral equation Z ½uðxÞTx Pðx,yÞPðx,yÞTuðxÞ dGðxÞ cðyÞuðyÞ Z @O Z þ Rðx,yÞuðxÞ dOðxÞ ¼ Pðx,yÞf ðxÞ dOðxÞ ð10Þ O

O

Identity (10) can be used for formulating either a BDIE or a BDIDE, with respect to u and its derivatives. Let us consider the two forms below. 2.1. Boundary-domain integral equation (BDIE) Substituting the boundary conditions (2) and (3) into (10), introducing a new variable t(x)¼Tu(x) for the unknown flux on

M.A. AL-Jawary, L.C. Wrobel / Engineering Analysis with Boundary Elements 35 (2011) 1279–1287

@D O and using Eq. (10) at y A O[@O reduce the BVP (1)–(3) to the following boundary-domain integral equation (BDIE) for u(x) at x A O [ @N O and t(x) at x A @D O, Z Z uðxÞTx Pðx,yÞ dGðxÞ þ Pðx,yÞtðxÞ dGðxÞ c0 ðyÞuðyÞ @N O @D O Z ð11Þ þ Rðx,yÞuðxÞ dOðxÞ ¼ C0 ðyÞ, y A O [ @O O

CðyÞ :¼

Z

ð12Þ

Pðx,yÞtðxÞ dGðxÞ þ

@N O

Z

¼ C ðxi Þ xi A O [ @O,

Pðx,yÞf ðxÞ dOðxÞ

ð13Þ

O

where c0 ðyÞ is ( 0 if y A @D O c0 ðyÞ ¼ cðyÞ if y A O [ @N O

xj A @D O

Kij uðxj Þ,

0

i ¼ 1, . . . ,J, no sum in i

where C ðx Þ is calculated from Eq. (12), and Z Z Cðxi Þ ¼ uðxÞTx Pðx,xi Þ dGðxÞ Pðx,xi ÞtðxÞ dGðxÞ Z

ð14Þ

2.2. Boundary-domain integro-differential equation (BDIDE) Using another approach, we can substitute the boundary conditions (2) and (3) into (10) but leave T as a differential flux operator acting on u on the Dirichlet boundary @D O and use the following boundary-domain integro-differential equation (BDIDE) at y A O [ @N O, Z Z uðxÞTx Pðx,yÞ dGðxÞ þ Pðx,yÞTuðxÞ dGðxÞ cðyÞuðyÞ @ O @D O Z N ð15Þ þ Rðx,yÞuðxÞ dOðxÞ ¼ CðyÞ, yA O [ @N O O

where CðyÞ is given by Eq. (13). As we will see below, this approach can lead, after discretisation, to a system with a reduced number of linear algebraic equations. Eqs. (11) and (15) will lead, after discretisation, to fully populated systems of linear algebraic equations.

3. Discretisation of the BDIE/BDIDE 3.1. Discretisation of the BDIE Let us discretise the domain O into a mesh of triangular elements Tk, k ¼ 1,2, . . . ,N, Th \ Tm ¼ |, ha m. Let J be the total number of nodes xi, i ¼ 1, . . . ,J, at the vertices of triangles, from which there are JD nodes on @D O. To obtain a system of linear algebraic equations from the BDIE (11), by the collocation method, we collocate at the nodes xi, i ¼ 1, . . . ,J and substitute an interpolation of u(x) of the form ( X fkj ðxÞ if x,xj A T k uðxÞ  uðxj ÞFj ðxÞ, Fj ðxÞ ¼ ð16Þ 0 otherwise o 3x j

where o j is the support of Fj ðxÞ, which consists of all triangular elements that have xj as a vertex; fkj ðxÞ are the shape functions localised on an element Tk, and associated with the node xj. For the triangular elements, fkj ðxÞ can be chosen as piecewise linear functions. We can also use an interpolation of tðxÞ ¼ ðTuÞðxj Þ along boundary nodes belonging to o ðxj Þ \ @D O X tðxÞ ¼ tðxj Þvj ðxÞ, x A o ðxj Þ \ @D O ð17Þ

ð18Þ

i

o ðxi Þ\@D O

As the last term on the left-hand side of Eq. (11) includes the unknown values of u over the whole domain O, this BDIE does not lead to a BIE as in the case when the parametrix is a fundamental solution.

xj A o ðxj Þ\@D O

X

0

xj A @D O

uðxÞTx Pðx,yÞ dGðxÞ @ O

ZD 

Here, vj ðxÞ are boundary shape functions, taken now as constant. Therefore, vj(x) will be equal to 1 at xj A o ðxj Þ \ @D O and vj ðxÞ ¼ 0 if xj 2 = o ðxj Þ \ @D O. Substituting the interpolations (16) and (17) in BDIE (11) and applying the collocation method, we arrive at the following system of J linear algebraic equations for J unknowns uðxj Þ, xj A O [ @N O and tðxj Þ ¼ ðTuÞðxj Þ, xj A @D O, X X Kij uðxj Þ þ Q 0ij tðxj Þ c0 ðxi Þuðxi Þ þ xj A O[@N O

C0 ðyÞ :¼ ½co ðyÞcðyÞuðyÞ þ CðyÞ

1281

o ðxi Þ\@N O

f ðxÞPðx,xi Þ dOðxÞ

þ

ð19Þ

oðxi Þ\O

Z

Kij ¼

oðxi Þ\O

Q 0ij ¼

Fj ðxÞRðx,xi Þ dOðxÞ

Z o ðxi Þ\@D O

Z o ðxi Þ\@N O

Fj ðxÞTx Pðx,xi Þ dGðxÞ

Pðx,xi Þvj ðxÞ dGðxÞ

ð20Þ

ð21Þ

3.2. Discretisation of the BDIDE To obtain a system of linear algebraic equations from the BDIDE (15) by the collocation method, we collocate at the nodes xi , i ¼ 1, . . . J, arriving at a system of JJD algebraic equations for JJD unknowns uðxj Þ, xj A O [ @N O. Substituting interpolation formulae (16) into the BDIDE (15) leads to the following system: X X K 0ij uðxj Þ ¼ Cðxi Þ K 0ij uðxj Þ, cðxi Þuðxi Þ þ xj A O[@N O

xj A @D O

xi A O [ @N O, no sum in i where Z K 0ij ¼

oðxi Þ\O

Fj ðxÞRðx,xi Þ dOðxÞ þ

Z 

o ðxi Þ\@N O

ð22Þ Z o ðxi Þ\@D O

Pðx,xi ÞT Fj ðxÞ dGðxÞ

Fj ðxÞTx Pðx,xi Þ dGðxÞ

ð23Þ

and Cðxi Þ is given in Eq. (19).

4. Assembling the system matrix and right-hand side for BDIE/BDIDE Let us start with a Laplace equation with a mesh of eight boundary elements, nine nodes and eight triangular cells, as shown in Fig. 1. For the BDIE method, the system of algebraic equations resulting from Eq. (18) has two unknown variables t and u, i.e. t in Dirichlet boundaries and u in Neumann boundaries, in addition to interior nodes. In this paper, we present two implementations using mixed boundary elements with linear variation of u and constant t, to avoid the discontinuity of t at corner points. In the first case, the collocation nodes for calculating t in Dirichlet boundaries are taken at the mid point of the boundary elements, while in the second case the collocation points are at the end nodes.

1282

M.A. AL-Jawary, L.C. Wrobel / Engineering Analysis with Boundary Elements 35 (2011) 1279–1287

u

u

6

3

9 e5

e7 e6

e8

5 t

2

u

8 t

2 e3

e1

e4

e2 1

7 u 4

u

u

Fig. 1. Simple mesh.

Therefore, for the first case, the system Ax ¼b is given by 2 3 2 3 A11 A12 A13 t A @D O 6 6 A21 D þ A22 7 A23 5 , x ¼ 4 u A @N O 7 A¼4 5 , A31 A32 I þ A33 77 u A O 71 2 6 b¼4

cvec þb11þ b12 þ b13 b21 þ b22 þ b23 b31 þ b32 þ b33

3 7 5 71

where ½A1144 , ½A2124 , ½A3114 are the integrals in Eq. (21) with the collocation nodes belonging to @D O, @N O and interior nodes in O and integration nodes xj belonging to @D O, respectively. Also, ½A1242 , ½A2222 , ½A3212 are the second integral in Eq. (20) with the collocation nodes xi belonging to @D O, @N O and interior nodes in O and integration nodes xj belonging to @N O, respectively. Moreover, ½A1341 , ½A2321 , ½A3311 are the second integral in Eq. (20) with the collocation nodes xi belonging to @D O, @N O and interior nodes in O and integration nodes xj belonging to O, respectively. In a general mesh the dimensions of the matrix are ½AJ2J2 , and right-hand side ½bJ21 , where J is the total number of nodes. 0 The matrix D in this simple case is given by D ¼ ½0:5 0 0:5 22 , and the matrix I is the identity matrix, in this case just equal to 1. The right-hand side can be assembled in the same way as matrix A, where ½b1141 , ½b2121 , ½b3111 are the first integral in Eq. (19) with the collocation nodes belonging to @D O, @N O and interior nodes in O, respectively. Also, ½b1241 , ½b2221 , ½b3211 are the second integral in Eq. (19) with the collocation nodes belonging to @D O, @N O and interior nodes in O, respectively. P Moreover, ½b1341 , ½b2321 and ½b3311 are equal to  xj A @D O Kij uðxj Þ (with only the second integral in Eq. (20), since R¼0 for the Laplace equation, therefore, the first integral disappears), with the collocation nodes belonging to @D O, @N O and interior nodes in O, respectively. In addition, ½cvec41 is a vector equal to cðxi Þuðxi Þ, with the collocation nodes belonging to @D O and the values of c given in Eq. (6). In the second case, where the collocation nodes for calculating t in Dirichlet boundaries are taken at the end points of the elements, the system is given by: Ax ¼b, where ½Amn and m Zn. This system can be solved in the least squares sense by solving the following system [19], AT Ax ¼ AT b, where for our simple mesh the system Ax ¼b is given by 2 3 2 3 A11 A12 A13 t A @D O 6 A21 D þ A22 7 6 A23 5 , x ¼ 4 u A @N O 7 A¼4 5 , A31

A32

I þ A33

97

uAO

71

cvec þ b11 þ b12 þ b13

6 b¼4

b21 þb22 þb23 b31 þb32 þb33

3 7 5 91

where ½A1164 , ½A2124 , ½A3114 , ½A1262 , ½A2222 , ½A3212 , ½A1361 , ½A2321 , ½A3311 and the matrices D, I are the same as before. By applying the least squares technique, the final system will be Cx¼d: ½C77 ¼ ½AT 79 ½A97 , and ½d71 ¼ ½AT 79 ½b91 . Also, the righthand side b can be calculated like in the previous case but with different dimensions for the sub-vectors, i.e., ½b1161 , ½b1261 , ½b1361 , ½cvec61 and the other sub-vectors are the same as before. For the BDIDE method, the system of algebraic equations in Eq. (22) has only one unknown variable u, i.e. u in Neumann parts in addition to interior nodes. In this case, the assembling of matrix A and vector b is much easier than in the BDIE, i.e. by just adding the sub-matrices or sub-vectors which have the same dimension J  J for matrix Afull or vector bfull. The matrix Afull and vector bfull only have coefficients on positions xi, xj A O [ @N O and zero elsewhere. So, for the simple mesh, we can construct matrix A and vector b from Afull, bfull, respectively, and the system Ax¼b is given by: A ¼ ½ A1 þ A2 þA3 kk , x ¼ ½uuAA@NOO k1 , b ¼ ½ b1 þ b2 þb3 þ b4 k1 , where k ¼ Jr, J is the total number of nodes and r is the number of nodes on Dirichlet boundaries. The matrix A1 in this simple case is 2 3 0:5 0 0 6 7 A1 ¼ 4 0 1 0 5 0

0

0:5

33

the matrix ½A233 is the second integral in Eq. (23) and ½A333 is the third integral in Eq. (23). Finally, the right-hand side b can be calculated by adding the sub-vectors b1, b2, b3, b4, where b1, b2 are the first and the second integrals in Eq. (19), respectively. The vector b3 and b4 are defined as X b3 ¼  K1ij uðxj Þ xj A @D O

b4 ¼ 

X

K2ij uðxj Þ

xj A @D O

where Z K1 ¼ o ðxi Þ\@D O

Pðx,xi ÞT Fj ðxÞ dGðxÞ

Z K2 ¼  o ðxi Þ\@N O

Fj ðxÞTx Pðx,xi Þ dGðxÞ

Remark. In order to assemble the system of algebraic equations for the Poisson equation, where the third integral in Eq. (19) appears, one can follow the same steps as before for both BDIE and BDIDE in addition to the domain integral on the right-hand side b. However, to solve Eqs. (1)–(3), which have been re-formulated to BDIE or BDIDE in Eqs. (18)–(22), an extra domain integral appears (first integral in Eqs. (20) and (23)). Therefore, extra sub-matrices will be added to the matrix A for both BDIE and BDIDE. 5. Numerical results In this section, we shall examine some test examples to assess the performance of the BDIE/BDIDE formulations. To verify the convergence of the methods, we applied the methods to some test problems on square and circular domains, for which an exact analytical solution, uexact, is available. Also, the relative error was

M.A. AL-Jawary, L.C. Wrobel / Engineering Analysis with Boundary Elements 35 (2011) 1279–1287

tðxÞ ¼ 2ðx21 þ x22 Þðx1 n1 ðxÞ þx2 n2 ðxÞÞ

calculated as j

rðJÞ ¼

j

max1 r j r J juapprox ðx Þuexact ðx Þj max1 r j r J juexact ðxj Þj

where uapprox is the numerical solution. The following numerical tests were performed. 5.1. Laplace’s equation with mixed boundary conditions The starting point for testing the BDIE/BDIDE formulations is to consider Laplace’s equation with mixed boundary conditions. In this case there is no domain integral, i.e. f¼0 and a(x)¼1 in our original BVP (1)–(3), @2 u @2 u r u :¼ Du ¼ 2 þ 2 ¼ 0 @x1 @x2

1283

for x1 ¼ 1 or x1 ¼ 2, 1 r x2 r 2

The exact solution for this problem is uexact ðxÞ ¼ x21 þ x22 , x A O . 5.5. Test 3 [10] Circular domain O ¼ fðx1 ,x2 Þ : ðx1 1:5Þ2 þ ðx2 1:5Þ2 r0:25g, aðxÞ ¼ x21 þ x22 , f ðxÞ ¼ 8ðx21 þx22 Þ for x A O , with boundary conditions: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi uðxÞ ¼ x21 þ x22 for x2 ¼ 0:25ðx1 1:5Þ2 þ 1:5 tðxÞ ¼ 2ðx21 þx22 Þðx1 n1 ðxÞ þ x2 n2 ðxÞÞ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi for x2 ¼  0:25ðx1 1:5Þ2 þ 1:5

The exact solution for this problem is uexact ðxÞ ¼ x21 þ x22 , x A O .

2

5.6. Test 4

This simple test involves a square domain O ¼ fðx1 ,x2 Þ : 2 r x1 ,x2 r 3g, aðxÞ ¼ 1, f ðxÞ ¼ 0 for x A O , with boundary conditions: u ¼ x1 þ 2

for x2 ¼ 2, 2 r x1 r 3

u ¼ x1 þ 3

for x2 ¼ 3, 2 r x1 r 3

tðxÞ ¼ n1 ðxÞ þ n2 ðxÞ

for x1 ¼ 2 or x1 ¼ 3, 2 rx2 r3

The exact solution for this problem is uexact ðxÞ ¼ x1 þx2 , x A O . 5.2. Poisson’s equation with mixed boundary conditions The next test considers Poisson’s equation, in which case there is a domain integral coming from f a 0; we still consider a(x)¼1 in the original BVP (1)–(3), and assume a square domain O ¼ fðx1 ,x2 Þ : 2 rx1 ,x2 r 3g, a(x)¼1, f(x) ¼4 for x A O , with boundary conditions,

Circular domain O ¼ fðx1 ,x2 Þ : x21 þx22 r1g, aðxÞ ¼ x1 þx2 þ4, f ðxÞ ¼ 2 for x A O , with boundary conditions: qffiffiffiffiffiffiffiffiffiffiffiffi uðxÞ ¼ x1 þ x2 for x2 ¼ 1x21 tðxÞ ¼ ðx1 þ x2 þ 4Þðn1 ðxÞ þ n2 ðxÞÞ

The exact solution for this problem is uexact ðxÞ ¼ x1 þ x2 , x A O . 5.7. Test 5 Square domain O ¼ fðx1 ,x2 Þ : 2 rx1 ,x2 r 3g, aðxÞ ¼ expðx1 þ x2 Þ, f ðxÞ ¼ 2expðx1 þ x2 Þ for x A O , with boundary conditions: uðxÞ ¼ x1 þ 2

for x2 ¼ 2, 2 rx1 r3 for x2 ¼ 3, 2 rx1 r3

uðxÞ ¼ 4 þx21

for x2 ¼ 2, 2 rx1 r 3

uðxÞ ¼ x1 þ 3

uðxÞ ¼ 9 þx21

for x2 ¼ 3, 2 rx1 r 3

tðxÞ ¼ expðx1 þx2 Þðn1 ðxÞ þ n2 ðxÞÞ

tðxÞ ¼ 2ðx1 n1 ðxÞ þx2 n2 ðxÞÞ

for x1 ¼ 2 or x1 ¼ 3, 2 r x2 r 3

The exact solution for this problem is uexact ðxÞ ¼ x21 þx22 , x A O . The next series of tests with variable coefficients involve simple square or circular geometries with increasing degree of complexity of the variation of both the material parameter coefficients and the body force term f. The exact solutions of the problems range from linear to cubic, and will be used to verify the convergence of the numerical solutions.

Square domain O ¼ fðx1 ,x2 Þ : 2 r x1 ,x2 r3g, aðxÞ ¼ 2ðx1 þx2 Þ, f ðxÞ ¼ 4 for x A O , with boundary conditions: uðxÞ ¼ 2 þx1

for x2 ¼ 2, 2 rx1 r 3

uðxÞ ¼ 3 þx1

for x2 ¼ 3, 2 rx1 r 3

tðxÞ ¼ 2ðx1 þx2 Þðn1 ðxÞ þ n2 ðxÞÞ

for x1 ¼ 2 or x1 ¼ 3, 2 rx2 r3

The exact solution for this problem is uexact ðxÞ ¼ x1 þx2 , x A O . 5.4. Test 2 [10] Square domain O ¼ fðx1 ,x2 Þ : 1r x1 ,x2 r 2g, aðxÞ ¼ x21 þ x22 , f ðxÞ ¼ 8ðx21 þ x22 Þ for x A O , with boundary conditions: uðxÞ ¼ 1 þx21

for x2 ¼ 1, 1 rx1 r 2

uðxÞ ¼ 4 þx21

for x2 ¼ 2, 1 rx1 r 2

for x1 ¼ 2 or x1 ¼ 3, 2r x2 r3

The exact solution for this problem is uexact ðxÞ ¼ x1 þ x2 , x A O . 5.8. Test 6 Square domain O ¼ fðx1 ,x2 Þ : 2 rx1 ,x2 r 3g, aðxÞ ¼ expðx1 þ x2 Þ, f ðxÞ ¼ expðx1 þ x2 Þð6x1 þ 3x21 þ 6x2 þ 3x22 Þ for x A O , with boundary conditions: uðxÞ ¼ 8 þ x31 uðxÞ ¼ 27 þ x31

5.3. Test 1

qffiffiffiffiffiffiffiffiffiffiffiffi for x2 ¼  1x21

for x2 ¼ 2, 2 rx1 r3 for x2 ¼ 3, 2 rx1 r3

tðxÞ ¼ expðx1 þ x2 Þð3x21 n1 ðxÞ þ 3x22 n2 ðxÞÞ

for x1 ¼ 2 or x1 ¼ 3, 2 r x2 r 3

The exact solution for this problem is uexact ðxÞ ¼ x31 þ x32 , x A O . Computer programs were developed by using Matlab. The graphs of relative error have the number of nodes on the horizontal axis and the relative error on the vertical axis. 5.9. Numerical results for BDIE The results in Figs. 2 and 3 are for the Laplace and Poisson tests, respectively. In these tests, the domain is square and the exact solutions are linear and quadratic, respectively. Since we are using linear basis functions, there is no interpolation error for Laplace’s equation test, but there is interpolation error for Poisson’s equation as the exact solution is quadratic. In addition, other errors come either from discretisation of the domain into triangles or from calculating the boundary and domain integrals numerically; very good results and high rates of convergence are

1284

M.A. AL-Jawary, L.C. Wrobel / Engineering Analysis with Boundary Elements 35 (2011) 1279–1287

10–3

Mid−node collocation End−node collocation

collocation is still better than mid-node collocation, with similar rates of convergence to that of Figs. 1 and 2.

5.10. Numerical results for BDIDE

10–4

This subsection presents the solutions obtained for Laplace and Poisson equations in addition to tests 1–6 by using the BDIDE method (Figs. 10–17). The results in Figs. 2–17 demonstrate that both the BDIE and BDIDE methods are able to generate accurate solutions with high rates of convergence for the BVP (1)–(3). When comparing the solutions obtained by using both methods, it can be seen that the BDIDE method produced better results for the Laplace equation and for tests 1 and 3–5. However, more accurate results were obtained for the BDIE method for the Poisson equation and for tests 2 and 6. A possible explanation is the approximation of the flux t in the BDIDE method using linear basis functions for u living

10–5

10–6 101

102

103

104

Fig. 2. Relative error of the approximate solutions for Laplace’s equation with mid-node and end-node collocation; when J ¼ 1089, the relative errors are r  2:48  105 and 1.01  10  6, respectively.

10–3

Mid−node collocation End−node collocation

10–3

Mid−node collocation End−node collocation

10–4

10–4 10–5

10–5 10–6 101

10–6

101

102

103

104

Fig. 4. Relative error of the approximate solutions for Test 1 with mid-node and end-node collocation; when J ¼1089, the relative errors are r  1:62  104 and 1.38  10  6, respectively.

102

103

104

Fig. 3. Relative error of the approximate solutions for Poisson’s equation with mid-node and end-node collocation; when J ¼ 1089, the relative errors are r  4:64  105 and 4.40  10  6, respectively.

obtained for both mid-node and end-node collocation methods. However, the results using end-node collocation are better than mid-node. Tests 1 and 2 in Figs. 4 and 5 analyse problems with variable coefficients, so there is one more domain integral coming from the remainder, i.e. R a0. Therefore, there are discretisation and numerical integration errors for test 1. Also, there is interpolation error for test 2, as the exact solution is quadratic. End-node collocation results are better than mid-node collocation for both tests, with high rates of convergence which are similar to those in Figs. 2 and 3. In addition, for the circular domain in tests 3 and 4 (Figs. 6 and 7), an extra approximation error has been added which comes from approximating the boundary curve by polygons. It can be seen that mid-node collocation gives slightly better results than end-node collocation in this case; a possible reason is the approximate calculation of the value of c at the end nodes. The exact solutions for tests 5 and 6 are linear and cubic, respectively. It can be seen from Figs. 8 and 9 that using end-node

10–2

Mid−node collocation End−node collocation

10–3

10–4

10–5 101

102

103

104

Fig. 5. Relative error of the approximate solutions for Test 2 with mid-node and end-node collocation; when J¼ 1089, the relative errors are r  5:00  104 and 1.40  10  5, respectively.

M.A. AL-Jawary, L.C. Wrobel / Engineering Analysis with Boundary Elements 35 (2011) 1279–1287

10–2

10−2

Mid−node collocation End−node collocation

1285

Mid−node collocation End−node collocation

10−3

10–3 10−4

10–4 101

102

103

104

Fig. 6. Relative error of the approximate solutions for Test 3 with mid-node and end-node collocation; when J ¼3715, the relative errors are r  2:63  104 and 4.01  10  4, respectively.

10–2

Mid−node collocation End−node collocation

10−5 101

102

103

104

Fig. 9. Relative error of the approximate solutions for Test 6 with mid-node and end-node collocation; when J ¼1089, the relative errors are r  9:00  104 and 2.26  10  5, respectively.

10−5

10–3 10−6

10–4

10–5 101

102

103

104

Fig. 7. Relative error of the approximate solutions for Test 4 with mid-node and end-node collocation; when J¼ 3715, the relative errors are r  1:00  105 and 2.00  10  5, respectively.

10−2

Mid−node collocation End−node collocation

10−7 101

102

103

104

Fig. 10. Relative error of the approximate solutions for Laplace’s equation; when J¼ 1089, r  5:80  107 .

10−2

10−3

10−4

10−3

10−5

10−6 101

102

103

104

Fig. 8. Relative error of the approximate solutions for Test 5 with mid-node and end-node collocation; when J ¼1089, the relative errors are r  8:00  104 and 2.70  10  6, respectively.

10−4 101

102

103

104

Fig. 11. Relative error of the approximate solutions for Poisson’s equation; when J¼ 1089, r  1:00  104 .

1286

M.A. AL-Jawary, L.C. Wrobel / Engineering Analysis with Boundary Elements 35 (2011) 1279–1287

10−5

10−3

10−4

10−5

10−6

10−6

10−7 101

102

103

104

Fig. 12. Relative error of the approximate solutions for Test 1; when J¼ 1089, r  6:10  107 .

10−7 101

10−5

10−3

10−6

102

103

104

Fig. 13. Relative error of the approximate solutions for Test 2; when J¼ 1089, r  1:00  104 .

10−2

103

104

Fig. 15. Relative error of the approximate solutions for Test 4; when J ¼3715, r  5:00  107 .

10−2

10−4 101

102

10−7 101

102

103

104

Fig. 16. Relative error of the approximate solutions for Test 5; when J¼ 1089, r  6:80  107 .

10−2

10−3 10−3 10−4

10−5 101

102

103

104

Fig. 14. Relative error of the approximate solutions for Test 3; when J¼ 3715, r  4:56  105 .

10−4 101

102

103

104

Fig. 17. Relative error of the approximate solutions for Test 6; when J¼ 1089, r  2:00  104 .

M.A. AL-Jawary, L.C. Wrobel / Engineering Analysis with Boundary Elements 35 (2011) 1279–1287

on triangles; thus, T Fj ðxÞ is constant within each triangle. These approximations are appropriate for the Laplace equation and for tests 1, 4 and 5 as the solution to these tests are all linear, while the solution to the Poisson equation and tests 2, 3 and 6 is quadratic or cubic. The accuracy of the BDIE for tests 3 and 4 is also reduced by the approximation of the flux t at the boundary nodes, as there is a slight flux discontinuity at these points which is avoided in the BDIDE method.

6. Concluding remarks In this paper, the BDIE and BDIDE methods are developed and implemented for solving two-dimensional second-order linear elliptic mixed problems with variable coefficients. Convergence studies with mesh refinement show that the present methods possess excellent rates of convergence. The integrals in formulae (19)–(21) and (23) have a weak singularity. To calculate the boundary integrals we used a standard Gaussian quadrature rule. For the domain integrals, we have implemented a Gaussian quadrature rule with Duffy transformation by mapping the triangles into squares and eliminating the weak singularity. The following remarks apply to the present approach:

 A parametrix (Levi function), which is available for equations with variable coefficients, is used as a test function.

 The values of the unknown variables are obtained accurately with the present methods.

 Unlike in the standard BEM, the unknown function u is







approximated using linear basis functions living on triangles for both BDIE and BDIDE methods, allowing to obtain the values of u at interior points directly. As both u and t along the boundary are calculated in the BDIE method, we implemented mixed boundary elements with linear u and constant t to avoid the discontinuities of t at corner points. In this case, collocation was tested at the mid and end points of each boundary element. It was shown that end-node collocation generally provides higher accuracy than mid-node collocation. The only boundary variable in the BDIDE method is u along Neumann boundaries, thus, there is no need for collocation along Dirichlet boundaries. Thus, the problem caused by the discontinuity of the normal derivative at corner points is avoided. The generation, assembly and solution of the system of linear equations for the BDIE method are more complicated and thus take longer than that for the BDIDE method.

1287

References [1] Ang WT. A beginner’s course in boundary element methods. Boca Raton, USA: Universal Publishers; 2007. [2] Ang K-C. Introducing the boundary element method with MATLAB. International Journal of Mathematical Education in Science and Technology 2008;39: 505–19. [3] Brebbia CA, Telles JCF, Wrobel LC. Boundary element techniques. Berlin: Springer; 1984. [4] Paris F, Canas J. Boundary element method fundamentals and applications. New York: Oxford University Press; 1997. [5] Clements DL. Fundamental solutions for second order linear elliptic partial differential equations. Computational Mechanics 1998;22:26–31. [6] Ang WT, Kusuma J, Clements DL. A boundary element method for a second order elliptic partial differential equation with variable coefficients. Engineering Analysis with Boundary Elements 1996;18:311–6. [7] Sladek V, Sladek J, Zhang Ch. Domain element local integral equation method for potential problems in anisotropic and functionally graded materials. Computational Mechanics 2005;37:78–85. [8] Brunton I. Solving variable coefficient partial differential equations using the boundary element method. PhD thesis. University of Auckland, New Zealand; 1996. [9] Mikhailov SE. Localized boundary-domain integral formulations for problems with variable coefficients. Engineering Analysis with Boundary Elements 2002;26:681–90. [10] Mikhailov SE, Nakhova IS. Mesh-based numerical implementation of the localized boundary-domain integral equation method to a variable-coefficient Neumann problem. Journal of Engineering Mathematics 2005;51: 251–9. [11] Skerget L, Rek Z. Boundary-domain integral method using a velocity–vorticity formulation. Engineering Analysis with Boundary Elements 1995;15:359–70. [12] Skerget L, Hribersek M, Kuhn G. Computational fluid dynamics by boundarydomain integral method. International Journal for Numerical Methods in Engineering 1999;46:1291–311. [13] Jecl R, Skerget L, Petresin E. Boundary-domain integral method for transport phenomena in porous media. International Journal for Numerical Methods in Fluids 2001;35:39–54. [14] Popov V, Power H. The DRM-MD integral equation method: an efficient approach for the numerical solution of domain-dominant problems. International Journal for Numerical Methods in Engineering 1999;44:327–53. [15] Florez W, Power H, Chejne F. Multi-domain dual reciprocity BEM approach for the Navier–Stokes system of equations. Communications in Numerical Methods in Engineering 2000;16:671–81. [16] Peratta A, Popov V. A new scheme for numerical modelling of flow and transport processes in 3D fractured porous media. Advances in Water Resources 2006;29:42–61. [17] Katsikadelis JT. The analog equation method-A powerful BEM-based solution technique for solving linear and nonlinear engineering problems. In: Boundary element method XVI. Southampton: Computational Mechanics Publications; 1994. p. 167. [18] Katsikadelis JT, Nerantzaki MS. The boundary element method for nonlinear problems. Engineering Analysis with Boundary Elements 1999;23:365–73. [19] Hoitinga W. Direct minimization of equation residuals in least squares hp-finite element methods: a direct and iterative solution method. Internal Report. Delft University of Technology, The Netherlands; 2004.