- Email: [email protected]
Solution of potential problems using combinations of the regular and derivative boundary integral equations
Solution of potential problems using combinations of the regular and derivative boundary integral equations Marc S. Ingber Department USA
of Mechanical
Engineering,
University
of New Mexico,
Albuquerque,
NM,
Thomas J. Rudolphi Department USA
of Engineering
Science
and Mechanics,
towa
State
University,
Ames,
IA,
The feasibility of using combinutions of the boundary integrul equation (BIE) and the normal derivative of the boundary integral equation (DBZE) is investigated for the two-dimensional Laplace equation. By using the combinations of these two equations it is possible to derive Fredholm integral equations of either the first or second kind regardless of the boundary conditions. Although the Fredholm equations of the second kind are well conditioned, in general they provide less uccurate interior results than the Fredholm equations of the first kind und the classical direct boundary element method (DBEM). On the other hand, the Fredholm equations of the first kind cun provide poor solutions in regions close to the boundary where large gradients exist in the boundury flux. Keywords: boundary element tion, hypersingular integral
method,
derivative
Introduction The direct boundary element method (DBEM) has proven to be a valuable tool in the analysis of a wide class of boundary value problems. The method has certain advantages over some of the popular domain methods in that only the boundary needs to be discretized, thus significantly reducing the order of the resulting algebraic equations. Further, in many cases the resulting algebraic equations are better conditioned for the DBEM than for the domain methods. The conditioning of the boundary element equations is strongly dependent on the Green’s function used in formulating the integral equation, the approximation functions used for the source densities, the geometry of the boundary, the boundary discretization, and the boundary conditions.1-3 If the boundary conditions are entirely of the Neumann type, then the integral equation is a Fredholm equation of the second kind, and the associated linear algebra statement is generally very well condiAddress reprint requests to Dr. Ingber at the Department of Mechanical Engineering, University of New Mexico, Albuquerque, NM 87131, USA. Received
536
5 September
Appl.
1989; accepted
Math. Modelling,
20 March
1990
1990, Vol. 14, October
boundary
integral
equation,
Fredholm
integral
equa-
tioned. If the boundary conditions are entirely of the Dirichlet type, then the integral equation is a Fredholm equation of the first kind, and the conditioning of the associated linear algebra statement is worse than for the Neumann problem. Of course, there is a spectrum of problems between these two extremes with mixed boundary conditions. In practice, it is sometimes necessary to discretize the domain in such a way that some of the collocation points lie close together in order to resolve a singularity in the potential or to accurately approximate highly curved geometries. The set of linear equations can become ill-conditioned for these types of problems, and the accompanying finite arithmetic errors can severely degrade the accuracy of the numerical solution. Motivated by the authors’ previous work in grid optimization,4-h the original intent of this research was to formulate a hybrid boundary element method using combinations of the boundary integral equation (BIE) and the normal derivative of the boundary integral equation (DBIE) that would minimize the conditioning of the associated linear equations. In this hybrid method, which we designate HBEMl, the BIE is written at collocation nodes along the boundary where Neumann boundary conditions are prescribed, and the DBIE is written at collocation nodes along the boundary where
0 1990 Butterworth-Heinemann
Combinations
of regular and derivative
Dirichlet boundary conditions are prescribed. Although the HBEMl enhanced the diagonal dominance of the resulting linear equations and consequently the matrix conditioning, in most example problems considered, the method yielded less accurate interior results than the DBEM. Because of this, an alternative approach to the HBEMl was conceived in which the BIE is written at collocation nodes along the boundary where Dirichlet boundary conditions are prescribed and the DBIE is written at collocation nodes along the boundary where Neumann conditions are prescribed. This alternative approach is designated HBEMZ. The HBEMl yields a system of Fredholm equations of the second kind, while the HBEM2 yields a system of Fredholm equations of the first kind. For mixed problems the DBEM yields a system of equations that is neither entirely of the first or second kind. Several mathematical studies have been performed to determine the properties of boundary element methods yielding Fredholm integral equations of the first or second kind.3.7-‘X Nevertheless, the mathematical theory is far from complete, especially for the more popular point collocation methods. Further, for the most part these studies have not considered problems with mixed boundary conditions. One of the purposes of this paper is to expose some possible misconceptions about Fredholm equations of the first and second kinds through numerical experiment for problems that are either singular or nearly singular. Since the DBIE is hypersingular and exists only at points where the source densities are Holder continuously differentiable,14 isoparametric Overhauser spline elements are employed in the HBEMl, the HBEM2, and the DBEM. The Overhauser spline elements possess properties that make them outperform some of the traditional Lagrangian elements.” In the following sections we formulate the hybrid methods, outline a method of evaluating the DBIE, and look at several example problems. We show the feasibility of combining the regular and derivative boundary integral equations.
Theoretical formulation Let 4(p) be harmonic in the domain R C R2. By following the standard arguments of the DBEM,16 4(p) can be represented in terms of a boundary integral as
rl(PM(P) =
~G(P, 4 --gy 4(q) 9
boundary
equations:
M. S. lngber and T. J. Rudolphi
fOif[#fl 27&E R/ln(p) = the internal angle between the two tangents on either side ofp if p E r We refer to equation (1) as the boundary integral equation (BIE). The normal derivative of the BIE can be written in terms of the hypersingular integral equation asI $(p)
P
=
aG(p, q) w(q)
d2G(P, q)
J[ r
an Pan4 "(Y)-T%T
Y
P
fi I
(2) where a/an, denotes the derivative in the direction of the outward normal to I’ at the field point p. We refer to equation (2) as the derivative boundary integral equation (DBIE). Since the DBIE exists only at points where the approximations for 4(p) and a4(p)ldn,, are Holder continuously differentiable,14 the geometry must be smooth, and hence the coefficient of p)ldn,, (unlike the coefficient q(p) in the BIE) is always 7~at these points. We discretize both the BIE (equation (1)) and the DBIE (equation (2)) by dividing the boundary I into boundary elements and approximate &q) and a+(q)/&, with polynomial interpolates. To maintain the necessary smoothness of 4(q) and a&)/an,and at the same time reduce geometric errors, isoparametric Overhauser spline elements are employed. The formulation of Ortiz et ~1.‘~ is implemented here with one modification. In general, the parametric cubic curve c(t) is the following linear blend of two parabolas, p(r) and q(s) (see Figure I): c(t) = (1 - t)p(r) + tq(s)
(3)
In this way the first derivative continuity is maintained at the nodes p2 and p3. The modification made here to Ortiz et ul.‘s formulation occurs at corners or kinks in the geometry. Since the appropriate smoothness for the normal derivative cannot be maintained at these points, semidiscontinuous elements are employed on either side of a corner. This requires the use of five different sets of shape functions for the five types of elements shown in Figure 2, Element type 3 splines together two continuous elements. Element types 2 and 4 spline together a continuous and a semidiscontinuous element (the difference being the relative lo-
(1) 41dr-(cd
a+(q) - G(P, q)r
where I is the boundary of a, G(p, q) is the function given by G(p, q) = log Ip - 41, a/an, the derivative in the direction of the outward to the boundary at the source point q, and defined by
Green’s denotes normal q(p) is
PZ
Figure 1.
Parametric
P,
cubic curve c(t) for the Overhauser
spline
element
Appl.
Math.
Modelling,
1990,
Vol.
14, October
537
Combinations
of regular and derivative
boundary
equations:
LLY 4 ox
1 w!
Figure 2.
tinuous approach, when using the same number of nodes as the discontinuous approach, yields a higher-order method. Within each boundary element, 4(q) and a+(s)/&, are approximated as
5
2
Definition
0
Geometric
X
Collocation
sketch of the Overhauser
Al. S. lngber and T. J. Rudolphi
nodes nodes
spline element ~(4W,
types
=
5
(4)
#Si(d
i=l
cation of the two elements). Finally, element types 1 and 5 are automatically continuous at their middle node and require no spline at all. The Overhauser shape functions for element types 2, 3, and 4 are given in the appendix. An alternative approach would be to use discontinuous boundary elements, which would require only one set of shape functions. Although the initial work associated with developing five sets of shape functions is tedious, the continuous and semidiscon-
T&P) =
2 5 I [ Wi(y,~4 - &‘S;(q)G(p,
r=l
i=l
,
(8)
[A&+%1= [&yMjI
where +j and 4; represent the values of c$(q) and d+(q)/&,, respectively, at the global nodes. For the HBEMl, upon rearrangement (assembly) of equation (8) into a coefficient matrix multiplying a vector of unknowns equal to a known (load) vector, the 7~from the left-hand side of equations (6) and (7) will lie along the main diagonal of the coefficient matrix, thus improving the diagonal dominance of the coefficient matrix. In this way the conditioning of the linear equations
dn,dn,
VP. V,G(p, q) = -V;G(p, Appl.
= ,g $‘S;(q)
Math.
(5)
where 4: and &’ represent the values of 4(q) and t@(q)ldn,, respectively, at local node i within the eth boundary element r,; S; represents the Overhauser spline shape functions; and ne is either 3 or 4, depending on the element type. Inserting equations (4) and (5) into equations (1) and (2) yields
q)] dl-
is minimized. On the other hand, for the HBEM2, upon assembly the 7~from the left-hand side of equations (6) and (7) is absorbed into the load vector. The numerical evaluation of the hypersingular integral in equation (2) is discussed in the next section. Numerical evaluation integral equation
of the derivative
Modelling,
fundamental
solution for the Laplace equation,
q) = 0
1990,
Vol. 14, October
boundary
The kernel contained in the DBIE (equation (2)) given by aG(p, q)ldn, is regular and can be integrated by using standard Gaussian quadrature formulas. However, the kernel contained in the DBIE given by a*G(p, q)ldn,dn, is hypersingular and cannot be integrated in the ordinary sense. These types of kernels have been encountered in boundary element formulations for acoustic problems*8.1’ and for crack-opening problems.20~21 Through the use of simple solutions to the governing differential equation, Rudolphi has shown how the hypersingular term can be converted to three completely regular integrals. Alternatively, Stallybrass23 converts the hypersingular term to a single Cauchy principal value integral, and we follow this approach here. First, we note the identity
(n, . n,>V, . V,G(p, q) + (np X n,). [VP X v,G(p, q)l - nq. v, X [n, X V,G(p, q)l
Since G(p, q) is the free-space
538
r,
c
where N is the total number of boundary elements. In the DBEM, equation (6) is collocated at all the nodes within the boundary elements. In the HBEMl, equation (6) is collocated at nodes within the boundary elements where Neumann conditions are specified, and equation (7) is collocated at nodes within the boundary elements where Dirichlet conditions are specified. In the HBEM2, equation (6) is collocated at nodes within the boundary elements where Dirichlet conditions are specified, and equation (7) is collocated at nodes within the boundary elements where Neumann conditions are specified. For all methods we generate a system of linear equations that can formally be written as
~*G(P,4 =
+I&) an,
(9)
we have (10)
Combinations VP x V,GCp,
4
=
of regular and derivative
- 0,
x V,Gtp,
d
boundary
equations:
M. S. lngber and T. J. Rudolphi (11)
= 0
Further, $@)V, x [n, x v,G(~,q)l Therefore, /
44d
using equations
= V, x [&dn,
9) dr = - j- 12,. v q lan
1‘
(12)
(9)-(12), we obtain
@G(P, an
x V,G(p,q)l - V,d4d x In, x V,G(p,q)l
+
v, x [Ndn,
x
V,G(p,
dl dr
I nq. V,+(q) x 1% x V,G(p, 41 dr
(13)
1‘
Upon application of the divergence theorem the first integral on the right-hand side of equation (13) is seen to be zero. Hence using a vector identity for the scalar triple product on the second integral on the right-hand side of equation (13), we obtain ~*G(P,
/ I-
4)
$47) an an P
Y
a-=I ’n, x
V,G(p,
411.
In, x V,#q)l
dr
(14)
I‘
In the derivation of equation (14) we have implicitly assumed that $(q) is Holder continuously differentiable in order to apply certain of the vector operations. Although the hypersingular integral can now be interpreted in the Cauchy principal value sense, the integral on the right-hand side of equation (14) is still cumbersome because it contains a tangential derivative of 4(q), which can be evaluated only by numerical differentiation. We employ a device of Meyer rt a1.19 to circumvent this problem. From the preceding development we note that I
s
&9)&/ . v, x [n, x V,G(p, y)] dr = - / [n, x Vq44q)l. In, x
V,G(p,
(15)
41 d-
I
If we allowed $(q) = 6(p) = constant to be taken as a special case, the right-hand side of equation (15) would be zero, implying that
I-
+(p)n, . V, x [n, x V,G(p,
0
(16)
Hence we can write a2G(p,
/ 1‘
an,an
4)
444)dr = /w(p)
Y
- d44)~, . v, x b,? x V,G(P, 41 dr
(17)
1
in which the right-hand side of the equation is integrable in the Cauchy principal value sense. We use this form in all our numerical computations. In an alternative fashion, one can arrive at equation (17) by subtracting the term 4(p) from 4(q) in the integrand on the left-hand side of equation (13) and then adding back the appropriate term with 4(p) outside the integral on the added-back term. Then taking, in equation (2), the special case of $(p) = constant, that is, the analog of the rigid body motion argument of elastostatics, one concludes that the added-back term must be zero and hence the result of equation (17).
Numerical examples We consider three numerical examples in this section. Two of the examples have smooth solutions, while one example has a singularity at a reentrant corner. However, even though one of the example problems has an analytic exact solution, its domain is an elongated
ellipse causing singularlike behavior toward the ends of the major axis. The geometry for the first example is the unit square. Boundary conditions are consistent with the exact solution given by 6(x, y) = cash (x) sin (y)
(18)
We consider the problem in which Dirichlet conditions are specified along the sides x = 0 and y = 0 and Neumann conditions are specified along the sides x = 1 and y = 1. The results generated by the DBEM, HBEMl, and HBEM2 using ten uniform Overhauser spline elements per side are shown and compared to the analytic solution in Table I. All three methods yield excellent results, and there is very little to distinguish between the methods. As a convergence test, the maximum absolute error within the domain as calculated from the results at 81 evenly spaced points in the interior versus the number of nodes (degrees of freedom) is shown in Figure 3 for the three methods. The DBEM and the HBEMl provide essentially the same maxi-
Appl. Math. Modelling,
1990, Vol. 14, October
539
Combinations
of regular and derivative
boundary
equations:
Table 1. Comparison of the calculatedand analyticvalues of the potential4 at selected points within the domain for examplel.Numbersequence: DBEM,HBEMl,HBEM2, analytic
M. S. lngber and T. J. Rudolphi
The geometry for the second example is given by the L-shaped region shown in Figure 4. In this example the boundary conditions are consistent with the exact solution given by
X Y
0.9
0.7
0.5
0.3
0.1
+(r, 0) = r2’3 sin (28/3)
0.1
0.3
0.5
0.7
0.9
0.7872 0.7872 0.7872 0.7872
0.8186 0.8186 0.8186 0.8188
0.8827 0.8827 0.8830 0.8833
0.9822 0.9823 0.9828 0.9832
1.1209 1.1209 1.1222 1.1226
0.6473 0.6473 0.6474 0.6474
0.6731 0.6731 0.6732 0.6734
0.7259 0.7259 0.7262 0.7264
0.8077 0.8078 0.8083 0.8086
0.9221 0.9223 0.9229 0.9232
0.4817 0.4816 0.4818 0.4818
0.5008 0.5009 0.5010 0.5011
0.5401 0.5402 0.5403 0.5406
0.6010 0.6012 0.6014 0.6018
0.6861 0.6864 0.6867 0.6871
0.2969 0.2968 0.2969 0.2970
0.3086 0.3087 0.3087 0.3089
0.3328 0.3329 0.3329 0.3332
0.3702 0.3706 0.3705 0.3709
0.4225 0.4231 0.4230 0.4235
0.1001 0.1001 0.1001 0.1003
0.1041 0.1041 0.1041 0.1044
0.1122 0.1124 0.1123 0.1126
0.1248 0.1251 0.1249 0.1253
0.1419 0.1429 0.1422 0.1431
Degrees
Figure 3. The maximum a function of the number
where r and 0 are the polar coordinates. Dirichlet boundary conditions are specified on the sides adjacent to the reentrant corner, and Neumann conditions are specified elsewhere. A singularity of order r-1’3 exists in the normal derivative at the origin. To resolve the singularity, we pack nodes close to the origin. In practice, this is sometimes done by using a grid optimization scheme. Here, we use the geometric stretching parameter A, given by the ratio of the lengths of two adjacent elements, on the two sides next to the reentrant corner and evenly spaced elements on the other sides. As an example, the case for A = 0.8 and 13 elements per side is shown in Figure 4. For this particular problem the results are more difficult to interpret. We focus attention on the case in which 20 Overhauser spline elements are placed on each of the six sides. The interior results for the DBEM, HBEMl, and HBEM2 in the first quadrant are shown in Table 2 for a stretching parameter of 0.9. This particular stretching parameter is chosen because it was determined to yield the most accurate results for the DBEM. Table 2 shows that the HBEM2 and the DBEM provided very good results. However, for this problem the results provided by the HBEMl are unacceptable. The maximum absolute error within the domain for the three methods as a function of the stretching parameter A as calculated from the results at 243 evenly spaced points in the interior is shown in Figure 5. The figure shows the general trend that appeared in the first example problem, which was that the HBEM2 pro-
of Freedom
absolute error within the domain of nodes for example 1
as
mum absolute errors, while the HBEM2 has a somewhat smaller maximum absolute error. Figure 3 indicates that the convergence rate is superlinear for all three methods. It is interesting to compare the condition numbers for the three methods. With 84 degrees of freedom the DBEM has a condition number of 153, the HBEMl has a condition number of 81, and the HBEM2 has a condition number of 4588.
540
(19)
Appl. Math. Modelling,
1990, Vol. 14, October
Figure 4.
The discretizedL-shaped domain for example 2
Combinations
of regular and derivative
Table 2. Comparison of the calculated and analytic values of the potential I#Jat selected points within the domain for example 2 computed with a stretchingfactorof A = 0.9. Number sequence: DBEM,HBEMl,HBEMZ,analytic
boundary
equations:
M. S. lngber and T. J. Rudolphi
Table 3. Comparison of the calculatedvalues, computed by the HBEMl with three different stretching factors,and analytic values of the potential4 at selected points within the domain for example 2. Number sequence: A = 0.7.A = 0.8.A = 0.9, analytic
X
Y 0.5
0.4
0.3
0.2
0.1
0.1
0.2
0.3
0.4
0.5
0.5064 0.5498 0.5059 0.5061
0.4721 0.5143 0.4716 0.4718
0.4429 0.4842 0.4424 0.4426
0.4181 0.4590 0.4176 0.4179
0.3971 0.4380 0.3966 0.3969
0.4286 0.4686 0.4281 0.4283
0.3937 0.4323 0.3933 0.3935
0.3653 0.4028 0.3649 0.3651
0.3422 0.3790 0.3418 0.3420
0.3232 0.3598 0.3228 0.3230
0.3436 0.3798 0.3432 0.3434
0.3088 0.3433 0.3085 0.3087
0.2824 0.3156 0.2821 0.2823
0.2621 0.2944 0.2618 0.2620
0.2461 0.2780 0.2458 0.2460
0.2480 0.2802 0.2478 0.2479
0.2155 0.2455 0.2153 0.2154
0.1936 0.2220 0.1934 0.1935
0.1779 0.2053 0.1778 0.1779
0.1661 0.1928 0.1659 0.1661
0.1357 0.1633 0.1356 0.1357
0.1121 0.1369 0.1119 0.1121
0.0988 0.1219 0.0987 0.0988
0.0900 0.1120 0.0899 0.0901
0.0837 0.1049 0.0836 0.0837
X V
0.1
0.2
0.3
0.4
0.5
0.5
0.4335 0.5059 0.5498 0.5061
0.4007 0.4715 0.5143 0.4718
0.3722 0.4421 0.4842 0.4426
0.3474 0.4172 0.4590 0.4179
0.3255 0.3960 0.4380 0.3969
0.4
0.3627 0.4284 0.4686 0.4283
0.3297 0.3934 0.4323 0.3935
0.3022 0.3648 0.4028 0.3651
0.2792 0.3415 0.3790 0.3420
0.2596 0.3223 0.3598 0.3230
0.3
0.2856 0.3438 0.3798 0.3434
0.2530 0.3089 0.3433 0.3087
0.2277 0.2823 0.3156 0.2823
0.2077 0.2617 0.2944 0.2620
0.1913 0.2455 0.2780 0.2460
0.2
0.1992 0.2488 0.2802 0.2479
0.1689 0.2160 0.2455 0.2154
0.1480 0.1937 0.2220 0.1935
0.1327 0.1778 0.2053 0.1779
0.1207 0.1657 0.1928 0.1661
0.1
0.0977 0.1372 0.1633 0.1357
0.0757 0.1129 0.1369 0.1121
0.0631 0.0992 0.1219 0.0988
0.0545 0.0901 0.1120 0.0901
0.0481 0.0834 0.1049 0.0837
Table 4. The condition numbers associated with the DBEM, HBEMl, and HBEMZ as a function of the stretchingparameter in examole 2 Stretching parameter, A 1.0
0.9 0.8 0.7 0.6
06
0.7
0.8
Stretching
Parameter
0.9
1
A
Figure 5. The maximum absolute error within the domain a function of the stretchingparameter for example 2
as
vided the most accurate results followed by the DBEM and then the HBEMl. However, in this example the spread in accuracies of the three methods is wider. The minimum absolute error for the HBEM2 occurred with no stretching, while the minimum absolute error for the HBEMl occurred for A = 0.8. It is interesting to note that in a very narrow range of the stretching parameter it is possible to obtain good solutions with the
Condition number DBEM 138 329 2229 28,889 6,235,200
HBEMl 79 23 63 390 2609
HBEM2 46,964 67,646 399,680 4,720,OOO 71,856,OOO
HBEMI. The interior results for the HBEMI in the first quadrant for X = 0.7, A = 0.8, and h = 0.9 are shown in Tublr 3. Only for A = 0.8 are the solutions acceptable. The condition numbers for the three methods as a function of the stretching parameter are shown in Table 4. As was expected, the HBEMI provided much smaller condition numbers than did the other two methods. Even with a very large condition number for the HBEM2 at a stretching of A = 0.6 the method still provided very good results. (The results were calculated by using single precision on a DECstation 3 100.) It is also interesting to see how well the various methods predict the normal derivative along the boundary near the singularity. These results are shown in Table 5. Although the HBEMl was better able to resolve the singularity in the normal derivative close to the singularity, this advantage was quickly lost at points away from the reentrant corner. The domain for the third example problem is given by the ellipse of aspect ratio 5 as shown in Figure 6.
APPl.
Math. Modelling,
1990, Vol. 14, October
541
Combinations
of regular and derivative
boundary
Table 5. Comparison of the calculated and analytic values of the flux a+lan at points along the boundary near the reentrant corner for examole 2
equations:
M. S. lngber and T. J. Rudolph
Table 7. Comparison of the calculated and analytic values of the potential 4 at selected interior points for example 3
@ian X
DBEM
Y
0.0051 0.0154 0.0325 0.0515 0.0726
0.0000 0.0000 0.0000 0.0000 0.0000
- 4.8235 -2.5818 -2.1468 -1.7912 - 1.6025
HBEMl
HBEM2
Analytic
-
-4.8184 - 2.5791 -2.1445 - 1.7892 - 1.6007
-
3.8665 2.6809 2.0898 1.7924 1.5985
elliptic domain
in the
3.9936 2.8147 2.2007 1.8898 1.6865
Figure 6. The portion of the discretized first quadrant for example 3
Table 6. Comparison of the calculated and analytic values of the flux d$ddnat points along the boundary in the first quadrant for example 3 adan X
Y
DBEM
HBEMI
5.0000 4.8496 4.6617 4.4267 4.1331 3.7660 3.3071 2.7335 2.0165 1.1203 0.0000
0.0000 0.2440 0.3616 0.4649 0.5628 0.6578 0.7500 0.8373 0.9151 0.9746 1 .oooo
35.2293 8.4780 6.5075 10.2323 8.7184 7.6508 6.1798 4.4945 2.5517 0.3313 - 2.1724
9.3677 13.5452 12.0922 10.7773 9.4217 7.9854 6.4198 4.6802 2.7242 0.4419 -2.0813
We assume Dirichlet boundary conditions with the exact solution given by cb(X,Y) = X2 - y2 + 2xy - 1
Analytic 10.0000
13.5534 12.2480 10.8868 9.4982 8.0345 6.4458 4.6852 2.7092 0.4836 - 2.0000
associated (20)
For this problem the HBEM2 is identical to the DBEM. A geometric stretching parameter is again chosen to place smaller elements near the ends of the major axis to resolve the large curvature and rapid variations in the normal derivative of #J in those regions. With stretching, the isoparametric Overhauser spline elements are able to approximate accurately the curved geometry. For the grid containing 40 isoparametric Overhauser spline elements as shown in Figure 6 the DBEM had a condition number of 174, while the HBEMI had a condition number of 5. The solutions for the flux (a@&) at points along the boundary in the first quadrant are shown in Tuble 6. The HBEMl is seen to yield uniformly more accurate results for the surface fluxes. The DBEM yields extremely poor solutions near the ends of the major axis with errors in the flux of more than 300%. The interior solutions for the two methods are compared to the analytic results in Tuble 7. The
542
Appl.
Math.
Modelling,
1990,
Vol. 14, October
X
Y
DBEM
HBEMI
Analytic
0.5444 0.5444 0.5444 0.5444 0.5444 0.5444 1.6333 1.6333 1.6333 1.6333 2.7222 2.7222 2.7222 2.7222 3.8111 3.8111 3.8111 3.8111 4.9000 4.9000 4.9000 4.9500
- 0.9000 - 0.7000 -0.1000 0.1000 0.7000 0.9000 - 0.5000 - 0.3000 0.3000 0.5000 - 0.7000 -0.1000 0.1000 0.7000 - 0.5000 - 0.3000 0.3000 0.5000 -0.1000 0.0000 0.1000 0.0000
- 2.3644 - 1.9574 - 0.8240 - 0.6050 - 0.4265 -0.5111 - 0.2252 0.5934 2.5643 3.0597 2.1060 5.8778 6.9727 9.7708 9.5121 11.2278 15.8306 17.1824 23.8131 25.0066 25.7550 25.9407
- 2.4852 - 2.0772 - 0.9454 - 0.7282 - 0.5586 - 0.6479 -0.3180 0.4947 2.4494 2.9401 2.0457 5.7810 6.8660 9.6426 9.4379 11.1174 15.6629 17.0160 22.0811 23.0331 23.8617 23.5199
- 2.4936 - 1.9558 - 0.8225 - 0.6047 -0.4314 - 0.5336 -0.2156 0.5978 2.5578 3.0511 2.1093 5.8560 6.9449 9.7316 9.4635 11.1479 15.7212 17.0857 22.0200 23.0100 23.9800 23.5025
HBEMl yields more accurate results near the ends of the major axis. It is in this region that the DBEM yielded very poor surface fluxes. However, despite the fact that the HBEMl yielded uniformly more accurate surface fluxes, away from the ends of the major axis the DBEM is more accurate. This result is surprising in light of the calculated surface fluxes and indicates that the HBEM2 (DBEM in this example) tends to smooth boundary errors as one moves into the interior analogous to a numerical St. Venant effect. In fact, the largest percent error in the interior at the selected points shown in Table 7 is 7.5% for the DBEM and over 20% for the HBEMl . Conclusions The feasibility of using combinations of the BIE and the DBIE has been demonstrated in the example problems. The three examples were chosen to show the relative merits of the classical DBEM compared to alternative methods for generating Fredholm integral equations of the first and second kinds for a smooth problem, for a singular problem, and for a “nearly” singular problem. For the smooth problem the two hybrid methods as well as the classical DBEM were shown to yield very good results, the results provided by the HBEM2 being slightly more accurate. For the singular problem with the reentrant corner the HBEMI provided poor solutions except for one value of the stretching parameter. Even in cases with extreme element stretching, the HBEM2 provided the most accurate results for the singular problem despite possessing very large condition numbers. The third example is perhaps the most perplexing. In this example the DBEM provided better solutions in a large portion of the domain despite yielding uni-
Combinations
of regular
and derivative
formly less accurate surface fluxes. Near the ends of the major axis the HBEMl was more accurate. In hindsight, one might argue that the DBIE should be written at nodes near the ends of the major axis and the BIE at all other nodes. Nevertheless, as a general statement, it appears that the HBEM2, which yields a system of Fredholm integral equations of the first kind, provides the most accurate interior solutions of the three methods. The general methodology discussed for generating systems of Fredholm integral equations of the first and second kinds is applicable to a much wider class of differential equations, and we assert our conclusions only for the two-dimensional Laplace equation, which generally does not yield excessive condition numbers because of the logarithmic fundamental solution. For other fundamental solutions it is quite possible that finite arithmetic errors and small errors in determining the boundary data could become critical in determining the overall accuracy of the solution.
boundary 14 15
I6 I7
18
19
20
21
22
References 1
2 3
4
5
6
7
8
Vable. M. Making the boundary element method less sensitive to changes or errors in the input data. Internar. J. Numer. Mefhods Engrg. 1987, 24, 1533-1540 Goldberg, M. A. Solution Methodsfor Integral Equations: Theory and Applications. Plenum Press, New York, 1978 Niessner, H. Significance of kernel singularities for the numerical solution of Fredholm integral equations. Bounda~ Elements IX, Vol. I, ed. C. A. Brebbia. Springer-Verlag, Berlin, 1987, pp. 213-227 Ingber, M. S. and Mitra, A. K. Grid optimization for the boundary element method. Internat. 1. Numer. Method Engrg. 1986, 23, 2121-2136 Mitra, A. K. and Ingber, M. S. The numerical solution of Stokes flow in a domain with re-entrant boundaries by the boundary element method. Internat. .I. Numer. Method Huids 1988, 8, 327-338 Ingber, M. S. and Mitra, A. K. Grid redistribution based on measurable error indicators for the direct boundary element method. Engrg. Analysis, in press Atkinson, K. E. A Survey of Numerical Methods for the Solution of Fredholm Integral Equations of the Second Kind. SIAM, Philadelphia, 1976 Giroire, J. and Nedelec, J. C. Numerical solution of an exterior Neumann problem using a double layer potential. Math. Comp. 1978, 32, 973-990 Hsiao, G. C., Koop, P. and Wendland, W. L. A Galerkin collocation method for some integral equations of the first kind. Computing 1980, 25, 89-130 Arnold, D. N. and Wendland, W. L. The convergence of spline collocation for strongly elliptic equations on curves. Numer. M&h. 1985, 47, 313-341 Wendland, W. L. and Christiansen, S. On the condition number of the influence matrix belonging to some first kind integral equations with logarithmic kernel. Appl. Anal. 1986.21, 175-183 Chandler, G. A. and Graham, I. G. Product integration-collocation methods for noncompact integral operator equations. Math. Comp. 1988, 50, 125-138 Costabel, M., Ervin, V. J. and Stephan, E. P. On the convergence of collocation methods for Symm’s integral equation on open curves. Math. Camp. 1988, 51, 167-179
23
equations:
M. S. lngber
and T. J. Rudolphi
Gtinter, N. M. Potential Theory and Its Applications to Basic Problems of Mathematical Physics. Ungar, New York, I967 Ortiz. J. C.. Walters. H. G. and Gipson. G. S. Development of Overhauser splines as boundary ‘elements. Boundary Elements IX, Vol. 1, ed. C. A. Brebbia. Springer-Verlag, Berlin, 1987, pp. 401-407 Brebbia, C. A., Telles, J. C. F. and Wrobel, L. C. Boandary EIement Techniques. Springer-Verlag, Berlin, 1984 Ingber, M. S. and Mitra, A. K. The evaluation of the normal derivative along the boundary in the direct boundary element method. Appl. Math. Modelling 1989, 13, 32-40 Burton, A. J. and Miller, G. F. The application of integral equation methods to the solution of some exterior boundaryvalue problems. Proc. Roy. Sot. London Ser. A 1971, 323, 201-210 Meyer, W. L., Bell, W. A.. Zinn, B. T. and Stallybrass, M. P. Boundary integral solutions ofthree dimensional acoustic radiation problems. J. Sound Vibr. 1978, .59(2), 245-262 Ioakimidis, N. I. A new singular integral equation for the classical crack problem in plane and antiplane elasticity. In/ernat. .I. Fracture- 1983, 21, 115-122 Budreck. D. E. and Achenbach. J. E. Scattering from threedimensional planar cracks by the boundary integral equation method. J. Appl. Mech. 1988, 55, 405-412 Rudolphi. T. J. The use of simple solutions in the regularization of hypersingular boundary integral equations. Comp. Math. Appl., in press Stallybrass. M. P. On a pointwise variational principle for the approximate solution of linear boundary value problems. J. Math. Mech. 1967. 16, 1247-1286
Appendix The shape functions for element types 2,3, and 4 (Fig14~ 2) are given below. Their derivations are straightforward using equation (3) to blend together the appropriate continuous and semidiscontinuous shape functions. Element type 2: N:(t) = -9P110 N#t)
+ 9P15 - 9t/lO
= 2t’ - 7t2/2 + t/2 + 1
A’s(t) = - St’/5 + 11t2/5 + 2tt5 N:(t)
= t3/2 - 1’12
Element type 3: N_:(t) = - t’/2 + t2 - t/2 N:(t)
= 3t3/2 - V/2
N:(t)
= -3t3/2
N:(t)
= t’l2 - t2/2
+ 1
+ 2t’ + t/2
Element type 4: N:(t)
= - t3/2 + t2 - t/2
N!?(t) = 8t3/5 - 13t2/5 + 1 Aqt)
= -2P
+ 5t2/2 + t/2
N;(t) = 9t3/10 - 9t2110
Appl.
Math. Modelling,
1990, Vol. 14, October
543
of Mechanical
Engineering,
University
of New Mexico,
Albuquerque,
NM,
Thomas J. Rudolphi Department USA
of Engineering
Science
and Mechanics,
towa
State
University,
Ames,
IA,
The feasibility of using combinutions of the boundary integrul equation (BIE) and the normal derivative of the boundary integral equation (DBZE) is investigated for the two-dimensional Laplace equation. By using the combinations of these two equations it is possible to derive Fredholm integral equations of either the first or second kind regardless of the boundary conditions. Although the Fredholm equations of the second kind are well conditioned, in general they provide less uccurate interior results than the Fredholm equations of the first kind und the classical direct boundary element method (DBEM). On the other hand, the Fredholm equations of the first kind cun provide poor solutions in regions close to the boundary where large gradients exist in the boundury flux. Keywords: boundary element tion, hypersingular integral
method,
derivative
Introduction The direct boundary element method (DBEM) has proven to be a valuable tool in the analysis of a wide class of boundary value problems. The method has certain advantages over some of the popular domain methods in that only the boundary needs to be discretized, thus significantly reducing the order of the resulting algebraic equations. Further, in many cases the resulting algebraic equations are better conditioned for the DBEM than for the domain methods. The conditioning of the boundary element equations is strongly dependent on the Green’s function used in formulating the integral equation, the approximation functions used for the source densities, the geometry of the boundary, the boundary discretization, and the boundary conditions.1-3 If the boundary conditions are entirely of the Neumann type, then the integral equation is a Fredholm equation of the second kind, and the associated linear algebra statement is generally very well condiAddress reprint requests to Dr. Ingber at the Department of Mechanical Engineering, University of New Mexico, Albuquerque, NM 87131, USA. Received
536
5 September
Appl.
1989; accepted
Math. Modelling,
20 March
1990
1990, Vol. 14, October
boundary
integral
equation,
Fredholm
integral
equa-
tioned. If the boundary conditions are entirely of the Dirichlet type, then the integral equation is a Fredholm equation of the first kind, and the conditioning of the associated linear algebra statement is worse than for the Neumann problem. Of course, there is a spectrum of problems between these two extremes with mixed boundary conditions. In practice, it is sometimes necessary to discretize the domain in such a way that some of the collocation points lie close together in order to resolve a singularity in the potential or to accurately approximate highly curved geometries. The set of linear equations can become ill-conditioned for these types of problems, and the accompanying finite arithmetic errors can severely degrade the accuracy of the numerical solution. Motivated by the authors’ previous work in grid optimization,4-h the original intent of this research was to formulate a hybrid boundary element method using combinations of the boundary integral equation (BIE) and the normal derivative of the boundary integral equation (DBIE) that would minimize the conditioning of the associated linear equations. In this hybrid method, which we designate HBEMl, the BIE is written at collocation nodes along the boundary where Neumann boundary conditions are prescribed, and the DBIE is written at collocation nodes along the boundary where
0 1990 Butterworth-Heinemann
Combinations
of regular and derivative
Dirichlet boundary conditions are prescribed. Although the HBEMl enhanced the diagonal dominance of the resulting linear equations and consequently the matrix conditioning, in most example problems considered, the method yielded less accurate interior results than the DBEM. Because of this, an alternative approach to the HBEMl was conceived in which the BIE is written at collocation nodes along the boundary where Dirichlet boundary conditions are prescribed and the DBIE is written at collocation nodes along the boundary where Neumann conditions are prescribed. This alternative approach is designated HBEMZ. The HBEMl yields a system of Fredholm equations of the second kind, while the HBEM2 yields a system of Fredholm equations of the first kind. For mixed problems the DBEM yields a system of equations that is neither entirely of the first or second kind. Several mathematical studies have been performed to determine the properties of boundary element methods yielding Fredholm integral equations of the first or second kind.3.7-‘X Nevertheless, the mathematical theory is far from complete, especially for the more popular point collocation methods. Further, for the most part these studies have not considered problems with mixed boundary conditions. One of the purposes of this paper is to expose some possible misconceptions about Fredholm equations of the first and second kinds through numerical experiment for problems that are either singular or nearly singular. Since the DBIE is hypersingular and exists only at points where the source densities are Holder continuously differentiable,14 isoparametric Overhauser spline elements are employed in the HBEMl, the HBEM2, and the DBEM. The Overhauser spline elements possess properties that make them outperform some of the traditional Lagrangian elements.” In the following sections we formulate the hybrid methods, outline a method of evaluating the DBIE, and look at several example problems. We show the feasibility of combining the regular and derivative boundary integral equations.
Theoretical formulation Let 4(p) be harmonic in the domain R C R2. By following the standard arguments of the DBEM,16 4(p) can be represented in terms of a boundary integral as
rl(PM(P) =
~G(P, 4 --gy 4(q) 9
boundary
equations:
M. S. lngber and T. J. Rudolphi
fOif[#fl 27&E R/ln(p) = the internal angle between the two tangents on either side ofp if p E r We refer to equation (1) as the boundary integral equation (BIE). The normal derivative of the BIE can be written in terms of the hypersingular integral equation asI $(p)
P
=
aG(p, q) w(q)
d2G(P, q)
J[ r
an Pan4 "(Y)-T%T
Y
P
fi I
(2) where a/an, denotes the derivative in the direction of the outward normal to I’ at the field point p. We refer to equation (2) as the derivative boundary integral equation (DBIE). Since the DBIE exists only at points where the approximations for 4(p) and a4(p)ldn,, are Holder continuously differentiable,14 the geometry must be smooth, and hence the coefficient of p)ldn,, (unlike the coefficient q(p) in the BIE) is always 7~at these points. We discretize both the BIE (equation (1)) and the DBIE (equation (2)) by dividing the boundary I into boundary elements and approximate &q) and a+(q)/&, with polynomial interpolates. To maintain the necessary smoothness of 4(q) and a&)/an,and at the same time reduce geometric errors, isoparametric Overhauser spline elements are employed. The formulation of Ortiz et ~1.‘~ is implemented here with one modification. In general, the parametric cubic curve c(t) is the following linear blend of two parabolas, p(r) and q(s) (see Figure I): c(t) = (1 - t)p(r) + tq(s)
(3)
In this way the first derivative continuity is maintained at the nodes p2 and p3. The modification made here to Ortiz et ul.‘s formulation occurs at corners or kinks in the geometry. Since the appropriate smoothness for the normal derivative cannot be maintained at these points, semidiscontinuous elements are employed on either side of a corner. This requires the use of five different sets of shape functions for the five types of elements shown in Figure 2, Element type 3 splines together two continuous elements. Element types 2 and 4 spline together a continuous and a semidiscontinuous element (the difference being the relative lo-
(1) 41dr-(cd
a+(q) - G(P, q)r
where I is the boundary of a, G(p, q) is the function given by G(p, q) = log Ip - 41, a/an, the derivative in the direction of the outward to the boundary at the source point q, and defined by
Green’s denotes normal q(p) is
PZ
Figure 1.
Parametric
P,
cubic curve c(t) for the Overhauser
spline
element
Appl.
Math.
Modelling,
1990,
Vol.
14, October
537
Combinations
of regular and derivative
boundary
equations:
LLY 4 ox
1 w!
Figure 2.
tinuous approach, when using the same number of nodes as the discontinuous approach, yields a higher-order method. Within each boundary element, 4(q) and a+(s)/&, are approximated as
5
2
Definition
0
Geometric
X
Collocation
sketch of the Overhauser
Al. S. lngber and T. J. Rudolphi
nodes nodes
spline element ~(4W,
types
=
5
(4)
#Si(d
i=l
cation of the two elements). Finally, element types 1 and 5 are automatically continuous at their middle node and require no spline at all. The Overhauser shape functions for element types 2, 3, and 4 are given in the appendix. An alternative approach would be to use discontinuous boundary elements, which would require only one set of shape functions. Although the initial work associated with developing five sets of shape functions is tedious, the continuous and semidiscon-
T&P) =
2 5 I [ Wi(y,~4 - &‘S;(q)G(p,
r=l
i=l
,
(8)
[A&+%1= [&yMjI
where +j and 4; represent the values of c$(q) and d+(q)/&,, respectively, at the global nodes. For the HBEMl, upon rearrangement (assembly) of equation (8) into a coefficient matrix multiplying a vector of unknowns equal to a known (load) vector, the 7~from the left-hand side of equations (6) and (7) will lie along the main diagonal of the coefficient matrix, thus improving the diagonal dominance of the coefficient matrix. In this way the conditioning of the linear equations
dn,dn,
VP. V,G(p, q) = -V;G(p, Appl.
= ,g $‘S;(q)
Math.
(5)
where 4: and &’ represent the values of 4(q) and t@(q)ldn,, respectively, at local node i within the eth boundary element r,; S; represents the Overhauser spline shape functions; and ne is either 3 or 4, depending on the element type. Inserting equations (4) and (5) into equations (1) and (2) yields
q)] dl-
is minimized. On the other hand, for the HBEM2, upon assembly the 7~from the left-hand side of equations (6) and (7) is absorbed into the load vector. The numerical evaluation of the hypersingular integral in equation (2) is discussed in the next section. Numerical evaluation integral equation
of the derivative
Modelling,
fundamental
solution for the Laplace equation,
q) = 0
1990,
Vol. 14, October
boundary
The kernel contained in the DBIE (equation (2)) given by aG(p, q)ldn, is regular and can be integrated by using standard Gaussian quadrature formulas. However, the kernel contained in the DBIE given by a*G(p, q)ldn,dn, is hypersingular and cannot be integrated in the ordinary sense. These types of kernels have been encountered in boundary element formulations for acoustic problems*8.1’ and for crack-opening problems.20~21 Through the use of simple solutions to the governing differential equation, Rudolphi has shown how the hypersingular term can be converted to three completely regular integrals. Alternatively, Stallybrass23 converts the hypersingular term to a single Cauchy principal value integral, and we follow this approach here. First, we note the identity
(n, . n,>V, . V,G(p, q) + (np X n,). [VP X v,G(p, q)l - nq. v, X [n, X V,G(p, q)l
Since G(p, q) is the free-space
538
r,
c
where N is the total number of boundary elements. In the DBEM, equation (6) is collocated at all the nodes within the boundary elements. In the HBEMl, equation (6) is collocated at nodes within the boundary elements where Neumann conditions are specified, and equation (7) is collocated at nodes within the boundary elements where Dirichlet conditions are specified. In the HBEM2, equation (6) is collocated at nodes within the boundary elements where Dirichlet conditions are specified, and equation (7) is collocated at nodes within the boundary elements where Neumann conditions are specified. For all methods we generate a system of linear equations that can formally be written as
~*G(P,4 =
+I&) an,
(9)
we have (10)
Combinations VP x V,GCp,
4
=
of regular and derivative
- 0,
x V,Gtp,
d
boundary
equations:
M. S. lngber and T. J. Rudolphi (11)
= 0
Further, $@)V, x [n, x v,G(~,q)l Therefore, /
44d
using equations
= V, x [&dn,
9) dr = - j- 12,. v q lan
1‘
(12)
(9)-(12), we obtain
@G(P, an
x V,G(p,q)l - V,d4d x In, x V,G(p,q)l
+
v, x [Ndn,
x
V,G(p,
dl dr
I nq. V,+(q) x 1% x V,G(p, 41 dr
(13)
1‘
Upon application of the divergence theorem the first integral on the right-hand side of equation (13) is seen to be zero. Hence using a vector identity for the scalar triple product on the second integral on the right-hand side of equation (13), we obtain ~*G(P,
/ I-
4)
$47) an an P
Y
a-=I ’n, x
V,G(p,
411.
In, x V,#q)l
dr
(14)
I‘
In the derivation of equation (14) we have implicitly assumed that $(q) is Holder continuously differentiable in order to apply certain of the vector operations. Although the hypersingular integral can now be interpreted in the Cauchy principal value sense, the integral on the right-hand side of equation (14) is still cumbersome because it contains a tangential derivative of 4(q), which can be evaluated only by numerical differentiation. We employ a device of Meyer rt a1.19 to circumvent this problem. From the preceding development we note that I
s
&9)&/ . v, x [n, x V,G(p, y)] dr = - / [n, x Vq44q)l. In, x
V,G(p,
(15)
41 d-
I
If we allowed $(q) = 6(p) = constant to be taken as a special case, the right-hand side of equation (15) would be zero, implying that
I-
+(p)n, . V, x [n, x V,G(p,
0
(16)
Hence we can write a2G(p,
/ 1‘
an,an
4)
444)dr = /w(p)
Y
- d44)~, . v, x b,? x V,G(P, 41 dr
(17)
1
in which the right-hand side of the equation is integrable in the Cauchy principal value sense. We use this form in all our numerical computations. In an alternative fashion, one can arrive at equation (17) by subtracting the term 4(p) from 4(q) in the integrand on the left-hand side of equation (13) and then adding back the appropriate term with 4(p) outside the integral on the added-back term. Then taking, in equation (2), the special case of $(p) = constant, that is, the analog of the rigid body motion argument of elastostatics, one concludes that the added-back term must be zero and hence the result of equation (17).
Numerical examples We consider three numerical examples in this section. Two of the examples have smooth solutions, while one example has a singularity at a reentrant corner. However, even though one of the example problems has an analytic exact solution, its domain is an elongated
ellipse causing singularlike behavior toward the ends of the major axis. The geometry for the first example is the unit square. Boundary conditions are consistent with the exact solution given by 6(x, y) = cash (x) sin (y)
(18)
We consider the problem in which Dirichlet conditions are specified along the sides x = 0 and y = 0 and Neumann conditions are specified along the sides x = 1 and y = 1. The results generated by the DBEM, HBEMl, and HBEM2 using ten uniform Overhauser spline elements per side are shown and compared to the analytic solution in Table I. All three methods yield excellent results, and there is very little to distinguish between the methods. As a convergence test, the maximum absolute error within the domain as calculated from the results at 81 evenly spaced points in the interior versus the number of nodes (degrees of freedom) is shown in Figure 3 for the three methods. The DBEM and the HBEMl provide essentially the same maxi-
Appl. Math. Modelling,
1990, Vol. 14, October
539
Combinations
of regular and derivative
boundary
equations:
Table 1. Comparison of the calculatedand analyticvalues of the potential4 at selected points within the domain for examplel.Numbersequence: DBEM,HBEMl,HBEM2, analytic
M. S. lngber and T. J. Rudolphi
The geometry for the second example is given by the L-shaped region shown in Figure 4. In this example the boundary conditions are consistent with the exact solution given by
X Y
0.9
0.7
0.5
0.3
0.1
+(r, 0) = r2’3 sin (28/3)
0.1
0.3
0.5
0.7
0.9
0.7872 0.7872 0.7872 0.7872
0.8186 0.8186 0.8186 0.8188
0.8827 0.8827 0.8830 0.8833
0.9822 0.9823 0.9828 0.9832
1.1209 1.1209 1.1222 1.1226
0.6473 0.6473 0.6474 0.6474
0.6731 0.6731 0.6732 0.6734
0.7259 0.7259 0.7262 0.7264
0.8077 0.8078 0.8083 0.8086
0.9221 0.9223 0.9229 0.9232
0.4817 0.4816 0.4818 0.4818
0.5008 0.5009 0.5010 0.5011
0.5401 0.5402 0.5403 0.5406
0.6010 0.6012 0.6014 0.6018
0.6861 0.6864 0.6867 0.6871
0.2969 0.2968 0.2969 0.2970
0.3086 0.3087 0.3087 0.3089
0.3328 0.3329 0.3329 0.3332
0.3702 0.3706 0.3705 0.3709
0.4225 0.4231 0.4230 0.4235
0.1001 0.1001 0.1001 0.1003
0.1041 0.1041 0.1041 0.1044
0.1122 0.1124 0.1123 0.1126
0.1248 0.1251 0.1249 0.1253
0.1419 0.1429 0.1422 0.1431
Degrees
Figure 3. The maximum a function of the number
where r and 0 are the polar coordinates. Dirichlet boundary conditions are specified on the sides adjacent to the reentrant corner, and Neumann conditions are specified elsewhere. A singularity of order r-1’3 exists in the normal derivative at the origin. To resolve the singularity, we pack nodes close to the origin. In practice, this is sometimes done by using a grid optimization scheme. Here, we use the geometric stretching parameter A, given by the ratio of the lengths of two adjacent elements, on the two sides next to the reentrant corner and evenly spaced elements on the other sides. As an example, the case for A = 0.8 and 13 elements per side is shown in Figure 4. For this particular problem the results are more difficult to interpret. We focus attention on the case in which 20 Overhauser spline elements are placed on each of the six sides. The interior results for the DBEM, HBEMl, and HBEM2 in the first quadrant are shown in Table 2 for a stretching parameter of 0.9. This particular stretching parameter is chosen because it was determined to yield the most accurate results for the DBEM. Table 2 shows that the HBEM2 and the DBEM provided very good results. However, for this problem the results provided by the HBEMl are unacceptable. The maximum absolute error within the domain for the three methods as a function of the stretching parameter A as calculated from the results at 243 evenly spaced points in the interior is shown in Figure 5. The figure shows the general trend that appeared in the first example problem, which was that the HBEM2 pro-
of Freedom
absolute error within the domain of nodes for example 1
as
mum absolute errors, while the HBEM2 has a somewhat smaller maximum absolute error. Figure 3 indicates that the convergence rate is superlinear for all three methods. It is interesting to compare the condition numbers for the three methods. With 84 degrees of freedom the DBEM has a condition number of 153, the HBEMl has a condition number of 81, and the HBEM2 has a condition number of 4588.
540
(19)
Appl. Math. Modelling,
1990, Vol. 14, October
Figure 4.
The discretizedL-shaped domain for example 2
Combinations
of regular and derivative
Table 2. Comparison of the calculated and analytic values of the potential I#Jat selected points within the domain for example 2 computed with a stretchingfactorof A = 0.9. Number sequence: DBEM,HBEMl,HBEMZ,analytic
boundary
equations:
M. S. lngber and T. J. Rudolphi
Table 3. Comparison of the calculatedvalues, computed by the HBEMl with three different stretching factors,and analytic values of the potential4 at selected points within the domain for example 2. Number sequence: A = 0.7.A = 0.8.A = 0.9, analytic
X
Y 0.5
0.4
0.3
0.2
0.1
0.1
0.2
0.3
0.4
0.5
0.5064 0.5498 0.5059 0.5061
0.4721 0.5143 0.4716 0.4718
0.4429 0.4842 0.4424 0.4426
0.4181 0.4590 0.4176 0.4179
0.3971 0.4380 0.3966 0.3969
0.4286 0.4686 0.4281 0.4283
0.3937 0.4323 0.3933 0.3935
0.3653 0.4028 0.3649 0.3651
0.3422 0.3790 0.3418 0.3420
0.3232 0.3598 0.3228 0.3230
0.3436 0.3798 0.3432 0.3434
0.3088 0.3433 0.3085 0.3087
0.2824 0.3156 0.2821 0.2823
0.2621 0.2944 0.2618 0.2620
0.2461 0.2780 0.2458 0.2460
0.2480 0.2802 0.2478 0.2479
0.2155 0.2455 0.2153 0.2154
0.1936 0.2220 0.1934 0.1935
0.1779 0.2053 0.1778 0.1779
0.1661 0.1928 0.1659 0.1661
0.1357 0.1633 0.1356 0.1357
0.1121 0.1369 0.1119 0.1121
0.0988 0.1219 0.0987 0.0988
0.0900 0.1120 0.0899 0.0901
0.0837 0.1049 0.0836 0.0837
X V
0.1
0.2
0.3
0.4
0.5
0.5
0.4335 0.5059 0.5498 0.5061
0.4007 0.4715 0.5143 0.4718
0.3722 0.4421 0.4842 0.4426
0.3474 0.4172 0.4590 0.4179
0.3255 0.3960 0.4380 0.3969
0.4
0.3627 0.4284 0.4686 0.4283
0.3297 0.3934 0.4323 0.3935
0.3022 0.3648 0.4028 0.3651
0.2792 0.3415 0.3790 0.3420
0.2596 0.3223 0.3598 0.3230
0.3
0.2856 0.3438 0.3798 0.3434
0.2530 0.3089 0.3433 0.3087
0.2277 0.2823 0.3156 0.2823
0.2077 0.2617 0.2944 0.2620
0.1913 0.2455 0.2780 0.2460
0.2
0.1992 0.2488 0.2802 0.2479
0.1689 0.2160 0.2455 0.2154
0.1480 0.1937 0.2220 0.1935
0.1327 0.1778 0.2053 0.1779
0.1207 0.1657 0.1928 0.1661
0.1
0.0977 0.1372 0.1633 0.1357
0.0757 0.1129 0.1369 0.1121
0.0631 0.0992 0.1219 0.0988
0.0545 0.0901 0.1120 0.0901
0.0481 0.0834 0.1049 0.0837
Table 4. The condition numbers associated with the DBEM, HBEMl, and HBEMZ as a function of the stretchingparameter in examole 2 Stretching parameter, A 1.0
0.9 0.8 0.7 0.6
06
0.7
0.8
Stretching
Parameter
0.9
1
A
Figure 5. The maximum absolute error within the domain a function of the stretchingparameter for example 2
as
vided the most accurate results followed by the DBEM and then the HBEMl. However, in this example the spread in accuracies of the three methods is wider. The minimum absolute error for the HBEM2 occurred with no stretching, while the minimum absolute error for the HBEMl occurred for A = 0.8. It is interesting to note that in a very narrow range of the stretching parameter it is possible to obtain good solutions with the
Condition number DBEM 138 329 2229 28,889 6,235,200
HBEMl 79 23 63 390 2609
HBEM2 46,964 67,646 399,680 4,720,OOO 71,856,OOO
HBEMI. The interior results for the HBEMI in the first quadrant for X = 0.7, A = 0.8, and h = 0.9 are shown in Tublr 3. Only for A = 0.8 are the solutions acceptable. The condition numbers for the three methods as a function of the stretching parameter are shown in Table 4. As was expected, the HBEMI provided much smaller condition numbers than did the other two methods. Even with a very large condition number for the HBEM2 at a stretching of A = 0.6 the method still provided very good results. (The results were calculated by using single precision on a DECstation 3 100.) It is also interesting to see how well the various methods predict the normal derivative along the boundary near the singularity. These results are shown in Table 5. Although the HBEMl was better able to resolve the singularity in the normal derivative close to the singularity, this advantage was quickly lost at points away from the reentrant corner. The domain for the third example problem is given by the ellipse of aspect ratio 5 as shown in Figure 6.
APPl.
Math. Modelling,
1990, Vol. 14, October
541
Combinations
of regular and derivative
boundary
Table 5. Comparison of the calculated and analytic values of the flux a+lan at points along the boundary near the reentrant corner for examole 2
equations:
M. S. lngber and T. J. Rudolph
Table 7. Comparison of the calculated and analytic values of the potential 4 at selected interior points for example 3
@ian X
DBEM
Y
0.0051 0.0154 0.0325 0.0515 0.0726
0.0000 0.0000 0.0000 0.0000 0.0000
- 4.8235 -2.5818 -2.1468 -1.7912 - 1.6025
HBEMl
HBEM2
Analytic
-
-4.8184 - 2.5791 -2.1445 - 1.7892 - 1.6007
-
3.8665 2.6809 2.0898 1.7924 1.5985
elliptic domain
in the
3.9936 2.8147 2.2007 1.8898 1.6865
Figure 6. The portion of the discretized first quadrant for example 3
Table 6. Comparison of the calculated and analytic values of the flux d$ddnat points along the boundary in the first quadrant for example 3 adan X
Y
DBEM
HBEMI
5.0000 4.8496 4.6617 4.4267 4.1331 3.7660 3.3071 2.7335 2.0165 1.1203 0.0000
0.0000 0.2440 0.3616 0.4649 0.5628 0.6578 0.7500 0.8373 0.9151 0.9746 1 .oooo
35.2293 8.4780 6.5075 10.2323 8.7184 7.6508 6.1798 4.4945 2.5517 0.3313 - 2.1724
9.3677 13.5452 12.0922 10.7773 9.4217 7.9854 6.4198 4.6802 2.7242 0.4419 -2.0813
We assume Dirichlet boundary conditions with the exact solution given by cb(X,Y) = X2 - y2 + 2xy - 1
Analytic 10.0000
13.5534 12.2480 10.8868 9.4982 8.0345 6.4458 4.6852 2.7092 0.4836 - 2.0000
associated (20)
For this problem the HBEM2 is identical to the DBEM. A geometric stretching parameter is again chosen to place smaller elements near the ends of the major axis to resolve the large curvature and rapid variations in the normal derivative of #J in those regions. With stretching, the isoparametric Overhauser spline elements are able to approximate accurately the curved geometry. For the grid containing 40 isoparametric Overhauser spline elements as shown in Figure 6 the DBEM had a condition number of 174, while the HBEMI had a condition number of 5. The solutions for the flux (a@&) at points along the boundary in the first quadrant are shown in Tuble 6. The HBEMl is seen to yield uniformly more accurate results for the surface fluxes. The DBEM yields extremely poor solutions near the ends of the major axis with errors in the flux of more than 300%. The interior solutions for the two methods are compared to the analytic results in Tuble 7. The
542
Appl.
Math.
Modelling,
1990,
Vol. 14, October
X
Y
DBEM
HBEMI
Analytic
0.5444 0.5444 0.5444 0.5444 0.5444 0.5444 1.6333 1.6333 1.6333 1.6333 2.7222 2.7222 2.7222 2.7222 3.8111 3.8111 3.8111 3.8111 4.9000 4.9000 4.9000 4.9500
- 0.9000 - 0.7000 -0.1000 0.1000 0.7000 0.9000 - 0.5000 - 0.3000 0.3000 0.5000 - 0.7000 -0.1000 0.1000 0.7000 - 0.5000 - 0.3000 0.3000 0.5000 -0.1000 0.0000 0.1000 0.0000
- 2.3644 - 1.9574 - 0.8240 - 0.6050 - 0.4265 -0.5111 - 0.2252 0.5934 2.5643 3.0597 2.1060 5.8778 6.9727 9.7708 9.5121 11.2278 15.8306 17.1824 23.8131 25.0066 25.7550 25.9407
- 2.4852 - 2.0772 - 0.9454 - 0.7282 - 0.5586 - 0.6479 -0.3180 0.4947 2.4494 2.9401 2.0457 5.7810 6.8660 9.6426 9.4379 11.1174 15.6629 17.0160 22.0811 23.0331 23.8617 23.5199
- 2.4936 - 1.9558 - 0.8225 - 0.6047 -0.4314 - 0.5336 -0.2156 0.5978 2.5578 3.0511 2.1093 5.8560 6.9449 9.7316 9.4635 11.1479 15.7212 17.0857 22.0200 23.0100 23.9800 23.5025
HBEMl yields more accurate results near the ends of the major axis. It is in this region that the DBEM yielded very poor surface fluxes. However, despite the fact that the HBEMl yielded uniformly more accurate surface fluxes, away from the ends of the major axis the DBEM is more accurate. This result is surprising in light of the calculated surface fluxes and indicates that the HBEM2 (DBEM in this example) tends to smooth boundary errors as one moves into the interior analogous to a numerical St. Venant effect. In fact, the largest percent error in the interior at the selected points shown in Table 7 is 7.5% for the DBEM and over 20% for the HBEMl . Conclusions The feasibility of using combinations of the BIE and the DBIE has been demonstrated in the example problems. The three examples were chosen to show the relative merits of the classical DBEM compared to alternative methods for generating Fredholm integral equations of the first and second kinds for a smooth problem, for a singular problem, and for a “nearly” singular problem. For the smooth problem the two hybrid methods as well as the classical DBEM were shown to yield very good results, the results provided by the HBEM2 being slightly more accurate. For the singular problem with the reentrant corner the HBEMI provided poor solutions except for one value of the stretching parameter. Even in cases with extreme element stretching, the HBEM2 provided the most accurate results for the singular problem despite possessing very large condition numbers. The third example is perhaps the most perplexing. In this example the DBEM provided better solutions in a large portion of the domain despite yielding uni-
Combinations
of regular
and derivative
formly less accurate surface fluxes. Near the ends of the major axis the HBEMl was more accurate. In hindsight, one might argue that the DBIE should be written at nodes near the ends of the major axis and the BIE at all other nodes. Nevertheless, as a general statement, it appears that the HBEM2, which yields a system of Fredholm integral equations of the first kind, provides the most accurate interior solutions of the three methods. The general methodology discussed for generating systems of Fredholm integral equations of the first and second kinds is applicable to a much wider class of differential equations, and we assert our conclusions only for the two-dimensional Laplace equation, which generally does not yield excessive condition numbers because of the logarithmic fundamental solution. For other fundamental solutions it is quite possible that finite arithmetic errors and small errors in determining the boundary data could become critical in determining the overall accuracy of the solution.
boundary 14 15
I6 I7
18
19
20
21
22
References 1
2 3
4
5
6
7
8
Vable. M. Making the boundary element method less sensitive to changes or errors in the input data. Internar. J. Numer. Mefhods Engrg. 1987, 24, 1533-1540 Goldberg, M. A. Solution Methodsfor Integral Equations: Theory and Applications. Plenum Press, New York, 1978 Niessner, H. Significance of kernel singularities for the numerical solution of Fredholm integral equations. Bounda~ Elements IX, Vol. I, ed. C. A. Brebbia. Springer-Verlag, Berlin, 1987, pp. 213-227 Ingber, M. S. and Mitra, A. K. Grid optimization for the boundary element method. Internat. 1. Numer. Method Engrg. 1986, 23, 2121-2136 Mitra, A. K. and Ingber, M. S. The numerical solution of Stokes flow in a domain with re-entrant boundaries by the boundary element method. Internat. .I. Numer. Method Huids 1988, 8, 327-338 Ingber, M. S. and Mitra, A. K. Grid redistribution based on measurable error indicators for the direct boundary element method. Engrg. Analysis, in press Atkinson, K. E. A Survey of Numerical Methods for the Solution of Fredholm Integral Equations of the Second Kind. SIAM, Philadelphia, 1976 Giroire, J. and Nedelec, J. C. Numerical solution of an exterior Neumann problem using a double layer potential. Math. Comp. 1978, 32, 973-990 Hsiao, G. C., Koop, P. and Wendland, W. L. A Galerkin collocation method for some integral equations of the first kind. Computing 1980, 25, 89-130 Arnold, D. N. and Wendland, W. L. The convergence of spline collocation for strongly elliptic equations on curves. Numer. M&h. 1985, 47, 313-341 Wendland, W. L. and Christiansen, S. On the condition number of the influence matrix belonging to some first kind integral equations with logarithmic kernel. Appl. Anal. 1986.21, 175-183 Chandler, G. A. and Graham, I. G. Product integration-collocation methods for noncompact integral operator equations. Math. Comp. 1988, 50, 125-138 Costabel, M., Ervin, V. J. and Stephan, E. P. On the convergence of collocation methods for Symm’s integral equation on open curves. Math. Camp. 1988, 51, 167-179
23
equations:
M. S. lngber
and T. J. Rudolphi
Gtinter, N. M. Potential Theory and Its Applications to Basic Problems of Mathematical Physics. Ungar, New York, I967 Ortiz. J. C.. Walters. H. G. and Gipson. G. S. Development of Overhauser splines as boundary ‘elements. Boundary Elements IX, Vol. 1, ed. C. A. Brebbia. Springer-Verlag, Berlin, 1987, pp. 401-407 Brebbia, C. A., Telles, J. C. F. and Wrobel, L. C. Boandary EIement Techniques. Springer-Verlag, Berlin, 1984 Ingber, M. S. and Mitra, A. K. The evaluation of the normal derivative along the boundary in the direct boundary element method. Appl. Math. Modelling 1989, 13, 32-40 Burton, A. J. and Miller, G. F. The application of integral equation methods to the solution of some exterior boundaryvalue problems. Proc. Roy. Sot. London Ser. A 1971, 323, 201-210 Meyer, W. L., Bell, W. A.. Zinn, B. T. and Stallybrass, M. P. Boundary integral solutions ofthree dimensional acoustic radiation problems. J. Sound Vibr. 1978, .59(2), 245-262 Ioakimidis, N. I. A new singular integral equation for the classical crack problem in plane and antiplane elasticity. In/ernat. .I. Fracture- 1983, 21, 115-122 Budreck. D. E. and Achenbach. J. E. Scattering from threedimensional planar cracks by the boundary integral equation method. J. Appl. Mech. 1988, 55, 405-412 Rudolphi. T. J. The use of simple solutions in the regularization of hypersingular boundary integral equations. Comp. Math. Appl., in press Stallybrass. M. P. On a pointwise variational principle for the approximate solution of linear boundary value problems. J. Math. Mech. 1967. 16, 1247-1286
Appendix The shape functions for element types 2,3, and 4 (Fig14~ 2) are given below. Their derivations are straightforward using equation (3) to blend together the appropriate continuous and semidiscontinuous shape functions. Element type 2: N:(t) = -9P110 N#t)
+ 9P15 - 9t/lO
= 2t’ - 7t2/2 + t/2 + 1
A’s(t) = - St’/5 + 11t2/5 + 2tt5 N:(t)
= t3/2 - 1’12
Element type 3: N_:(t) = - t’/2 + t2 - t/2 N:(t)
= 3t3/2 - V/2
N:(t)
= -3t3/2
N:(t)
= t’l2 - t2/2
+ 1
+ 2t’ + t/2
Element type 4: N:(t)
= - t3/2 + t2 - t/2
N!?(t) = 8t3/5 - 13t2/5 + 1 Aqt)
= -2P
+ 5t2/2 + t/2
N;(t) = 9t3/10 - 9t2110
Appl.
Math. Modelling,
1990, Vol. 14, October
543