A boundary integral equation method for the two-dimensional diffusion equation subject to a non-local condition

Engineering Analysis with Boundary Elements 25 (2001) 1±6

www.elsevier.com/locate/enganabound

A boundary integral equation method for the two-dimensional diffusion equation subject to a non-local condition W.T. Ang* Computational and Mathematical Sciences Group, Faculty of Information Technology, Universiti Malaysia Sarawak, 94300 Kota Samarahan, Malaysia Received 24 January 2000; revised 26 October 2000; accepted 26 October 2000

Abstract A boundary integral equation method is proposed for the numerical solution of the two-dimensional diffusion equation subject to a nonlocal condition. The non-local condition is in the form of a double integral giving the speci®cation of mass in a region which is a subset of the solution domain. A speci®c test problem is solved using the method. q 2001 Elsevier Science Ltd. All rights reserved. Keywords: Boundary element method; Diffusion equation; Laplace transform; Non-local condition

1. Introduction In non-dimensionalized form, the partial differential equation that governs two-dimensional linear and isotropic diffusion processes is given by 22 u 22 u 2u : 1 ˆ 2 2 2t 2x 2y

…1†

A class of problems of practical interest is to solve Eq. (1) for the unknown function u…x; y; t† for time t $ 0 in a twodimensional region R (on the 0xy plane) subject to the initial and boundary conditions u…x; y; 0† ˆ f …x; y† u…x; y; t† ˆ g…x; y; t†

