Procedia Computer Science
Available online at www.sciencedirect.com
Procedia Computer Science 1 (2010) 2583–2591 Procedia Computer Science 00 (2009) 000–000
www.elsevier.com/locate/procedia www.elsevier.com/locate/procedia
International Conference on Computational Science (ICCS 2010)
The Singular Function Boundary Integral Method for Laplacian problems with boundary singularities in two and three–dimensions Christos Xenophontos*, Evgenia Christodoulou, Georgios Georgiou Department of Mathematics and Statistics, University of Cyprus, PO BOX 20537, Nicosia 1678, CYPRUS
Abstract We present a Singular Function Boundary Integral Method (SFBIM) for solving elliptic problems with a boundary singularity. In this method the solution is approximated by the leading terms of the asymptotic solution expansion, which exists near the singular point and is known for many benchmark problems. The unknowns to be calculated are the singular coefficients, i.e. the coefficients in the asymptotic expansion, also called (generalized) stress intensity factors. The discretized Galerkin equations are reduced to boundary integrals by means of Green’s theorem and the Dirichlet boundary conditions are weakly enforced by means of Lagrange multipliers, the values of which are introduced as additional unknowns in the resulting linear system. The method is described for two–dimensional Laplacian problems for which the analysis establishes exponential rates of convergence as the number of terms in the asymptotic expansion is increased. We also discuss the extension of the method to three–dimensional Laplacian problems with exhibits edge singularities. c 2010 Published by Elsevier Ltd.
Keywords: boundary singularities; Lagrange multipliers; stress intensity factors; boundary approximation methods
1. Introduction In the past few decades several methods for treating elliptic boundary value problems with boundary singularities, have been proposed. Among them one finds the so-called hybrid methods [1] which incorporate, directly or indirectly, the form of the local asymptotic expansion for the solution u, in the approximation scheme. This expansion is known in many occasions and has the following form: f
u (r ,T )
¦ D W r ,T , j
j
(1)
j 1
where (r, ș) are polar coordinates centered at the singular point, D j \ are the singular coefficients and W j r ,T
P
r j f j (T ) ,
* Corresponding author. Tel.: +357-228-92610; fax: +357-228-92601. E-mail address:
[email protected].
c 2010 Published by Elsevier Ltd. 1877-0509 doi:10.1016/j.procs.2010.04.292
(2)
C. Xenophontos et. al./ Procedia Computer Science 00 (2010) 000–000
2584
C. Xenophontos et al. / Procedia Computer Science 1 (2010) 2583–2591
are the singular functions, with ȝj, fj being the eigenvalues and eigenfunctions of the problem which are uniquely determined by the geometry and the boundary conditions along the boundaries sharing the singularity. It is important to note that the singular functions Wj satisfy the governing partial differential equation (PDE) and the boundary conditions (BCs) along the boundary parts sharing the singularity. Knowledge of the singular coefficients is of great importance in many engineering fields, such as fracture mechanics. Many methods have been proposed in the literature for their effective and efficient approximation, including high order p/hp finite element methods with post-processing [2, 3] and Trefftz methods [1]. In the former, the solution is first approximated on a refined grid designed especially to capture the singularity and the coefficients are obtained by an extraction formula which uses the computed solution. Methods that do not require any postprocessing and/or include information about the exact solution in the approximation scheme, such as Trefftz methods, are more attractive if the approximation of the coefficients is the main objective. The SFBIM, which was developed by Georgiou and co-workers [4] for Laplacian problems with a boundary singularity, also has this trait. In the SFBIM the solution is approximated by the leading terms of the expansion (1), i.e. N
uN
¦D
N j
W j r ,T ,
(3)
j 1
where D Nj are the approximate singular coefficients. The method has been successfully applied to a number of benchmark problems [5, 6]. In [7] it was shown theoretically that its convergence rate is exponential as N o f. The method has also been extended to biharmonic problems in two-dimensions arising from solid and fluid mechanics [8, 9]. The main advantages of the SFBIM are: x The dimension of the problem is reduced by one, leading to considerable computational savings x The singular coefficients are calculated directly, hence avoiding the need for post-processing In this work we present the method and its properties for a model two–dimensional Laplacian problem with a boundary singularity, including numerical experiments for illustration purposes. We then discuss its extensions to three-dimensional Laplacian problems with an edge singularity.
2. The SFBIM
2.1. Description of the method For concreteness, we consider the following Laplacian problem (see also Figure 1 below): Find u such that
2u
0 in : ,
wu wn
0 on S1 , u
0 on S 2 , u
f on S3 ,
Fig. 1. A planar Laplacian problem with a boundary singularity at O.
wu wn
g on S4 .
(4)
C. Xenophontos et. al./ Procedia Computer Science 00 (2010) 000–000
2585
C. Xenophontos et al. / Procedia Computer Science 1 (2010) 2583–2591
It is assumed that the data of the problem, f and g, are smooth and such that no other singularities arise. Now, by applying Galerkin’s principle we have
³³
:
W j 2uN
0 , j 1,..., N .
A double application of Green’s second identity, reduces the above integral to
³
w:
wW j wu N uN u N 2W j wn ³w: wn ³³:
Wj
0 , j 1,..., N .
Taking into account the boundary conditions of the problem, as well as the fact that the singular functions satisfy the governing PDE, we further obtain
³
S3
uN
wW j
wW j
³ uN
wn
wn
S4
³ Wi S3
wu N wn
³
S4
Wi g , j 1,..., N .
(5)
Letting OM { wu N / wn S and assuming that 3
M
¦O )
OM
A
, OA \ ,
A
(6)
A 1
where ) A , A 1,..., M are piecewise polynomials of degree p defined on a subdivision of S3 in elements of mesh width h; M is then proportional to 1/h. In the case when the ) A ’s are the typical (Lagrange based) basis functions used in the Finite Element Method, we have that the constants OA in (6) are the nodal values of OM . We then impose the boundary condition on S3 weakly, by means of Lagrange multipliers, i.e. we require that
³
S3
0 , A 1,..., M .
) A (u N f )
Equations (5) and (7) yield the following block system of (N + M)× (N + M) equations G G ª K1 K 2 º ªD º ª F1 º «K T 0 » « G» « G » , ¬ 2 ¼ N ¬O ¼ N ¬ F2 ¼ G G K
G
in which D
N 1
N N
T
G
ª¬D ,..., D º¼ , O
X
>O1 ,..., OM @
,
> K1 @i , j ³ §¨ W j ©
wWi · § wWi · Wj , i, j 1,..., N wn ¸¹ S³4 ¨© wn ¸¹
> K 2 @i ,A
A
G ª F1 º ¬ ¼i G ª¬ F2 º¼ A
³W ) i
(8)
b
T
S3
(7)
, i 1,..., N , A 1,..., M
S3
(9)
³W g , i i
1,..., N
S4
³ f)
A
, A 1,..., M
S3
Solving the system given by eq. (8) will produce the approximate singular coefficients as well as the discrete Lagrange multipliers. We emphasize that all integrals appearing above are one dimensional (thus the dimension of the problem is reduced by one) and are carried out on portions of the domain away from the singularity; therefore they can be evaluated by standard techniques, such as Gaussian quadrature. We should also mention that in the above system, the coefficient matrix K becomes singular when M > N, so the number of (discrete) Lagrange multipliers has to be smaller than the number of singular functions.
C. Xenophontos et. al./ Procedia Computer Science 00 (2010) 000–000
2586
C. Xenophontos et al. / Procedia Computer Science 1 (2010) 2583–2591
2.2. Error Analysis The SFBIM has been analyzed in the context of Laplacian problems and in this subsection we summarize the main results from [7]. Let H k (:) denote the usual Sobolev space of functions on ȍ whose (generalized) partial derivatives of order 0, 1, …, k are square integrable, and let || · ||k,ȍ denote the associated norms. We then define
^
`
H *1 { w H 1 : : w | S2 0 , H 1/ 2 (w:) { ^w H 1 : : w | w: H 0 (w:)` , and H –1/ 2 (w:) as the closure of H 0 (w:) with respect to the norm
w
³
w:
sup
1/ 2, w:
wv
v 1/ 2,w:
vH 1/2 ( w: )
.
Now, for any function w such that f
¦D W
w
j
j
(10)
,
j 1
wN rN , where
we can always write w
f
N
¦D W
wN
j
j
¦DW
, rN
j
j
(11)
.
j N 1
j 1
Under the assumption that there exist positive constants C1, C2 and ȕ(0, 1) such that
rN d C1 E N
,
wrN d C2 N E N , wn
(12)
it was shown in [7] that u uN where V1N
span ^W j `
N j 1
1, :
wu OM wn
1/ 2, S3
° wu d C ® infN u v 1,: infM P P V2 wn °¯vV1
½° ¾, ° 1, : ¿
(13)
and the space V2M is defined as follows: Let S3 be divided into sections īi , i = 1, …, n
such that n
S3
**
i
, hi
*i , h
max hi . 1d i d n
i 1
Then, with p ( I ) the set of polynomials of degree p on I, we set V2M
^O : O |
*i
`
p *i , i 1,..., n .
(14)
(We have N dim(V1N ), M M ( p, n) dim(V2M ) .) Using (13) it was also shown in [7] that if wu / wn H k (w:) for some k 1, then there exists a positive constant C, independent of N and M, such that
u uN
1, :
wu OM wn
dC 1/ 2, S3
where m = min{k, p + 1}. Moreover, since a j a Nj d u u N
^
`
N E N hm p k ,
0, :
(15)
, it follows that
a j a Nj d C E N , which shows that the method approximates the singular coefficients at an exponential rate as N o f.
(16)
C. Xenophontos et. al./ Procedia Computer Science 00 (2010) 000–000
2587
C. Xenophontos et al. / Procedia Computer Science 1 (2010) 2583–2591
3. Numerical results We now present the results of numerical computations for the test problem depicted graphically in Figure 2 below.
T DS wu wT
S3 0
O
S1
f (T )
u 2u
S2
0
T
u 0
0
Fig. 2. Test Laplacian problem on a sector.
The function f is taken as f T T T 2 / 2DS , with Įʌ being the angle of the circular sector. The local solution is given by f 2 j 1 . P (17) u ¦ D j r j sin P jT , P j 2D j 1 Since this problem can be solved analytically, we have that the exact singular coefficients are given by 16D
Dj
S 2R
Pj
2 j 1 3
(18)
,
where R is the radius of the sector. The numerical results that follow correspond to Į = R = 1 and our goal is to illustrate the method and its convergence rate. In figure 3 we show the percentage relative error in the solution u versus N, in a semi-log scale, for different values of the polynomial degree p used to approximate the Lagrange multipliers. We see that, independently of p, the method converges at an exponential rate, as predicted by equation (13). Error in the approximation of the solution u
1
10
p=1 p=2 p=3 100 u ||u - uSFBIM||1,: /||u||1,:
0
10
-1
10
-2
10
-3
10
6
8
10
12
14
16
ND
Fig. 3. Effect of the order p on the convergence of the solution u.
18
20
C. Xenophontos et. al./ Procedia Computer Science 00 (2010) 000–000
2588
C. Xenophontos et al. / Procedia Computer Science 1 (2010) 2583–2591
In figure 4 we show the percentage relative error in the first four singular coefficients versus N, in a semi-log scale, for the case when p = 1 (other values of p gave similar results). The exponential convergence, as predicted by equation (14), is again clearly visible. Error in jth coefficient, p = 1
2
10
j j j j
0
100 u |D j - D SFBIM | / |D j| j
10
= = = =
1 2 3 4
-2
10
-4
10
-6
10
-8
10
-10
10
10
15
20
25
30
35
ND
Fig. 4. Error in the first four singular coefficients for p = 1.
Finally, in figure 5 we show the error in wu / wn versus M, in a log-log scale, when wu / wn on the curved side of our domain is discretized by piece-wise polynomials of (fixed) degree p on a mesh with width h, i.e. M ~ p/h. The error estimate (13) states that the convergence rate is algebraic of order p, and indeed this is what figure 5 shows. Error in the approximation of the Lagrange multipliers
1
10
100 u max|O(Tk ) - Oh(Tk )|/|O(Tk )|
p=1 p=2 p=3
0
10
slope | - 1 -1
10
slope | - 3
-2
10
0
10
1
10 NO
Fig. 5. Error in wu / wn for different values of p.
slope | - 2 2
10
C. Xenophontos et. al./ Procedia Computer Science 00 (2010) 000–000
2589
C. Xenophontos et al. / Procedia Computer Science 1 (2010) 2583–2591
4. Extension to three-dimensions In this section we discuss the extension of the SFBIM to three-dimensional Laplacian problems, an area that, to our knowledge, still possesses several important open questions. In particular we consider the following: Find u such that 2u 0, in u g1 on u g 2 on u g3 on wu 0 on wn
½ ° ° °° ¾ ° ° S 4 S5 ° °¿
: S1 S2 S3
(19)
where the functions gi , i = 1, 2, 3 are given and the domain ȍ is shown in figure 6 below.
Fig. 6. Domain ȍ for the 3-D Laplace equation given by (19).
A singularity will arise along the edge AB, and we assume that the given data is smooth enough such that no other singularities are present. The difference between two and three-dimensions is substantial, since now the local solution is given by
§ ° ¨ ° Pj ¨ ®r f j T ¨ D j x3 1° ¨ ¨ ° © ¯
f
u
¦ j
§ 1 · ·½ r 2i ¨ ¸ ¸ ° © 4 ¹ ¸° , ¾ i n 1 n( P j n) ¸ ° ¸ ¸° ¹¿ i
f
2i
i 1
3
w ¦ wx D 2i
j
x3
(20)
where fj and ȝj are the eigenfunctions and eigenvalues, respectively, of the two-dimensional problem (posed on the face S1). Moreover, the singular coefficients Įj are no longer constants but now they are functions of the third coordinate x3; for this reason they are called Edge Stress Intensity Functions (ESIFs) [10, 11]. Nevertheless, they are known to be smooth functions [11], hence they can be approximated by polynomials, of say degree N, as N 1
D Nj x3
¦a
k 1 jk x3
(21)
k 1
where the coefficients a jk must be determined for each j. Once again, the solution u is approximated by the leading
C. Xenophontos et. al./ Procedia Computer Science 00 (2010) 000–000
2590
C. Xenophontos et al. / Procedia Computer Science 1 (2010) 2583–2591
terms of the asymptotic expansion (20), as in the two-dimensional case: N
¦ rP
uN
j
j 1
° ° f j T ®D Nj x3 ° ° ¯
f
2i
w ¦ wx D 2i
N j
3
i 1
x3
i § 1· ½ r 2i ¨ ¸ ° © 4¹ °. ¾ in 1 n( P j n) ° ° ¿
(22)
We note that the infinite sum in (22) will terminate after a finite number of terms, since D Nj x3 is a polynomial. Consequently, the number of unknowns in the above expression is N u (N + 1). To determine them we weigh the governing PDE by
Wˆ j l (r , T , x3 )
r
Pj
° ° f j T ® E lj x3 ° ° ¯
f
2i
i 1
3
w ¦ wx E 2i
l j
x3
i § 1· ½ r 2i ¨ ¸ ° © 4¹ °, ¾ in 1 n( P j n) ° ° ¿
(23)
where the functions E lj x3 are at our disposal – we may choose them as polynomials. It is easy to show that Wˆ j l satisfies the governing PDE and the boundary conditions on either side of the singular edge, independently of the choice for E lj x3 . As a result we have
³³³ Wˆ u l j
:
2
N
0 , j , l 1,..., N .
(24)
Using Green’s second identity twice, we get
³³³
:
u N 2Wˆ jl
wu Wˆ jl N w: wn
³³
³³
w:
uN
wWˆ jl wn
0 , j , l 1,..., N .
(25)
Since the functions Wˆ j l satisfy the PDE and the boundary conditions along S4 and S5, we have § wu wWˆ jl · ¨ N Wˆ jl u N ¸ S1 S 2 S3 ¨ wn wn ¸ © ¹
³³
0 , j , l 1,..., N ,
(26)
which is the analog of our 2-D method (cf. eq. (5)). The next step is to represent wu N / wn by a Lagrange multiplier function OMi on each face Si , i = 1, 2, 3 (expanded as a sum of polynomial basis functions ) A ) and to impose the Dirichlet conditions weakly, viz. M
¦O )
OMi
i A
A
, OAi \, i 1, 2,3 ,
(27)
A 1
³³
Si
) A gi u N
0 , A 1,..., M .
Then, eq. (26) becomes 3
§ wWˆ jl · i ˆl ¨ OM ¸ W j uN Si ¨ wn ¸ © ¹
¦ ³³ i 1
0 , j , l 1,..., N ,
(28)
C. Xenophontos et. al./ Procedia Computer Science 00 (2010) 000–000
C. Xenophontos et al. / Procedia Computer Science 1 (2010) 2583–2591
2591
G G G G G G which along with (27) yield a linear system analogous to (8), with unknown vector [a1 , a2 ,..., a N , O 1 , O 2 , O 3 ]T , where G G G a1 [a1,1 , a1,2 ,..., a1, N 1 ], a2 [a2,1 , a2,2 ,..., a2, N 1 ],..., aN [aN ,1 , aN ,2 ,..., aN , N 1 ] , G
G
1 O 1 [O11 , O21 ,..., OM ], O 2
G [O12 , O22 ,..., OM2 ], O 3
3 [O13 , O23 ,..., OM ].
Solving the system will produce the coefficients for the functions D Nj x3 (cf. (21)), as well as for the Lagrange i multipliers OM (cf. (27)).
Currently we are investigating possible choices for the polynomials E lj x3 and ) A . The results of this study and the implementation of the method will be reported in a future communication.
Acknowledgements EC was partially supported by the Cyprus Research Foundation, ȆǼȃǼȀ/ǼȃǿȈȋ/0308/014.
References 1. Z. C. Li, Combined Methods for Elliptic Equations with Singularities, Kluwer, The Netherlands, 1998. 2. I. Babuška and A. Miller, Int. J. Numer. Meth. Engg. 20 (1984) 1111. 3. B. Szabó and Z. Yosibash, Int. J. Numer. Meth. Engg. 39 (1996) 409. 4. G. Georgiou, L. Olson and Y. S. Smyrlis, Comm. Numer. Meth. Engg. 12 (1996) 127. 5. M. Elliotis, G. Georgiou and C. Xenophontos, Comm. Numer. Meth. Engg. 18 (2002) 213. 6. M. Elliotis, G. Georgiou and C. Xenophontos, Appl. Math. Comput. 169 (2005) 485. 7. C. Xenophontos, M. Elliotis and G. Georgiou, SIAM J. Sci. Comp. 28 (2006) 517. 8. M. Elliotis, G. Georgiou and C. Xenophontos, Eng. Anal. Bound. Elem. 30 (2006) 100. 9. M. Elliotis, G. Georgiou and C. Xenophontos, Int. J. Numer. Meth. Fluids 48 (2005) 1001. 10. Z. Yosibash, R. Actis and B. Szabó, Int. J. Numer. Meth. Fluids 53 (2002) 225. 11. Z. Yosibash, N. Omer, M. Costabel and M. Dauge, Int. J. Fracture 136 (2005) 37.