Analytical integration of singular kernels in symmetric boundary element analysis of Kirchhoff plates

Mirko  Mazza

An analytical integration method for evaluating the singular integrals arising in the construction of symmetric boundary element models is proposed, referring to the analysis of Kirchhoff plates. Kernels involved in the symmetric boundary formulation of Kirchhoff plates exhibit singularities up to O(1/r4). The proposed technique temporarily deactivates the singularities by using a limit approach and carries out the integration on contiguous boundary elements in order to smooth the singularities of the kernels through sufficiently continuous shape functions. The paper contains a brief presentation of the boundary integral formulation of the Kirchhoff plate, with the explicit expressions of the fundamental solutions involved and the derivation of the symmetric boundary element model. A detailed description of the proposed integration technique is presented, which discusses the main steps of the analytical procedure and provides the final results. The technique is illustrated through the discussion of a series of relevant cases, including high singularities and referring to both collinear and corner supports

INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING Int. J. Numer. Meth. Engng 2008; 76:127–155 Published online 18 March 2008 in Wiley InterScience (www.interscience.wiley.com). DOI: 10.1002/nme.2292 Analytical integration of singular kernels in symmetric boundary element analysis of Kirchhoff plates M. Mazza∗, † , L. Leonetti and M. Aristodemo Dipartimento di Modellistica per l’ingegneria, Università della Calabria, 87036, Rende, Cosenza, Italy SUMMARY An analytical integration method for evaluating the singular integrals arising in the construction of symmetric boundary element models is proposed, referring to the analysis of Kirchhoff plates. Kernels involved in the symmetric boundary formulation of Kirchhoff plates exhibit singularities up to O(1/r 4 ). The proposed technique temporarily deactivates the singularities by using a limit approach and carries out the integration on contiguous boundary elements in order to smooth the singularities of the kernels through sufficiently continuous shape functions. The paper contains a brief presentation of the boundary integral formulation of the Kirchhoff plate, with the explicit expressions of the fundamental solutions involved and the derivation of the symmetric boundary element model. A detailed description of the proposed integration technique is presented, which discusses the main steps of the analytical procedure and provides the final results. The technique is illustrated through the discussion of a series of relevant cases, including high singularities and referring to both collinear and corner supports. Copyright q 2008 John Wiley & Sons, Ltd. Received 23 February 2007; Revised 8 November 2007; Accepted 5 December 2007 KEY WORDS: boundary element method; symmetric formulation; Kirchhoff plates; analytical integration; singular kernels 1. INTRODUCTION The usual formulation of the boundary element method follows the collocation approach. The boundary integral equations are obtained weighting the differential equations with known singular fields and transferring the domain integrals to the boundary by the Gauss theorem [1]. Enforcing these equations at the collocation points and arranging the integral coefficients on the basis of the boundary conditions, non-symmetric systems are obtained. Their entries are computed from single ∗ Correspondence to: M. Mazza, Dipartimento di Modellistica per l’ingegneria, Università della Calabria, 87036, Rende, Cosenza, Italy. † E-mail: mirko.mazza@labmec.unical.it Contract/grant sponsor: Publishing Arts Research Council; contract/grant number: 98-1846389 Copyright q 2008 John Wiley & Sons, Ltd. 128 M. MAZZA, L. LEONETTI AND M. ARISTODEMO integrals of the product between the fundamental solution and the shape functions approximating the boundary field. Boundary formulations providing symmetric systems appear attractive in order to improve approximations and to control the convergence, having evident advantages in the analysis of dynamic problems and when boundary and finite elements have to be coupled. These formulations have been developed using different approaches. Variational formulations [2, 3], derivations from the Betti theorem [2] and weighting procedures of the standard boundary integral equations [4] have been proposed. All these approaches employ fundamental solutions associated with both static and kinematic sources and interpolate the auxiliary and actual boundary fields with the same shape functions, following the Galerkin procedure. The system entries result from a double boundary integration of the product between the shape functions and the fundamental solution. The presence of kernels associated with kinematic sources introduces higher singularities requiring a special treatment. To this end various analytical integration techniques have been proposed. The limit approach [5], the finite part of a divergent integral [6, 7] and the integration by parts [3] are examples of these. The analysis of Kirchhoff plates represents a demanding structural context for testing models based on boundary discretizations. The difficulties arise mainly from the high singularities involved in the kernels of the boundary integral equations. Many boundary models have been developed within the collocation approach, starting from the early works by Jaswon and Maiti [8], Forbes and Robinson [9], Altiero and Sikarskie [10], Bezine and Gamby [11], Hartmann and Zotemantel [12] and Stern [13]. Collocation models are usually based on boundary equations involving transversal displacement and normal slope, with kernels exhibiting singularities up to O(1/r 2 ). Different choices of boundary equations, involving the integral description of the bending moment and the Kirchhoff shear, have been proposed by Knöpke and by Frangi and Guiggiani, in conjunction with a regularization process using simple solutions [14] or a direct approach [15, 16]. Few symmetric boundary element models have been developed for Kirchhoff plates. The first works were by Hartmann et al. [4], Tottenham [17], Giroire and Nedelec [18] and Nazaret [19]. More recently Frangi and Bonnet introduced a symmetric boundary element model based on regularized boundary integral equations descending from the stationarity of an augmented potential energy functional [3]. A symmetric boundary integral formulation using boundary equations associated with distributed static and kinematic sources and corner sources has been developed by the authors [20]. This formulation gives rise to boundary integrals with kernels exhibiting orders of singularity ranging from O(lnr ) to O(1/r 4 ). The evaluation of singular integrals plays a crucial role for the actual construction of symmetric boundary models for Kirchhoff plates. The paper is devoted to this subject, proposing an analytical integration technique based on a limit approach [20, 21]. The method considers overlapping integration domains temporarily separated for deactivating the singularities of the kernels. In order to smooth these singularities, the integration is carried out on a pair of contiguous boundary elements forming a support and using sufficiently continuous shape functions. A limit process restores the true overlapping configurations. Starting from the differential form of the Kirchhoff plate problem, the boundary integral equations associated with static and kinematic sources are presented. All the fundamental solutions involved in these equations are provided in explicit form, illustrating out the reciprocity relationships between the kernels and the different orders of singularities. Afterwards the symmetric boundary integral formulation and the corresponding discrete model are briefly described. A complete description of the proposed integration technique is presented for the most relevant Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2008; 76:127–155 DOI: 10.1002/nme ANALYTICAL INTEGRATION OF SINGULAR KERNELS 129 cases, discussing the main steps of the analytical procedure and providing the final results for both collinear or including corner supports. 2. BOUNDARY INTEGRAL FORMULATION OF KIRCHHOFF PLATES The bending of a thin plate, defined over the domain delimited by the boundary and subjected to the transversal load p, is governed by the differential form: D w = p (1) where D is the flexural rigidity, is the Laplace operator and w is the transversal displacement [22]. Further boundary variables are the normal slope , the bending moment m and the Kirchhoff shear t: = w,n (2) m = −D(w,nn + w, tt ) (3) t = −D(w,nnn +(2−)w, ntt ) (4) where is the Poisson ratio and the comma denotes the differentiation with respect to the outward normal n or to the tangent t direction at the boundary points. The boundary fields meet the conditions t = t¯ w = w̄ or = ¯ or m = m̄ (5) (6) At the corners the twisting moment m t = −D(1−)w,nt (7) gives rise to the reactions ( j +) R ( j) = m t ( j +) ( j −) −m t (8) ( j −) and mark, respectively, the sides forward and backward to the where the superscripts corner. The Kirchhoff plate problem can also be described through boundary integral equations obtained weighting the field equation (1) by a function w ∗ and using the Gauss theorem to transfer the domain integrals onto the boundary [1]. Assuming as weighting functions the fundamental solutions w ∗F , associated with a unit force F, and wC∗ , associated with a unit couple C, the following transversal displacement and normal slope boundary integral equations are found [12]: nc ∗( j) kw + (t F∗ w +m ∗F ) d+ R F w ( j) = Copyright q (t w ∗F +m ∗F ) d+ 2008 John Wiley & Sons, Ltd. j=1 nc j=1 R ( j) ∗ ( j) wF + pw ∗F d (9) Int. J. Numer. Meth. Engng 2008; 76:127–155 DOI: 10.1002/nme 130 M. MAZZA, L. LEONETTI AND M. ARISTODEMO ∗ ) d+ (tC∗ w +m C ∗ (t wC∗ +m C ) d+ k+ = nc ∗ ( j) RC j=1 nc R ( j) j=1 w ( j) ∗ ( j) wC + p wC∗ d (10) where the factor k distinguishes sources located on the boundary (k = 12 ) from sources inside the domain (k = 1) and n c is the number of corners. The starred quantities are derived from w ∗F and wC∗ on the basis of relationships (2)–(4) and (7). Their expressions will be reported in the following section. It is worth noting that the presence of the corner reactions R ( j) implies the point variables w ( j) that are present in the boundary integral formulation. These point variables are known or unknown according to the boundary conditions assigned for the edges merging at the corner. In order to develop symmetric boundary models, further boundary integral equations have to be considered. Assuming as test functions the fundamental solutions corresponding to a unit slope discontinuity , a unit tangent rotation discontinuity t and a unit transversal displacement discontinuity w, the following bending moment, twisting moment and pure shear, q, boundary integral equations can be introduced [20]: nc ∗( j) ∗ km + (tw +m∗ ) d+ R ( j) w = km t + ∗ (t w +m ∗ ) d+ = (t∗ w w +m ∗w ) d+ j=1 ∗ (t w +m ∗t ) d+ t (tw∗ w +m∗w ) d+ nc ∗ (tw +m∗t ) d+ t kq + = j=1 ∗( j) R w ( j) − nc j=1 nc j=1 nc j=1 nc j=1 ∗ pw d (11) ∗ pw d t (12) p w∗ w d (13) ∗ R ( j) w t ∗ R w ( j) − t ∗( j) R ( j) ww ∗ ( j) Rw w ( j) − It is worth noting that Equation (12) expressed for the ends of two sides merging at the corner j provides a boundary integral description of the corner reaction R ( j) , whereas its tangent derivative added to Equation (13) provides the boundary integral description of the Kirchhoff shear. 3. FUNDAMENTAL SOLUTIONS The boundary integral equations presented in Section 2 involve fundamental solutions representing the effects produced by unit actions on the unbounded plate. Once the fundamental solution and Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2008; 76:127–155 DOI: 10.1002/nme ANALYTICAL INTEGRATION OF SINGULAR KERNELS 131 the derived quantities corresponding to a unit point force have been defined, the kernels associated with different actions can be obtained by using the Betti theorem. The equilibrium of a circular neighbor around a transversal unit force F, acting at the source point x, provides the transversal displacement w ∗F at the field point y: w ∗F [x, y, −, −] = 1 2 r ln |r | 8D (14) where r = y − x. In order to explain the reciprocity relationships involving the fundamental solutions, a specific notation has been introduced. The four quantities in square brackets denote the dependance from the vectors x and y and, if present, on the normal vectors n and at these points. Successive differentiations of w ∗F with respect to the normal and tangent directions, and , at the point y (Figure 1), produce the quantities ∗F , m ∗F , t F∗ and m ∗t F : 1 {rr (1+2 ln |r |)} 8D 1 m ∗F [x, y, −, ] = − {2(1+) ln(r )+(3+)r2 +(1+3)r2 } 8 1 {2r +(1−)r (r2 −r2 )} t F∗ [x, y, −, ] = − 4r 1 m ∗t F [x, y, −, ] = − {(1−)rr } 4 ∗F [x, y, −, ] = (15) (16) (17) (18) where r and r are directional derivatives defined by the expressions r = r1 1 +r2 2 (19) r = r1 1 +r2 2 = −r1 2 +r2 1 (20) Figure 1. Source and field points. Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2008; 76:127–155 DOI: 10.1002/nme 132 M. MAZZA, L. LEONETTI AND M. ARISTODEMO containing the components of the unit vectors , and the derivatives ri = (yi − xi )/r of the distance vector r with respect to the cartesian coordinates yi where i = 1, 2 [1]. The expressions obtained for the source F can be used to evaluate other fundamental solutions by the Betti theorem. For instance, the quantity wC∗ is derived from the quantity ∗F exchanging the points y and x and vectors n and [1]: wC∗ [x, y, n, −] = ∗F [y, x, , −] (21) Using formulae (2)–(4) and (7) the remainder quantities associated with the source C can be obtained ∗ C [x, y, n, −] = 1 {2(rrn +rt r ) ln |r |+3rrn +rt r } 8D ∗ [x, y, n, −] = − mC tC∗ [x, y, n, −] = − m ∗t C [x, y, n, −] = − (22) 1 {(1+)rn +2(1−)rrrt } 4r (23) 1 {(3−−2(1−)r2 )(rrt −rrn )+4(1−)r2rrt } 4r 2 (24) 1− 2 2 (r −r )rt 4r (25) containing the new directional derivatives rn = −r1 n 1 −r2 n 2 (26) rt = r1 n 2 −r2 n 1 (27) Following a similar path, the expressions related to the source are evaluated [14]: ∗ w [x, y, n, ] = m ∗F [y, x, , n] (28) ∗ ∗ [x, y, n, ] = m C {y, x, , n] (29) m ∗ [x, y, n, ] = − D {(1−2 )(r2 −r2 )+2(1−)(r2 +r2 )(rn2 −rt2 ) 4r 2 + 4(1−)2rrrn rt } ∗ t [x, y, n, ] = (30) D {(2(1−)(3−−2(1−)r2 )((rn2 −rt2 )r −2rrn rt ) 4r 3 + 8(1−)rr (r (rt2 +rn2 )−(1−)rrn rt ) m ∗t [x, y, n, ] = Copyright q + 4(1−)r ((r2 −r2 )(rt2 +rn2 )−2(1−)rrrn rt )} (31) D (1−){(1−)(r2 −r2 )rn rt −2rr (rt2 +rn2 )} 2r 2 (32) 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2008; 76:127–155 DOI: 10.1002/nme ANALYTICAL INTEGRATION OF SINGULAR KERNELS 133 The relationships associated with the source w are [21] w∗ w [x, y, n, ] = t F∗ [y, x, , n] (33) ∗w [x, y, n, ] = tC∗ [y, x, , n] (34) ∗ m ∗w [x, y, n, ] = t [y, x, , n] (35) t∗ w [x, y, n, ] = D {3(−1+)(rn r3 (−2+(−5+3)rn2 +(13−7)rt2 ) 2r 4 +r2rt (6+(27−21)rn2 +(−7+5)rt2 )r −rn r (−6+(−5+7)rn2 +3(3−5)rt2 )r2 −rt (2−(1+5)rn2 +(1+)rt2 )r3 ))} m ∗t w [x, y, n, ] = − (36) D {(−1+)(r2rt (2+(−1+)(rt2 −5rn2 ) 2r 3 − 4rn rr (−1+(−1+)(rn2 −2rt2 )) +rt r2 (−2+(−1+)(5rn2 −rt2 ))} (37) Finally, the fundamental solutions for the source t are expressed by the following relationships [20]: w∗ t [x, y, −, ] = m ∗t F [y, x, n, −] (38) ∗t [x, y, n, ] = m ∗t C [y, x, , n] (39) m ∗t [x, y, n, ] = m ∗t [y, x, , n] (40) ∗ t [x, y, n, ] = m ∗t w [y, x, , n] t (41) m ∗t t [x, y, −, ] = − D(−1+)2 {−4rn rt rr +(rn2 −rt2 )(r2 −r2 )} 4r 2 (42) All the kernels involved in the analysis, except for w ∗F , ∗F , m ∗t F and wC∗ , turn out to be singular functions of the distance r . A summary of the order of singularity characterizing the different kernels is shown in Table I, which also reports the existing reciprocity relationships. Table I shows that the kernels associated with kinematic sources exhibit higher order singularities. In order to exploit the reciprocity relationships connecting the kernels, the procedure for constructing symmetric boundary models has to employ all the presented boundary integral equations. Consequently, the treatment of more singular kernels makes the computation of the boundary coefficients involved in symmetric models heavy. Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2008; 76:127–155 DOI: 10.1002/nme 134 M. MAZZA, L. LEONETTI AND M. ARISTODEMO Table I. Orders of singularity and reciprocity relationships. F C w t w∗ ∗ t∗ m∗ m ∗t r 2 ln(r ) r ln(r ) ln(r ) 1/r 1/r 2 1/r 4 ln(r ) 1/r 1/r 3 1/r 2 — 1/r 1/r 3 1/r 2 1/r 2 (reciprocal) 4. BOUNDARY FORMULATION FOR SYMMETRIC MODELS A boundary formulation for Kirchhoff plates suitable for symmetric models can be introduced considering two states associated with their reciprocal work. The first is the actual state, constituted by the boundary fields w, , m, t and the corner variables w ( j) and R ( j) . The second is an auxiliary ˆ ŵ and the corner contributions R̂ ( j) and ŵ ( j) . The state described by the boundary fields tˆ, m̂, , Betti theorem links these states by the integral equation: (w tˆ + m̂)d+ nc R̂ ( j) w ( j) = j=1 nc ˆ d+ R ( j) ŵ ( j) (t ŵ +m ) (43) j=1 which can be written in the symbolic form: gT f d = 0 (44) where the vector g g = {tˆ m̂ − ŵ − ˆ R̂ ( j) − ŵ ( j) }T , j = 1, . . . , n c (45) collects the variables of the auxiliary boundary state, constituted by functions and corner values, while the vector f f = {w t m w ( j) R ( j) }T , j = 1, . . . , n c (46) contains the actual boundary variables. The actual state can be represented by the integral equations introduced in Section 2: f[x] = ∗ K [x, y] f[y] d y + Copyright q 2008 John Wiley & Sons, Ltd. w∗ [x, y] p[y] d (47) Int. J. Numer. Meth. Engng 2008; 76:127–155 DOI: 10.1002/nme 135 ANALYTICAL INTEGRATION OF SINGULAR KERNELS where the matrix K∗ collects all the involved kernels, including the terms corresponding to the point variables associated with the corners −t F∗ [x, y] ⎢ ∗ ⎢ −t [x, y] ⎢ C ⎢ ⎢ t ∗ [x, y] ⎢ w K∗ = ⎢ ⎢ ∗ ⎢ t [x, y] ⎢ ⎢ ∗ i ⎢ −t [x , y] ⎣ F ⎡ ∗ [x i , y] t t −m ∗F [x, y] w∗F [x, y] ∗F [x, y] −R ∗F [x, y j ] ∗ −m C [x, y] wC∗ [x, y] ∗ C [x, y] −RC∗ [x, y j ] m ∗w [x, y] −w∗ w [x, y] −∗w [x, y] R∗ w [x, y j ] m ∗ [x, y] ∗ −w [x, y] −∗ [x, y] ∗ R [x, y j ] −m ∗F [x i , y] w∗F [x i , y] ∗F [x i , y] −R ∗F [x i , y j ] m ∗t [x i , y] ∗ [x i , y] −w t −∗t [x i , y] ∗ [x i , y j ] R t w∗F [x, y j ] ⎤ ⎥ wC∗ [x, y j ] ⎥ ⎥ ⎥ ∗ j ⎥ −ww [x, y ] ⎥ ⎥ ⎥ ∗ −w [x, y j ] ⎥ ⎥ ⎥ w∗F [x i , y j ] ⎥ ⎦ (48) ∗ [x i , y j ] −w t and the vector w∗ coincides with the third column of this matrix. It is worth noting that the kernels associated with the tangent rotation discontinuity t correspond to a pair of point sources applied to the ends of the two sides merging at the corner j. By introducing Equation (47) into expression (44), the following integral relationship is derived: (49) gT [x] K∗ [x, y]f[y] d y dx = − gT [x] w∗ [x, y] p[y] d dx where the subscripts appearing in dx and d y specify the integration variables. This equation represents the starting point for constructing a symmetric boundary element model. 5. SYMMETRIC BOUNDARY SYSTEM The symmetric system can be generated by introducing a boundary element discretization, approximating the boundary fields in the form f = W[y]f̄ (50) g = W[x]ḡ (51) where f̄ and ḡ are the vectors that collect the vector of parameters used to interpolate the boundary fields, printed in boldface, and the corner values: f̄ = {w h t m w ( j) R ( j) }T , j = 1, . . . , n c ḡ = {t̂ m̂ − ŵ − ĥ R̂ ( j) − ŵ ( j) }T , j = 1, . . . , n c (52) (53) and W represents the row matrix that contains the shape functions, identical for auxiliary and actual fields associated in the work sense, and the unit values corresponding to the corner variables. Introducing these interpolations, Equation (49) can be rewritten in the discrete form: ḡT Kf̄ = ḡT p (54) The symmetric system is extracted from the above equation on the basis of the prescribed boundary conditions. The coefficients of the system matrix are taken from the columns of K corresponding to the unknown part of vector f, whereas the other columns and the vector p form Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2008; 76:127–155 DOI: 10.1002/nme 136 M. MAZZA, L. LEONETTI AND M. ARISTODEMO the known term vector. This assembling procedure, based on the use of reciprocal kernels, ensures symmetric coefficients. In order to classify the types of integrals involved in the boundary system, matrix K and vector p can be partitioned as follows: K= A(n×n) B(n×n c ) C(n c ×n) D(n c ×n c ) , p= l(n×1) m(n c ×1) (55) where the dimensions of the sub-matrices are described by the number of interpolation parameters n and corners n c . The entries of these matrices are formed by contributions belonging to the following types: A: [x]k ∗ [x, y][y] d y dx (56) x B: y [x]k ∗ [x, y j ] dx (57) [y]k ∗ [x i , y] d y (58) x C: y D : k ∗ [x i , y j ] l: p[y] k ∗ [x, y] d dx x m: (59) (60) p[y]k ∗ [x i , y] d (61) where the functions [x] and [y], grouped in the matrix W, interpolate the source and field distributions and k ∗ denotes the generic kernel. The coefficients of the boundary system are derived from the above types by selecting the appropriate shape functions and the kernel on the basis of the boundary conditions. The evaluation of such integrals poses different problems in relation to the order of singularity characterizing the fundamental solution k ∗ . In particular, the integrals (60) and (61), owing to the singularity of k ∗ being less than or equal to O(1/r ), can be evaluated easily by numerical techniques or transferring them to the boundary. Standard numerical quadratures can be used for evaluating the integrals (56), (57) and (58) when the source and the field points are separate, whereas specific integration techniques have to be used for overlapping integration domains or in the case of a source (field) point located on the integration domain y (x ). The remainder of this paper deals specifically with the evaluation of these integrals. As can be shown in Section 6.2, when corner contributions have to be evaluated, singularities of the type (59) are eliminated by the presence of opposite contributions arising from integrals (56), (57) and (58). 6. ANALYTICAL INTEGRATION OF SINGULAR CONTRIBUTIONS Various techniques have been proposed for integrating singular kernels. Some procedures transform the boundary integral equations through integrations by parts [23] or Stokes’ theorem [24] in order Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2008; 76:127–155 DOI: 10.1002/nme ANALYTICAL INTEGRATION OF SINGULAR KERNELS 137 to reduce the original order of singularity. The concept of the finite part of a divergent integral, introduced by Hadamard [6], has also been used to regularize the integral cutting off the boundary around the source point where a first-order Taylor expansion of the singular kernel is assumed [7]. A different approach considers the overlapping situation as the limit configuration of separate integration domains [5]. Few proposals regard the integration of the highly singular kernels arising in the context of symmetric boundary element formulation of Kirchhoff plates [3]. Here an analytical technique based on the limit approach is presented for these problems [20, 21]. The features of the proposed limit integration process are the following. A boundary coefficient constituted by double integrals is evaluated carrying out the inner integration over two contiguous boundary elements forming a support and moving the source point at a distance from the integration domain. This procedure allows the continuity of the shape functions [y], making the kernel temporarily regular, to be exploited. Once the integration has been completed, the true configuration is restored by taking the limit → 0. Singular terms located at the ends of the support are neutralized by the zero values of the hat shape functions at these points. Moreover some of the singular terms, located where the boundary elements merge, are neutralized adding the contributions of these boundary elements in order to exploit the inter-element continuity of the shape functions [y] and their derivatives. A similar procedure is used to evaluate the outer integral. In this case the shape function [x] and the inner integration extremes are defined in terms of a local abscissa s, which is displaced at a distance from the outer integration extremes. A limit for → 0 restores the true integration domain. The computation of a single integral corresponds to one step (inner or outer) of the double integration. The order of continuity of the hat interpolation functions has to be sufficient for a continuous description of the boundary variables and for smoothing the singularities affecting the kernels. To this end the integration of the kernels work-conjugated with the boundary variables , m and t can be made using the following interpolation laws: (i) [ ] = [ ]¯ , m[ ] = m [ ]m̄ (i) , t[ ] = t [ ]t¯(i) (62) where , m and t are the linear shape functions defined over the support of the node i (Figure 2): [ ] = m [ ] = t [ ] = 1+ , ∈ [−1, 0] 1− , ∈ [0, +1] (63) In order to ensure the continuity conditions at corners and the smoothing effect on singular kernels, the transversal displacement is interpolated separating the contributions of the transversal displacement parameter ū and the tangential rotation parameter ¯ t : w[ ] = w [ ]w̄ (i) = u [ ]ū (i) +t [ ]¯t (i) (64) Figure 2. Linear shape functions over contiguous boundary elements. Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2008; 76:127–155 DOI: 10.1002/nme 138 M. MAZZA, L. LEONETTI AND M. ARISTODEMO where the interpolations u and t are cubic functions defined over the support of the node i (Figure 3): u [ ] = t [ ] = 1−3 2 −2 3 , 1−3 2 +2 , +2 2 + 3 , ∈ [−1, 0] −2 2 + 3 , ∈ [0, +1] 3 ∈ [−1, 0] (65) ∈ [0, +1] (66) The evaluation of all the coefficients containing singular kernels requires the analysis of an extended set of cases characterized by different shape functions, kernels and geometrical situations. The following subsections contain a detailed description of the proposed integration procedure for some significant cases selected to illustrate the computational steps involved in the analysis. These boundary coefficients involve the highest singularities and refer to partially and totally overlapping domains for both the cases of collinear supports and corner situations. 6.1. Integration over collinear supports This situation occurs when the supports of the shape functions are collinear and partially or totally overlapping. As a first case, the evaluation of the coefficient on a boundary free of supports is considered. In particular, the coefficient corresponding to the source interpolation parameter û (i) , (i+1) associated with the node i, and to the field interpolation parameter t , associated with the subsequent node i +1, is obtained developing the double integral (i) (i+1) C[û , t t∗ w [x, y]t [y] dy dx (67) w [x] ]= (i) (i+1) The source support (i) is constituted by the pair of elements a and b merging at the node i and the field support (i+1) by the parts b and c merging at the node i +1 (Figure 4). The above Figure 3. Cubic shape functions over contiguous boundary elements. Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2008; 76:127–155 DOI: 10.1002/nme 139 ANALYTICAL INTEGRATION OF SINGULAR KERNELS boundary coefficient is given by the sum of the following four contributions: (i+1) C[û (i) , t ] = Iab + Iac + Ibb + Ibc = (a) u [x] a + b (b) t∗ w [x, y]t [y] dy dx + b (b) u [x] b (a) u [x] a t∗ w [x, y](b) t [y] dy dx + b (c) t∗ w [x, y]t [y] dy dx c (b) u [x] c t∗ w [x, y](c) t [y] dy dx (68) where the superscripts specify the domains of the shape functions. The analytical integration is used to evaluate the contribution Ibb , corresponding to overlapping elements, and the terms Iab and Ibc over adjacent elements, whereas a standard numerical quadrature can be used to evaluate the contribution Iac , which refers to separate elements. The inner integration starts by moving the source support away from the field support in order to temporarily deactivate the singularity of the kernel. A local reference system ( , ) is introduced over the elements b and c of the field distribution, assuming its origin O varies with the local abscissa s used to describe the source distribution over the elements a and b of the source distribution (Figure 4). Making use of these reference systems, the contributions Ibb , Iab and Ibc are rewritten in the form s2 2 (b) (a) [s]· lim t∗ w [ , ]t [ ] d ds (69) Iab = u →0 s1 Ibb = Ibc = s2 s1 s2 s1 (b) u [s]· lim →0 (b) u [s]· lim →0 1 2 (b) (70) (c) (71) t∗ w [ , ]t [ ] d ds 1 2 t∗ w [ , ]t [ ] d ds 1 (i+1) Figure 4. Shape functions associated with the source û (i) and the field t Copyright q 2008 John Wiley & Sons, Ltd. interpolation parameters. Int. J. Numer. Meth. Engng 2008; 76:127–155 DOI: 10.1002/nme 140 M. MAZZA, L. LEONETTI AND M. ARISTODEMO The shape functions t [ ] over the elements b and c are expressed in terms of the ends of the inner integration domains by the formulas: (b) t [ ] = (c) t [ ] = ( − 1 )2 ( − ( 1 − 2 )2 2) 1, 2 (72) ( − 1 )( − 2 )2 ( 1 − 2 )2 (73) The fundamental solution t∗ w [ , ] takes the form t∗ w [ , ] = D 2( 2 +2 )5 {3(−3+2+2 ) 6 −5(1−6+52 ) 4 2 + 5(11−18+72 ) 2 4 +(−7+10−32 )6 } (74) obtained from expression (36) writing the distance r , their derivatives r1 and r2 and the outward normals n and as follows: 2 +2 , r= r1 = , r r2 = − , r n = = (0, −1) (75) Once the integration on has been carried out, the contributions Iab , Ibb and Ibc become s2 5 (a) 2h Iab = u [s]· lim f h [ 1 , 2 , ] ds (76) →0 h=0 (b) u [s]· lim →0 5 s1 Ibb = Ibc = s2 s1 s2 s1 (b) u [s]· lim →0 2h fh [ 1, 2 , ] ds (77) 2 , ] ds (78) h=0 5 2h fˆh [ 1 , h=0 The unique terms of the above sums, which remain after the limit on , are the functions f0[ 1, 2 , ] = k1 2 2 3 3 4 ( 3 −6 2 +3 2 2 2 2 1 2 +2 2 +3 1 2 ln(( 1 + )/( 2 + ))) 1 2 1 1 2 ( 1 − 2 )2 ( 21 +2 )3 ( 22 +2 )2 (79) with k1 = D (−1+)(3+) 4 (80) and the function fˆ0 [ 1 , 2 , ], obtained from f 0 [ 1 , 2 , ], exchanging 1 with 2 while the terms corresponding to h = 1, . . . , 5 disappear once the limit → 0 is performed to restore the true location of the source support. The outer integration requires the inner integration ends 1 and 2 be expressed in terms of the local abscissa s by using the quantities summarized in Table II, where the quantities to be used in Iab (s<0) and in Ibb , Ibc (s>0) are identified. Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2008; 76:127–155 DOI: 10.1002/nme 141 ANALYTICAL INTEGRATION OF SINGULAR KERNELS Table II. Inner integration ends 1 2 1 and 2 in terms of the local abscissa s. Iab Ibb Ibc −s ℓb −s −s ℓb −s ℓb −s ℓb +ℓc −s The integrals (69), (70) and (71) become Iab = s2 s1 Ibb = Ibc = s2 s1 s2 s1 (a) u [s]g1 [s] ds (81) (b) u [s]g1 [s] ds (82) (b) u [s]g2 [s] ds (83) where the shape functions u [s] are defined as (a) u [s] = (b) u [s] = (ℓa −2s)(ℓa +s)2 ℓa3 (ℓb −s)2 (ℓb +2s) ℓ3b (s<0) (84) (s>0) (85) and the functions g1 [s] and g2 [s] have the form g1 [s] = k1 g2 [s] = k1 6 6 1 4 ℓb −s + − ln − − s (ℓb −s)2 s(ℓb −s) ℓb s ℓ2b (86) 2 6 6 1 ℓb +ℓc −s + + ln − ℓb −s (ℓb −s)2 ℓc (ℓb −s) (ℓb −s)(ℓb +ℓc −s) ℓ2c (87) Various singular coefficients appear in the quantities (86) and (87) defining the contributions Iab , Ibb and Ibc . Some of them, like 1/(ℓb −s)2 involved in Ibb and Ibc , disappear when these integrals, which have the same shape functions u(b) [s], are added together. The weakly singular terms are neutralized by the integration process. The integration of the remainder singular terms is performed by deactivating the singularities with a temporary change of the integration ends s1 and s2 of a quantity , as detailed in Table III. The integrals Iab , Ibb and Ibc are evaluated restoring the true domain by the limit → 0: k1 k1 Iab + Ibb = lim g̃1 []− 3 ln() + g̃2 []+ 3 ln() (88) →0 2ℓb 2ℓb Ibc = lim g̃3 [] →0 Copyright q 2008 John Wiley & Sons, Ltd. (89) Int. J. Numer. Meth. Engng 2008; 76:127–155 DOI: 10.1002/nme 142 M. MAZZA, L. LEONETTI AND M. ARISTODEMO Table III. Outer integration ends s1 and s2 in terms of the quantity . s1 s2 Iab Ibb Ibc −ℓa − ℓb − ℓb − where g̃i are regular functions of . Summing the contributions, the explicit expression of the required coefficient is obtained (i+1) C[û (i) , t ] = Iab + Ibb + Ibc = D 8 ℓa3 ℓ3b ℓ2c {−(−1 + )(3 + )(2ℓa3 ℓb ℓ2c (3ℓa +2ℓb ) ln(ℓa ) − 2ℓ4b (−3ℓa3 +ℓb ℓ3c ) ln(ℓb )+ℓb ℓc (6ℓa3 ℓ2b +ℓa ℓ2b ℓc (ℓa −2ℓb )−2ℓa3 ℓ2c + 2ℓc (−3ℓa4 −2ℓa3 ℓb +ℓ4b ) ln(ℓa +ℓb ))−2ℓa3 ℓ4c ln(ℓc ) + 2ℓa3 (−3ℓ4b +ℓ4c ) ln(ℓb +ℓc ))} (90) which has to be completed by adding the value resulting from the numerical integration of the non-singular term Iac . A slightly different situation arises when the source and the field parameters are located at the same node. As an example, referring to the evaluation of the boundary coefficient corresponding to the source parameter û (i) and to the field parameter u (i) (Figure 5): (i) (i) C[û , u ] = t∗ w [x, y]u [y] dy dx (91) u [x] (i) (i) the subdivision of the integration domains, merging at the same node i, gives the contributions C[û (i) , u (i) ] = Iaa + Iab + Iba + Ibb (a) ∗ t [x, y] [y] dy dx + [x] (a) = u u w + a a a b (b) u [x] (a) u [x] a t∗ w [x, y](a) u [y] dy dx + b t∗ w [x, y](b) u [y] dy dx b (b) u [x] b t∗ w [x, y](b) u [y] dy dx (92) where all the contributions have to be handled with the limit approach. Introducing the local reference system ( , ) and the local abscissa s, the above double integrals become s2 2 Iaa = (a) [s]· lim t∗ w [ , ] (a) (93) u u [ ] d ds →0 s1 Iab = Copyright q s2 s1 (a) u [s]· lim 2008 John Wiley & Sons, Ltd. →0 1 2 t∗ w [ , ] (b) u [ ] d ds (94) 1 Int. J. Numer. Meth. Engng 2008; 76:127–155 DOI: 10.1002/nme ANALYTICAL INTEGRATION OF SINGULAR KERNELS 143 Figure 5. Shape functions associated with the source û (i) and the field u (i) interpolation parameters. Iba = Ibb = s2 s1 s2 s1 (b) u [s]· lim →0 (b) u [s]· lim →0 2 t∗ w [ , ] (a) u [ ] d ds (95) t∗ w [ , ] (b) u [ ] d ds (96) 1 2 1 (b) where the shape functions (a) u [ ] and u [ ] are defined as (a) u [ ]= (b) u [ ]= ( − 1) 2 (2 ( ( − 2) + 1− 2 1 −3 2 ) 3 ) 2 (2 (− −3 1 + 3 + 1 2) 2) (97) (98) After the inner integration, the limit for → 0 is taken and the integration ends 1 and 2 are expressed in terms of the local abscissa s, as detailed in Table IV. Moreover the shape functions describing the distribution of the auxiliary state are expressed in the form: (a) u [s] = (b) u [s] = (ℓa −2s)(ℓa +s)2 ℓa3 (ℓb −s)2 (ℓb +2s) ℓ3b (s<0) (99) (s>0) (100) and the outer integration ends are expressed in terms of the quantity , which moves the singularities away from the ends of the integration domains (Table V). Functions g1 [s] and g2 [s] contain various singular terms. Some of them are eliminated adding the contributions, which present the same shape functions. Weakly singular terms are regularized Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2008; 76:127–155 DOI: 10.1002/nme 144 M. MAZZA, L. LEONETTI AND M. ARISTODEMO Table IV. Inner integration ends 1 2 1 and 2 in terms of the local abscissa s. Iaa Iab Iba Ibb −ℓa −s −s −s ℓb −s −ℓa −s −s −s ℓb −s Table V. Outer integration ends s1 and s2 in terms of the quantity . s1 s2 Iaa Iab Iba Ibb −ℓa + − −ℓa − ℓb ℓb − by the integration process. The remainder singularities disappear adding together the results of the integration over the contiguous domains. Taking the limit for → 0 for all the four contributions, the process furnishes the result C[û (i) , u (i) ] = Iaa + Iab + Iba + Ibb D = 4ℓa3 ℓ3b {−3(−1+)(3+)(ℓa +ℓb )(ℓa ℓb (ℓa +ℓb ) +4ℓ3b ln(ℓb )+4ℓa3 ln(ℓa )−4(ℓa3 +ℓ3b ) ln(ℓa +ℓb ))} (101) 6.2. Integration over supports including corners In the presence of supports of shape functions including corner points, the evaluation of a boundary coefficient can involve one or more types of kernels because the continuity requirements at the corner can link different types of variables on the merging sides. For the sake of simplicity, the following discussion will be restricted to corners between orthogonal elements. In order to explain the role of the continuity conditions at the corner, at first the evaluation of the (i) coefficient corresponding to the source ˆ and the field (i) interpolation parameters (Figure 6) is analyzed. In particular, the continuity links the normal rotation on the side a with the tangent rotation on the side b. This case arises when both the sides merging at the corner are hinged. The expression of this boundary coefficient can be split as shown before (i) C[ˆ , (i) ] = Iaa + Iab + Iba + Ibb (a) (a) m ∗ [x, y] [y] dy dx + [x] = + Copyright q a a a (b) b t [x] (a) [x] a (a) m ∗w [x, y] [y] dy dx + 2008 John Wiley & Sons, Ltd. b (b) (b) ∗ t [x, y] t [y] dy dx b t [x] b (b) t∗ w [x, y] t [y] dy dx (102) Int. J. Numer. Meth. Engng 2008; 76:127–155 DOI: 10.1002/nme 145 ANALYTICAL INTEGRATION OF SINGULAR KERNELS (i) Figure 6. Shape functions associated with the source ˆ and field (i) interpolation parameters. Introducing the local reference systems and separating the boundary configurations which imply singularities, the integrals appearing in (102) can be expressed as s2 2 (a) (a) Iaa = [s]· lim m ∗ [ , ] [ ] d ds (103) →0 s1 Iab = Iba = Ibb = s2 s1 s2 s1 s2 s1 (a) [s] (b) t [s] 2 1 (b) ∗ [ , ] t [ ] d ds t (104) 1 2 (a) m ∗w [ , ] [ ] d ds (105) (106) 1 (b) t [s]· lim →0 2 (b) t∗ w [ , ] t [ ] d ds 1 where the shape functions u [ ] are defined by the relationships: (a) u [ ]=− (b) u [ ]= ( − 1) ( 1 − 2) (107) ( − 1 )( − 2 )2 ( 1 − 2 )2 (108) It is worth noting that, in contrast to the case of collinear domains, the limit for → 0 regards only the contributions of the overlapping parts of the integration domains. As discussed previously, the inner integration ends 1 and 2 are expressed in terms of the local abscissa s by using the relationships summarized in Table VI which also shows the dependence of on the variable s (Figure 7). Using the relationships of Table VI, the integral contributions can be rewritten as s2 s2 Iaa = (a) (109) (a) [s] f [s] ds, I = 1 ab u [s] f 2 [s] ds u s1 s1 Iba = s2 s1 Copyright q (b) u [s] f 3 [s] ds, 2008 John Wiley & Sons, Ltd. Ibb = s2 s1 (b) u [s] f 4 [s] ds (110) Int. J. Numer. Meth. Engng 2008; 76:127–155 DOI: 10.1002/nme 146 M. MAZZA, L. LEONETTI AND M. ARISTODEMO where the functions gi resulting from the inner integrals have the form k2 k2 ln(−ℓa −s)2 g1 [s] = k1 − + ln s 2ℓa s2 9s(1+3) s2 k2 3s(−1+21) 9(1+3) + + g2 [s] = k1 + s ln 2 2 2 2 s 2ℓb 2ℓb ℓb +s 2(ℓ2b +s 2 ) 3s 3 (1+11) 96 ℓb s2 s4 + 2 2 − 2 − 3 arctan 2 2 2 s 2ℓb (ℓb +s ) (ℓb +s ) ℓb ℓb (111) (112) k2 2ℓ2 (1−) k2 g3 [s] = k1 − 2 + 2 2 + 2a 2 2 s ℓa +s (ℓa +s ) Table VI. Integration ends 1 2 1, 2 (113) and local variable in terms of the abscissa s. Iaa Iab Iba Ibb −ℓa −s −s — 0 ℓb −s −ℓa 0 s −s ℓb −s — Figure 7. Different locations of the reference system , and of the abscissa s. Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2008; 76:127–155 DOI: 10.1002/nme 147 ANALYTICAL INTEGRATION OF SINGULAR KERNELS 33ℓb s 2 −40s 3 (9ℓ2b +25ℓb s −40s 2 ) k2 g4 [s] = k1 2 +k2 +k 2 s 2ℓ3b s 2 2ℓ3b (ℓb −s) s2 3k2 (3ℓb −8s) ln + (ℓb −s)2 ℓ3b (114) with k1 = −1+ , 4 k2 = 3+ (115) It is worth noting that some of the singular contributions arising from the inner integration are eliminated by the continuity at the corner of the shape functions approximating the boundary fields. Afterwards, the outer integrals can be evaluated expressing the shape functions u [s] as (a) u [s] = ℓa +s ℓa (b) u [s] = s(ℓb −s)2 ℓ2b (116) (s<0) (s > 0) (117) and defining the outer integration ends in terms of the previously introduced distance as summarized in Table VII. Taking the limit for → 0 for all the contributions, the following final result is attained: (i) C[ˆ , (i) ] = Iaa + Iab + Iba + Ibb ℓb D (−1+)(32ℓa ℓb (1+) arctan = 2 ℓa 16ℓb +4(−ℓ2b (3+)+ℓa2 (2+6)) ln(ℓa )−ℓ2b (−7+3+4(3+) ln(ℓb )) −4(−ℓ2b (3+)+ℓa2 (1+3)) ln(ℓ2b +ℓ2b )) (118) When the boundary conditions allow the transversal displacement of the corner, the evaluation of the boundary coefficients requires the computation of contributions associated with both distributed and point sources. In these situations the double boundary integrals appear together with single boundary integrals and point terms. The case presenting the greatest number of different kernels Table VII. Integration ends s1 , s2 in terms of the quantity . s1 s2 Copyright q Iaa Iab Iba −ℓa + − −ℓa − ℓb ℓb − 2008 John Wiley & Sons, Ltd. Ibb Int. J. Numer. Meth. Engng 2008; 76:127–155 DOI: 10.1002/nme 148 M. MAZZA, L. LEONETTI AND M. ARISTODEMO Figure 8. Shape functions associated with the source û (i) and the field u (i) interpolation parameters. arises when the coefficient corresponding to the source û (i) and the field u (i) parameters is considered (Figure 8): u [x]R∗ w [x, y (i) ] dx t∗ w [x, y]u [y] dy dx + u [x] C[û (i) , u (i) ] = + (i) (i) (i) (i) ∗ ∗ t [x (i) , y] u [y] dy + R [x (i) , y (i) ] t t (119) where the kernels exhibit orders of singularity ranging from O(1/r 2 ) to O(1/r 4 ). This case arises when both the sides merging at the corner are free of supports. Splitting the overlapping source and field supports (i) into the boundary elements a and b, expression (119) can be rewritten in the form: C[û (i) , u (i) ] = Iaa + Iab + Iba + Ibb + Ja + Jb + Ha + Hb + K (a) (a) ∗ = (a) [x] t [x, y] [y] dy dx + [x] u u u w a + + + a b a a u(b) [x] a a t∗ w [x, y] (a) u [y] dy dx + ∗ (i) (a) u [x]Rw [x, y ] dx + ∗ t [x (i) , y] (a) u [y] dy + t b b t∗ w [x, y] u(b) [y] dy dx b b u(b) [x] b t∗ w [x, y] (b) u [y] dy dx ∗ (i) (b) u [x]Rw [x, y ] dx ∗ ∗ t [x (i) , y] u(b) [y] dy + R [x (i) , y (i) ] t t (120) In the above expression, the double integrals represent distributed contributions produced by distributed sources (Iaa , Iab , Iba , Ibb ), the single integrals are distributed contributions associated with point sources (Ja , Jb ) or points contributions associated with distributed sources (Ha , Hb ). The point quantity K is constituted by the twisting moments on the two sides merging at a corner Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2008; 76:127–155 DOI: 10.1002/nme 149 ANALYTICAL INTEGRATION OF SINGULAR KERNELS produced by the point tangent rotation discontinuity applied on the sides merging at the same corner. The evaluation of the double boundary integrals in (120) is carried out introducing the reference system ( , ) over each boundary element where the field distribution is considered and using the abscissa s over the elements where the source distribution is located: Iaa = Iab = Iba = Ibb = s2 s1 s2 s1 s2 s1 s2 s1 (a) u [s]· lim →0 (a) u [s] (b) u [s] 2 2 t∗ w [ , ] (a) u [ ] d ds (121) 1 t∗ w [ , ] (b) u [ ] d ds (122) t∗ w [ , ] (a) u [ ] d ds (123) (124) 1 2 1 (b) u [s]· lim →0 2 t∗ w [ , ] (b) u [ ] d ds 1 In these expressions the shape functions are defined as (a) u [ ]= (b) u [ ]= (a) u [s] = (b) u [s] = ( − 2 (2 1 −3 2 ) 3 ) (125) + 2 −3 1 ) (− 1 + 2 )3 (126) 1) ( ( − 2) + 1− 2 2 (2 (ℓa −2s)(ℓa +s)2 ℓa3 (ℓb +2s)(ℓb −s)2 ℓ3b (s<0) (127) (s>0) (128) and the integration ends 1 , 2 and s1 , s2 are defined in terms of the local abscissa s and of the distance , respectively, using the relationships specified in Table VIII. Table VIII. Expressions of 2, and s1 , s2 . Iaa Iab Iba Ibb 1 2 −ℓa −s −s — 0 ℓb −s −ℓa 0 s −s ℓb −s — s1 s2 −ℓa + − −ℓa − ℓb ℓb − Copyright q 1, 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2008; 76:127–155 DOI: 10.1002/nme 150 M. MAZZA, L. LEONETTI AND M. ARISTODEMO The single boundary integrals Ha , Hb are expressed in the local reference system ( , ) in the form: 2 2 (a) (+) ∗ (−) ∗ Ha = Ha − Ha = (a) (129) u [ ] t (+) [ , ] d − u [ ] t (−) [ , ] d t 1 (+) Hb = Hb (−) − Hb = 2 ∗ (b) u [ ]t (+) t 1 t 1 [ , ] d − 2 ∗ (b) u [ ]t (−) t 1 [ , ] d (130) where the superscripts (−) and (+) distinguish the point sources located on the sides before and after the corner. It is worth noting that the singularities of the kernels are temporarily deactivated defining the integration ends 1 , 2 and the variable in terms of a quantity as detailed in Table IX. This expedient separates the distributed fields on the two sides merging at the corner from the point sources. A similar procedure can be used to deal with the single boundary integrals Ja , Jb , represented in terms of the local abscissa s: s2 s2 ∗(−) ∗(+) Ja = Ja(+) − Ja(−) = (a) (131) (a) [s]m [s] ds − u [s]m t w [s] ds u t w s1 s1 (+) Jb = Jb (−) − Jb = s2 s1 ∗(+) (b) u [s]m t w [s] ds − s2 s1 ∗(−) (b) u [s]m t w [s] ds (132) where the superscripts (−) and (+) denote the twisting moments acting on the sides before and after the corner. After the integrals have been computed, the integration ends s1 , s2 are expressed in terms of the distance , as specified in Table X. The computation of the point contribution K only involves the evaluation of the kernels m ∗t t , expressed by the local abscissa s. In order to deactivate the singularities, the point tangent rotation sources of are moved away from the corner, giving K = (m ∗ (−) ∗ (+) ∗ (−) ∗ (+) (i) −(m (−) [s = −]−m (−) [s = −])(i) (+) [s = ]) (+) [s = ]−m t t t t t t t t Table IX. Expressions of (+) Ha 1 2 −ℓa + 0 1, 2 (133) and in terms of the distance . (+) Hb 0 ℓb (−) (−) Ha Hb −ℓa 0 − ℓb − 0 Table X. Expressions of s1 and s2 in terms of the distance . (−) Ja s1 s2 Copyright q −ℓa + − 2008 John Wiley & Sons, Ltd. (+) Ja −ℓa − (−) Jb (+) Jb ℓb ℓb − Int. J. Numer. Meth. Engng 2008; 76:127–155 DOI: 10.1002/nme ANALYTICAL INTEGRATION OF SINGULAR KERNELS 151 The superscripts (·)(+) and (·)(−) are used to specify both the locations of the point sources t and the locations of the point effects m t with respect to the corner (Figure 8). Grouping all the contributions described above, the singular contributions disappear. The limit on the variable restores the real situation, leading to the following final result: C[û (i) , u (i) ] = lim (Iaa + Iab + Iba + Ibb + Ja + Jb + Ha + Hb + K ) →0 = D 4ℓa3 ℓ3b 3(−1+)(ℓa ℓb (ℓa2 +ℓ2b )(1+11)−4ℓa4 (−1+) + 2(−1+)(4(ℓa4 −ℓ4b ) arctan − (ℓa2 +ℓ2b ) ln(ℓa2 +ℓ2b ))) ℓa +ℓa ℓb (2ℓa2 ln(ℓa )+2ℓ2b ln(ℓb ) ℓb (134) 6.3. Remainder cases In the previous subsections four examples explain how the proposed integration method works for some singular kernels. The cases, which have been selected, refer to the highest singularities, including the presence of corner contributions. The development of a computer code for the symmetric boundary element analysis of Kirchhoff plates, with various boundary conditions, requires the availability of the results for an extended set of cases, which cannot be provided here. A complete collection of the analytical results for the analysis of rectangular plates can be found in the technical report [25], which contains all the integral coefficients, together with the shape functions and the kernels involved. This report can be downloaded at the URL address: www.labmec.unical.it/pubblicazioni/collana.php. 7. TEST OF THE INTEGRATION TECHNIQUE In this section the effectiveness of the proposed integration technique is checked by means of a numerical test. Owing to the difficulty in finding analytical results for kernels exhibiting the high orders of singularity involved in the symmetric boundary element analysis of Kirchhoff plates, the double integration of an O(1/r 2 ) singular kernel for a potential problem has been considered. For this problem Bonnet and Guiggiani provide the numerical results obtained by a direct approach in comparison with the analytical solution [26]. The case of totally overlapping supports, constituted by non-collinear elements, is considered. The boundary coefficient is defined by the integral: (i) (i) ˆ C[ g ∗ [x, y] [y] dy dx (135) [x] , ]= (i) (i) where the kernel g ∗ is g ∗ [x, y, n, ] = Copyright q 2008 John Wiley & Sons, Ltd. 1 (n 1 1 +n 2 2 −2rn r ) 2r 2 (136) Int. J. Numer. Meth. Engng 2008; 76:127–155 DOI: 10.1002/nme 152 M. MAZZA, L. LEONETTI AND M. ARISTODEMO Figure 9. Local reference system ( , ) and abscissa s for the contribution Iba . and the shape functions [x] and [y] are associated with the source interpolation parameter ˆ (i) and to the field interpolation parameter (i) located at the same node i. According to the previous notation, the boundary integral (135) is split into four terms (i) C[ ˆ , (i) ] = Iaa + Iab + Iba + Ibb s2 2 (a) (a) [s] lim = g ∗ [ , ] [ ] d ds + →0 s1 + s2 s1 (b) [s] (a) s1 1 2 s2 (a) g ∗ [ , ] [ ] d ds + 1 s2 s1 (b) [s] [s] lim →0 2 (b) g ∗ [ , ] [ ] d ds 1 2 (b) g ∗ [ ] [ ] d ds 1 (137) where the following linear shape functions are used ( − 2) ( 1 − 2) (a) [ ]= (b) [ ]=− ( − 1) ( 1 − 2) s ℓa s (b) [s] = 1+ ℓb (a) (138) [s] = 1− (139) (s<0) (140) (s>0) (141) The inner integration bounds 1 , 2 and the local variable in terms of the distance s and the angle between the boundary elements (Figure 9) are detailed in Table XI which also defines the outer integration bounds. Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2008; 76:127–155 DOI: 10.1002/nme 153 ANALYTICAL INTEGRATION OF SINGULAR KERNELS Table XI. Inner and outer integration ends. Iaa ( = 0) Iab ( > 0) 2 −s −ℓa −s → 0 −s cos ℓb −s cos −s sin s1 s2 −ℓa + − −ℓa − 1 Ibb ( = 0) Iba ( <0) −ℓa −s cos −s cos s sin −s ℓb −s → 0 ℓb ℓb − Table XII. Integrals for different angles . ℓa /ℓb = 1 = /4 =0 = /2 = 3/4 −1 −1.772588722239 −1.772588722239 −1 −1 −1.651597320576 −1.651597320576 −1 −1 −1.263943507354 −1.263943507354 −1 −1 −0.509441809342 −0.509441809342 −1 (i) C[ ˆ , (i) ] −5.545177444479 −5.303194641153 −4.527887014709 −3.018883618685 Iaa Iab Iba Ibb −1 −3.513809508233 −0.741220785993 −1 −1 −3.422460201748 −0.649871479508 −1 −1 −3.144582412480 −0.371993690240 −1 −1 −2.677817314297 0.094771407942 −1 −6.255030294227 −6.072331681256 −5.516576102721 −4.58304590635 Iaa Iab Iba Ibb (i) C[ ˆ , (i) ] Carrying out the integrations and the limits, the following results are found: Iaa = Ibb = 1 4 (142) ℓa +ℓb cos ℓb +ℓa cos 1 2 2 +ℓb arctan 2 sin ℓa arctan Iab = 8ℓa ℓb ℓa sin ℓb sin + ln(ℓa2 +2 cos ℓa ℓb +ℓ2b )(2ℓa ℓb +cos (ℓa2 +ℓ2b )) − ln ℓ2b (ℓ2b cos +2ℓa ℓb )−2ℓa ℓb −ℓa2 cos ln(ℓa2 )+(2 −) sin (ℓa2 +ℓ2b ) Iba = Iab [ℓa ↔ ℓb ] (143) (144) The numerical values of these integrals, multiplied by 4, are shown in Table XII, which refers to some specific geometric data. These values coincide with the results derived by the analytical formula used in [26] to compare the results obtained by the direct approach. Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2008; 76:127–155 DOI: 10.1002/nme 154 M. MAZZA, L. LEONETTI AND M. ARISTODEMO 8. CONCLUSIONS Symmetric boundary element models may be a preferable alternative to standard collocation boundary models. For this reason, it is interesting to analyze the computational aspects encountered in their actual construction. Indeed, their formulation appears to be heavier and contains kernels having a less compact form and higher singularities. These features have limited the development of symmetric boundary models. One of the crucial computational steps in the construction of these models is represented by the accurate evaluation of the boundary coefficients. In particular, the accuracy in evaluating the self-coefficients, corresponding to overlapping distributions of sources and boundary variables, is a delicate issue. In this work an analytical integration method able to evaluate all the types of self-coefficients, using simple mathematical tools, has been presented. The integration procedure has been described referring to the symmetric boundary element analysis of Kirchhoff plates. This context has been chosen in order to explore the computational problems arising in severe cases, characterized by singularities of up to O(1/r 4 ). As discussed in the previous sections, the proposed method proves to be effective for integrating the coefficients containing the highest singularities in both the cases of collinear and corner domains. In contrast to other techniques, usually based on the extraction of the singular part, the proposed method evaluates the boundary coefficients through a fully analytical procedure. No explicit regularization process, such as transformation by parts, is carried out, whereas the smoothing effect of the interpolation functions and their continuity are exploited. In the proposed method the self-coefficients are evaluated introducing a limit process to make regular the kernels temporarily and developing the integration on a support constituted by two contiguous boundary elements where the continuity of the shape functions is exploited. This choice allows terms that would be singular in the case of separate evaluation to be removed. The remainder singularities are neutralized exploiting the zeroing of the shape functions and their derivatives at the ends of the support. Owing to the variety of possible boundary conditions, the complete set of coefficients is rather extensive. The examples, explaining the details of the proposed integration procedure, refer to higher singularities and include corner situations. It is worth noting that the integration method has been described with reference to rectangular plates but could be used for generic polygonal plates. In this case the presence of the parametric description of the orientation of the boundary elements forming the integration domains increases the number of kernels involved in the continuity condition and reduces the compactness of the analytical results. REFERENCES 1. Hartmann F. Introductions to Boundary Elements—Theory and Applications. Springer: Berlin, 1989. 2. Sirtori S, Maier G, Novati G, Miccoli S. A Galerkin symmetric boundary-element method in elasticity: formulation and implementation. International Journal for Numerical Methods in Engineering 1992; 35:255–282. 3. Frangi A, Bonnet M. A Galerkin symmetric and direct BIE method for Kirchhoff elastic plates: formulation and implementation. International Journal for Numerical Methods in Engineering 1998; 41:337–369. 4. Hartmann F, Katz C, Protopsaltis B. Boundary elements and symmetry. Ingenieur Archiv 1985; 55:440–449. 5. Holzer S. How to deal with hypersingular integrals in the symmetric BEM. Communications in Numerical Methods in Engineering 1993; 9:219–232. 6. Hadamard J. Lectures on Cauchy’s Problem in Linear Partial Differential Equations. Yale University Press: New Haven, CT, 1923. Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2008; 76:127–155 DOI: 10.1002/nme ANALYTICAL INTEGRATION OF SINGULAR KERNELS 155 7. Carini A, Diligenti M, Maranesi P, Zanelli M. Analytical integration for two-dimensional elastic analysis by the symmetric Galerkin boundary element methods. Computational Mechanics 1999; 23:308–323. 8. Jaswon MA, Maiti M. An integral equation formulation for plate bending problems. Journal of Engineering Mathematics 1968; 2:83–93. 9. Forbes DJ, Robinson AR. Numerical analysis of elastic plates and shallow shells by an integral equation method. University of Illinois Structural Research Series Report 345/, 1969. 10. Altiero NJ, Sikarskie DL. A boundary integral method applied to plates of arbitrary plane forms. Computers and Structures 1978; 9:163–168. 11. Bézine G, Gamby D. Boundary integral formulation for plate flexure with arbitrary boundary conditions. Mechanics Research Communications 1978; 5:197–206. 12. Hartmann F, Zotemantel R. The direct boundary element method in plate bending. International Journal for Numerical Methods in Engineering 1986; 23:75–93. 13. Stern M. A general boundary integral formulation for the numerical solution of plate bending problems. International Journal of Solids and Structures 1979; 15:769–782. 14. Knöpke B. The hypersingular integral equation for the bending moments m x x , m x y and m yy of the Kirchhoff plate. Computational Mechanics 1994; 15:19–30. 15. Frangi A, Guiggiani M. A direct approach for boundary integral equations with high-order singularities. International Journal for Numerical Methods in Engineering 2000; 49:871–898. 16. Frangi A, Guiggiani M. Boundary element analysis of Kirchhoff plates with direct evaluation of hypersingular integrals. International Journal for Numerical Methods in Engineering 1999; 46:1845–1863. 17. Tottenham H. The boundary element method for plates and shells. In Developments in Boundary Element Methods, vol. I, Banerjee PK, Butterfield R (eds), 1979. 18. Giroire J, Nedelec JC. A new system of boundary integral equations for plates with free edges. Mathematical Methods in the Applied Sciences 1995; 8:755–772. 19. Nazaret C. Equations intègrals de frontiere pour des problèmes polygonales à bord libre. Comptes Rendus de l’ Academie des Sciences Paris Sèrie I 1996; 322:989–994. 20. Mazza M, Leonetti L, Aristodemo M. Symmetric boundary element analysis of Kirchhoff plates. Proceedings ECCOMAS 2004, Finland, Jyväskylä, 2004. 21. Mazza M, Aristodemo M. Integration of hypersingular kernels in symmetric BEM analysis of Kirchhoff plates. Proceedings Sixth International Conference on Computational Structures Technology, Czech Republic, Prague, 2002. 22. Timoshenko SP, Woinowsky-Krieger S. Theory of Plates and Shells. McGraw-Hill: New York, 1984 23. Frangi A. A new regularized BE formulation for Kirchhoff plates. European Journal of Mechanics – A/Solids 1996; 15:915–931. 24. Bonnet M. A regularized Galerkin symmetric BIE formulation for mixed 3D elastic boundary values problems. Boundary Elements Abstract and Newsletters 1993; 4(109):113. 25. Mazza M, Leonetti L, Aristodemo M. A table of results from analytical integration of singular kernels arising in the symmetric boundary element analysis of Kirchhoff plates. LabMec Technical Report n.49, 2007. (Available from: www.labmec.unical.it/pubblicazioni/collana.php.) 26. Bonnet M, Giuggiani M. Direct evaluation of double singular integrals and new free terms in 2D (symmetric) Galerkin BEM. Computer Methods in Applied Mechanics and Engineering 2003; 192:2565–2596. Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2008; 76:127–155 DOI: 10.1002/nme

RELATED PAPERS

RELATED TOPICS

Log In

Analytical integration of singular kernels in symmetric boundary element analysis of Kirchhoff plates

Analytical integration of singular kernels in symmetric boundary element analysis of Kirchhoff plates

Related Papers

RELATED PAPERS

RELATED TOPICS