for …x; y† [ R; for …x; y† [ C1 and t $ 0;

u…x; y; t† ˆ h…x; y†z…t† 2 ‰u…x; y; t†Š ˆ k…x; y; t† 2n

for …x; y† [ C2 and t $ 0;

…2† …3† …4†

for …x; y† [ C3 and t $ 0; …5†

and the non-local (integral) condition ZZ u…x; y; t† dx dy ˆ m…t† for t $ 0: S

…6†

where f, g, h, k and m are known and suitably prescribed functions, z is an unknown function to be determined, the region R is bounded by a simple closed curve C, the curves * Fax: 160-82-672301. E-mail address: wtang@mailhost.®t.unimas.my (W.T. Ang).

C1, C2 and C3 are non-intersecting and such that C1 < C2 < C3 ˆ C; S is a given subregion of R that is independent of time t and is bounded by a simple closed curve D given by D ˆ C2 < C4 ; the open curve C4 lies completely in the interior of R, and 2u=2n ˆ n´7u; n is the unit normal vector on C pointing away from R. From a physical standpoint, Eq. (6) speci®es the total amount of mass of the diffusing quantity u (or the total amount of heat energy, in the case of heat diffusion) which the region S can possess at any time t. Condition (4) with z…t† being unknown implies that the concentration of the diffusing quantity (or the temperature) on some part of the boundary must be controlled in a certain way in order that the region S carries the speci®ed amount of mass (or heat energy). For a sketch of the geometry of the problem, refer to Fig. 1. The class of problems de®ned by Eqs. (1)±(6) arises in many practical applications in heat transfer, control theory, thermoelasticity and medical sciences. A speci®c application which involves the use of the absorption of light to measure the concentration of a diffusing chemical is described by Noye and Dehghan [6]. Special cases of the problem, such as R being a rectangular region, have been solved directly by many researchers, e.g. Gumel et al. [5], Noye and Dehghan [6,7], Noye et al. [8], Cannon et al. [3], and Wang and Lin [11] using the ®nite-difference methods. With the exception of Ref. [6], the case S ˆ R; i.e. C 4 ˆ C1 < C3 ; was studied in all the references just cited. The present paper makes use of a boundary integral equation method (BIEM) for the numerical solution of Eqs. (1)±(6) in the Laplace transform (LT) space. The physical solution is recovered by using the Stehfest

0955-7997/01/$ - see front matter q 2001 Elsevier Science Ltd. All rights reserved. PII: S 0955-799 7(00)00068-0

2

W.T. Ang / Engineering Analysis with Boundary Elements 25 (2001) 1±6

the LT space is then to solve 22 U 22 U 1 2 pU ˆ 2f …x; y† 2x2 2y2

…8†

for U…x; y; p† ˆ L{u…x; y; t†; t ! p} (for …x; y† [ R and a suitably selected value of p) subject to for …x; y† [ C1 ;

U…x; y; p† ˆ G…x; y; p† U…x; y; p† ˆ h…x; y†F…p† Fig. 1. A sketch of the geometry of the problem. The solution domain is bounded by a simple closed curve C, which consists of three distinct parts, C1, C2 and C3, over which different boundary conditions are given. The speci®cation of mass (non-local condition) is given in the shaded region bounded by C2 < C4 : The curve C4 lies in the interior of the solution domain.

2 ‰U…x; y; p†Š ˆ K…x; y; p† 2n and ZZ S

algorithm [10] for the numerical inversion of the Laplace transformation. This approach should provide a useful and interesting alternative to the existing ®nite-difference methods. It can be easily implemented on the computer for solution domains of arbitrary shape. The BIEM for solving Eq. (1) with either u or 2u/2n completely speci®ed at each and every point on the boundary of the solution domain is well established but continues to be a research subject of considerable interest, see e.g. Rizzo and Shippy [9] and Chen et al. [4]. The application of the BIEM to diffusion problems with non-local conditions, such as the one de®ned by Eqs. (1)±(6), would lead to a formulation involving unknowns at nodal points throughout the physical solution domain, if the non-local conditions are not properly treated. In the present paper, the LT of the non-local condition (6) is recast into a form which involves the LT of u at points …x; y† on only the boundary C. Consequently, the linear algebraic equations in the BIEM do not involve unknowns at points in the interior of the physical domain R.

We shall reformulate the problem described in Section 1 by taking the LT of all functions and equations involved with respect to the time parameter t. To this end, let us de®ne the LT operator L on a function r…x; y; t† …t $ 0† by def

L{r…x; y; t†; t ! p} ˆ

Z1 0

r…x; y; t† exp…2pt† dt

…7†

where p is the LT parameter. In the present paper, we shall assume that p is real and positive. Applying L on Eqs. (1)±(6), we ®nd that the problem in

for …x; y† [ C3 ;

U…x; y; p† dx dy ˆ M…p†;

…10† …11†

…12†

where G…x; y; p† ˆ L{g…x; y; t†; t ! p}; F…p† ˆ L{z…t†; t ! p}; K…x; y; p† ˆ L{k…x; y; t†; t ! p} and M…p† ˆ L{m…t†; t ! p}: Notice that F…p† is an unknown function to be determined. Using Eq. (8) and the divergence theorem, we ®nd that Eq. (12) can be rewritten as ZZ Z 2 ‰U…x; y; p†Š ds ˆ pM…p† 2 f …x; y† dx dy: C 2 < C4 2n S …13† Notice that the right hand side of Eq. (13) is known and the left hand side containing an unknown function is expressed in terms of a boundary integral. To facilitate the task of solving Eqs. (8)±(11) together with Eq. (13), we make the substitution U…x; y; p† ˆ Upart …x; y; p† 1 V…x; y; p†;

…14†

where Upart …x; y; p† is any particular solution of Eq. (8). According to Atkinson [2], a particular solution of Eq. (8) is given by Upart …x; y; p† ˆ 2

2. Formulation in LT space

for …x; y† [ C 2 ;

…9†

ZZ T

f …j; h†G…x; y; j; h; p† dj dh;

…15†

where G…x; y; j; h; p† is given by Eq. (23) (in Section 3) and T is a two-dimensional region which may be chosen to assume any speci®c geometry convenient for computing numerically the double integral as long as R # T: In the speci®c example considered in Section 5, we take T ˆ R because R is a simple square region. With the substitution in Eq. (14), Eqs. (8)±(11) and (13) become: 2 2V 22 V 1 2 pV ˆ 0 2x2 2y2

…16†

W.T. Ang / Engineering Analysis with Boundary Elements 25 (2001) 1±6

V…x; y; p† ˆ G…x; y; p† 2 Upart …x; y; p†

for …x; y† [ C1 ;

V…x; y; p† < V …k†

3

for …x; y† [ C …k† …k ˆ 1; 2; ¼; N†; …25†

…17† V…x; y; p† ˆ h…x; y†F…p† 2 U part …x; y; p†

for …x; y† [ C2 ; …18†

2 2 ‰V…x; y; p†Š ˆ K…x; y; p† 2 ‰U …x; y; p†Š 2n 2n part

…19†

for …x; y† [ C3 ; Z C2 < C 4

ZZ 2 ‰V…x; y; p†Š ds ˆ pM…p† 2 p f …x; y† dx dy 2n S 2

Z C 2 < C4

2 ‰U …x; y; p†Š ds; 2n part …20†

Since Upart satis®es Eq. (8), the integral over C2 < C4 on the right hand side of Eq. (20) can be transformed back to a domain integral over S in order to rewrite Eq. (20) in the form Z 2 ‰V…x; y; p†Š ds C2 < C4 2n ˆ pM…p† 2 p

ZZ S

Upart …x; y; p† dx dy;

…21†

if one wishes to avoid evaluating 2‰Upart …x; y; p†Š=2n when dealing with the non-local condition. 3. BIEM For …j; h† [ R < C; the standard boundary integral solution of the homogeneous modi®ed Helmholtz equation (16) is given by (see, e.g. Rizzo and Shippy [9]) Z n 2 V…x; y; p† ‰G…x; y; j; h; p†Š l…j; h†V…j; h; p† ˆ 2n C o 2 2 G…x; y; j; h; p† ‰V…x; y; p†Š ds; …22† 2n where l…j; h† ˆ 1 if …j; h† [ R and 0 , l…j; h† , 1 if …j; h† [ C and 1 p K 0 …r p† G…x; y; j; h; p† ˆ 2 …23† 2p p where r ˆ …x 2 j† 2 1 …y 2 h†2 and K0 is the modi®ed Bessel function of the second kind and of order zero. The boundary C is discretized by putting N closely packed points on it. Two consecutive points are joined by straight line segments or elements C …1† ; C…2† ; ¼; C …N21† and C …N† : We make the following approximations: C < C …1† < C …2† < ¼ < C…N21† < C …N† ;

…24†

2 ‰V…x; y; p†Š < W …k† 2n

for …x; y† [ C …k† …k ˆ 1; 2; ¼; N†; …26†

where V …k† and W …k† …k ˆ 1; 2; ¼; N† are constants to be determined. With Eqs. (24)±(26), by letting …j; h† be the midpoint of C …i† given by …j …i† ; h…i† †; Eq. (22) can be rewritten approximately as N  Z 1 …i† X 2 V ˆ ‰G…x; y; j …i† ; h…i† ; p†Š ds V …k† …k† 2 2n C kˆ1  Z G…x; y; j …i† ; h…i† ; p† ds 2 W …k†

…27†

C …k†

for i ˆ 1; 2; ¼; N: In Eq. (27), in taking l…j …i† ; h…i† † ˆ 1=2; we assume that …j …i† ; h…i† † lies on a smooth part of C. Now if F…p† in Eq. (18) is given, then from Eqs. (17)± (19) either V …k† or W …k† (not both) is known over C …k† : It follows that Eq. (27) constitutes a system of N linear algebraic equations in N unknowns. Since F…p† is not known, additional equations are needed to complete the system of linear algebraic equations. The additional equations are obtained as outlined below. Discretize the curve C4 into J straight line segments denoted by C …N11† ; C…N12† ; ¼; C…N1J21† and C…N1J† and make the approximation 2‰V…x; y; p†Š=2n < W …N1k† for …x; y† [ C …N1k† …k ˆ 1; 2; ¼; J†;

…28†

where W …N1k† are constants to be determined and the normal derivative of V on C…N1k† (or C4) is in the direction pointing away from S. If the boundary elements obtained by discretizing C2 are labelled C …1† ; C …2† ; ¼; C…Q21† and C…Q† …Q , N† then Eq. (21) may be used to give the approximation Q X kˆ1

W …k† L…k† 1

J X

W …N1k† L…N1k†

kˆ1

ˆ pM…p† 2 p

ZZ S

Upart …x; y; p† dx dy;

…29†

where L…k† is the length of the straight line C …k† : Although Eq. (29) caters for the unknown parameter F…p†; it contains an additional J unknown constants W …N1k† …k ˆ 1; 2; ¼; J†: More equations are required. Applying the standard boundary integral equation in Eq. (22) for the regions S and R\S; respectively, and letting

4

W.T. Ang / Engineering Analysis with Boundary Elements 25 (2001) 1±6

…j; h† be the midpoint of C …N1i† …i ˆ 1; 2; ¼; J†; we obtain Q  Z X 2 1 …N1i† ‰G…x; y; j …N1i† ; h…N1i† ; p†Š ds ˆ V …k† 2V …k† 2n C kˆ1  Z G…x; y; j …N1i† ; h…N1i† ; p† ds 2 W …k†

V…x; y; p† at any interior point …j; h† approximately via V…j; h; p† <

kˆ1

2 ‰G…x; y; j …N1i† ; h…N1i† ; p†Š ds  …N 1 k† 2n C  Z …N1k† …N1i† …N1i† G…x; y; j ; h ; p† ds 2W C …N 1 k†

for i ˆ 1; 2; ¼; J;

(30)

and

(

N X

V …N1i† ˆ

V …k†

Z

2W 2

J X

2 ‰G…x; y; j …N1i† ; h…N1i† ; p†Š ds 2n

C …k†

kˆQ 1 1 …k†

)

Z

G…x; y; j

C …k†

( V

…N1k†

2W …N1k†

…N1i†

Z C …N 1 k†

kˆ1

Z C …N 1 k†

;h

…N1i†

; p† ds

2 ‰G…x; y; j …N1i† ; h…N1i† ; p†Š ds 2n )

2 W …k†

kˆ1

W …k† L…k† ˆ pM…p† 2 p

u…x; y; t† <

R

Upart …x; y; p† dx dy;

Z C…k†

2 ‰G…x; y; j; h; p†Š ds 2n )

G…x; y; j; h; p† ds :

  2P ln…2† X n ln…2† c n U x; y; ; t nˆ1 t

  2P ln…2† X n ln…2† c F ; z…t† < t nˆ1 n t

…33†

…34†

where P is a positive integer and c n ˆ …21†n1P min…n;P† X

mP …2m†! ; …P 2 m†!m!…m 2 1†!…n 2 m†!…2m 2 n†! …35†

(31)

ZZ

C …k†

The physical solution u…x; y; t† and z…t† can be recovered approximately from U…x; y; p† and F…p† by using an LT inversion technique. According to the Stehfest algorithm [10], which is nowadays increasingly used in applied mechanics for the numerical inversion of LT (e.g. Ang [1]), we obtain

mˆ‰…n 1 1†=2Š

where V …N1i† is the approximate constant value of V…x; y† over C…N1i† and …j …N1i† ; h…N1i† † is the midpoint of C…N1i† : Thus, after employing the boundary conditions (17)±(19), the system of linear algebraic equations (27)±(31) can be solved for the unknowns which include F…p†; V …i† and W …i† for i ˆ 1; 2; ¼; Q and also for i ˆ N 1 1; N 1 2; ¼; N 1 J; and either V …i† or W …i† for i ˆ Q 1 1; Q 1 2; ¼; N: Notice that if the boundary condition (17) or (19) holds over C …i† then V …i† or W …i† is known, respectively. This is why either V …i† or W …i† (not both) for i ˆ Q 1 1; Q 1 2; ¼; N are unknowns. (Boundary elements where Eq. (17) or (19) holds are labelled C …Q11† ; C …Q12† ; ¼; C …N21† and C…N† :† The derivation above is for S , R: What about the special case where S ˆ R? For this special case, Eqs. (29)±(31) can be replaced by just a single equation, i.e.

Z

4. LT inversion



G…x; y; j …N1i† ; h…N1i† ; p† ds

for i ˆ 1; 2; ¼; J;

N X

V …k†

J  X V …N1k†

Z

1 2

(

kˆ1

C …k†

1

N X

where ‰rŠ denotes the integer part of the real number r. As numerical techniques for inverting LT are highly susceptible to round-off errors, P cannot be selected to be as large as we like. On the other hand, choosing P to be too small may yield numerical results of lower accuracy. The optimum choice of P depends on the arithmetical precision of the computer (Stehfest [10]). Perhaps the best way to choose the optimum P is through testing the computer code of Eq. (34) on inverting known LTs of some elementary test functions. For the test problem in Section 5, programming in double precision using the Fortran 77 language on a Pentium III 450 MHz (64 Mb SDRAM) desktop computer, we ®nd that to achieve satisfactory results, the range of P which may be used is from 4 to 6, i.e. between 8 and 12 terms may be used in the approximation (34). With better machines, it may be possible to use a higher value of P, e.g. P ˆ 10; as reported by Stehfest [10].

…32†

and the unknowns are F…p†; V …i† and W …i† for i ˆ 1; 2; ¼; Q and either V …i† or W …i† for i ˆ Q 1 1; Q 1 2; ¼; N: Once V …k† and W …k† are determined, we can compute

5. A test problem For testing the BIEM, we construct a problem involving a square region R given by 0 , x , 1; 0 , y , 1: The region

W.T. Ang / Engineering Analysis with Boundary Elements 25 (2001) 1±6

ZZ S

u…x; y; t† dx dy ˆ a exp 2

5

p2 t 2

! 1 b;

…42†

where  Z1 Zx…1 2 x†  p …x 1 y† dy dx < 0:126 280 295 9; aˆ sin 2 0 0     Z1 Zx…1 2 x† px py bˆ exp sin dy dx 2 2 0 0 < 5:941 796 306 £ 1022 :

Fig. 2. A graphical comparison of the numerical and exact z…t†:

S where the speci®cation of mass is given is S ˆ {…x; y† : y , x…1 2 x†; x $ 0; y $ 0}:

…36†

The initial±boundary conditions and the non-local condition are:       p px py ‰x 1 yŠ 1 exp sin u…x; y; 0† ˆ sin 2 2 2 …37† for …x; y† [ R; " u…0; y; t† ˆ 1 1 exp 2

2

pt 2

!#

 sin

py 2

 …38†

for 0 , y , 1; !   p2 t p sin ‰1 1 yŠ u…1; y; t† ˆ exp 2 2 2     p py sin 1 exp 2 2

…39†

u…x; 1; t† ˆ exp 2

p2 t 2

!

 sin





p px ‰x 1 1Š 1 exp 2 2

 …40†

for 0 , x , 1;  px z…t† u…x; 0; t† ˆ sin 2



for 0 , x , 1;

…41†

Table 1 Numerical results for u…x; y; 1=10† at selected points …x; y† in the interior of the solution domain Point (x,y)

BIEM

Exact

(0.50,0.50) (0.30,0.80) (0.20,0.30) (0.70,0.10) (0.10,0.10) (0.90,0.80)

2.170 2.129 1.059 1.067 0.3745 4.194

2.161 2.127 1.053 1.050 0.3717 4.187

The exact solution for the speci®c problem above is given by !   p2 t p sin ‰x 1 yŠ u…x; y; t† ˆ exp 2 2 2     px py 1 exp sin …44† 2 2 with

! p2 t : z…t† ˆ exp 2 2

…45†

For the purpose of executing the BIEM, following Eq. (15), we take a particular solution of Eq. (8) to be given by " ! Z1 Z1 p ‰ j 1 hŠ sin Upart …x; y; p† ˆ 2 2 0 0 ! !# pj ph sin 1 exp 2 2 £ G…x; y; j; h; p† dj dh

for 0 , y , 1;

(43)

…46†

which can be evaluated using a numerical quadrature. The integrand in Eq. (46) behaves like ln(r) for small r, where r is the distance separating …x; y† and …j; h†: From numerical experiments, we found that integrals with such integrands may be evaluated with reasonably good accuracy by merely using a simple midpoint rule, provided that the domain of integration is ®nely discretized into many square cells. Thus, the midpoint rule is applied to compute Eq. (46). Discretizing the boundary of the square region into 80 equal length elements (each of length 0.05 units) and the curve y ˆ x…1 2 x† for x [ ‰0; 1Š into 20 elements with endpoints given by …k=20; k‰20 2 kŠ=400† for k ˆ 0; 1; 2; ¼; 20; we apply the BIEM to solve the speci®c problem described above in the LT space. We use P ˆ 4 in Eq. (34) to recover the solution in the physical space. The numerical results for the function z…t† thus obtained for t [ ‰0; 1=2Š are graphed and compared in Fig. 2 with the exact solution given by Eq. (45). As can be seen in Fig. 2, there is a reasonably good agreement between the numerical and the exact values of z. For t ˆ 1=10; we also compute u

6

W.T. Ang / Engineering Analysis with Boundary Elements 25 (2001) 1±6

numerically at selected points …x; y† in the interior of the square (solution) domain. The numerical values are compared and found to agree well with the exact ones given by Eq. (44). Refer to Table 1. When the number of elements used in the computation is doubled or trebled, an improvement in the accuracy of the numerical results may be observed. 6. Conclusions The problem of solving the two-dimensional diffusion equation subject to a non-local condition in the form of an integral taken over part of the solution domain is formulated in the LT space. A simple BIEM is then proposed for the numerical solution of the problem in the LT space. The physical solution is recovered through the use of a numerical technique for inverting LT. Such a method of solution is used to solve a speci®c problem which has a known exact solution. The numerical results obtained agree well with the exact solution. In the proposed method, the non-local condition in the LT space is recast in such a way that the LT of the unknown function u (in the diffusion equation) appears only in the integrand of a boundary integral. Although domain integrals do appear in the integral formulation, their integrands do not contain unknown functions. Thus, the manipulation of unknown data is carried out on only the boundary of the solution domain, as in a typical BIEM. This makes the proposed method easier to implement on the computer, compared with other techniques such as the ®nite-difference and the ®nite-element methods. It is particularly suitable for

problems having solution domains that are irregular in shape. The execution speed of the BIEM can be greatly improved if the computation for the different values of the LT parameter p is carried out in parallel using a computer with multiple processors. References [1] Ang WT. A crack in an anisotropic layered material under the action of impact loading. ASME J Appl Mech 1988;55:122±5. [2] Atkinson KE. The numerical evaluation of particular solutions for Poisson's equation. IMA J Numer Anal 1985;5:319±38. [3] Cannon JR, Lin Y, Matheson AL. The solution of the diffusion equation in two-space variables subject to the speci®cation of mass. J Appl Anal 1993;50:1±19. [4] Chen CS, Golberg MA, Hon YC. The method of fundamental solutions and quasi-Monte-Carlo method for diffusion equations. Int J Numer Methods Engng 1998;43:1421±35. [5] Gumel AB, Ang WT, Twizell EH. Ef®cient parallel algorithm for the two-dimensional diffusion equation subject to speci®cation of mass. Int J Comput Math 1997;64:153±63. [6] Noye BJ, Dehghan M. Explicit solution of two-dimensional diffusion subject to speci®cation of mass. Math Comput Simulation 1994;37:37±45. [7] Noye BJ, Dehghan M. A time-splitting ®nite difference method for two-dimensional diffusion with an integral condition. Commun Numer Methods Engng 1994;10:649±60. [8] Noye BJ, Dehghan M, van der Hoek J. Explicit ®nite difference methods for two-dimensional diffusion with a non-local boundary condition. Int J Engng Sci 1994;32:1829±34. [9] Rizzo F, Shippy DJA. A method of solution for certain problems of transient heat conduction. AIAA J 1970;8:2004±9. [10] Stehfest H. Numerical inversion of the Laplace transform. Commun ACM 1970;13:47±9 (see also p. 624). [11] Wang S, Lin Y. A numerical method for the diffusion equation with nonlocal boundary speci®cations. Int J Engng Sci 1990;28:543±6.