- Email: [email protected]
Numerical solution of the generalized Laplace equation by coupling the boundary element method and the perturbation method
Numerical solution of the generalized Laplace equation by coupling the boundary element method and the perturbation method R. Rangogni ENEL-DSR-CRIS. Via Ornato 90/14, Milan. Italy (Received January 1986)
The boundary element method (BEM) has, in general, some advantages with respect to domain methods, in so far as no internal discretization of the domain is required. Such an advantage, however, is lost if the BEM is used directly for the numerical solution of the generalized Laplace equation (GLE). This paper demonstrates that the GLE can also be dealt with advantageously by coupling the BEM with the perturbation method (PM). The technique involves transforming the starting equation (GLE) into an equation without partial first derivatives and then solving a few problems by BEM. The method requires an internal network but the unknowns are only on the boundary. The procedure is applied to numerically solve two test problems with known analytical solutions. Key words: mathematical models, generalized Laplace equation, boundary element method The boundary element method (BEM) has become a powerful and alternative method to the well known methods of finite differences (FDM) and finite elements (FEM). This is due mainly to the simplicity in the arrangement of the input data, to a lesser number of unknowns to be dealt with and to its applicability to problems defined on infinite domains. The application range of BEM, however, is confined to simple differential equations owing to the difficulty of finding a fundamental solution (FS) for a general differential operator. The applications of BEM mainly concern constant coefficient differential equations such as the Laplace, Poisson, Helmholtz, wave, heat equations etc. It is well known that the Laplace equation, the first one to which BEM was applied, is a particular case of a more general equation called the generalized Laplace equation (GLE).I Unfortunately, it is not advantageous to apply BEM
266
Appl. Math. Modelling, 1986, Vol. 10, August
directly to GLE since its integral formulation requires the computation of both domain and boundary integrals and the solution of a larger algebraic system of equations because the unknowns are on the boundary and in the domain. BEM has been used successfully2 to obtain an approximate numerical solution of the GLE. This involves transforming the GLE into an equivalent differential equation in which the first partial derivatives have disappeared, thus rendering BEM applicable. The transformed equation (GLET), however, has non-constant coefficients so that it would be difficult, in general, to find its FS which, in any case, would be rather complicated. The GLET in reference 2 is approximated, under suitable conditions, by a constant coefficient equation having a known FS, to which BEM can be applied in the usual fashion. Whenever the condition is not satisfied, the technique'- is not applicable. This paper presents a different approach for solving 0307-904X/86/04284-10/$03.00 © 1986 Butterworth & Co. (Publishers) Ltd
Solution of the generalized Laplace equation: R. Rangogni the G L E T when the technique in reference 2 cannot be applied. This goal is achieved by coupling BEM with a perturbation method (PM) that implies solving more than one problem, each of which has only boundary unknowns. The technique requires an internal network but can be rather coarse and, in any case, does not increase the number of unknown values. The numerical values of the original unknown function can be easily obtained from the numerical solution of the GLET. The procedure is presented for problems defined in twodimensional domains; the extension to three-dimensional problems is straightforward. To show how the technique works, two test problems are solved numerically and the results compared with the ones obtained using the technique in reference 2 and with a known analytical solution.
u = [a(x,y)p-gl = G, Ou
1 Oa
On
2a On
on Si
u + (a)-~g:
on $2
Ou ( l Oa ]u - - + g3-0,, ~Tn/ = 0
(6)
onS3
Equation (5) is suitable for applying BEM but there is still the difficulty of finding its FS. In reference 2, a parameter R is defined such that: maxf R = minf
on D
(7)
If R is close to unity, then equation (5) is approximated by the following equation: V2U + kZu = 0
Transformation of the governing equation Many physical problems lead to a GLE~: V.[a(x,y)VV(x,y)] = 0
on
D
(1)
where V. and V denote the divergence and gradient operators. Associated with equation (1) are the general boundary conditions: V(x,y) = gt(x,y)
on S,
avian = gz(x,y)
on $2 S~ u $2 u S3= S
aV/an =g3(x,y)V
on S 3
(2)
Equation (1), for the unknown function V = V(x,y), is supposed to be defined on a domain Dc_R e with a boundary S. The known function a = a(x,y) describes the physical properties of the domain D (i.e. it may be the electrical, thermal, hydraulic conductivity etc.), and it is supposed to be greater than zero on D. In this hypothesis on function a, equation (1) can be written in the form: V2V + 1-VV.Va = 0
(3)
a
which shows clearly that for constant a it reduces to the simpler Laplace equation. It is shown in reference (2) that the application of BEM to equation (3) is not advantageous, because a domain integral appears with the unknown function V in the domain integral that increases the number of unknown values to be searched for. The substitution: V(x,y) = u(x,y)[a(x,y)] -~-
(4)
however, leads equation (3) to the equation: V2u +f(x,y)u = 0 Ira[ 2 f ( x , y ) = 4a 2
V2a 2a
(5)
(8)
where k 2 is the average value o f f in D. Equation (8) can be easily solved by BEM and the values of function V can be obtained using equation (4). When parameter R is not close to unity, however, (i.e. if it is greater than, e.g, R > 5) the accuracy of the results decreases drastically so that the technique is now not applicable. Here, the problem is approached differently, using a perturbation method in order to obtain a better accuracy also for the case when R is larger than 5.
Perturbation method Instead of equation (5), the following family of equations is considered: 0~
V2llq-Ef(lg)=O
(9)
or, in other words, equation (5) is thought of as a perturbation of the Laplace equation; for e = 0, equation (9) becomes the Laplace equation while for e = 1, equation (5) is obtained. The solution of equation (9) is sought in the form: lg = bl0 "4- ELl I q- g21.12 4- . . . = ~ U j g j=0
(10)
j
Substituting function (10) and its partial derivatives into equation (9) and ordering according to the powers of e, gives: V-'Uo + e(V2ul + fuo) + e2(V2u2+ ful) + . . . .
0 (11)
Equation (11) is satisfied for all values of e e [0,1] if problems (12) and (13), Pi (i~>0), are satisfied (it is assumed that the boundary condition will be satisfied by the first term u. in expression (10)): V 2Zlo
=
Ilo ~
G I
Ouo
On 0u0 -
-
-
on D
0
on Sl
10a Uo + (a)~gz 2a On (
1 0a -
)uo = O
on S 2
P0
02)
on $3
with boundary conditions:
Appl. Math. Modelling, 1986,Vol. 10,August
267
Solution of the generalized Laplace equation: R. Rangogni V2Ul = --f(Uo)
on D
UI=0
on S t
OU1
On Oul
on S 2 P l
-0
a,, +
(
1 Oa~
=0
(13)
on $3
Problem Po is the Laplace problem which does not present any difficulty for solution by BEM, while problems Pj (j~> 1) are Poisson problems, with a known source term whenever problem Pi-I is solved. The FS for all problems Pi 0" ~>0) is:
G(P,Q)=log(r)
r=[P-Q[
PeD
Numerical
Q~D (14)
Using equations (14) and the second Green identity, it is possible to rewrite each problem Pi 0"I> 1) in the following integral form:
au,(P) =
( / Ouj u OG'~ds JsIG foGf(ui-l)dD(15) J~n)
M
k=l
k=l
f o r j = 1,qlj = OuJOn, where ~b0k, ~lk, q'0kand 6k denote, respectively:
o% son
•
6k =
I
e2(%)
max I Vi - 9"i] i × 100 maxV D
--
F~i(Vi - 15"i2}X maxV
(18)
100
D Subscript i extends over points Pi --- (xi,yi) c D U S where a numerical value is computed. In the two problems, parameter R is very high so the technique of reference 2 is not applicable. To compare the two procedures, however, the errors given by equations (18) obtained using the reference 2 technique are shown between brackets ( ) for both problems. For all problems, only solutions u0 and ut are computed. The computer times refer to a Vax 750.
Example 1
V2V+Ivv.Va=O
on
D - - {(x,y):
0
0
a(x,y) = (x + 1)(y + 2) V(x -= 0,y) = 0
(17)
GfuodD
(19)
V(x = 5,y) = 5 OV On
--=0
y=0,5;0~
D k
In equations (17), se is a parameter to describe a point on Sk. Integrals ~b0~, q'0k and qqk are computed analytically, while 6k is computed using nine Gauss points on each Dk. It is worth noting that the domain integral, in equations (15) or (16), does not contain the unknown function u, so the number of unknowns is limited by the number of boundary points (it does not exceed N).
268
e=(%)-
a
(16)
s[~n dS
In this section, two problems, both with known analytical solutions and solved numerically in reference 2, are solved numerically using the procedure presented in the previous paragraphs. To measure the error between the analytical solution V and the numerical solution V, two percentage errors are defined as:
The solution of the following problem is required:
M
au,(P) = Z q,kchO~-- Uk( Ot,k -- qhk) -- Uk+t~O,k-- ~ 6k
01k
examples
--
where n denotes the inward normal vector on S and a is the angle in P described by the vector r = I PQ[ when Q extends over boundary S. Equation (15) is as usual discretized in the BEM technique 3-6 for the numerical treatment, subdividing boundary S into subelements Sk (k= 1. . . . . N) and approximating functions uj and OuJOn, on each element Sk, by polynomials. The domain integral is computed by dividing domain D into subdomains Dk (k = 1. . . . . M) and approximating function uj-t by suitable shape functions. For example, in the numerical problems which will be solved, function ui is chosen to be linear on each S k , function OuJOn is chosen to be constant on each element Sk and function uj_l is approximated by a second order polynomial on each Dk. This choice of ui and OuJOn is suggested by the fact that linear functions have constant partial derivatives. Different approximations can be found in the literature. Using these approximations, the discretized form of equation (15) is:
&Ok=fGdS s~
When problems Pi are solved for a particular value o f j the solution u is obtained by writing expression (10) with e = 1 and, in practice, it turns out to be sufficient to compute only a few terms of the series (in the examples only u0 and Ul are computed). In conclusion, it is important to note that PM, here applied to G L E T equation (5), can be applied directly to equation (3); this would avoid computation of the second partial derivatives of function a that, from a mathematical point of view, could not exist. The presentation of the procedure, however, was carried out with reference to equation (5) for a comparison with the technique given in reference 2.
Appl. Math. Modelling, 1986,Vol. 10,August
The analytical solution of problem (19) is:
V(x,y) = ~
5
log (1 + x)
while f u n c t i o n f a n d parameter R are given by:
f(x,y) = ~
+ (2 +
(20)
Solution of the generalized Laplace equation: R. Rangogni '
'
'
'
5
'
op,
W
7
IF;--" o
o o
x
x
b
a
21- \',
Figure 1 N e t w o r k s u s e d in e x a m p l e 1
The boundary S is discretized by 40 linear equal elements. To check the influence of the internal discretization two internal meshes are used: the first very coarse and the second finer as shown in Figure 1. The results are checked on the 40 boundary points as well as on the internal points (five for the coarse mesh and 56 for the finer one). The errors given by equations (18) are: e~, = 2.34% e2
=
Mesh a
1.76%
e~ = 2.27%
(e= = 7.4%) Mesh b
(22)
(e2 = 17.9%)
e2 = 2.46% The maximum error e= has about the same value for the two computational networks considered. This explains the fact that e_, for network b (96 points) is larger than e, for network a (45 points). The computing time is 37.8s using mesh b against 30.5 s using the reference 2 technique. The improvement in the results should be noted.
Example 2 The electrical potential is required of a cylindrical capacitor with non-homogeneous dielectric. Mathematically, the problem is defined by:
1 V2V + - V V . V a = O a
on
D ~ {(x,y): 1 < r < 5 , x > O , y > O }
a(x,y) = 1 + r V= 1
r= 1
V=5
r=5
OV
r 2=x2 + y2 (23)
fory=O,l<~x<~5;x=O, an 0 <~y <~5 The analytical solution of problem (23) is: 4 [ 2r \ V(x,y) = 1 -I log ~/3ilog~]--~r)___
(24)
The boundary S is discretized with eight equal linear elements on the rectilinear sides and 10 equal elements on the curved sides. For this problem, function f and parameter R are:
1{ 1
t
Z
vn
3
4
Figure 2 Errors e= a n d e 2 a g a i n s t n u m b e r o f c o m p u t e d s o l u t i o n
The errors obtained with this problem, using 36 boundary points and 43 internal points are: e= = 0.35%
(2.68%)
e2 = 0.12%
(12.2%)
(26)
The run time is 28.87 s (21.84). This second example was also solved by applying PM directly to equation (3) in order to compare the two approaches. The number of boundary and internal points is always the same and the results are (using just V0 and Vl): e= = 1.2%
e2 = 1.8%
(27)
for a run time of 31 s. Essentially, the same numerical technique is applied to both equations (3) and (5). In the computation of domain integrals, however, the approximation of the term fui_ ~ is higher by one than that of term VaVu/_ I. This would entail relatively larger errors for the second approach. Further work should clarify this point. To give an idea of the speed of convergence the calculus was finally extended as far as solution V4 and the results are shown in Figure 2 (the run time was 95 s). The same solution was obtained by integrating on each Dk using only four Gauss points; the results were the same and the run time was 59.4 s.
Conclusions
-0
f(x,y)=~ 2(l~r) 2
0
r(l+r
R=19.3
(25)
A procedure is presented in this paper for numerically solving the G L E . In particular, this procedure should cover the case where the accuracy of technique given in reference 2 is not enough. The two techniques should be complementary to each other. The present procedure requires an internal mesh but, as shown in example 1, can be rather coarse, reducing the effort in preparing the input data. It is important to stress that the presence of the internal network does not increase the number of unknown values. The execution time is greater with respect to reference 2 using the same number of boundary points but the errors are drastically reduced. With regard to run time, it has to be noted that after computing ut (or Vt) for adding a new solution, it is not necessary to construct and factorize the matrix
Appl. Math. Modelling, 1986,Vo1.10,August
269
Solution of the generalized Laplace equation: R. Rangogni because, owing to the fact that for it> 1 all problems P~ have the same boundary conditions, it remains the same. It is adequate to compute and factorize the new known term. An important feature of the procedure presented in this paper is that it can be extended to other partial differential equations and also to some nonlinear ones.
References 1 Morse, P. M. and Feshbach, H. "Methods of theoretical physics', vols I-II, McGraw-Hill, New York, 1953
270 Appl. Math. Modelling, 1986,Vol. 10,August
2 3 4 5
6
Rangogni, R. and Occhi, R. "Numerical solution of the generalized Laplace equation by the boundary clement method', submitted to Appl. Math. Modelling Jaswon, M. A. and Symm, G. T. 'Integral equation methods in potential theory and elastostatics', Academic Press, London, 1977 Brebbia. C. A. and Walker, S. "Boundary element techniques in engineering', Butterworth and Co., London, 1980 Di Monaco, A. and Rangogni, R. 'Boundary element method: processing of the source term of the Poisson equation by means of boundary integral only', Proc. 4th Int. conJ~ finite elements in water resources, Hannover, West Germany, June 1982 (eds Holtz, K. P. et al.), Springer Verlag, Berlin Rangogni, R. 'The solution of the ngn-homogeneus Hehnholtz equation by means of the boundary element method', Appl. Math. Modelling, 1984, 8 (6), 442--444
The boundary element method (BEM) has, in general, some advantages with respect to domain methods, in so far as no internal discretization of the domain is required. Such an advantage, however, is lost if the BEM is used directly for the numerical solution of the generalized Laplace equation (GLE). This paper demonstrates that the GLE can also be dealt with advantageously by coupling the BEM with the perturbation method (PM). The technique involves transforming the starting equation (GLE) into an equation without partial first derivatives and then solving a few problems by BEM. The method requires an internal network but the unknowns are only on the boundary. The procedure is applied to numerically solve two test problems with known analytical solutions. Key words: mathematical models, generalized Laplace equation, boundary element method The boundary element method (BEM) has become a powerful and alternative method to the well known methods of finite differences (FDM) and finite elements (FEM). This is due mainly to the simplicity in the arrangement of the input data, to a lesser number of unknowns to be dealt with and to its applicability to problems defined on infinite domains. The application range of BEM, however, is confined to simple differential equations owing to the difficulty of finding a fundamental solution (FS) for a general differential operator. The applications of BEM mainly concern constant coefficient differential equations such as the Laplace, Poisson, Helmholtz, wave, heat equations etc. It is well known that the Laplace equation, the first one to which BEM was applied, is a particular case of a more general equation called the generalized Laplace equation (GLE).I Unfortunately, it is not advantageous to apply BEM
266
Appl. Math. Modelling, 1986, Vol. 10, August
directly to GLE since its integral formulation requires the computation of both domain and boundary integrals and the solution of a larger algebraic system of equations because the unknowns are on the boundary and in the domain. BEM has been used successfully2 to obtain an approximate numerical solution of the GLE. This involves transforming the GLE into an equivalent differential equation in which the first partial derivatives have disappeared, thus rendering BEM applicable. The transformed equation (GLET), however, has non-constant coefficients so that it would be difficult, in general, to find its FS which, in any case, would be rather complicated. The GLET in reference 2 is approximated, under suitable conditions, by a constant coefficient equation having a known FS, to which BEM can be applied in the usual fashion. Whenever the condition is not satisfied, the technique'- is not applicable. This paper presents a different approach for solving 0307-904X/86/04284-10/$03.00 © 1986 Butterworth & Co. (Publishers) Ltd
Solution of the generalized Laplace equation: R. Rangogni the G L E T when the technique in reference 2 cannot be applied. This goal is achieved by coupling BEM with a perturbation method (PM) that implies solving more than one problem, each of which has only boundary unknowns. The technique requires an internal network but can be rather coarse and, in any case, does not increase the number of unknown values. The numerical values of the original unknown function can be easily obtained from the numerical solution of the GLET. The procedure is presented for problems defined in twodimensional domains; the extension to three-dimensional problems is straightforward. To show how the technique works, two test problems are solved numerically and the results compared with the ones obtained using the technique in reference 2 and with a known analytical solution.
u = [a(x,y)p-gl = G, Ou
1 Oa
On
2a On
on Si
u + (a)-~g:
on $2
Ou ( l Oa ]u - - + g3-0,, ~Tn/ = 0
(6)
onS3
Equation (5) is suitable for applying BEM but there is still the difficulty of finding its FS. In reference 2, a parameter R is defined such that: maxf R = minf
on D
(7)
If R is close to unity, then equation (5) is approximated by the following equation: V2U + kZu = 0
Transformation of the governing equation Many physical problems lead to a GLE~: V.[a(x,y)VV(x,y)] = 0
on
D
(1)
where V. and V denote the divergence and gradient operators. Associated with equation (1) are the general boundary conditions: V(x,y) = gt(x,y)
on S,
avian = gz(x,y)
on $2 S~ u $2 u S3= S
aV/an =g3(x,y)V
on S 3
(2)
Equation (1), for the unknown function V = V(x,y), is supposed to be defined on a domain Dc_R e with a boundary S. The known function a = a(x,y) describes the physical properties of the domain D (i.e. it may be the electrical, thermal, hydraulic conductivity etc.), and it is supposed to be greater than zero on D. In this hypothesis on function a, equation (1) can be written in the form: V2V + 1-VV.Va = 0
(3)
a
which shows clearly that for constant a it reduces to the simpler Laplace equation. It is shown in reference (2) that the application of BEM to equation (3) is not advantageous, because a domain integral appears with the unknown function V in the domain integral that increases the number of unknown values to be searched for. The substitution: V(x,y) = u(x,y)[a(x,y)] -~-
(4)
however, leads equation (3) to the equation: V2u +f(x,y)u = 0 Ira[ 2 f ( x , y ) = 4a 2
V2a 2a
(5)
(8)
where k 2 is the average value o f f in D. Equation (8) can be easily solved by BEM and the values of function V can be obtained using equation (4). When parameter R is not close to unity, however, (i.e. if it is greater than, e.g, R > 5) the accuracy of the results decreases drastically so that the technique is now not applicable. Here, the problem is approached differently, using a perturbation method in order to obtain a better accuracy also for the case when R is larger than 5.
Perturbation method Instead of equation (5), the following family of equations is considered: 0~
V2llq-Ef(lg)=O
(9)
or, in other words, equation (5) is thought of as a perturbation of the Laplace equation; for e = 0, equation (9) becomes the Laplace equation while for e = 1, equation (5) is obtained. The solution of equation (9) is sought in the form: lg = bl0 "4- ELl I q- g21.12 4- . . . = ~ U j g j=0
(10)
j
Substituting function (10) and its partial derivatives into equation (9) and ordering according to the powers of e, gives: V-'Uo + e(V2ul + fuo) + e2(V2u2+ ful) + . . . .
0 (11)
Equation (11) is satisfied for all values of e e [0,1] if problems (12) and (13), Pi (i~>0), are satisfied (it is assumed that the boundary condition will be satisfied by the first term u. in expression (10)): V 2Zlo
=
Ilo ~
G I
Ouo
On 0u0 -
-
-
on D
0
on Sl
10a Uo + (a)~gz 2a On (
1 0a -
)uo = O
on S 2
P0
02)
on $3
with boundary conditions:
Appl. Math. Modelling, 1986,Vol. 10,August
267
Solution of the generalized Laplace equation: R. Rangogni V2Ul = --f(Uo)
on D
UI=0
on S t
OU1
On Oul
on S 2 P l
-0
a,, +
(
1 Oa~
=0
(13)
on $3
Problem Po is the Laplace problem which does not present any difficulty for solution by BEM, while problems Pj (j~> 1) are Poisson problems, with a known source term whenever problem Pi-I is solved. The FS for all problems Pi 0" ~>0) is:
G(P,Q)=log(r)
r=[P-Q[
PeD
Numerical
Q~D (14)
Using equations (14) and the second Green identity, it is possible to rewrite each problem Pi 0"I> 1) in the following integral form:
au,(P) =
( / Ouj u OG'~ds JsIG foGf(ui-l)dD(15) J~n)
M
k=l
k=l
f o r j = 1,qlj = OuJOn, where ~b0k, ~lk, q'0kand 6k denote, respectively:
o% son
•
6k =
I
e2(%)
max I Vi - 9"i] i × 100 maxV D
--
F~i(Vi - 15"i2}X maxV
(18)
100
D Subscript i extends over points Pi --- (xi,yi) c D U S where a numerical value is computed. In the two problems, parameter R is very high so the technique of reference 2 is not applicable. To compare the two procedures, however, the errors given by equations (18) obtained using the reference 2 technique are shown between brackets ( ) for both problems. For all problems, only solutions u0 and ut are computed. The computer times refer to a Vax 750.
Example 1
V2V+Ivv.Va=O
on
D - - {(x,y):
0
0
a(x,y) = (x + 1)(y + 2) V(x -= 0,y) = 0
(17)
GfuodD
(19)
V(x = 5,y) = 5 OV On
--=0
y=0,5;0~
D k
In equations (17), se is a parameter to describe a point on Sk. Integrals ~b0~, q'0k and qqk are computed analytically, while 6k is computed using nine Gauss points on each Dk. It is worth noting that the domain integral, in equations (15) or (16), does not contain the unknown function u, so the number of unknowns is limited by the number of boundary points (it does not exceed N).
268
e=(%)-
a
(16)
s[~n dS
In this section, two problems, both with known analytical solutions and solved numerically in reference 2, are solved numerically using the procedure presented in the previous paragraphs. To measure the error between the analytical solution V and the numerical solution V, two percentage errors are defined as:
The solution of the following problem is required:
M
au,(P) = Z q,kchO~-- Uk( Ot,k -- qhk) -- Uk+t~O,k-- ~ 6k
01k
examples
--
where n denotes the inward normal vector on S and a is the angle in P described by the vector r = I PQ[ when Q extends over boundary S. Equation (15) is as usual discretized in the BEM technique 3-6 for the numerical treatment, subdividing boundary S into subelements Sk (k= 1. . . . . N) and approximating functions uj and OuJOn, on each element Sk, by polynomials. The domain integral is computed by dividing domain D into subdomains Dk (k = 1. . . . . M) and approximating function uj-t by suitable shape functions. For example, in the numerical problems which will be solved, function ui is chosen to be linear on each S k , function OuJOn is chosen to be constant on each element Sk and function uj_l is approximated by a second order polynomial on each Dk. This choice of ui and OuJOn is suggested by the fact that linear functions have constant partial derivatives. Different approximations can be found in the literature. Using these approximations, the discretized form of equation (15) is:
&Ok=fGdS s~
When problems Pi are solved for a particular value o f j the solution u is obtained by writing expression (10) with e = 1 and, in practice, it turns out to be sufficient to compute only a few terms of the series (in the examples only u0 and Ul are computed). In conclusion, it is important to note that PM, here applied to G L E T equation (5), can be applied directly to equation (3); this would avoid computation of the second partial derivatives of function a that, from a mathematical point of view, could not exist. The presentation of the procedure, however, was carried out with reference to equation (5) for a comparison with the technique given in reference 2.
Appl. Math. Modelling, 1986,Vol. 10,August
The analytical solution of problem (19) is:
V(x,y) = ~
5
log (1 + x)
while f u n c t i o n f a n d parameter R are given by:
f(x,y) = ~
+ (2 +
(20)
Solution of the generalized Laplace equation: R. Rangogni '
'
'
'
5
'
op,
W
7
IF;--" o
o o
x
x
b
a
21- \',
Figure 1 N e t w o r k s u s e d in e x a m p l e 1
The boundary S is discretized by 40 linear equal elements. To check the influence of the internal discretization two internal meshes are used: the first very coarse and the second finer as shown in Figure 1. The results are checked on the 40 boundary points as well as on the internal points (five for the coarse mesh and 56 for the finer one). The errors given by equations (18) are: e~, = 2.34% e2
=
Mesh a
1.76%
e~ = 2.27%
(e= = 7.4%) Mesh b
(22)
(e2 = 17.9%)
e2 = 2.46% The maximum error e= has about the same value for the two computational networks considered. This explains the fact that e_, for network b (96 points) is larger than e, for network a (45 points). The computing time is 37.8s using mesh b against 30.5 s using the reference 2 technique. The improvement in the results should be noted.
Example 2 The electrical potential is required of a cylindrical capacitor with non-homogeneous dielectric. Mathematically, the problem is defined by:
1 V2V + - V V . V a = O a
on
D ~ {(x,y): 1 < r < 5 , x > O , y > O }
a(x,y) = 1 + r V= 1
r= 1
V=5
r=5
OV
r 2=x2 + y2 (23)
fory=O,l<~x<~5;x=O, an 0 <~y <~5 The analytical solution of problem (23) is: 4 [ 2r \ V(x,y) = 1 -I log ~/3ilog~]--~r)___
(24)
The boundary S is discretized with eight equal linear elements on the rectilinear sides and 10 equal elements on the curved sides. For this problem, function f and parameter R are:
1{ 1
t
Z
vn
3
4
Figure 2 Errors e= a n d e 2 a g a i n s t n u m b e r o f c o m p u t e d s o l u t i o n
The errors obtained with this problem, using 36 boundary points and 43 internal points are: e= = 0.35%
(2.68%)
e2 = 0.12%
(12.2%)
(26)
The run time is 28.87 s (21.84). This second example was also solved by applying PM directly to equation (3) in order to compare the two approaches. The number of boundary and internal points is always the same and the results are (using just V0 and Vl): e= = 1.2%
e2 = 1.8%
(27)
for a run time of 31 s. Essentially, the same numerical technique is applied to both equations (3) and (5). In the computation of domain integrals, however, the approximation of the term fui_ ~ is higher by one than that of term VaVu/_ I. This would entail relatively larger errors for the second approach. Further work should clarify this point. To give an idea of the speed of convergence the calculus was finally extended as far as solution V4 and the results are shown in Figure 2 (the run time was 95 s). The same solution was obtained by integrating on each Dk using only four Gauss points; the results were the same and the run time was 59.4 s.
Conclusions
-0
f(x,y)=~ 2(l~r) 2
0
r(l+r
R=19.3
(25)
A procedure is presented in this paper for numerically solving the G L E . In particular, this procedure should cover the case where the accuracy of technique given in reference 2 is not enough. The two techniques should be complementary to each other. The present procedure requires an internal mesh but, as shown in example 1, can be rather coarse, reducing the effort in preparing the input data. It is important to stress that the presence of the internal network does not increase the number of unknown values. The execution time is greater with respect to reference 2 using the same number of boundary points but the errors are drastically reduced. With regard to run time, it has to be noted that after computing ut (or Vt) for adding a new solution, it is not necessary to construct and factorize the matrix
Appl. Math. Modelling, 1986,Vo1.10,August
269
Solution of the generalized Laplace equation: R. Rangogni because, owing to the fact that for it> 1 all problems P~ have the same boundary conditions, it remains the same. It is adequate to compute and factorize the new known term. An important feature of the procedure presented in this paper is that it can be extended to other partial differential equations and also to some nonlinear ones.
References 1 Morse, P. M. and Feshbach, H. "Methods of theoretical physics', vols I-II, McGraw-Hill, New York, 1953
270 Appl. Math. Modelling, 1986,Vol. 10,August
2 3 4 5
6
Rangogni, R. and Occhi, R. "Numerical solution of the generalized Laplace equation by the boundary clement method', submitted to Appl. Math. Modelling Jaswon, M. A. and Symm, G. T. 'Integral equation methods in potential theory and elastostatics', Academic Press, London, 1977 Brebbia. C. A. and Walker, S. "Boundary element techniques in engineering', Butterworth and Co., London, 1980 Di Monaco, A. and Rangogni, R. 'Boundary element method: processing of the source term of the Poisson equation by means of boundary integral only', Proc. 4th Int. conJ~ finite elements in water resources, Hannover, West Germany, June 1982 (eds Holtz, K. P. et al.), Springer Verlag, Berlin Rangogni, R. 'The solution of the ngn-homogeneus Hehnholtz equation by means of the boundary element method', Appl. Math. Modelling, 1984, 8 (6), 442--444