- Email: [email protected]
A new h-adaptive refinement scheme for the boundary element method using local reanalysis
N O R T H = HOt,LAND
A New
h-Adaptive
Boundary
Element
Refinement Method
Scheme Using
for the
Local
Reanalysis
A. C h a r a f i
School of Mathematical Studies University of Portsmouth Portsmouth POI 2EG, UK and L. C. W r o b e l
Department of Mechanical Engineering Brunel University Middlesex UB8 3PH, UK Transmitted by M. A. Golberg
ABSTRACT This paper presents a new way of estimating the discretization error in the boundary element method (BEM) based on the concept of local reanalysis (LR). The principle of LR is implemented for the particular case of two-dimensional potential problems, in an a t t e m p t to find a technique to calculate a cheaplycomputable predicted solution to use in the evaluation of the quality of the boundary element results. The error estimators are defined from the difference between the BEM and LR solutions, which measures the "sensitivity" of the boundary element results to a local change in the mesh. An h-adaptive method using LR is developed and different numerical experiments are carried out to test its performance. (~) Elsevier Science Inc., 1997
1.
INTRODUCTION
A n y a d a p t i v e s t r a t e g y has to answer t h e t h r e e basic questions of when, where, a n d how t o refine. R e f i n e m e n t will b e n e c e s s a r y w h e n t h e s o l u t i o n o b t a i n e d f r o m t h e n u m e r i c a l a n a l y s i s is n o t satisfactory. T h e a c c u r a c y of t h e s o l u t i o n m a y b e assessed b y use of e r r o r e s t i m a t o r s , b u t e n g i n e e r i n g c o n s i d e r a t i o n s m a y also b e used t o t e s t t h e r e l i a b i l i t y of t h e results. O n c e it h a s b e e n j u d g e d n e c e s s a r y to refine, error i n d i c a t o r s a r e u s e d t o show
APPLIED MATHEMATICS AND COMPUTATION82:239-271 (1997) (~) Elsevier Science Inc., 1997 655 Avenue of the Americas, New York, NY 10010
0096-3003/97/$17.00 PII S0096-3003(96)00028-8
240
A. CHARAFI AND L. C. WROBEL
where the refinement should occur. To the question of how to proceed with the refinement there are basically three schemes. • Nodal redistribution or r-type. • Element subdivision or h-type. • Local variation of the order of polynomial approximation or p-type. Hybrid methods, combining the merits (and demerits) of any of these schemes, are also possible. The hp-method, combining h- and p-adaptive schemes, is the most common. The solution of an adaptive problem is similar for all the various methods and obeys the following algorithm. 1. From an initial coarse mesh, a first approximation to the solution is computed. 2. The accuracy of the solution is assessed using suitable error estimators. 3. The mesh is refined with local error indicators guiding where the refinement should occur. 4. From the new mesh an improved solution is calculated. 5. The procedure (1) to (4) is repeated until the error estimator indicates that no further refinements are necessary or a specified stopping criteria is satisfied. In the h-method, the order of the polynomial approximation inside the element is kept constant, and the mesh is refined locally by subdividing the elements with large errors. This results in a graded mesh with a concentration of smaller elements near sources of errors. The first attempt on developing an h-version of BEM was by Rank [1] for a simple Galerkin boundary element method with constant elements. He presented the first theoretical results on a posteriori error estimation for the h-adaptive BEM, that were later extended to more general integral operators as well as to collocation BEM by Wendland and Yu [2]. Rank [1, 3], based on theoretical results of a posteriori error estimation, proposed the use of the H i - n o r m of the element residual to guide the hrefinement. The scheme was implemented for two-dimensional potential problems. Katz [4] developed an h-version of adaptive BEM to obtain the warping function of beams, for which the Galerkin method was employed to approximate the integral equation. Rencis and Mullen [5, 6] developed an h-adaptive boundary element for two-dimensional elasticity problems, in which they used an error indicator that predicts the number of uniform element subdivisions needed to obtain
Adaptive Refinement for Boundary Element Method
241
a pre-assigned accuracy. Although their approach was based on an oversimplification of an a priori error estimation result, and limited to constant elements, it nevertheless highlights some potentials of the h-refinement in improving the quality of the solution. Rencis and Jong [7, 8] extended the concept of Zienkiewicz and Zhu [9] for the F E M to develop a predicted error in a direct collocation BEM for plane and axisymmetric problems governed by the Laplace equation. Guiggiani [10, 11] has introduced non-residual-type error indicators by calculating the difference in the solutions obtained from two different sets of collocation points. He has applied his approach to two-dimensional elasticity problems with quadratic elements. Kamiya and Kawaguchi [12, 13] developed an h-version of BEM using the concept of "sample point error analysis" to solve two-dimensional potential and linear elastostatics problems. Abe [14, 15] proposed an error indicator and an evaluation of the nodal error, in the h-adaptive BEM, based on a new residue calculated at collocation points only. The nodal error was evaluated from the approximated residue by solving a matrix equation and the proposed strategy was applied to solve two-dimensional potential problems. The h-version has proved to be satisfactory especially when dealing with problems with singularities but the non-uniform element spacing obtained may lead to possible problems with the numerical integration and the conditioning of the system of equations, unless a sophisticated numerical integration scheme is used.
2.
T H E P R I N C I P L E O F LOCAL REANALYSIS
In this approach, starting from the assumption that a local change in the mesh will affect the solution only if the latter is not accurate enough, we introduce the concept of local reanalysis to measure the sensitivity of the solution once the mesh has been locally modified. In what follows, the concept of local reanalysis is explained. Consider the following boundary value problem in the two-dimensional domain ~ and its boundary F Au = 0 u = ~ q= ~
in on F~ on F~
(1)
where u and q = (Ou/On) are the potential and the flux, respectively, and n is the unit outward normal vector of F.
242
A. CHARAFI AND L. C. WROBEL
The boundary integral equation (BIE) of this two-dimensional potential problem for an arbitrary source point x i on the boundary is given by
f u ~ OE2 xi
c(xi)u(x') + Jr ( )-~n(
- ()d~ = ~r q(~)E2(x' - ~)d,~,
(2)
where E2 is the fundamental solution of the Laplace equation. After discretization of the boundary into L boundary elements, (2) becomes
c(x')u(x') +
u(()--~-~n ix - ()d( = e=l
e
~.q(()E2(x' - ()d~
(3)
e=l
After an approximation of the variation of u and q within each element using interpolation functions and nodal values, this equation is then applied for i = 1 , . . . , n, where n is the number of collocation points on the boundary F. The quality of the boundary element solution is strongly dependent on the discretization employed (number of elements and degree of polynomial approximation). It is this discretization error that is studied here, although there are other sources of error, such as round-off errors, integration errors, and the improper representation of curved boundaries a n d / o r of some boundary conditions. After a first approximation is obtained using the BEM, we need to test its quality. This is achieved by using LR over each element in turn, to decide whether it needs to be refined or not. This step consists of using the BEM, on the selected element for reanalysis, to solve the integral equation obtained when assuming that the boundary values on all the other elements are already known from a previous full BEM analysis. Suppose that an element Fe has been selected for the reanalysis. Equation (3) holds only at collocation points used to get the BEM solution. LR consists of using the values of (u, q) on the portion r - Fe from the initial solution, namely (~, q~, and solving the following BIE for the unknowns (~, q~ on element Fe
. OE2 ,(x - ~)d~ = b(x) + f c(x)~(x) + a(x) + f r ~(~)--0~n
~(~)E2(x - ~)d~ (4)
JFe
e
where a(=)
=
-
(5)
Adaptive Refinement for Boundary Element Method
243
f / #
\
/
x
/ /
\
I I
t
Original n o d a l
I
points
\
!
I I
/
\\ \ \
~/
, /'
Additional nodal
points
FIG. 1. Local discretization for local reanalysis. and
b(x) = E [ ~(~)E2(x - ~)d~ j~e 4r
(6)
To solve (4), the element Fe is subdivided into two subelements of the same type as those used to get the initial BEM solution, which leads to the discretized equation ne
ne
c , ~ + as + E H, kuk = b, + E k=l
aikq'k
(7)
k----1
where n~ represents the number of collocation points on F~ (if quadratic elements are used n~ = 5, see Figure 1), and G~k and H~k are influence coefficients obtained by integrating the fundamental solution and its normal derivative, respectively, over the two subelements. After applying the prescribed boundary conditions on F~ and rearranging the coefficients H~k and Gik accordingly, a linear system of equations is formed [AI{X} = {Y} + {b} - {a}
(8)
which is very small compared to the one obtained when doing the BEM analysis on the whole boundary, since the degrees of freedom are only those on the portion F~ of the boundary F. Hence, local reanalysis is relatively
244
A. CHARAFI AND L. C. WROBEL
r-l I
I I,
!/
n,
i
[A]
I I
',
I I
I
I
I
:
[×!/ /
/
Effect of additional collocation points
I .... Global Analysis Local Analysis
FIG. 2. Matrix illustrationof localreanalysis.
quick to perform. Figure 2 shows how the local system of equations for the local reanalysis is obtained from the global one. Computing a new solution using the principle of local reanalysis, and comparing it to the BEM solution, provides a quick feedback on the quality of the results, enabling to build an error estimator---substituting the exact solution by the one obtained from the local reanalysis--which provides a measure of how the solution is affected by a local change in the mesh. This idea has been first suggested by Casale [16] in his dissertation thesis. Its use as an error estimator has already proved reliable in the case of potential problems, as shown by Chara~i et al. [17]. With regards to the computational efficiency of the LR procedure, while it is true that, for the new collocation points, new integrals have to be evaluated on all boundary elements, most of these integrals can be reused in the next solution step if the new collocation points axe retained. Thus, very little additional work is wasted in building the error estimates, and the scheme works very efficiently. 3.
STUDY OF T H E LOCAL REANALYSIS SOLUTION
To have a good assessment of the solution obtained from the local reanalysis procedure, one has to minimize errors arising from the approximation of the geometry and boundary conditions, as well as those caused by the numerical approximation of the integrals involved in BEM. Straight-line
Adaptive Refinement for Boundary Element Method
245
FIG. 3. Possible boundary conditions at a corner for potential problems. Double node
, Collocation p o i n t Fro. 4. Definition of quarter-point elements. geometries are the simplest to treat. However, they introduce corner points which have to be appropriately dealt with. These numerical considerations are described first, before analyzing the local reanalysis procedure.
Numerical considerations Geometry and boundary conditions
As linear elements are used for the numerical t r e a t m e n t of the boundary integral equation, the examples are chosen so t h a t the geometry is correctly represented (straight lines) and the b o u n d a r y conditions exactly applied (linear variation at most).
Treatment of corners In the B E M solution of a model defined by straight lines, corners have to be appropriately dealt with. There are three types of possible boundary conditions in the vicinity of a corner as shown in Figure 3. Because of the continuity of the potential, the only unknown at the corner is the potential u in the first case, and the flux q (on the side where u is prescribed) in the second. But because the flux can be discontinuous, it is necessary in the third case to solve for the flux before and after the corner. This is achieved by use of non-nodal collocation elements presented by Marques [lS]. T h e basic idea behind this type of elements is to define two identical geometrical nodes at the corner point. To avoid obtaining a singular system of equations, the collocation points are located inside the elements at the quarter-points, as explained in Figure 4. Special integration scheme To minimize the effect of the error from numerical integration, and to focus on the discretization error generated by
246
A. CHARAFI AND L. C. WROBEL
the collocation method, a suitable integration scheme proposed by Telles [19] is used. It consists of an adaptive Gaussian quadrature scheme using a third degree polynomial transformation which bunches the integration points closer to the collocation points, and thus allows a smaller number of integration points for the same accuracy. Furthermore, Telles' transformation gives better results than if Gaussian quadrature is applied directly, and permits the computation of terms with singularities of order O(ln(r)). Numerical experiments Numerical experiments are carried out in order to compare the "predicted" solution obtained from the local reanalysis and the "next" boundary element solution obtained from a uniformly refined model. By "next" boundary element solution we mean the solution obtained when subdividing all the elements. The aim of this numerical investigation is to assess whether local reanalysis can (firstly) predict the "behavior" of the solution for the next step, and (secondly) provide information about areas of large errors.
EXAMPLE 1. A rectangle subject to Robin boundary conditions. Description
This example deals with a steady-state heat conduction in an exposed exterior column of rectangular cross-section, studied in detail by Wrobel [20]. A simplified model is presented in Figure 5, where one half of the column is studied, subjected to boundary conditions of the Robin type, q = au + b. The values of the coefficients a and b are given by a=-0.5 a=0 a a=-6
and and and and
b=50alongABC, b=0atD, b linear along CD and DE, b=0alongEFG.
C
D
E
|
I
l
I,
8
I
/,=36 FIG. 5. Initialdiscretizationfor Example I.
(9)
F =7
Adaptive Refinement for Boundary Element Method
247
84.0 75.2
~
u
~ Initial BEM s o l u t i o n 0 Predicted LR s o l u t i o n
• ~-~I~%
66.4
o_ 8, N. t
"
,olutlo.
57.8 48.8 .,..$ c 40.0 0 O.
31.2 22.4 13.6 4.8 -4.0
B
C
FIG. 6. Comparison of B E M
D
E
F
and L R solutions for Example I.
Results To study the behavior of the "predicted" solution obtained from local reanalysis, a comparison is made with the boundary element solution obtained from both the initial mesh described in Figure 5, and the "next" refined mesh. By "next" mesh, we mean the mesh obtained by subdividing into two all the elements of the initial one. The three solutions on side [BF] are presented in Figure 6. The results of the "predicted" solution from LR and the "next" solution from BEM are in good agreement, which suggests that the LR procedure can be thought of as a "quick" way of obtaining an approximation to the "next" BEM solution on the whole boundary. EXAMPLE 2. Heat conduction in a transformer coil.
Description In this example the heat conduction in a transformer coil [7] was considered as shown in Figure 7. The initial discretization employed 15 boundary elements. Results The BEM solutions obtained from the initial discretization (15 elements) of Figure 7 and from the mesh obtained by subdividing into two all the elements of the initial mesh, as well as the solution obtained from the LR procedure applied to the initial mesh, are presented in Figures 8-9.
248
A. C H A R A F I
L
u=0
A N D L. C. W R O B E L
K
!-
l~-V
u=0
~0
q=5~---
'~
"I
.
I
D I q=~-~_ C
AT
u--O
l-
'
i
s
I
B
.~J
J
J
FIG. 7. Geometry and boundary conditions for Example 2.
6O 54
o Initial BEM solution I t o Predicted LR solution o Next BEM solution J
48 42 v
36 30 O
24 18 12 6 0
B
C
D
E
F
-
-
•
s
_
G
H
I
J
K
FIG. 8. T h e potential u from B E M and LB. solutions for E x a m p l e 2.
Adaptive Refinement/or Boundary Element Method
249
1
--3
j
K
L
a
B
FIG. 9. T h e flux q from B E M a n d L K solutions for E x a m p l e 2.
When solving for the potential u (smooth solution), see Figure 8, the "predicted" solution from LR does not approximate the "next" BEM solution as well as in the previous example, especially between the points E and H. However, the LR procedure has correctly predicted the "behavior" of the solution, and this should be sufficient to indicate the need to refine this area. When solving for the flux q, see Figure 9, the solution obtained from LR still predicts the "behavior" of the BEM solution , and the difference between the "predicted" and "next" solutions are more apparent in the neighbourhood of points K and B. EXAMPLE 3. Singular flux problem
Description This example is a square domain ABCD of length 10 and with boundary conditions specified by u = 0 u = 10 q= 0
along AB at C and linear along BC along CDA.
The flux q is singular at the point C.
(10)
250
A. CHARAFI AND L. C. WROBEL 10
9 8
i
\,.
7
~
6
0
o Initial BEN s o l u t i o n "
~
~
Predicted LR solution
o
NextBEg
solution
\
4
3
~
1
~
0
C
~...~
s
.
D
A
FIQ. 10. The potential u from BEM and LR solutions for Example 3.
Results The difference in Figures 10-11 between the initial BEM solution on one hand, and the next BEM as well as the"predicted" LR solutions on the other hand, strongly indicates the need for more refinement of the mesh on elements [Be], [CO] and IDA]. Moreover, Figure 11 shows that the LR has predicted the discontinuity of the flux q at point B and its singularity at point C (side [BC]).
Conclusions from the numerical experiments It is clear from the study of the previous examples, and the analysis of the corresponding results, that the LR solution succeeds in giving an indication of the behavior of the BEM solution. It also has the ability to predict the existence of singularities. It should not, however, be used as an attempt to improve the solution itself as the LR measures mainly how "sensitive" the BEM solution is to a local change in the mesh. 4.
AN h-ADAPTIVE T E C H N I Q U E W I T H LOCAL REANALYSIS
After performing the local reanalysis over each element of the boundary F we get a solution ~ and ~', considered as the predicted solution. The exact
Adaptive Refinement for Boundar~d Element Method
251
3.0 2.5
1~ InitlalBEM solution o Predicted LR~mlution
.
o NextsE~~ i u u o n
Z.O
2
//~/-
/
1,5
o"
/
1.0
f
v
>¢
0.5 _
Z
0.0
-0.5 -1.0 -1.5 -20
"'0
A
B
C
FIG. 11. The flux q from BEM and LR solutions for Example 3. error for the potential (resp. flux) is approximated by e~ ~- ~ - ~
(resp. eq ~_ q - q~).
(11)
T h e basic (although empirical) idea is t h a t the difference between these two solutions, the one obtained from the boundary element discretization and the one obtained from the local reanalysis, provides some information about where the refinement should occur. T h e "error" measured is more a kind of "variation" or "sensitivity" of the solution due to a local change in the mesh. It should be noted t h a t the solution from L R cannot be regarded as more accurate t h a t the solution from BEM. T h e basic idea is t h a t the two solutions will be "more different" in those parts of the boundary where the discretization is poorer. g~'For no77r~8
To obtain a measure of the error (or rather of its estimate), we found it convenient to compare its H ° - n o r m on each b o u n d a r y element. We define the H°-local norms for the potential u and its associated error e~ on an element i as
1'' i
/12/ i
252
A. CHARAFI AND L. C. WROBEL
J
//// --(
1
! I
3
2
FIG. 12. Representation of the solutions ~ and ~ on a linear element. Because linear elements are used, the quantities [[u[[i and [[eu[[i can be evaluated analytically. On an element Fi the solutions ~ and ~ are represented as shown in Figure 12. E v a l u a t i o n o f [[ulli We have
~(~/)=
~'I+
- -
~2
-1_
(13)
therefore, r
Ilull2 = - ../r~u2(x)dx
(14)
(15) l'[ (~2- ~1)~ ] - 4 § + (~' + ~2)2 '
(16)
where li represents the length of element i. Evaluation of We have
[[e.lh
{ (1
- ~ 1 + (1 + ~)~3
~(~)
=
-
~)~3 + ~ 2
-l~y
(17)
-1<~0 0<~<1,
(18)
therefore, e~'(~) = { air/2 + bl~ + Cl a2~7 2 -b b2~] -}- c2
AdaptiveRefinementfor BoundaryElementMethod
253
where
al :
[U2 -- E1 ]2 T -~-~1 -- u3
bl : 2 [ ~ - F U l - U Z ]
[ ~1+~2 2
~3]
(20)
u2
(21)
a2---- - - - - 7 - +
u3 -- u2
(22)
b2= 2 [^~
^ "4-~2 +u3 -u2] [Ul 2
C1 =
ul +2 u2
(19)
and
C2 :
u 1 +2~ 2
~3]
(23)
u3
(24)
The value of the local error norm is then given by
l~ [(al +a2) Ile~,ll~ = ~" 3
+ (~-~)
+(Cl +C2)]
(25)
The H°-global norms for the potential and its associated error for a discretization of L elements are derived from the local ones by
[ L
~1/2
[ L
,~1/2
The H°-local and global norms and error norms for the flux are defined in exactly the same way as for the potential.
Error indicators As the error norms defined have little physical significance and can be misleading, we introduce nondimensional error percentage values. Following the work by Zienkiewicz and Zhu [9], we define the global error percentage for both the potential and the flux respectively, by
1 •,-~, = k, ll~.ll 2 + ile,~ll2 ]
x loo
and
Tq =
i1~12 + ileqll ~
× lOO (27)
where [[~[12 + [[e~l[ 2 and [[q~[2 + Heq[[2 approximate the exact potential norm Ilull2 and the exact flux norm Ilqll2, respectively.
254
A. CHARAFI AND L. C. WROBEL
Assuming a uniform distribution of the error on the boundary F, the relative error percentage on the element Fi for both the potential and the flux respectively, is defined by 1
ii~--~i+ ]le~]12]
1
x 100
and
r 4 = ~H~I2 +
11%112] × 100.
(28)
Computational strategy After solving the problem with the initial data (geometry and boundary conditions), the elements are analyzed one by one. O n each element, the percentage error obtained from the difference between the LR and BEM solutions is computed to decide whether to keep the refinement or not. The new mesh is then defined, refining the elements where the error is not within the imposed tolerance, and the problem is solved as if it is a new one. The self-adaptive mesh refinement strategy is controlled by specifying for both the potential and the flux on all elements that ~_<~
and
Tq<~
(29)
where ~ is a user-specified tolerance. It is important to notice that ~ is a local requirement. The refinement procedure can be terminated by different means. • By specifying an error bound for the global error percentage of both the potential and the flux, if one is interested in the global accuracy. However, it should be kept in mind that a good overall accuracy could well mask large errors in small parts of the boundary. • By storage considerations, that is, by specifying the maximum number of degrees of freedom the user is prepared to include in the analysis. • By time limitation considerations specifying, for example, a maximum number of resolutions. The details of the programming procedure of the error estimation and the adaptive mesh refinement strategy are given in Figure 13.
Numerical examples The h-adaptive BEM with LR is tested on five examples. The first three examples are those presented in Section 3, the fourth example is a sheet pile problem with a singularity and the fifth example is an L-shape domain along which the solution is smooth. In each case, the final potential and flux solutions obtained from the adaptive procedure are compared to a "reference" solution obtained by solving the same problem with a uniform
Adaptive Refinement for Boundary Element Method
255
[ Define geometry and boundary conditions
[
•] Perform standard B E M analysis i
l Compute the norms for u and q : [[~[[2 and
I1~112
L
] Loop on the elements ]
Local reanalysis : , Subdivide element and compute new H and G for new nodes , Form small system of equations AX = B , Solve the small system of equations , Compute the local norms of the errors eu and eq : Ile.ll~ and Ile, ll~
I End of loop] l Compute the global error norms :
I1=,,11= = r.~ Ile,,ll~
II~qll 2 = ~i Ileqll~
and
1
lu Compute the global errors percentage r,, and rq :
= ( ~ ) ~
× ~oo , ~ = ( ~ ) ~
× ~oo
l TGet relativepercentage errorson the elements :
[
•
{
r~=
~
½
Xll,~ll%lle,,ll~ /
~ ~nn
. . . . .
r ~-
q -
f
~
½
~l1411,+ll~dl, I
x 100
1 [ For aU elements r ~ < ~ and r ~ _ _ f ? J Yes I No 1
~~or
!
elementsi such that
: r~ > ~ or ¢~ > keep the subdivision
FIQ. 13. Flow chart of/z-adaptive BEM with local reanalysis.
256
A. CHARAFI AND L. C. WROBEL
mesh of a large number of elements enough to capture the behavior of the exact solution. EXAMPLE 1. This is the same Example 1 defined in Section 3. The model was analyzed starting from the coarse initial mesh of 6 elements, shown in Figure 5, which is sufficient to describe the geometry and boundary conditions. The adaptive scheme was tested with a specified local tolerance ~ = 0.04%. The adaptive solution was compared to a "reference" solution obtained by solving the same problem with a uniform mesh of 100 elements. Table 1 shows the decreasing values of the percentage error for the potential through the different iterations. The final mesh presented in Figure 14 was obtained with 50 elements. Figure 15 shows t h a t the distribution of the nodes in the final mesh is in good agreement with the "reference" solution. The same example was reanalyzed using a uniform mesh of 50 elements, for which the global percentage error was T~ = 0.0651% and the m a x i m u m local percentage error was z max = 0.2545%. With the adaptive procedure, the value of the global percentage error has decreased to T~ = 0.0248%, and the m a x i m u m local percentage error has decreased to Tmax = 0.0126%. TABLE 1 GLOBAL THE
PERCENTAGE
ADAPTIVE
Iteration
ERROP~S
PROCESS
FOR
THROUGH EXAMPLE
1
N u m b e r of e l e m e n t s
~'u [%]
6 12 22 38 48 50
3.2296 0.7753 0.2075 0.0580 0.0267 0.0248
0 1 2 3 4 5
C
D
E -v
A -
vA
A I ......................T....... T..... G FIG. 14. F i n a l a d a p t i v e m e s h for E x a m p l e 1.
Adaptive Refinement for Boundary Element Method
257
80 72 • Adaptive s o l u t i o n J - - Reference s o l u t i o n J
64 56 48 .,.d
Q) O
4O 32 24 16 8
0 B
C
D
E
F
Boundary length FIG. 15. F i n a l solution for E x a m p l e 1.
EXAMPLE 2. The geometry and boundary conditions for this example are the same as in Example 2 in Section 3. The adaptive scheme was started with the minimum number of elements sufficient to describe the geometry and boundary conditions, that is, 12 elements. The local tolerance imposed was ~ - 0.25% and the values of the global percentage error for both the potential and the flux, for each iteration, are presented in Figure 16, which indicates the convergence of the adaptive process. The final mesh from the present adaptive procedure, presented in Figure 17, was obtained with the global percentage tolerances ~'u = 0.2285% and rq = 0.1386%. The adaptive solution, for both the potential and the flux, is compared to a "reference" solution obtained from a quasi-uniform mesh of 216 elements. By quasi-uniform is meant a mesh in which the ratio of two adjacent elements is equal to either 1 or 1.5. Figure 18 shows that the two potential solutions are in good agreement, with the distribution of nodes capturing the curvature of the solution as well as the convex-concave changes in the curve. The flux solution obtained from the final adaptive mesh in Figure 19 is also in good agreement with the "reference" solution. Solving the same problem using linear elements and an adaptive scheme with the "sample point error analysis", Kamiya and Kawaguchi [12] obtained the mesh presented in Figure 20. It is not always easy to compare
258
A. C H A R A F I A N D L. C. W R O B E L
14.0 12.8 I o Global potential error Tu I I ~ Global flux error Tq ]
It.2 O I. t~
9.8
@ 8.4 U0 c- 7.0 5,6 4.2 O
2.6 1.4 0.0
10
19
28
37
48
Number
55
64
73
82
of e l e m e n t s
FIG. 16. Convergence graph for Example 2.
L.
K
G N u m b e r of e l e m e n t s : E
F
C B FIG. 17. Final adaptive mesh for Example 2.
9t
100
Adaptive Refinement for Boundary Element Method
259
60 54 48
42
"~
30
0
24
O.,
18
12
B
C
D
E
F
G
H
I
J
K
Boundary length FiG. 18. Final potential solution for Example 2. two different adaptive schemes as their respective error indicators and estimators may measure different entities (residual error, interpolation-type error, sensitivity, etc...). However, one can always compare the final meshes with respect to a "reference" solution obtained from a very fine mesh, or indeed the exact solution when it is available. ComParing the mesh obtained from the present adaptive scheme with t h a t obtained by Kamiya and Kawaguchi [12] in Figure 20, one notices that the present scheme has produced a graded mesh towards K along the side [KL], contrary to the mesh produced by Kamiya and Kawaguchi [12]. A close look at the "reference" solution in Figure 19 shows the strong variation of the flux in the neighbourhood of K (side [KL]). Overall, the final meshes obtained by the two methods are qualitatively similar. EXAMPLE 3. This is the same Example 3 defined in Section 3. To model this problem, all that is needed is to define the 4 nodes of the square. The adaptive procedure automatically increases the n u m f e r of elements until the specified tolerance of ~ = 0.125% has been achieved or a stopping criteria reached. Table 2 shows the values of the percentage errors for both the potential and the flux for each iteration. Note the oscillation of the percentage error
260
A. C H A R A F I A N D L. C. W R O B E L
0 -I -2
p~ v ,--4
r-~
-3 -3 -4 -5 -6 -7 -8
K
L
I
L
A Boundary
B
length
FIG. 19. Final flux solution for Example 2.
L
K
H! . . . . . . . . . .
Number of elements: 96
_G
I F
D t ...... ......
A
e e 2 s s ~ ; : s ; : : s s
t
C
B
FIG. 20. Final adaptive mesh for Example 2 by Kamiya and Kawaguchi.
Adaptive Refinement for Boundary Element Method
261
TABLE 2 GLOBAL PERCENTAGE ERRORS THROUGH THE ADAPTIVE PROCESS FOR EXAMPLE 3 Iteration 0 1 2 3 4 5 6
N u m b e r of elements
~u [%]
Tq [%]
4 8 16 30 50 63 75
15.4175 5.2777 1.8933 0.6739 0.2395 0.1007 0.0740
27.2363 22.0559 34.4350 10.7786 7.6949 5.4770 3.8964
N u m b e r of e l e m e n t s : 75
D
A
FIG. 21. Final adaptive mesh for E x a m p l e 3.
for the flux at iteration 2, due to the high value of the local percentage error for the flux near the point C (side [PC]), caused by the singularity of the flux solution at that point. The final mesh is represented in Figure 21. The final solution from the adaptive scheme is compared to the "reference" solution obtained from a uniform mesh of 320 elements. Figure 22 shows that the nodal partition captures very well the behavior of the potential solution with more nodes towards C (side [CD]) and less towards A (side [DA]). The fluxes from the adaptive and the "reference" solutions are in good agreement as shown in Figure 23. The singularity at the point C (side [BC]) is very well captured by the graded mesh towards point C.
262
A. C H A R A F I A N D L. C. W R O B E L 10
\
[* Adaptivesolution ] ~. ReferencesolutionI
0
i
D Boundary
length
FIG. 22. Final potential solution for Example 3. 5.00
4.34
* Adaptive solution m Reference solution
3.68 3.02
j
,---- 2.36 Ce
1.70 r,..
1.04
0.38 f
-0.28 -0.94 -
I
1.60
A
B Boundary
length
FIG. 23. Final flux solution for Example 3.
Adaptive Refinement for Boundary Element Method
A
263
E D,
15
C
B L I-
18
-I
FIG. 24. Initial model for Example 4. EXAMPLE 4. This example is a square A B C D E representing a sheet pile problem, see Figure 24. The boundary conditions are given by q u q u
= = -=
0 0 0 157
along along along along
ABC, CD, DE, EA.
(30)
This example is particularly interesting because it presents a singularity at the end of the cut-wall DE (the flux is singular at point D, side [CD]). The tolerance specified was Y -- 0.25% and the results obtained are shown in Figures 25-27. Note t h a t the final mesh, with 51 linear elements, is more refined (as expected) in the vicinity of point D due to the singular behavior of the flux. As the error indicators are obtained by computing the integral of the difference between the BEM and LR solutions over the element considered, their value is expected to be high near the singularity. This is what caused the oscillation at the beginning of the adaptive process for the global percentage error Tq, see Table 3. This influence, however, decreases as the length of the element near the singularity gets smaller. T h e high value of Tq is due to the large influence of the local error in the vicinity of point D (on the other parts of the boundary the local errors are very small), and also to the increasing gap between the predicted solution from the local reanalysis and the B E M solution at point D, as shown in Figure 28.
264
A. CHARAFI AND L. C. WROBEL
E
A
D
Number
=
=
of elements:
=
=
|
51
i
C
B FIG. 25. F i n a l a d a p t i v e m e s h for E x a m p l e 4.
162.0
145.8 129.6 113.4 v
97.2 81.0 q) 64.8 O 48.6
\
32.4 16.2 0.0
A
8 B
"4 C
Boundary length FIG. 26. Final potential solution for Example 4.
D
Adaptive Refinement for Boundary Element Method
265
TABLE 3 G L O B A L P E R C E N T A G E E R R O R S T H R O U G H THE A D A P T I V E PROCESS FOR E X A M P L E 4
Iteration
Number of elements
r u [%]
Tq [%]
5 10 19 33 41 47 51
6.7605 2.1247 0.6175 0.2605 0.1967 0.1775 0.1724
42.9869 35.6410 37.1137 32.1726 30.5618 29.8507 29.1731
0 1 2 3 4 5 6
400
36
380
320
33 • --
Adaptive solution Reference So ution
30
280
27
240
24
~ 200 X
x
mo
r,-- Is
120
15
80
12
40
9 g
D
E
A
FIG. 27. Final flux solution for E x a m p l e 4.
The oscillation of the flux solution near point D (side [CD]), is a numerical phenomenon that occurs sometimes near singularities when using B E M with linear elements. It is not linked to the adaptive process, as both adaptive and "reference" solutions show the same behavior, see Figure 27. This oscillation did not occur in Example 3, which has a singularity at point C (side [BC]), leading to think that it could be linked to the order of the singularity. EXAMPLE 5. This example is an L-shape domain with the boundary conditions specified in Figure 29. The "reference" solution was obtained from a uniform mesh of 320 elements for which the global tolerances were Tu = 0.029% and ~'a = 0.0008%. The adaptive procedure was tested with a specified tolerance of Y = 0.1%, and starting with an initial mesh of 6 elements which are sufficient to describe the geometry and boundary conditions.
266
A. CHARAFI AND L. C. WROBEI u
J'
',
'
"
,
"
.
' •
.
.
7
'
.
,
•
I
I
",
I
'
: BEMsolution I P~edicted solution
-60 -120 ~ C~ -160 ~J
-240
x ~=~
,~ 0
-300
\
-360
-420
2
~ -460 -540 -600
2
2
3
3
4
5
5
6
6
N um ber of r e s o l u t i o n s
FIG. 28. Flux variation at the point D for different resolutions for Example 4.
E
D
u=10
•
q=O
._!___
C i
u=O
5
q=0
F
•
/
# I
..... I_ i
/
, i i
, I
/
/
'
," -" / I , , "
," '
i I.,£2
10
--',
FIG. 29. Geometry and boundary conditions for Example 5.
7
Adaptive Refinement .for Boundary Element Method
267
TABLE 4 GLOBAL PERCENTAGE ERRORS THROUGH THE ADAPTIVE PROCESS FOR EXAMPLE 5
Iteration 0 1 2 3 4
Number of elements
Tu [%]
~'q [%]
6 12 24 40 60
8.5704 2.2188 0.6352 0.2218 0.0993
11.3012 4.6408 0.3436 0.0938 0.0759
C
B
N u m b e r of e l e m e n t s : 60
F
h FIG. 30. Final adaptive mesh for Example 5.
T a b l e 4 shows t h e decreasing values of t h e global p e r c e n t a g e errors, for b o t h t h e p o t e n t i a l a n d t h e flux, t h r o u g h t h e a d a p t i v e process. F i g u r e 30 shows t h e final n o d a l d i s t r i b u t i o n . N o t e t h e g r a d e d m e s h t o w a r d s t h e corners C a n d F which reflects t h e b e h a v i o r of t h e p o t e n t i a l "reference" sol u t i o n as p r e s e n t e d in F i g u r e 31. T h e m a x i m u m local e r r o r p e r c e n t a g e for t h e p o t e n t i a l was ~max __ 0.4506% for t h e a d a p t i v e s o l u t i o n c o m p a r e d t o 7umax = 0.3591% for t h e "reference" solution. F i g u r e 32 shows t h a t t h e values of t h e flux v a r y in a n o n l i n e a r w a y b u t w i t h i n a v e r y s m a l l interval (form 0.75 t o 0.82 in a b s o l u t e values), a n d t h e a d a p t i v e p r o c e d u r e has c a p t u r e d t h i s b e h a v i o r . F i g u r e 32 also i n d i c a t e s t h a t t h e r e is still need for refinement t o w a r d s t h e ends of s e g m e n t s lAB] a n d [DE]. I n d e e d , t h e m a x i m u m local t o l e r a n c e for t h e flux was Tq ax = 0.2698% for t h e a d a p t i v e s o l u t i o n c o m p a r e d t o T~nax = 0.0062% for t h e "reference" solution.
268
A. CHARAFI AND L. C. WROBEL
10
9
I " Adaptive so]ution solution I
7 6 v
5
0J o
4
c~
3 2 I.
0
E
F
A
Boundary
B
length
FIG. 31. Final potential solution for Example 5.
0.02
-0.74, -0.78
0.81
Io Adaptivesolution ]
0.80
-0.76
0.80
-0.76
3
,,~
-0.77
0+719 0.78
-0,70 r... -0.79
0.77
-0.80
0.78
-0.80
0 .?mS
-0.81
0.75 0.74
-0.82
B
D
FIQ. 32. Final flux solution for Example 5.
C
D
Adaptive Refinement for Boundary Element Method
5.
269
CONCLUSION
A new adaptive BEM scheme has been developed, using error estimates based on the concept of LR, for the construction of near optimal computational models in two dimensions. The self-adaptive boundary element technique was developed based on an estimation of the discretization error in each element using the concept of LR. It consists of measuring the sensitivity of the boundary element solution to local changes in the mesh by comparing the solution obtained from the boundary element discretization with that obtained by solving a series of small linear systems of equations--3 x 3 for potential problems with linear elements. Numerical experiments with two-dimensional potential problems have shown that, although the LR solution is not always an improvement of the BEM solution (as can be seen in Section 3), it nevertheless indicates the behavior of the BEM solution and has the ability to predict the existence of singularities. Error indicators derived from the LR procedure measure how sensitive the numerical solution is to a local refinement in the mesh, yielding a mechanism to determine if the overall mesh is sufficiently accurate or if mesh refinement is needed. This was implemented in an h-version with linear elements for two-dimensional potential problems and has produced promising results by providing significant information about where to refine the mesh. It should be pointed out that the cost of the adaptive boundary element strategy described in this work is not always cheaper than the cost of other adaptive strategies available. However, while in many adaptive schemes-especially those based on residual-type error estimators--the auxiliary data computed to evaluate the error indicators and estimators is lost, all the new integrals computed for the evaluation of the LR solution can be re-used in the next solution step if the new collocation points are retained. Therefore, very little additional work is wasted in building the error estimates, which makes the scheme very efficient. REFERENCES 1 E. Rank, A posteriori error estimates and adaptive refinement for some boundary integral element methods, in Proc. Int. Conf. on Accuracy Estimates and Adaptive Refinements in Finite Elements Computations, Lisbon, 1984. 2 W . L . Wendland and De-hao Yu, Adaptive boundary element methods for strongly elliptic integral equations, Num. Math. 53:539-558 (1988).
270 3
4
5
6
7
8
9
10
11
A. CHARAFI AND L. C. WROBEL E. Rank, Adaptivity and accuracy estimation for finite element and boundary integral element methods, in Accuracy Estimates and Adaptive Refinements in Finite Element Computations, John Wiley & Sons, pp. 79-94, 1986. C. Katz, Self-adaptive boundary elements for the shear stress in beams, in Proc. 2nd BEM Conf., Computational Mechanics Publications, 1986, pp. 759-769. J . J . Rencis and R. L. Mullen, Solution of elasticity problems by a selfadaptive mesh refinement technique for boundary element computation, Int. J. Num. Meth. Eng. 23:1509-1527 (1986). J . J . Rencis and R. L. Mullen, A self-adaptive mesh refinement technique for boundary element solution of the Laplace equation, Comp. Mech. 3:309-319 (1988). J . J . Rencis and K. Y. Jong, A self adaptive h-refinement technique for the boundary element method, Comp. Meth. Appl. Mech. Eng. 73:295-316 (1989). J . J . Rencis and K. Y. Jong, An accuracy post-processor for boundary element analysis, in Proc. lOth BEM Conf., Computational Mechanics Publications, 1988, Vol. 1, pp. 187-203. O. C. Zienkiewicz and J. Z. Zhu, A simple error estimator and adaptive procedure for practical engineering analysis, Int. J. Num. Meth. Eng. 24:337357 (1987). M. Guiggiani, Non-residual type accuracy estimates and adaptive refinement in boundary element analysis, in Boundary Element Techniques: Applications in Engineering, Proc. 4th BETECH Conf., Computational Mechanics Publications, 1989, pp. 323-340. M. Guiggiani, Error indicators for adaptive mesh refinement in the boundary element method: A new approach, Int. J. Num. Meth. Eng. 29:1247-1269
(1990). 12 N. Kamiya and K. Kawaguchi, An attempt of boundary elements, in Proc. 12th Int. BEM Conf., Computational Mechanics Publications and SpringerVerlag 2:527-538 (1990). 13 N. Kamiya and K. Kawaguchi, Error analysis and adaptive refinement of boundary elements, in Reliability and Robustness of Engineering Software II, C. A. Brebbia and A. J. Ferrante (eds.), Computational Mechanics Publications and Elsevier Applied Science, 1991, pp. 187-201. 14 K. Abe, h-Adaptive boundary element method and evaluation of error, in Proc. IABEM Symposium, Springer-Verlag, 1993, pp. 1-10. 15 K. Abe, A new residue and nodal error evaluation in h-adaptive boundary element, Adv. Eng. Soft. 15:231-247 (1992). 16 M. Casale, The Integration of Geometric Modeling and Structural Analysis using Trimmed Patches. Ph.D. Thesis, University of California, Irvine, 1989. 17 A. Charafi, L. C. Wrobel, and R. Adey, An approach to h-adaptive boundary element methods using local reanalysis, in Proc. 7th BETECH Conf., Computational Mechanics Publications and Elsevier Applied Science, 1992, pp. 905-918.
Adaptive Refinement for Boundary Element Method
271
18 E. Marques, Coupling of the Finite Element and the Boundary Element Method: An Application to Potential Problems. M.Sc. Thesis, COPPE/UFRJ, Rio de Janeiro, 1986. (in Portuguese). 19 J . C . F . TeUes, A self-adaptive coordinate transformation for efficient numerical evaluation of general boundary element integrals, Int. J. Num. Meth. Eng. 24:959-973 (1987). 20 L.C. Wrobel, Potential problems, in Topics in Boundary Element Research, Vol. 3: Computational Aspects, Chapter 8, Springer-Verlag, Berlin, 1987.
A New
h-Adaptive
Boundary
Element
Refinement Method
Scheme Using
for the
Local
Reanalysis
A. C h a r a f i
School of Mathematical Studies University of Portsmouth Portsmouth POI 2EG, UK and L. C. W r o b e l
Department of Mechanical Engineering Brunel University Middlesex UB8 3PH, UK Transmitted by M. A. Golberg
ABSTRACT This paper presents a new way of estimating the discretization error in the boundary element method (BEM) based on the concept of local reanalysis (LR). The principle of LR is implemented for the particular case of two-dimensional potential problems, in an a t t e m p t to find a technique to calculate a cheaplycomputable predicted solution to use in the evaluation of the quality of the boundary element results. The error estimators are defined from the difference between the BEM and LR solutions, which measures the "sensitivity" of the boundary element results to a local change in the mesh. An h-adaptive method using LR is developed and different numerical experiments are carried out to test its performance. (~) Elsevier Science Inc., 1997
1.
INTRODUCTION
A n y a d a p t i v e s t r a t e g y has to answer t h e t h r e e basic questions of when, where, a n d how t o refine. R e f i n e m e n t will b e n e c e s s a r y w h e n t h e s o l u t i o n o b t a i n e d f r o m t h e n u m e r i c a l a n a l y s i s is n o t satisfactory. T h e a c c u r a c y of t h e s o l u t i o n m a y b e assessed b y use of e r r o r e s t i m a t o r s , b u t e n g i n e e r i n g c o n s i d e r a t i o n s m a y also b e used t o t e s t t h e r e l i a b i l i t y of t h e results. O n c e it h a s b e e n j u d g e d n e c e s s a r y to refine, error i n d i c a t o r s a r e u s e d t o show
APPLIED MATHEMATICS AND COMPUTATION82:239-271 (1997) (~) Elsevier Science Inc., 1997 655 Avenue of the Americas, New York, NY 10010
0096-3003/97/$17.00 PII S0096-3003(96)00028-8
240
A. CHARAFI AND L. C. WROBEL
where the refinement should occur. To the question of how to proceed with the refinement there are basically three schemes. • Nodal redistribution or r-type. • Element subdivision or h-type. • Local variation of the order of polynomial approximation or p-type. Hybrid methods, combining the merits (and demerits) of any of these schemes, are also possible. The hp-method, combining h- and p-adaptive schemes, is the most common. The solution of an adaptive problem is similar for all the various methods and obeys the following algorithm. 1. From an initial coarse mesh, a first approximation to the solution is computed. 2. The accuracy of the solution is assessed using suitable error estimators. 3. The mesh is refined with local error indicators guiding where the refinement should occur. 4. From the new mesh an improved solution is calculated. 5. The procedure (1) to (4) is repeated until the error estimator indicates that no further refinements are necessary or a specified stopping criteria is satisfied. In the h-method, the order of the polynomial approximation inside the element is kept constant, and the mesh is refined locally by subdividing the elements with large errors. This results in a graded mesh with a concentration of smaller elements near sources of errors. The first attempt on developing an h-version of BEM was by Rank [1] for a simple Galerkin boundary element method with constant elements. He presented the first theoretical results on a posteriori error estimation for the h-adaptive BEM, that were later extended to more general integral operators as well as to collocation BEM by Wendland and Yu [2]. Rank [1, 3], based on theoretical results of a posteriori error estimation, proposed the use of the H i - n o r m of the element residual to guide the hrefinement. The scheme was implemented for two-dimensional potential problems. Katz [4] developed an h-version of adaptive BEM to obtain the warping function of beams, for which the Galerkin method was employed to approximate the integral equation. Rencis and Mullen [5, 6] developed an h-adaptive boundary element for two-dimensional elasticity problems, in which they used an error indicator that predicts the number of uniform element subdivisions needed to obtain
Adaptive Refinement for Boundary Element Method
241
a pre-assigned accuracy. Although their approach was based on an oversimplification of an a priori error estimation result, and limited to constant elements, it nevertheless highlights some potentials of the h-refinement in improving the quality of the solution. Rencis and Jong [7, 8] extended the concept of Zienkiewicz and Zhu [9] for the F E M to develop a predicted error in a direct collocation BEM for plane and axisymmetric problems governed by the Laplace equation. Guiggiani [10, 11] has introduced non-residual-type error indicators by calculating the difference in the solutions obtained from two different sets of collocation points. He has applied his approach to two-dimensional elasticity problems with quadratic elements. Kamiya and Kawaguchi [12, 13] developed an h-version of BEM using the concept of "sample point error analysis" to solve two-dimensional potential and linear elastostatics problems. Abe [14, 15] proposed an error indicator and an evaluation of the nodal error, in the h-adaptive BEM, based on a new residue calculated at collocation points only. The nodal error was evaluated from the approximated residue by solving a matrix equation and the proposed strategy was applied to solve two-dimensional potential problems. The h-version has proved to be satisfactory especially when dealing with problems with singularities but the non-uniform element spacing obtained may lead to possible problems with the numerical integration and the conditioning of the system of equations, unless a sophisticated numerical integration scheme is used.
2.
T H E P R I N C I P L E O F LOCAL REANALYSIS
In this approach, starting from the assumption that a local change in the mesh will affect the solution only if the latter is not accurate enough, we introduce the concept of local reanalysis to measure the sensitivity of the solution once the mesh has been locally modified. In what follows, the concept of local reanalysis is explained. Consider the following boundary value problem in the two-dimensional domain ~ and its boundary F Au = 0 u = ~ q= ~
in on F~ on F~
(1)
where u and q = (Ou/On) are the potential and the flux, respectively, and n is the unit outward normal vector of F.
242
A. CHARAFI AND L. C. WROBEL
The boundary integral equation (BIE) of this two-dimensional potential problem for an arbitrary source point x i on the boundary is given by
f u ~ OE2 xi
c(xi)u(x') + Jr ( )-~n(
- ()d~ = ~r q(~)E2(x' - ~)d,~,
(2)
where E2 is the fundamental solution of the Laplace equation. After discretization of the boundary into L boundary elements, (2) becomes
c(x')u(x') +
u(()--~-~n ix - ()d( = e=l
e
~.q(()E2(x' - ()d~
(3)
e=l
After an approximation of the variation of u and q within each element using interpolation functions and nodal values, this equation is then applied for i = 1 , . . . , n, where n is the number of collocation points on the boundary F. The quality of the boundary element solution is strongly dependent on the discretization employed (number of elements and degree of polynomial approximation). It is this discretization error that is studied here, although there are other sources of error, such as round-off errors, integration errors, and the improper representation of curved boundaries a n d / o r of some boundary conditions. After a first approximation is obtained using the BEM, we need to test its quality. This is achieved by using LR over each element in turn, to decide whether it needs to be refined or not. This step consists of using the BEM, on the selected element for reanalysis, to solve the integral equation obtained when assuming that the boundary values on all the other elements are already known from a previous full BEM analysis. Suppose that an element Fe has been selected for the reanalysis. Equation (3) holds only at collocation points used to get the BEM solution. LR consists of using the values of (u, q) on the portion r - Fe from the initial solution, namely (~, q~, and solving the following BIE for the unknowns (~, q~ on element Fe
. OE2 ,(x - ~)d~ = b(x) + f c(x)~(x) + a(x) + f r ~(~)--0~n
~(~)E2(x - ~)d~ (4)
JFe
e
where a(=)
=
-
(5)
Adaptive Refinement for Boundary Element Method
243
f / #
\
/
x
/ /
\
I I
t
Original n o d a l
I
points
\
!
I I
/
\\ \ \
~/
, /'
Additional nodal
points
FIG. 1. Local discretization for local reanalysis. and
b(x) = E [ ~(~)E2(x - ~)d~ j~e 4r
(6)
To solve (4), the element Fe is subdivided into two subelements of the same type as those used to get the initial BEM solution, which leads to the discretized equation ne
ne
c , ~ + as + E H, kuk = b, + E k=l
aikq'k
(7)
k----1
where n~ represents the number of collocation points on F~ (if quadratic elements are used n~ = 5, see Figure 1), and G~k and H~k are influence coefficients obtained by integrating the fundamental solution and its normal derivative, respectively, over the two subelements. After applying the prescribed boundary conditions on F~ and rearranging the coefficients H~k and Gik accordingly, a linear system of equations is formed [AI{X} = {Y} + {b} - {a}
(8)
which is very small compared to the one obtained when doing the BEM analysis on the whole boundary, since the degrees of freedom are only those on the portion F~ of the boundary F. Hence, local reanalysis is relatively
244
A. CHARAFI AND L. C. WROBEL
r-l I
I I,
!/
n,
i
[A]
I I
',
I I
I
I
I
:
[×!/ /
/
Effect of additional collocation points
I .... Global Analysis Local Analysis
FIG. 2. Matrix illustrationof localreanalysis.
quick to perform. Figure 2 shows how the local system of equations for the local reanalysis is obtained from the global one. Computing a new solution using the principle of local reanalysis, and comparing it to the BEM solution, provides a quick feedback on the quality of the results, enabling to build an error estimator---substituting the exact solution by the one obtained from the local reanalysis--which provides a measure of how the solution is affected by a local change in the mesh. This idea has been first suggested by Casale [16] in his dissertation thesis. Its use as an error estimator has already proved reliable in the case of potential problems, as shown by Chara~i et al. [17]. With regards to the computational efficiency of the LR procedure, while it is true that, for the new collocation points, new integrals have to be evaluated on all boundary elements, most of these integrals can be reused in the next solution step if the new collocation points axe retained. Thus, very little additional work is wasted in building the error estimates, and the scheme works very efficiently. 3.
STUDY OF T H E LOCAL REANALYSIS SOLUTION
To have a good assessment of the solution obtained from the local reanalysis procedure, one has to minimize errors arising from the approximation of the geometry and boundary conditions, as well as those caused by the numerical approximation of the integrals involved in BEM. Straight-line
Adaptive Refinement for Boundary Element Method
245
FIG. 3. Possible boundary conditions at a corner for potential problems. Double node
, Collocation p o i n t Fro. 4. Definition of quarter-point elements. geometries are the simplest to treat. However, they introduce corner points which have to be appropriately dealt with. These numerical considerations are described first, before analyzing the local reanalysis procedure.
Numerical considerations Geometry and boundary conditions
As linear elements are used for the numerical t r e a t m e n t of the boundary integral equation, the examples are chosen so t h a t the geometry is correctly represented (straight lines) and the b o u n d a r y conditions exactly applied (linear variation at most).
Treatment of corners In the B E M solution of a model defined by straight lines, corners have to be appropriately dealt with. There are three types of possible boundary conditions in the vicinity of a corner as shown in Figure 3. Because of the continuity of the potential, the only unknown at the corner is the potential u in the first case, and the flux q (on the side where u is prescribed) in the second. But because the flux can be discontinuous, it is necessary in the third case to solve for the flux before and after the corner. This is achieved by use of non-nodal collocation elements presented by Marques [lS]. T h e basic idea behind this type of elements is to define two identical geometrical nodes at the corner point. To avoid obtaining a singular system of equations, the collocation points are located inside the elements at the quarter-points, as explained in Figure 4. Special integration scheme To minimize the effect of the error from numerical integration, and to focus on the discretization error generated by
246
A. CHARAFI AND L. C. WROBEL
the collocation method, a suitable integration scheme proposed by Telles [19] is used. It consists of an adaptive Gaussian quadrature scheme using a third degree polynomial transformation which bunches the integration points closer to the collocation points, and thus allows a smaller number of integration points for the same accuracy. Furthermore, Telles' transformation gives better results than if Gaussian quadrature is applied directly, and permits the computation of terms with singularities of order O(ln(r)). Numerical experiments Numerical experiments are carried out in order to compare the "predicted" solution obtained from the local reanalysis and the "next" boundary element solution obtained from a uniformly refined model. By "next" boundary element solution we mean the solution obtained when subdividing all the elements. The aim of this numerical investigation is to assess whether local reanalysis can (firstly) predict the "behavior" of the solution for the next step, and (secondly) provide information about areas of large errors.
EXAMPLE 1. A rectangle subject to Robin boundary conditions. Description
This example deals with a steady-state heat conduction in an exposed exterior column of rectangular cross-section, studied in detail by Wrobel [20]. A simplified model is presented in Figure 5, where one half of the column is studied, subjected to boundary conditions of the Robin type, q = au + b. The values of the coefficients a and b are given by a=-0.5 a=0 a a=-6
and and and and
b=50alongABC, b=0atD, b linear along CD and DE, b=0alongEFG.
C
D
E
|
I
l
I,
8
I
/,=36 FIG. 5. Initialdiscretizationfor Example I.
(9)
F =7
Adaptive Refinement for Boundary Element Method
247
84.0 75.2
~
u
~ Initial BEM s o l u t i o n 0 Predicted LR s o l u t i o n
• ~-~I~%
66.4
o_ 8, N. t
"
,olutlo.
57.8 48.8 .,..$ c 40.0 0 O.
31.2 22.4 13.6 4.8 -4.0
B
C
FIG. 6. Comparison of B E M
D
E
F
and L R solutions for Example I.
Results To study the behavior of the "predicted" solution obtained from local reanalysis, a comparison is made with the boundary element solution obtained from both the initial mesh described in Figure 5, and the "next" refined mesh. By "next" mesh, we mean the mesh obtained by subdividing into two all the elements of the initial one. The three solutions on side [BF] are presented in Figure 6. The results of the "predicted" solution from LR and the "next" solution from BEM are in good agreement, which suggests that the LR procedure can be thought of as a "quick" way of obtaining an approximation to the "next" BEM solution on the whole boundary. EXAMPLE 2. Heat conduction in a transformer coil.
Description In this example the heat conduction in a transformer coil [7] was considered as shown in Figure 7. The initial discretization employed 15 boundary elements. Results The BEM solutions obtained from the initial discretization (15 elements) of Figure 7 and from the mesh obtained by subdividing into two all the elements of the initial mesh, as well as the solution obtained from the LR procedure applied to the initial mesh, are presented in Figures 8-9.
248
A. C H A R A F I
L
u=0
A N D L. C. W R O B E L
K
!-
l~-V
u=0
~0
q=5~---
'~
"I
.
I
D I q=~-~_ C
AT
u--O
l-
'
i
s
I
B
.~J
J
J
FIG. 7. Geometry and boundary conditions for Example 2.
6O 54
o Initial BEM solution I t o Predicted LR solution o Next BEM solution J
48 42 v
36 30 O
24 18 12 6 0
B
C
D
E
F
-
-
•
s
_
G
H
I
J
K
FIG. 8. T h e potential u from B E M and LB. solutions for E x a m p l e 2.
Adaptive Refinement/or Boundary Element Method
249
1
--3
j
K
L
a
B
FIG. 9. T h e flux q from B E M a n d L K solutions for E x a m p l e 2.
When solving for the potential u (smooth solution), see Figure 8, the "predicted" solution from LR does not approximate the "next" BEM solution as well as in the previous example, especially between the points E and H. However, the LR procedure has correctly predicted the "behavior" of the solution, and this should be sufficient to indicate the need to refine this area. When solving for the flux q, see Figure 9, the solution obtained from LR still predicts the "behavior" of the BEM solution , and the difference between the "predicted" and "next" solutions are more apparent in the neighbourhood of points K and B. EXAMPLE 3. Singular flux problem
Description This example is a square domain ABCD of length 10 and with boundary conditions specified by u = 0 u = 10 q= 0
along AB at C and linear along BC along CDA.
The flux q is singular at the point C.
(10)
250
A. CHARAFI AND L. C. WROBEL 10
9 8
i
\,.
7
~
6
0
o Initial BEN s o l u t i o n "
~
~
Predicted LR solution
o
NextBEg
solution
\
4
3
~
1
~
0
C
~...~
s
.
D
A
FIQ. 10. The potential u from BEM and LR solutions for Example 3.
Results The difference in Figures 10-11 between the initial BEM solution on one hand, and the next BEM as well as the"predicted" LR solutions on the other hand, strongly indicates the need for more refinement of the mesh on elements [Be], [CO] and IDA]. Moreover, Figure 11 shows that the LR has predicted the discontinuity of the flux q at point B and its singularity at point C (side [BC]).
Conclusions from the numerical experiments It is clear from the study of the previous examples, and the analysis of the corresponding results, that the LR solution succeeds in giving an indication of the behavior of the BEM solution. It also has the ability to predict the existence of singularities. It should not, however, be used as an attempt to improve the solution itself as the LR measures mainly how "sensitive" the BEM solution is to a local change in the mesh. 4.
AN h-ADAPTIVE T E C H N I Q U E W I T H LOCAL REANALYSIS
After performing the local reanalysis over each element of the boundary F we get a solution ~ and ~', considered as the predicted solution. The exact
Adaptive Refinement for Boundar~d Element Method
251
3.0 2.5
1~ InitlalBEM solution o Predicted LR~mlution
.
o NextsE~~ i u u o n
Z.O
2
//~/-
/
1,5
o"
/
1.0
f
v
>¢
0.5 _
Z
0.0
-0.5 -1.0 -1.5 -20
"'0
A
B
C
FIG. 11. The flux q from BEM and LR solutions for Example 3. error for the potential (resp. flux) is approximated by e~ ~- ~ - ~
(resp. eq ~_ q - q~).
(11)
T h e basic (although empirical) idea is t h a t the difference between these two solutions, the one obtained from the boundary element discretization and the one obtained from the local reanalysis, provides some information about where the refinement should occur. T h e "error" measured is more a kind of "variation" or "sensitivity" of the solution due to a local change in the mesh. It should be noted t h a t the solution from L R cannot be regarded as more accurate t h a t the solution from BEM. T h e basic idea is t h a t the two solutions will be "more different" in those parts of the boundary where the discretization is poorer. g~'For no77r~8
To obtain a measure of the error (or rather of its estimate), we found it convenient to compare its H ° - n o r m on each b o u n d a r y element. We define the H°-local norms for the potential u and its associated error e~ on an element i as
1'' i
/12/ i
252
A. CHARAFI AND L. C. WROBEL
J
//// --(
1
! I
3
2
FIG. 12. Representation of the solutions ~ and ~ on a linear element. Because linear elements are used, the quantities [[u[[i and [[eu[[i can be evaluated analytically. On an element Fi the solutions ~ and ~ are represented as shown in Figure 12. E v a l u a t i o n o f [[ulli We have
~(~/)=
~'I+
- -
~2
-1_
(13)
therefore, r
Ilull2 = - ../r~u2(x)dx
(14)
(15) l'[ (~2- ~1)~ ] - 4 § + (~' + ~2)2 '
(16)
where li represents the length of element i. Evaluation of We have
[[e.lh
{ (1
- ~ 1 + (1 + ~)~3
~(~)
=
-
~)~3 + ~ 2
-l~y
(17)
-1<~0 0<~<1,
(18)
therefore, e~'(~) = { air/2 + bl~ + Cl a2~7 2 -b b2~] -}- c2
AdaptiveRefinementfor BoundaryElementMethod
253
where
al :
[U2 -- E1 ]2 T -~-~1 -- u3
bl : 2 [ ~ - F U l - U Z ]
[ ~1+~2 2
~3]
(20)
u2
(21)
a2---- - - - - 7 - +
u3 -- u2
(22)
b2= 2 [^~
^ "4-~2 +u3 -u2] [Ul 2
C1 =
ul +2 u2
(19)
and
C2 :
u 1 +2~ 2
~3]
(23)
u3
(24)
The value of the local error norm is then given by
l~ [(al +a2) Ile~,ll~ = ~" 3
+ (~-~)
+(Cl +C2)]
(25)
The H°-global norms for the potential and its associated error for a discretization of L elements are derived from the local ones by
[ L
~1/2
[ L
,~1/2
The H°-local and global norms and error norms for the flux are defined in exactly the same way as for the potential.
Error indicators As the error norms defined have little physical significance and can be misleading, we introduce nondimensional error percentage values. Following the work by Zienkiewicz and Zhu [9], we define the global error percentage for both the potential and the flux respectively, by
1 •,-~, = k, ll~.ll 2 + ile,~ll2 ]
x loo
and
Tq =
i1~12 + ileqll ~
× lOO (27)
where [[~[12 + [[e~l[ 2 and [[q~[2 + Heq[[2 approximate the exact potential norm Ilull2 and the exact flux norm Ilqll2, respectively.
254
A. CHARAFI AND L. C. WROBEL
Assuming a uniform distribution of the error on the boundary F, the relative error percentage on the element Fi for both the potential and the flux respectively, is defined by 1
ii~--~i+ ]le~]12]
1
x 100
and
r 4 = ~H~I2 +
11%112] × 100.
(28)
Computational strategy After solving the problem with the initial data (geometry and boundary conditions), the elements are analyzed one by one. O n each element, the percentage error obtained from the difference between the LR and BEM solutions is computed to decide whether to keep the refinement or not. The new mesh is then defined, refining the elements where the error is not within the imposed tolerance, and the problem is solved as if it is a new one. The self-adaptive mesh refinement strategy is controlled by specifying for both the potential and the flux on all elements that ~_<~
and
Tq<~
(29)
where ~ is a user-specified tolerance. It is important to notice that ~ is a local requirement. The refinement procedure can be terminated by different means. • By specifying an error bound for the global error percentage of both the potential and the flux, if one is interested in the global accuracy. However, it should be kept in mind that a good overall accuracy could well mask large errors in small parts of the boundary. • By storage considerations, that is, by specifying the maximum number of degrees of freedom the user is prepared to include in the analysis. • By time limitation considerations specifying, for example, a maximum number of resolutions. The details of the programming procedure of the error estimation and the adaptive mesh refinement strategy are given in Figure 13.
Numerical examples The h-adaptive BEM with LR is tested on five examples. The first three examples are those presented in Section 3, the fourth example is a sheet pile problem with a singularity and the fifth example is an L-shape domain along which the solution is smooth. In each case, the final potential and flux solutions obtained from the adaptive procedure are compared to a "reference" solution obtained by solving the same problem with a uniform
Adaptive Refinement for Boundary Element Method
255
[ Define geometry and boundary conditions
[
•] Perform standard B E M analysis i
l Compute the norms for u and q : [[~[[2 and
I1~112
L
] Loop on the elements ]
Local reanalysis : , Subdivide element and compute new H and G for new nodes , Form small system of equations AX = B , Solve the small system of equations , Compute the local norms of the errors eu and eq : Ile.ll~ and Ile, ll~
I End of loop] l Compute the global error norms :
I1=,,11= = r.~ Ile,,ll~
II~qll 2 = ~i Ileqll~
and
1
lu Compute the global errors percentage r,, and rq :
= ( ~ ) ~
× ~oo , ~ = ( ~ ) ~
× ~oo
l TGet relativepercentage errorson the elements :
[
•
{
r~=
~
½
Xll,~ll%lle,,ll~ /
~ ~nn
. . . . .
r ~-
q -
f
~
½
~l1411,+ll~dl, I
x 100
1 [ For aU elements r ~ < ~ and r ~ _ _ f ? J Yes I No 1
~~or
!
elementsi such that
: r~ > ~ or ¢~ > keep the subdivision
FIQ. 13. Flow chart of/z-adaptive BEM with local reanalysis.
256
A. CHARAFI AND L. C. WROBEL
mesh of a large number of elements enough to capture the behavior of the exact solution. EXAMPLE 1. This is the same Example 1 defined in Section 3. The model was analyzed starting from the coarse initial mesh of 6 elements, shown in Figure 5, which is sufficient to describe the geometry and boundary conditions. The adaptive scheme was tested with a specified local tolerance ~ = 0.04%. The adaptive solution was compared to a "reference" solution obtained by solving the same problem with a uniform mesh of 100 elements. Table 1 shows the decreasing values of the percentage error for the potential through the different iterations. The final mesh presented in Figure 14 was obtained with 50 elements. Figure 15 shows t h a t the distribution of the nodes in the final mesh is in good agreement with the "reference" solution. The same example was reanalyzed using a uniform mesh of 50 elements, for which the global percentage error was T~ = 0.0651% and the m a x i m u m local percentage error was z max = 0.2545%. With the adaptive procedure, the value of the global percentage error has decreased to T~ = 0.0248%, and the m a x i m u m local percentage error has decreased to Tmax = 0.0126%. TABLE 1 GLOBAL THE
PERCENTAGE
ADAPTIVE
Iteration
ERROP~S
PROCESS
FOR
THROUGH EXAMPLE
1
N u m b e r of e l e m e n t s
~'u [%]
6 12 22 38 48 50
3.2296 0.7753 0.2075 0.0580 0.0267 0.0248
0 1 2 3 4 5
C
D
E -v
A -
vA
A I ......................T....... T..... G FIG. 14. F i n a l a d a p t i v e m e s h for E x a m p l e 1.
Adaptive Refinement for Boundary Element Method
257
80 72 • Adaptive s o l u t i o n J - - Reference s o l u t i o n J
64 56 48 .,.d
Q) O
4O 32 24 16 8
0 B
C
D
E
F
Boundary length FIG. 15. F i n a l solution for E x a m p l e 1.
EXAMPLE 2. The geometry and boundary conditions for this example are the same as in Example 2 in Section 3. The adaptive scheme was started with the minimum number of elements sufficient to describe the geometry and boundary conditions, that is, 12 elements. The local tolerance imposed was ~ - 0.25% and the values of the global percentage error for both the potential and the flux, for each iteration, are presented in Figure 16, which indicates the convergence of the adaptive process. The final mesh from the present adaptive procedure, presented in Figure 17, was obtained with the global percentage tolerances ~'u = 0.2285% and rq = 0.1386%. The adaptive solution, for both the potential and the flux, is compared to a "reference" solution obtained from a quasi-uniform mesh of 216 elements. By quasi-uniform is meant a mesh in which the ratio of two adjacent elements is equal to either 1 or 1.5. Figure 18 shows that the two potential solutions are in good agreement, with the distribution of nodes capturing the curvature of the solution as well as the convex-concave changes in the curve. The flux solution obtained from the final adaptive mesh in Figure 19 is also in good agreement with the "reference" solution. Solving the same problem using linear elements and an adaptive scheme with the "sample point error analysis", Kamiya and Kawaguchi [12] obtained the mesh presented in Figure 20. It is not always easy to compare
258
A. C H A R A F I A N D L. C. W R O B E L
14.0 12.8 I o Global potential error Tu I I ~ Global flux error Tq ]
It.2 O I. t~
9.8
@ 8.4 U0 c- 7.0 5,6 4.2 O
2.6 1.4 0.0
10
19
28
37
48
Number
55
64
73
82
of e l e m e n t s
FIG. 16. Convergence graph for Example 2.
L.
K
G N u m b e r of e l e m e n t s : E
F
C B FIG. 17. Final adaptive mesh for Example 2.
9t
100
Adaptive Refinement for Boundary Element Method
259
60 54 48
42
"~
30
0
24
O.,
18
12
B
C
D
E
F
G
H
I
J
K
Boundary length FiG. 18. Final potential solution for Example 2. two different adaptive schemes as their respective error indicators and estimators may measure different entities (residual error, interpolation-type error, sensitivity, etc...). However, one can always compare the final meshes with respect to a "reference" solution obtained from a very fine mesh, or indeed the exact solution when it is available. ComParing the mesh obtained from the present adaptive scheme with t h a t obtained by Kamiya and Kawaguchi [12] in Figure 20, one notices that the present scheme has produced a graded mesh towards K along the side [KL], contrary to the mesh produced by Kamiya and Kawaguchi [12]. A close look at the "reference" solution in Figure 19 shows the strong variation of the flux in the neighbourhood of K (side [KL]). Overall, the final meshes obtained by the two methods are qualitatively similar. EXAMPLE 3. This is the same Example 3 defined in Section 3. To model this problem, all that is needed is to define the 4 nodes of the square. The adaptive procedure automatically increases the n u m f e r of elements until the specified tolerance of ~ = 0.125% has been achieved or a stopping criteria reached. Table 2 shows the values of the percentage errors for both the potential and the flux for each iteration. Note the oscillation of the percentage error
260
A. C H A R A F I A N D L. C. W R O B E L
0 -I -2
p~ v ,--4
r-~
-3 -3 -4 -5 -6 -7 -8
K
L
I
L
A Boundary
B
length
FIG. 19. Final flux solution for Example 2.
L
K
H! . . . . . . . . . .
Number of elements: 96
_G
I F
D t ...... ......
A
e e 2 s s ~ ; : s ; : : s s
t
C
B
FIG. 20. Final adaptive mesh for Example 2 by Kamiya and Kawaguchi.
Adaptive Refinement for Boundary Element Method
261
TABLE 2 GLOBAL PERCENTAGE ERRORS THROUGH THE ADAPTIVE PROCESS FOR EXAMPLE 3 Iteration 0 1 2 3 4 5 6
N u m b e r of elements
~u [%]
Tq [%]
4 8 16 30 50 63 75
15.4175 5.2777 1.8933 0.6739 0.2395 0.1007 0.0740
27.2363 22.0559 34.4350 10.7786 7.6949 5.4770 3.8964
N u m b e r of e l e m e n t s : 75
D
A
FIG. 21. Final adaptive mesh for E x a m p l e 3.
for the flux at iteration 2, due to the high value of the local percentage error for the flux near the point C (side [PC]), caused by the singularity of the flux solution at that point. The final mesh is represented in Figure 21. The final solution from the adaptive scheme is compared to the "reference" solution obtained from a uniform mesh of 320 elements. Figure 22 shows that the nodal partition captures very well the behavior of the potential solution with more nodes towards C (side [CD]) and less towards A (side [DA]). The fluxes from the adaptive and the "reference" solutions are in good agreement as shown in Figure 23. The singularity at the point C (side [BC]) is very well captured by the graded mesh towards point C.
262
A. C H A R A F I A N D L. C. W R O B E L 10
\
[* Adaptivesolution ] ~. ReferencesolutionI
0
i
D Boundary
length
FIG. 22. Final potential solution for Example 3. 5.00
4.34
* Adaptive solution m Reference solution
3.68 3.02
j
,---- 2.36 Ce
1.70 r,..
1.04
0.38 f
-0.28 -0.94 -
I
1.60
A
B Boundary
length
FIG. 23. Final flux solution for Example 3.
Adaptive Refinement for Boundary Element Method
A
263
E D,
15
C
B L I-
18
-I
FIG. 24. Initial model for Example 4. EXAMPLE 4. This example is a square A B C D E representing a sheet pile problem, see Figure 24. The boundary conditions are given by q u q u
= = -=
0 0 0 157
along along along along
ABC, CD, DE, EA.
(30)
This example is particularly interesting because it presents a singularity at the end of the cut-wall DE (the flux is singular at point D, side [CD]). The tolerance specified was Y -- 0.25% and the results obtained are shown in Figures 25-27. Note t h a t the final mesh, with 51 linear elements, is more refined (as expected) in the vicinity of point D due to the singular behavior of the flux. As the error indicators are obtained by computing the integral of the difference between the BEM and LR solutions over the element considered, their value is expected to be high near the singularity. This is what caused the oscillation at the beginning of the adaptive process for the global percentage error Tq, see Table 3. This influence, however, decreases as the length of the element near the singularity gets smaller. T h e high value of Tq is due to the large influence of the local error in the vicinity of point D (on the other parts of the boundary the local errors are very small), and also to the increasing gap between the predicted solution from the local reanalysis and the B E M solution at point D, as shown in Figure 28.
264
A. CHARAFI AND L. C. WROBEL
E
A
D
Number
=
=
of elements:
=
=
|
51
i
C
B FIG. 25. F i n a l a d a p t i v e m e s h for E x a m p l e 4.
162.0
145.8 129.6 113.4 v
97.2 81.0 q) 64.8 O 48.6
\
32.4 16.2 0.0
A
8 B
"4 C
Boundary length FIG. 26. Final potential solution for Example 4.
D
Adaptive Refinement for Boundary Element Method
265
TABLE 3 G L O B A L P E R C E N T A G E E R R O R S T H R O U G H THE A D A P T I V E PROCESS FOR E X A M P L E 4
Iteration
Number of elements
r u [%]
Tq [%]
5 10 19 33 41 47 51
6.7605 2.1247 0.6175 0.2605 0.1967 0.1775 0.1724
42.9869 35.6410 37.1137 32.1726 30.5618 29.8507 29.1731
0 1 2 3 4 5 6
400
36
380
320
33 • --
Adaptive solution Reference So ution
30
280
27
240
24
~ 200 X
x
mo
r,-- Is
120
15
80
12
40
9 g
D
E
A
FIG. 27. Final flux solution for E x a m p l e 4.
The oscillation of the flux solution near point D (side [CD]), is a numerical phenomenon that occurs sometimes near singularities when using B E M with linear elements. It is not linked to the adaptive process, as both adaptive and "reference" solutions show the same behavior, see Figure 27. This oscillation did not occur in Example 3, which has a singularity at point C (side [BC]), leading to think that it could be linked to the order of the singularity. EXAMPLE 5. This example is an L-shape domain with the boundary conditions specified in Figure 29. The "reference" solution was obtained from a uniform mesh of 320 elements for which the global tolerances were Tu = 0.029% and ~'a = 0.0008%. The adaptive procedure was tested with a specified tolerance of Y = 0.1%, and starting with an initial mesh of 6 elements which are sufficient to describe the geometry and boundary conditions.
266
A. CHARAFI AND L. C. WROBEI u
J'
',
'
"
,
"
.
' •
.
.
7
'
.
,
•
I
I
",
I
'
: BEMsolution I P~edicted solution
-60 -120 ~ C~ -160 ~J
-240
x ~=~
,~ 0
-300
\
-360
-420
2
~ -460 -540 -600
2
2
3
3
4
5
5
6
6
N um ber of r e s o l u t i o n s
FIG. 28. Flux variation at the point D for different resolutions for Example 4.
E
D
u=10
•
q=O
._!___
C i
u=O
5
q=0
F
•
/
# I
..... I_ i
/
, i i
, I
/
/
'
," -" / I , , "
," '
i I.,£2
10
--',
FIG. 29. Geometry and boundary conditions for Example 5.
7
Adaptive Refinement .for Boundary Element Method
267
TABLE 4 GLOBAL PERCENTAGE ERRORS THROUGH THE ADAPTIVE PROCESS FOR EXAMPLE 5
Iteration 0 1 2 3 4
Number of elements
Tu [%]
~'q [%]
6 12 24 40 60
8.5704 2.2188 0.6352 0.2218 0.0993
11.3012 4.6408 0.3436 0.0938 0.0759
C
B
N u m b e r of e l e m e n t s : 60
F
h FIG. 30. Final adaptive mesh for Example 5.
T a b l e 4 shows t h e decreasing values of t h e global p e r c e n t a g e errors, for b o t h t h e p o t e n t i a l a n d t h e flux, t h r o u g h t h e a d a p t i v e process. F i g u r e 30 shows t h e final n o d a l d i s t r i b u t i o n . N o t e t h e g r a d e d m e s h t o w a r d s t h e corners C a n d F which reflects t h e b e h a v i o r of t h e p o t e n t i a l "reference" sol u t i o n as p r e s e n t e d in F i g u r e 31. T h e m a x i m u m local e r r o r p e r c e n t a g e for t h e p o t e n t i a l was ~max __ 0.4506% for t h e a d a p t i v e s o l u t i o n c o m p a r e d t o 7umax = 0.3591% for t h e "reference" solution. F i g u r e 32 shows t h a t t h e values of t h e flux v a r y in a n o n l i n e a r w a y b u t w i t h i n a v e r y s m a l l interval (form 0.75 t o 0.82 in a b s o l u t e values), a n d t h e a d a p t i v e p r o c e d u r e has c a p t u r e d t h i s b e h a v i o r . F i g u r e 32 also i n d i c a t e s t h a t t h e r e is still need for refinement t o w a r d s t h e ends of s e g m e n t s lAB] a n d [DE]. I n d e e d , t h e m a x i m u m local t o l e r a n c e for t h e flux was Tq ax = 0.2698% for t h e a d a p t i v e s o l u t i o n c o m p a r e d t o T~nax = 0.0062% for t h e "reference" solution.
268
A. CHARAFI AND L. C. WROBEL
10
9
I " Adaptive so]ution solution I
7 6 v
5
0J o
4
c~
3 2 I.
0
E
F
A
Boundary
B
length
FIG. 31. Final potential solution for Example 5.
0.02
-0.74, -0.78
0.81
Io Adaptivesolution ]
0.80
-0.76
0.80
-0.76
3
,,~
-0.77
0+719 0.78
-0,70 r... -0.79
0.77
-0.80
0.78
-0.80
0 .?mS
-0.81
0.75 0.74
-0.82
B
D
FIQ. 32. Final flux solution for Example 5.
C
D
Adaptive Refinement for Boundary Element Method
5.
269
CONCLUSION
A new adaptive BEM scheme has been developed, using error estimates based on the concept of LR, for the construction of near optimal computational models in two dimensions. The self-adaptive boundary element technique was developed based on an estimation of the discretization error in each element using the concept of LR. It consists of measuring the sensitivity of the boundary element solution to local changes in the mesh by comparing the solution obtained from the boundary element discretization with that obtained by solving a series of small linear systems of equations--3 x 3 for potential problems with linear elements. Numerical experiments with two-dimensional potential problems have shown that, although the LR solution is not always an improvement of the BEM solution (as can be seen in Section 3), it nevertheless indicates the behavior of the BEM solution and has the ability to predict the existence of singularities. Error indicators derived from the LR procedure measure how sensitive the numerical solution is to a local refinement in the mesh, yielding a mechanism to determine if the overall mesh is sufficiently accurate or if mesh refinement is needed. This was implemented in an h-version with linear elements for two-dimensional potential problems and has produced promising results by providing significant information about where to refine the mesh. It should be pointed out that the cost of the adaptive boundary element strategy described in this work is not always cheaper than the cost of other adaptive strategies available. However, while in many adaptive schemes-especially those based on residual-type error estimators--the auxiliary data computed to evaluate the error indicators and estimators is lost, all the new integrals computed for the evaluation of the LR solution can be re-used in the next solution step if the new collocation points are retained. Therefore, very little additional work is wasted in building the error estimates, which makes the scheme very efficient. REFERENCES 1 E. Rank, A posteriori error estimates and adaptive refinement for some boundary integral element methods, in Proc. Int. Conf. on Accuracy Estimates and Adaptive Refinements in Finite Elements Computations, Lisbon, 1984. 2 W . L . Wendland and De-hao Yu, Adaptive boundary element methods for strongly elliptic integral equations, Num. Math. 53:539-558 (1988).
270 3
4
5
6
7
8
9
10
11
A. CHARAFI AND L. C. WROBEL E. Rank, Adaptivity and accuracy estimation for finite element and boundary integral element methods, in Accuracy Estimates and Adaptive Refinements in Finite Element Computations, John Wiley & Sons, pp. 79-94, 1986. C. Katz, Self-adaptive boundary elements for the shear stress in beams, in Proc. 2nd BEM Conf., Computational Mechanics Publications, 1986, pp. 759-769. J . J . Rencis and R. L. Mullen, Solution of elasticity problems by a selfadaptive mesh refinement technique for boundary element computation, Int. J. Num. Meth. Eng. 23:1509-1527 (1986). J . J . Rencis and R. L. Mullen, A self-adaptive mesh refinement technique for boundary element solution of the Laplace equation, Comp. Mech. 3:309-319 (1988). J . J . Rencis and K. Y. Jong, A self adaptive h-refinement technique for the boundary element method, Comp. Meth. Appl. Mech. Eng. 73:295-316 (1989). J . J . Rencis and K. Y. Jong, An accuracy post-processor for boundary element analysis, in Proc. lOth BEM Conf., Computational Mechanics Publications, 1988, Vol. 1, pp. 187-203. O. C. Zienkiewicz and J. Z. Zhu, A simple error estimator and adaptive procedure for practical engineering analysis, Int. J. Num. Meth. Eng. 24:337357 (1987). M. Guiggiani, Non-residual type accuracy estimates and adaptive refinement in boundary element analysis, in Boundary Element Techniques: Applications in Engineering, Proc. 4th BETECH Conf., Computational Mechanics Publications, 1989, pp. 323-340. M. Guiggiani, Error indicators for adaptive mesh refinement in the boundary element method: A new approach, Int. J. Num. Meth. Eng. 29:1247-1269
(1990). 12 N. Kamiya and K. Kawaguchi, An attempt of boundary elements, in Proc. 12th Int. BEM Conf., Computational Mechanics Publications and SpringerVerlag 2:527-538 (1990). 13 N. Kamiya and K. Kawaguchi, Error analysis and adaptive refinement of boundary elements, in Reliability and Robustness of Engineering Software II, C. A. Brebbia and A. J. Ferrante (eds.), Computational Mechanics Publications and Elsevier Applied Science, 1991, pp. 187-201. 14 K. Abe, h-Adaptive boundary element method and evaluation of error, in Proc. IABEM Symposium, Springer-Verlag, 1993, pp. 1-10. 15 K. Abe, A new residue and nodal error evaluation in h-adaptive boundary element, Adv. Eng. Soft. 15:231-247 (1992). 16 M. Casale, The Integration of Geometric Modeling and Structural Analysis using Trimmed Patches. Ph.D. Thesis, University of California, Irvine, 1989. 17 A. Charafi, L. C. Wrobel, and R. Adey, An approach to h-adaptive boundary element methods using local reanalysis, in Proc. 7th BETECH Conf., Computational Mechanics Publications and Elsevier Applied Science, 1992, pp. 905-918.
Adaptive Refinement for Boundary Element Method
271
18 E. Marques, Coupling of the Finite Element and the Boundary Element Method: An Application to Potential Problems. M.Sc. Thesis, COPPE/UFRJ, Rio de Janeiro, 1986. (in Portuguese). 19 J . C . F . TeUes, A self-adaptive coordinate transformation for efficient numerical evaluation of general boundary element integrals, Int. J. Num. Meth. Eng. 24:959-973 (1987). 20 L.C. Wrobel, Potential problems, in Topics in Boundary Element Research, Vol. 3: Computational Aspects, Chapter 8, Springer-Verlag, Berlin, 1987.