- Email: [email protected]
An improved direct method for evaluating hypersingular stress boundary integral equations in BEM
Engineering Analysis with Boundary Elements 61 (2015) 274–281
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
An improved direct method for evaluating hypersingular stress boundary integral equations in BEM Wei-Zhe Feng, Jian Liu, Xiao-Wei Gao n State Key Laboratory of Structural Analysis for Industrial Equipment, School of Aeronautics and Astronautics, Dalian University of Technology, Dalian 116024, China
art ic l e i nf o
a b s t r a c t
Article history: Received 1 April 2015 Received in revised form 3 August 2015 Accepted 6 August 2015
An efficient numerical method for evaluating all kinds of singular boundary integrals presented in Ref. Gao, 2002 [16] is improved by using a newly derived formulation for computing the spatial derivative of the global distance, which is inevitably used in the finite part of radial integrals. Based on this improvement, more accurate and stable results are obtained in evaluating singular integrals on curved boundary elements. Applying this improved singular integral method to evaluate hypersingular boundary integrals in the stress boundary integral equation, and a high order accuracy of stresses is achieved. Numerical examples are given to illustrate the correctness of the proposed method, and computational results show that the accurate stresses can be obtained even on a coarse mesh, compared with other conventional techniques. & 2015 Elsevier Ltd. All rights reserved.
Keywords: Boundary element method Hypersingular integrals Stress analysis Power series expansion
1. Introduction In the application of boundary element method to solve elasticity problems, basic unknown variables (boundary displacements and tractions) and internal stresses can be evaluated by directly solving the boundary integral equation. However, when the source point is located on the boundary of the problem, strongly singular and hypersingular integrals arise in the stress boundary integral equation. Therefore, the stresses over boundary nodes cannot be solved directly from the stress integral equation [1,2]. So far there are two major approaches for evaluating boundary stresses. One is the commonly used approach, known as “traction recovery method (TRM)” [1–4]. By differentiating shape functions of boundary elements and employing Hooke's law, boundary stresses can be determined by the tractions over the source point and displacements over the elements adjacent to the source point. Since derivatives of the shape functions are used in this method, continuity of stresses across the adjacent elements cannot be ensured, as a result, stresses over adjacent elements have to be averaged, meanwhile, the stress precision decreases. Another approach for computing boundary stresses is the method evaluating hypersingular stress integral equations directly. Since stresses are computed from the integral equation, they have the accuracy as high as that of displacements, provided that
n
Corresponding author. Tel.: þ 86 8470 6332. E-mail address: [email protected] (X.-W. Gao).
http://dx.doi.org/10.1016/j.enganabound.2015.08.002 0955-7997/& 2015 Elsevier Ltd. All rights reserved.
hypersingular integrals appearing in the stress integral equation can be evaluated accurately. Among the methods directly evaluating hypersingular stress integral equations, one is the regularizing method. Rudolphi [5] proposed two “simple solutions” for fundamental functions, and using these “simple solutions”, the normal derivative of the boundary integral equation can be regularized. Chati [6] applied this method to the calculation of hypersingular potential gradients and elastic stresses. Liu [7] established some integral identities for fundamental functions, and with these identities the conventional boundary integral equation and hypersingular boundary integral equation, can both be written in weakly singular forms. Besides the regularizing methods, Guiggiani [8,9] proposed a general algorithm for evaluating hypersingular curved surface integrals based on expansion of geometry quantities as Laurent series in the parameter plane of intrinsic coordinates. Huber [10] applied this algorithm to the evaluation of 3D boundary elastic stresses. Frangi and Guiggiani [11] extended the direct method by Guiggiani et al. for second-order problems to forth-order problems, like bending of thin elastic plates. However, the shape of the exclusion neighborhood is not assumed to be circular, and all singular integrals and free-term coefficients are evaluated in the general framework. And then the systematic direct approach was extended to problems with higher singularities up to r 4 [12]. Based on the similar direct approach, Frangi and Guiggiani [13] proposed a new general method for the evaluation of free-term coefficients in 3D hypersingular boundary integral equation for general cases of the source point being located at a non-smooth boundary point, and proofs of compatibility conditions were given.
W.-Z. Feng et al. / Engineering Analysis with Boundary Elements 61 (2015) 274–281
Chen and Hong [14] reviewed the application of dual boundary element methods, and summarized the regularizations of hypersingular integrals and divergent series in detail. Recently Chen and his coworkers [15] proposed dual integral equations for Ndimensional Laplace problems, and free terms of the dual boundary integral equations for the four-dimensional Laplace problems were derived by using the bump-contour approach. Recently Gao [16] proposed a method for handling supersingular boundary integrals with the order of singularities up to 6 by expanding the non-singular part of the integrand as a power series in the intrinsic plane. In this way, the singularities are condensed in the radial integral and removed analytically by expressing the non-singular part of the resulting integrand as a power series. The numerical results coincide well with results in Ref. [8], and the method can solving all kinds of singular integrals. But when expanding non-singular part of the integrand (a function of r ;i , and r ;i ¼ ðxi xpi Þ=r) as a power series in the intrinsic plane, an approximate way is taken for evaluating the derivative variable r ;i for the first collocation point. This approximation brings in errors and unstabilities in the radial integrals, especially when the source point is close to the counter of the element (e.g. when employing discontinuous elements). In this paper, exact analytical formulation for r ;i is derived when the collocation point is located on the source point. Based on this small improvement, more accurate and stable results are obtained in evaluating singular integrals on curved boundary elements. Applying this improved direct method to evaluate hypersingular integrals in the boundary stress integral equation, a high order accuracy of stresses is achieved. The first numerical example is given to demonstrate the effectiveness of the improved method in evaluating high order singular boundary integrals. The other two examples are given to illustrate the correctness of the proposed method in evaluating boundary stresses, and results show that accurate stress values can be obtained even on a coarse mesh, compared with conventional technique and FEM method.
2. Elasticity boundary integral equations The classic elasticity boundary integral equation has the form as follows [1,2]: Z Z cij ðxp Þuj ðxp Þ ¼ U ij ðx; xp Þt j ðxÞdΓðxÞ T ij ðx; xp Þuj ðxÞdΓðxÞ; ð1Þ Γ
evaluation. For the strongly singular integral including the kernel function Tij, rigid body motion method [1] can be used to determine the singular diagonal term by cumulating the regular off-diagonal items. Taking a derivative of Eq. (1) with respect to the coordinates of source point xp , and noticing that ∂ðnÞ=∂xpi ¼ ∂ðnÞ=∂xi , the displacement gradient integral equation can be derived as follows: Z Z cui;k ðxp Þ ¼ U ij;k ðx; xp Þt j ðxÞdΓ T ij;k ðx; xp Þuj ðxÞdΓ; ð4Þ Γ
Γ
where ∂U ij ðx; xp Þ ; ∂xk
ð5Þ
∂T ij ðx; xp Þ : ∂xk
ð6Þ
U ij;k ðx; xp Þ ¼ T ij;k ðx; xp Þ ¼
Generalized Hooke's law reads: σ ij ¼ Dijkl uk;l ;
ð7Þ
where
2ν δij δkl þ δik δjl þ δil δjk : Dijkl ¼ μ 1 2ν
ð8Þ
Substitute Eq. (4) into Eq. (7), yields the stress integral equation: Z Z cσ ij ðxp Þ ¼ U ijk ðx; xp Þt k ðxÞdΓðxÞ T ijk ðx; xp Þuk ðxÞdΓðxÞ; ð9Þ Γ
Γ
where U ijk ¼
1 1 Cðδki r ;j þ δkj r ;i δij r ;k Þ þ βr ;i r ;j r ;k ; 4παð1 νÞ r α
μ 1 βr;m nm ½Cδij r k þ νðδik r;j þ δjk r i Þ γ r;i r;j r;k 2παð1 νÞ r β þ βνðni r;j r;k þ nj r;i r;k Þ þ Cðβnk r;i r;j þ nj δik þni δjk Þ Dnk δij :
ð10Þ
T ijk ¼
ð11Þ
In Eqs. (10) and (11), C ¼ 1 2ν; D ¼ 1 4ν; γ ¼ β þ 2, α ¼ β 1, with β being the dimension of problems (β ¼ 2 for 2D, and β ¼ 3 for 3D problems). When the source point xp is located on the boundary, the two boundary integrals in Eq. (9) are strongly singular and hypersingular respectively.
Γ
where cij ¼ cδij . When source point xp is located inside the computational domain c¼ 1; when source point is located on the boundary, c depends on the geometry of the boundary. In Eq. (1), U ij ðx; xp Þ is the Kelvin fundamental solution for elasticity problems, x and xp stand for field point and source point respectively: 8 < 8πð11 νÞμ½ð3 4νÞδij ln 1r þ r ;i r ;j for 2D U ij ¼ ; ð2Þ : 16πð11 νÞμ r ð3 4νÞδij þ r ;i r ;j for 3D where μ is shear modulus, ν is Poisson's ratio, and δij represents the Kronecker delta. From U ij ðx; xp Þ, the expression of T ij ðx; xp Þ can be derived as follows: T ij ¼
275
1 ð1 2νÞðni r ;j nj r ;i Þ þ nm r m ½ð1 2νÞδij þ βr ;i r ;j : 4απð1 νÞr α ð3Þ
In Eq. (1) when source point is approaching to the boundary, the first and second integral on the right hand side are weakly and strongly singular respectively. For the weakly singular integral, logarithmic quadrature formulation [1] for 2D and element subdivision technique [1,17] for 3D problems can be used for precise
3. Discretization of boundary integrals In order to numerically evaluate the boundary integrals shown in Eq. (9), coordinates and physical variables in an element can be expressed in terms of their nodal values as follows: xi ðξÞ ¼
N node X
N α ðξÞ xαi ;
ð12Þ
N α ðξÞ uαi ;
ð13Þ
Nα ðξÞ t αi :
ð14Þ
α¼1
ui ðξÞ ¼
N node X α¼1
t i ðξÞ ¼
N node X α¼1
where N node is the number of element nodes; ξ represents the intrinsic coordinate (for a surface element it being understood that ðξÞ ¼ ðξ1 ; ξ2 Þ); Nα is the shape function of the α-th node [2], and xαi is the i-th component of coordinates at the α-th node. Eq. (12) serves as a mapping from the global coordinate system to the intrinsic coordinate system. After mapping curved element from the global coordinate system, it becomes a straight line or a square element in the coordinate system, see Figs. 1 and 2.
276
W.-Z. Feng et al. / Engineering Analysis with Boundary Elements 61 (2015) 274–281
4. Evaluation of high order singular boundary integrals
Fig. 1. Nodes over line element in intrinsic coordinate system: (a) linear and (b) quadratic.
For convenience, the two integrals included in (15) can be expressed as Z T ijk ðx; xp ÞN α ðξÞj J ε j I 1 ðxp Þ ¼ dΓðξÞ; ð17aÞ r β ðx; xp Þ Γξ Z I 2 ðxp Þ ¼
U ijk ðx; xp ÞN α ðξÞj J ε j dΓðξÞ: r α ðx; xp Þ Γξ
ð17bÞ
where U ijk ðx; xp Þ and T ijk ðx; xp Þ are regular part of the singular kernel functions U ijk ðx; xp Þ and T ijk ðx; xp Þ. Eq. (17a) is hypersingular, and Eq. (17b) is strongly singular. Considering the integration over a 3D curved element, Eq. (17a) becomes: Z 1 Z 1 T ijk ðx; xp ÞNα ðξÞ I 1 ðxp Þ ¼ J e dξdη: ð18Þ r β ðx; xp Þ 1 1 Fig. 2. Nodes over quadrilateral elements in intrinsic coordinate system: (a) 4noded linear element and (b) 8-noded quadratic element.
n
q
The integration region in Eq. (18) is a square plane, and therefore the radial integration method (RIM) [18] can be employed to transform it into a contour integral along four sides of the element (see Fig. 3) as follows: Z 1 ∂ρðx; xp Þ Fðx; xp ÞdLðxÞ; ð19Þ I 1 ðxp Þ ¼ p ∂n' L ρðx; x Þ where
q
Z
q
Fðx; xp Þ ¼ lim ε-0
p
p
p
Fig. 3. Integration plane in intrinsic coordinate system, five collocation points are distributed on the radial direction from source point p to field point q.
For e-th element, the element integration can be expressed as follows: Ie1 ðxp Þ ¼
N node X
uαi
Z
α¼1
Γe
T ijk ðx; xp ÞN α ðξÞdΓðxÞ ¼
N node X
uαi
Z
α¼1
Γξ
T ijk ðx; xp ÞN α ðξÞj J ε j dΓðξÞ;
ð15aÞ Ie2 ðxp Þ ¼
N node X α¼1
t αi
Z Γe
U ijk ðx; xp ÞN α ðξÞdΓðxÞ ¼
N node X α¼1
t αi
Z Γξ
U ijk ðx; xp ÞN α ðξÞj J ε j dΓðξÞ:
ρðx;xp Þ ρξ ðεÞ
T ijk ðx; xp ÞN α ðξÞ J e ρ dρ; r β ðx; xp Þ
ð20Þ
In Eq. (19), x represents the field point which takes value on the boundary of the square as shown in Fig. 3; n0 is the outward normal to the boundary of the integration square in the intrinsic coordinate system (ξ; η); L is the contour of the intrinsic square of the element; and ρξ ðεÞ is an infinitesimal radius of a local circle isolating the source point (ξp ; ηp ); and ρ is the local distance defined as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ρ ¼ ðξ ξP Þ2 þ ðη ηP Þ2 ; ð21Þ where ðξp ; ηp Þ is intrinsic coordinate of the source point. From Fig. 3, it can be seen that ∂ρðx; xp Þ ¼ ρξ n0ξ þ ρη n0η ; ∂n0
ð22Þ
where x ξxp ρ;ξ ¼ ξρðx;x p Þ ¼ cos θ
ρη ¼
ηx ηxp ¼ sin θ: ρðx; xp Þ
ð23Þ
ð15bÞ In Eq. (15) Γ e is the global integration region of element e; Γ ξ is a straight line from 1 to þ 1 for 2D and a square shown in Fig. 2 for 3D problems, j J ε j is the transformation Jacobian from the global coordinates to intrinsic coordinates: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi j J ε j ¼ ð∂x=∂ξÞ2 þ ð∂y=∂ξÞ2 for 2D : ð16Þ ∂x for 3D j J ε j ¼ ∂x ∂ξ ∂η where x is the coordinate vector, and ð∂x=∂ξÞ can be determined by differentiating Eq. (12). For regular integrals, i.e., when xp is not located on the element under integration, Eq. (15) can be evaluated using standard Gaussian quadrature [1,2]. For singular integrals, i.e., xp is located on the element under integration, Gauss quadrature will generate unacceptable results and particular integral method is required.
4.1. Expanding global distance r as series of local distance ρ In order to evaluate the singular radial integral, i.e. Eq. (20), the global distance r need to be expressed as a power series in the local distance ρ. To derive the series, the intrinsic coordinates ξ and η need to be expressed in terms of ρ first. Referring to Fig. 3, the following relationships are written: ( ξ ¼ ξp þ ρρ;ξ ¼ ξp þ ρ cos θ : ð24Þ η ¼ ηp þ ρρ;η ¼ ηp þ ρ sin θ The expansion for the following power series is proposed vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u M N uX X r n ¼ρ¼ Cnρ ¼ t Gm ρm ; ð25Þ ρ n¼0 m¼0
W.-Z. Feng et al. / Engineering Analysis with Boundary Elements 61 (2015) 274–281
where N and M are the orders of power series, C n and Gm are coefficients to be determined. To determine the coefficients C n and Gm we rewrite the intrinsic coordinate as follows: ξ ¼ ξp þξs ρ:
ð26Þ
The global distance r can be determined as follows vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi )2 u β (N u M uX X node uX t p ¼ ρt N l ðξðρÞÞxli xi Gm ρm ;
ð27Þ
i¼1
m¼0
l¼1
where β is the dimension of the problem and N node is the number of element nodes. Substituting shape function N l into Eq. (27), and expanding the left-hand side of Eq.(27) into a power series in ρ, the coefficients Gm and order M can be determined by comparing equal order terms of ρ on the left- and right-hand sides of Eq.(27). The expressions for coefficients Gm can be seen in [16]. After obtaining the coefficients Gm , we rewrite Eq.(25) in the following form: !2 N M X X n Cnρ ¼ Gm ρm : ð28Þ n¼0
n¼0
Expanding the left-hand side of the above equation and comparing the coefficients for the same powers of ρ on the two sides, the following recursive relationship can be obtained: Gn ¼
n X
CiCn i:
ð29Þ
i¼0
Solving this equation set for C n , we can obtain the following recursive relationship: pffiffiffiffiffiffi C 0 ¼ G0 ; ð30aÞ !
Cn ¼
n 1 X 1 Gn CiCn i : 2C 0 i¼1
ð30bÞ
4.2. Evaluation of singular radial integrals Substituting Eq. (25) into the radial integrals Eq. (20) yields Z ρðx;xp Þ FðρÞ dρ; ð31Þ Fðx; xp Þ ¼ lim β1 ε-0 ρ ðεÞ ρ ξ where FðρÞ is the regular part, i.e. T ijk ðx; xp Þ J e Nα ðξÞ FðρÞ ¼ : β ρ ðρÞ
ð32Þ
In order to evaluate the singular integral shown in Eq. (31), the non-singular part of the integrand is expressed as a power series in the local distanceρ, i.e. FðρÞ ¼
N X
Bn ρn ;
ð33Þ
n¼0
where N is the order of the series, and Bn the coefficients which are determined by collocating N þ 1 points over the integration region ρðξρ ; ξs Þ, see Fig. 3. Eqs. (24) and (12) are used to determine the intrinsic and global coordinates needed for calculating the quantities appearing in the right-hand side of Eq. (32). The coefficient for the first point (n ¼0) is computed using T ijk ðxp ; xp Þj J ε j Nα ðξÞ ðC 0 Þβ
¼ B0 :
ð34Þ
Other coefficients can be solved using ½RðNÞfBg ¼ fYg;
ð35Þ
277
where [R(N)] is a square matrix of order N, i.e. 2 3 1 1; ρ1 ; ⋯ ρN 1 6 1 7 6 1; ρ2 ; ⋯ ρN 7 2 7: ½R ¼ 6 6⋮ ⋮ ⋯ ⋮ 7 4 5 1 1; ρN ; ⋯ ρN N fBg and fYg are vectors as follows: 9 8h i 9 8 > > > Fðρ1 Þ Bð0Þ =ρ1 > 1 > > > > > > B > > > > h i > > > > > > > = < Fðρ Þ Bð0Þ =ρ > < B2 = 2 2 : ; fYg ¼ fBg ¼ > > > > > > > > > >⋮ >⋮ > h i > > > > > > : N; > > ð0Þ B > ; : FðρN Þ B =ρN >
ð36Þ
ð37Þ
Once the coefficients Bn are obtained from above equations, substituting Eq. (33) into (31) yields Z ρðξp ;ξs Þ N N X X Bn lim ρn β dρ ¼ Bn E n ; ð38Þ Fðx; xp Þ ¼ n¼0
ε-0
ρξ ðεÞ
n¼0
where 8 > > < n 1β þ 1 ρβ n 11 ðξ ;ξ Þ limρ β n1 1 ðεÞ p s ε-0 ξ En ¼ > > ln ρðξ ; ξ Þ lim ln ρ p s ξ ðεÞ : ε-0
ðn aβ 1Þ ðn ¼ β 1Þ
ð39Þ
To separate the finite parts from the limiting terms in Eq. (39), we need to express the local infinitesimal distance ρξ ðεÞ in terms of the global one, say ε, and infinite terms [8], i.e. lim ln ρξ ðεÞ ¼ ln H 0 þ Inf inite terms;
ð40aÞ
ε-0
lim ε-0
1 ρβξ n 1 ðεÞ
¼ H β n 1 þ Inf inite terms
where 8 H 0 ¼ 1=C 0 > > > > > H1 ¼ C 1 > > > < 2 H 2 ¼ 2C 2 C 1 ; > 3 > > H ¼ 3C 3C C þ C > 3 3 1 2 1 > > > 2 2 4 > : H 4 ¼ 4C 4 þ 4C 1 C 2 4C 1 C 3 2C 2 C 1 C i ¼ C i =C 0 :
ð0 rn rβ 2Þ:
ð40bÞ
ð41Þ
ð42Þ
For a physical problem, the boundary integrals should exist. This means that the infinite terms in Eq. (40) should be eliminated after considering contributions of all adjacent elements around the source point. The final results can be obtained as
8 1 1 > 0 r n rβ 2 > β n 1 Hβ n 1 ; n β þ 1 > ρ < : ð43Þ En ¼ ln ρ ln H 0 ; n ¼ β1 > > > ρn β þ 1 ; : n4β1 nβþ1
Now, Eqs. (19) and (38) are regularized expressions for evaluating 3D singular boundary integrals, and can be applied to evaluate the singular integrals included in Eq. (17). When expanding non-singular part of the integrand (a function of r ;i , and r ;i ¼ ðxi xpi Þ=r) as a power series in the radial finite part integral on the intrinsic plane, see Eqs. (33)–(37), the first collocation point is taken as the source point, see Eq. (34) and Fig. 3, the distance variable r exactly equals to 0, therefore program bug arises in evaluating r ;i by using the function r ;i ¼ ðxi xpi Þ=r. In Ref. [16], an approximate way is taken to bypass this program bug, the first collocation point (p1 in Fig. 4) shrinks in a little distance ρο in the radial direction on the intrinsic plane, keeps a little distance away from the source point to evaluate the variable r ;i ¼ ðxqi xpi Þ=r
278
W.-Z. Feng et al. / Engineering Analysis with Boundary Elements 61 (2015) 274–281
Substituting Eqs. (47) and (49) into Eq. (44) yields r ;i ¼
Fig. 4. Collocation points distributed on the intrinsic radial direction, and pðξp ; ηp Þ is the source point, qðξq ; ηq Þ is the field point, p1 p5 are the collocation points for power series expansion, where p1 is the first collocation point with a little distance shrinking in from pðξp ; ηp Þ.
successfully. Because the value of ρο is very small, this approximation is acceptable for general case. But when the source point is very close to the counter of the element (e.g. when employing discontinuous elements), errors and unstabilities brought in by this approximation are obvious. In the next section, the exact analytical formulation for computing the spatial derivative of the global distance r ;i is derived, so that the collocation point can be located on the source point, and the expanding process (Eqs. (33)–(37)) is accurately satisfied.
5. Formulation for computing the spatial derivative of global distance r ;i In boundary element method, most kernel functions include the spatial derivative r ;i , which can be expressed as p
r ;i ¼
∂r xi xi ; ¼ ∂xi r
ð44Þ
where xpi and xi are the coordinates of source point and field point, and r is the distance between source point xpi and field xi , with the qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi p 2 expression r ¼ ðxi xi Þ . When xi is approaching to xpi , the expression of Eq. (44) is regular and can be directly computed. But when xi is coincide with xpi , r ¼ 0, Eq. (44) will lead to indeterminacy. In this case, a particular treatment needs to be performed in following explanation. Appling Taylor expansion, it follows that: xi xpi ¼
∂xi 1 ∂2 xi ðξK ξpK Þ þ ðξ ξpK ÞðξJ ξpJ Þ þ ⋯: 2 ∂ξK ∂ξJ K ∂ξK
ð45Þ
where K and J take the value of 1 for 2D problems, and 1 to 2 for 3D problems, and the repeated subscripts imply summation. From Eq. (24), we can write ξK ¼ ξpK þ ρρ;K ;
ð46Þ
in which, ρ; K is determined by Eqs. (23) for 3D problems, and ρ; 1 is 1 or þ 1 for 2D problems. Substituting Eq. (46) into Eq. (45) yields xi xpi ¼ Ai ρ þ Bi ρ2 þ Oðρ3 Þ;
Ai þBi ρ : MðρÞ
ð51Þ
When the collocation point is located on the source point,ρ ¼ 0, so that the expression of r ;i is as follows (repeated subscripts imply summation): Ai ffi: r ;i ¼ pffiffiffiffiffiffiffiffi Aj Aj
ð52Þ
Eq. (52) is the exact formulation for computing r ;i for the first collocation point, and this is also suitable for 2D case, in 2D i problems Ai ðφÞ ¼ ∂x ∂ξ ξs (ξs takes the value of 1 or þ1). In real application, when evaluating the radial integral using power series expansion and analytical approach, r ;i involved in the integrand for the first collocation point is evaluated using Eq. (52), for other collocation points, using Eq.(44).
6. Numerical examples 6.1. Singular integral over a 11-degree curved surface This example is to demonstrate the effectiveness of the modification in evaluating integrands in radial integrals. The superR singular integral Iðxp Þ ¼ Γ ½3r ;3 ∂r=∂n n3 =ð4πr β ÞdΓ over a cylindrical surface with radius 10, formed by rotating a line with length of 2 for 111 as shown in Fig. 5. The source point b, with intrinsic coordinates of (ξ ¼ 0:7η ¼ 0:7) are computed. Table 1 gives the results for singular order β ¼ 5, with varying power series order N in Eq. (33). For comparison, results using the method in Ref. [19] are also listed, using a same power series order N. From Table 1 we can see that, results of the improved method is very stable compared with the original method in Ref. [16], and coincide well with the results in Ref. [19], it is obvious that the improvement of the original method is effective. 6.2. Singular integral over a 90-degree cylindrical surface Comparing with the well-established intrinsic Laurent series expansion approach, this example is adopted to examine the correctness and stabilities of the proposed method. The following singular integral is operated on a 90-degree cylinder surface element (radius¼ 1, length ¼2) as shown in Fig. 6
Z 1 ∂r Iðxp Þ ¼ 3r ;3 n3 dΓ: ð53Þ β 4πr ∂n Γ
ð47Þ
where Ai ¼
∂xi ρ ; ∂ξK ; K
ð48aÞ
Bi ¼
1 ∂2 xi ρ ρ : 2 ∂ξK ∂ξJ ; K ; J
ð48bÞ
The global distance r can be evaluated using the following expression (repeated subscripts imply summation): qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r ¼ ðxi xpi Þðxi xpi Þ ¼ MðρÞρ; ð49Þ where MðρÞ ¼
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Ai Ai þ 2Ai Bi ρ þBi Bi ρ2 :
ð50Þ
Fig. 5. 11-degree cylindrical surface element.
W.-Z. Feng et al. / Engineering Analysis with Boundary Elements 61 (2015) 274–281
This example was analyzed by Guiggiani et al. [8], in which analytical Laurent series expansion in intrinsic space is adopted to eliminate singularities for the case β ¼ 3. Their results are used here as comparison. Table 2 gives the results for source points a, b and c with the intrinsic coordinates (0, 0), (0.66, 0) and (0.66, 0.66), respectively. From Table 2, it can be seen that the current results are in good agreement with the results given by Guiggiani el al. in Ref. [8], demonstrating the correctness of the formulations in this paper. To examine the influence of N, the order of power series in Eq. (38) determined by the number of collocation points distributed in the local distance ρ, on the computed results, Table 3 lists the computed results with varying power series order N. And relative errors are calculated taking the numerical results of Ref. [8] as the standard values. It can be seen that the results is stable, and when
Table 1 Results of supersingular integral on curved boundary element for source point b (β ¼ 5) N
Improved method
Original method
Ref. [19]
4 5 6 7 8
2.4055649 2.4055663 2.4055629 2.4055633 2.4055634
2.4068766 2.4113997 2.4210097 2.4377839 2.4638292
2.4055632 2.4055646 2.4055612 2.4055616 2.4055617
z 2
c b a
1
y
x Fig. 6. A 90-degree curved element.
the order of power series N is equal or more than 6, accurate results can be achieved. The number of Gauss points Ng used to evaluate the line integral shown in Eq. (19) is automatically determined in code based on a criterion described in Ref. [2]. To investigate the dependence of the computational results on the number of Gauss points, computations using different Gauss points were performed and the results are listed in Table 4, which shows that remarkably stable results can be achieved. 6.3. 3D high-order singular integral over a rectangular region This example is concerned with the evaluation of the integral R4 R1 r Iðxp Þ ¼ 0 1 r;1β dxdy over a rectangular region with x ¼ f 1; 1g and y ¼ f0; 4g, see Fig. 7. Source point a with the intrinsic coordinatesð0:99; 0:5Þ is computed. Table 5 lists the computed results for singularity orders of β from 3 to 6. Analytical results are easy to obtain for this example and are also listed in Table 5. From Table 5 we can see that the numerical results remarkably coincide with the analytical results. 6.4. 3D thick cylinder under internal pressure This example is about a thick cylinder under internal press. The internal diameter of the model is 200, and external diameter 400. The value of internal pressure is P ¼ 10, material constants: shear modulus μ ¼ 400, Poisson's ratio ν ¼ 0:25. Due to symmetry, a quarter of the model is selected for computing. The height of the cylinder is Z¼400, see Fig. 8. Two end sections are constrained in Z direction. Fig. 8 shows the BEM model for computation, in which the internal and external surfaces are discretized into 10 continuous equally-spaced quadratic elements in hoop direction, each straight line over sections is discretized into 10 continuous equally-spaced elements, and along thickness 10 elements are also used. Totally, there are 1082 nodes and 600 boundary elements. Because the mesh is fine, “Traction Recovery Method” (TRM) has corresponding high accuracy, so we uses the results of “Traction Recovery Method” as comparison. Table 6 shows the results of shear stress τxy over internal hoop line of the middle cross-section. Fig. 9 shows the distribution of normal stress σ xx over the internal line of the end cross-section. Table 4 Computational results for different numbers of Gauss points (β¼ 3).
Table 2 Computed results for three points a, b and c (β¼ 3).
Ref. [8] Current
279
Point a
Point b
Point c
0.343807 0.343808
0.497099 0.497095
0.877214 0.877370
Ng
Point a
Point b
Point c
4 8 10 14 20
0.343802427838 0.343802427838 0.343802427832 0.343802427838 0.343802427838
0.49709467538 0.49709467538 0.49709467532 0.49709467538 0.49709467538
0.87737047821 0.87737047821 0.87737047823 0.87737047821 0.87737047821
Table 3 Computed results with varying order of power series (β ¼ 3). Order N
4 5 6 7 8
Point a
Point b
Point c
Value
Error (%)
Value
Error (%)
Value
Error (%)
0.340241 0.343924 0.343986 0.343817 0.343802
1.037E0 3.403E 2 5.206E 2 2.909E 3 1.454E 3
0.493786 0.497201 0.497267 0.497109 0.497095
6.665E 1 2.052E 2 3.380E 2 2.012E 3 8.050E 4
0.881048 0.881320 0.876125 0.876885 0.877370
4.370E 1 4.680E 1 1.241E 1 3.751E 2 1.778E 2
280
W.-Z. Feng et al. / Engineering Analysis with Boundary Elements 61 (2015) 274–281
4
y
xp
Fig. 9. Normal stress σ xx over internal hoop line.
x
2 Fig. 7. 8-noded rectangular element.
Fig. 10. Rectangular plate with a hole. Table 5 Computed results for varying order of singularity. β
Analytical
Current
3 4 5 6
156.048 6666.37 392699 2.66667E7
156.048 6666.37 392699.0 2.66667E7
Fig. 11. Quarter of the model to be analyzed.
Fig. 8. 3D BEM model of the thick-wall cylinder. Fig. 12. BEM model.
Table 6 Shear stresses τxy over internal hoop line. Angle (1)
TRM
SIEPPEM
Error (%)
90 81 72 63 54 45 36 27 18 9 0
1.6742 4.1202 7.8371 10.7869 12.6807 13.3333 12.6807 10.7869 7.8371 4.1202 1.6742
1.6669 4.1205 7.8376 10.7875 12.6815 13.3341 12.6815 10.7875 7.8376 4.1205 1.6669
0.4391 0.0054 0.0058 0.0058 0.0059 0.0058 0.0059 0.0058 0.0058 0.0054 0.4391
From Table 6 and Fig. 9 we can see that the approach proposed in this paper coincide well with the classic method TRM, that demonstrates the accuracy of this method. 6.5. A plate with a circular hole under an extension load This example is about a tension P applied to a rectangle plate including a hole. The thickness of the plate is 0.4, and other geometrical dimensions have been shown in Fig. 10. The applied tension load is P ¼ 10. The shear modulus is μ ¼ 400, and Poisson's ratio is ν ¼ 0:25. Due to symmetry, a quarter of the model is selected for computing, see Fig. 11. Fig. 12 shows the BEM model, in which surface of the hole is dispersed into 6 quadratic elements, and 2 in
W.-Z. Feng et al. / Engineering Analysis with Boundary Elements 61 (2015) 274–281
281
results even using a very coarse mesh, demonstrating the accuracy of the presented method.
7. Conclusions
Fig. 13. FEM model.
A new technique for computing the spatial derivative of the global distance r is derived in the paper, which is inevitably used in the power series expansion. And by applying the improved singular boundary integral evaluation method to the stress integral equation, boundary stresses can be directly evaluated from the hypersingular stress boundary integral equations. Numerical examples show that accurate results can be obtained using the proposed method.
Acknowledgments The authors gratefully acknowledge the National Natural Science Foundation of China (11172055, 11202045). References
Fig. 14. Normal stress σ yy over symmetric arc line (z ¼ 0.2).
Fig. 15. Normal stress σ xx over arc line (z¼ 0.4).
the thickness direction. There are total of 338 nodes and 112 boundary elements. To test the accuracy, the FEM results computed by commercial software ANSYS is also given based on a very fine mesh shown in Fig. 13. And another hypersingular integral method proposed in Ref. [19] is also given for comparison. Fig. 14 shows the normal stress σ yy in y direction over the middle arc of the curved surface and Fig. 15 shows the normal stress σ xx over the end arc of the curved surface. Form Figs. 14 and 15, it can be seen that the current results coincide well with FEM
[1] Brebbia CA, Dominguez J. Boundary Elements: An introductory course. London: Mc Graw-Hill Book Co.; 1992. [2] Gao XW, Davies TG. Boundary Element Programming in Mechanics. Cambridge,UK: Cambridge University Press; 2002. [3] Banerjee PK, Raveendra ST. Advanced boundary blement analysis of twodimensional and 3-dimensional problems of elastoplasticity. Int J Numer Methods Eng 1986;23(6):985–1002. [4] Kane JH. Boundary Element Analysis in Engineering Continuum Mechanics. Englewood Cliffs, NJ: Prentice Hall; 1994. [5] Rudolphi TJ. The use of simple solutions in the regularization of hypersingular boundary integral equations. Math Comput Model 1991;15(3):269–78. [6] Chati MK, Mukherjee S. Evaluation of gradients on the boundary using fully regularized hypersingular boundary integral equations. Acta Mech 1999;135 (1–2):41–55. [7] Liu YJ, Rudolphi TJ. Some identities for fundamental solutions and their applications to weakly-singular boundary element formulations. Eng Anal Bound Elem 1991;8(6):301–11. [8] Guiggiani M, Krishnasamy G, Rudolphi T, Rizzo F. A general algorithm for the numerical solution of hypersingular boundary integral equations. ASME J Appl Mech 1992;59:604. [9] Guiggiani M. Hypersingular formulation for boundary stress evaluation. Eng Anal Bound Elem 1994;13(2):169–79. [10] Huber O, Lang A, Kuhn G. Evaluation of the stress tensor in 3D elastostatics by direct solving of hypersingular integrals. Comput Mech 1993;12(1-2):39–50. [11] Frangi A, Guiggiani M. Boundary element analysis of Kirchhoff plates with direct evaluation of hypersingular integrals. Int J Numer Methods Eng 1999;46 (11):1845–63. [12] Frangi A, Guiggiani M. A direct approach for boundary integral equations with high-order singularities. Int J Numer Methods Eng 2000;49(7):871–98. [13] Frangi A, Guiggiani M. Free terms and compatibility conditions for 3D hypersingular boundary integral equations. Z Angew Math Mech 2001;81 (10):651–64. [14] Chen JT, Hong HK. Review of dual boundary element methods with emphasis on hypersingular integrals and divergent series. ASME Appl Mech Rev 1999;52 (1):17–33. [15] Chen JT, Huang WS, Lee JW, Hong HK. On the free terms of the dual BIE for Ndimensional Laplace problems. Eng Anal Bound Elem 2015;59:123–8. [16] Gao XW. An effective method for numerical evaluation of general 2D and 3D high order singular boundary integrals. Comput Methods Appl Mech Eng 2010;199(45):2856–64. [17] Lachat JC. A Further Development of the Boundary Integral Technique for Elastostatics. Ph.D. thesis. Southampton: University of Southampton; 1975. [18] Gao XW. The radial integration method for evaluation of domain integrals with boundary-only discretization. Eng Anal Bound Elem 2002;26(10):905–16. [19] Gao XW, Feng WZ, Yang K, Cui M. Projection plane method for evaluation of arbitrary high order singular boundary integrals. Eng Anal Bound Elem 2015;50:265–74.
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
An improved direct method for evaluating hypersingular stress boundary integral equations in BEM Wei-Zhe Feng, Jian Liu, Xiao-Wei Gao n State Key Laboratory of Structural Analysis for Industrial Equipment, School of Aeronautics and Astronautics, Dalian University of Technology, Dalian 116024, China
art ic l e i nf o
a b s t r a c t
Article history: Received 1 April 2015 Received in revised form 3 August 2015 Accepted 6 August 2015
An efficient numerical method for evaluating all kinds of singular boundary integrals presented in Ref. Gao, 2002 [16] is improved by using a newly derived formulation for computing the spatial derivative of the global distance, which is inevitably used in the finite part of radial integrals. Based on this improvement, more accurate and stable results are obtained in evaluating singular integrals on curved boundary elements. Applying this improved singular integral method to evaluate hypersingular boundary integrals in the stress boundary integral equation, and a high order accuracy of stresses is achieved. Numerical examples are given to illustrate the correctness of the proposed method, and computational results show that the accurate stresses can be obtained even on a coarse mesh, compared with other conventional techniques. & 2015 Elsevier Ltd. All rights reserved.
Keywords: Boundary element method Hypersingular integrals Stress analysis Power series expansion
1. Introduction In the application of boundary element method to solve elasticity problems, basic unknown variables (boundary displacements and tractions) and internal stresses can be evaluated by directly solving the boundary integral equation. However, when the source point is located on the boundary of the problem, strongly singular and hypersingular integrals arise in the stress boundary integral equation. Therefore, the stresses over boundary nodes cannot be solved directly from the stress integral equation [1,2]. So far there are two major approaches for evaluating boundary stresses. One is the commonly used approach, known as “traction recovery method (TRM)” [1–4]. By differentiating shape functions of boundary elements and employing Hooke's law, boundary stresses can be determined by the tractions over the source point and displacements over the elements adjacent to the source point. Since derivatives of the shape functions are used in this method, continuity of stresses across the adjacent elements cannot be ensured, as a result, stresses over adjacent elements have to be averaged, meanwhile, the stress precision decreases. Another approach for computing boundary stresses is the method evaluating hypersingular stress integral equations directly. Since stresses are computed from the integral equation, they have the accuracy as high as that of displacements, provided that
n
Corresponding author. Tel.: þ 86 8470 6332. E-mail address: [email protected] (X.-W. Gao).
http://dx.doi.org/10.1016/j.enganabound.2015.08.002 0955-7997/& 2015 Elsevier Ltd. All rights reserved.
hypersingular integrals appearing in the stress integral equation can be evaluated accurately. Among the methods directly evaluating hypersingular stress integral equations, one is the regularizing method. Rudolphi [5] proposed two “simple solutions” for fundamental functions, and using these “simple solutions”, the normal derivative of the boundary integral equation can be regularized. Chati [6] applied this method to the calculation of hypersingular potential gradients and elastic stresses. Liu [7] established some integral identities for fundamental functions, and with these identities the conventional boundary integral equation and hypersingular boundary integral equation, can both be written in weakly singular forms. Besides the regularizing methods, Guiggiani [8,9] proposed a general algorithm for evaluating hypersingular curved surface integrals based on expansion of geometry quantities as Laurent series in the parameter plane of intrinsic coordinates. Huber [10] applied this algorithm to the evaluation of 3D boundary elastic stresses. Frangi and Guiggiani [11] extended the direct method by Guiggiani et al. for second-order problems to forth-order problems, like bending of thin elastic plates. However, the shape of the exclusion neighborhood is not assumed to be circular, and all singular integrals and free-term coefficients are evaluated in the general framework. And then the systematic direct approach was extended to problems with higher singularities up to r 4 [12]. Based on the similar direct approach, Frangi and Guiggiani [13] proposed a new general method for the evaluation of free-term coefficients in 3D hypersingular boundary integral equation for general cases of the source point being located at a non-smooth boundary point, and proofs of compatibility conditions were given.
W.-Z. Feng et al. / Engineering Analysis with Boundary Elements 61 (2015) 274–281
Chen and Hong [14] reviewed the application of dual boundary element methods, and summarized the regularizations of hypersingular integrals and divergent series in detail. Recently Chen and his coworkers [15] proposed dual integral equations for Ndimensional Laplace problems, and free terms of the dual boundary integral equations for the four-dimensional Laplace problems were derived by using the bump-contour approach. Recently Gao [16] proposed a method for handling supersingular boundary integrals with the order of singularities up to 6 by expanding the non-singular part of the integrand as a power series in the intrinsic plane. In this way, the singularities are condensed in the radial integral and removed analytically by expressing the non-singular part of the resulting integrand as a power series. The numerical results coincide well with results in Ref. [8], and the method can solving all kinds of singular integrals. But when expanding non-singular part of the integrand (a function of r ;i , and r ;i ¼ ðxi xpi Þ=r) as a power series in the intrinsic plane, an approximate way is taken for evaluating the derivative variable r ;i for the first collocation point. This approximation brings in errors and unstabilities in the radial integrals, especially when the source point is close to the counter of the element (e.g. when employing discontinuous elements). In this paper, exact analytical formulation for r ;i is derived when the collocation point is located on the source point. Based on this small improvement, more accurate and stable results are obtained in evaluating singular integrals on curved boundary elements. Applying this improved direct method to evaluate hypersingular integrals in the boundary stress integral equation, a high order accuracy of stresses is achieved. The first numerical example is given to demonstrate the effectiveness of the improved method in evaluating high order singular boundary integrals. The other two examples are given to illustrate the correctness of the proposed method in evaluating boundary stresses, and results show that accurate stress values can be obtained even on a coarse mesh, compared with conventional technique and FEM method.
2. Elasticity boundary integral equations The classic elasticity boundary integral equation has the form as follows [1,2]: Z Z cij ðxp Þuj ðxp Þ ¼ U ij ðx; xp Þt j ðxÞdΓðxÞ T ij ðx; xp Þuj ðxÞdΓðxÞ; ð1Þ Γ
evaluation. For the strongly singular integral including the kernel function Tij, rigid body motion method [1] can be used to determine the singular diagonal term by cumulating the regular off-diagonal items. Taking a derivative of Eq. (1) with respect to the coordinates of source point xp , and noticing that ∂ðnÞ=∂xpi ¼ ∂ðnÞ=∂xi , the displacement gradient integral equation can be derived as follows: Z Z cui;k ðxp Þ ¼ U ij;k ðx; xp Þt j ðxÞdΓ T ij;k ðx; xp Þuj ðxÞdΓ; ð4Þ Γ
Γ
where ∂U ij ðx; xp Þ ; ∂xk
ð5Þ
∂T ij ðx; xp Þ : ∂xk
ð6Þ
U ij;k ðx; xp Þ ¼ T ij;k ðx; xp Þ ¼
Generalized Hooke's law reads: σ ij ¼ Dijkl uk;l ;
ð7Þ
where
2ν δij δkl þ δik δjl þ δil δjk : Dijkl ¼ μ 1 2ν
ð8Þ
Substitute Eq. (4) into Eq. (7), yields the stress integral equation: Z Z cσ ij ðxp Þ ¼ U ijk ðx; xp Þt k ðxÞdΓðxÞ T ijk ðx; xp Þuk ðxÞdΓðxÞ; ð9Þ Γ
Γ
where U ijk ¼
1 1 Cðδki r ;j þ δkj r ;i δij r ;k Þ þ βr ;i r ;j r ;k ; 4παð1 νÞ r α
μ 1 βr;m nm ½Cδij r k þ νðδik r;j þ δjk r i Þ γ r;i r;j r;k 2παð1 νÞ r β þ βνðni r;j r;k þ nj r;i r;k Þ þ Cðβnk r;i r;j þ nj δik þni δjk Þ Dnk δij :
ð10Þ
T ijk ¼
ð11Þ
In Eqs. (10) and (11), C ¼ 1 2ν; D ¼ 1 4ν; γ ¼ β þ 2, α ¼ β 1, with β being the dimension of problems (β ¼ 2 for 2D, and β ¼ 3 for 3D problems). When the source point xp is located on the boundary, the two boundary integrals in Eq. (9) are strongly singular and hypersingular respectively.
Γ
where cij ¼ cδij . When source point xp is located inside the computational domain c¼ 1; when source point is located on the boundary, c depends on the geometry of the boundary. In Eq. (1), U ij ðx; xp Þ is the Kelvin fundamental solution for elasticity problems, x and xp stand for field point and source point respectively: 8 < 8πð11 νÞμ½ð3 4νÞδij ln 1r þ r ;i r ;j for 2D U ij ¼ ; ð2Þ : 16πð11 νÞμ r ð3 4νÞδij þ r ;i r ;j for 3D where μ is shear modulus, ν is Poisson's ratio, and δij represents the Kronecker delta. From U ij ðx; xp Þ, the expression of T ij ðx; xp Þ can be derived as follows: T ij ¼
275
1 ð1 2νÞðni r ;j nj r ;i Þ þ nm r m ½ð1 2νÞδij þ βr ;i r ;j : 4απð1 νÞr α ð3Þ
In Eq. (1) when source point is approaching to the boundary, the first and second integral on the right hand side are weakly and strongly singular respectively. For the weakly singular integral, logarithmic quadrature formulation [1] for 2D and element subdivision technique [1,17] for 3D problems can be used for precise
3. Discretization of boundary integrals In order to numerically evaluate the boundary integrals shown in Eq. (9), coordinates and physical variables in an element can be expressed in terms of their nodal values as follows: xi ðξÞ ¼
N node X
N α ðξÞ xαi ;
ð12Þ
N α ðξÞ uαi ;
ð13Þ
Nα ðξÞ t αi :
ð14Þ
α¼1
ui ðξÞ ¼
N node X α¼1
t i ðξÞ ¼
N node X α¼1
where N node is the number of element nodes; ξ represents the intrinsic coordinate (for a surface element it being understood that ðξÞ ¼ ðξ1 ; ξ2 Þ); Nα is the shape function of the α-th node [2], and xαi is the i-th component of coordinates at the α-th node. Eq. (12) serves as a mapping from the global coordinate system to the intrinsic coordinate system. After mapping curved element from the global coordinate system, it becomes a straight line or a square element in the coordinate system, see Figs. 1 and 2.
276
W.-Z. Feng et al. / Engineering Analysis with Boundary Elements 61 (2015) 274–281
4. Evaluation of high order singular boundary integrals
Fig. 1. Nodes over line element in intrinsic coordinate system: (a) linear and (b) quadratic.
For convenience, the two integrals included in (15) can be expressed as Z T ijk ðx; xp ÞN α ðξÞj J ε j I 1 ðxp Þ ¼ dΓðξÞ; ð17aÞ r β ðx; xp Þ Γξ Z I 2 ðxp Þ ¼
U ijk ðx; xp ÞN α ðξÞj J ε j dΓðξÞ: r α ðx; xp Þ Γξ
ð17bÞ
where U ijk ðx; xp Þ and T ijk ðx; xp Þ are regular part of the singular kernel functions U ijk ðx; xp Þ and T ijk ðx; xp Þ. Eq. (17a) is hypersingular, and Eq. (17b) is strongly singular. Considering the integration over a 3D curved element, Eq. (17a) becomes: Z 1 Z 1 T ijk ðx; xp ÞNα ðξÞ I 1 ðxp Þ ¼ J e dξdη: ð18Þ r β ðx; xp Þ 1 1 Fig. 2. Nodes over quadrilateral elements in intrinsic coordinate system: (a) 4noded linear element and (b) 8-noded quadratic element.
n
q
The integration region in Eq. (18) is a square plane, and therefore the radial integration method (RIM) [18] can be employed to transform it into a contour integral along four sides of the element (see Fig. 3) as follows: Z 1 ∂ρðx; xp Þ Fðx; xp ÞdLðxÞ; ð19Þ I 1 ðxp Þ ¼ p ∂n' L ρðx; x Þ where
q
Z
q
Fðx; xp Þ ¼ lim ε-0
p
p
p
Fig. 3. Integration plane in intrinsic coordinate system, five collocation points are distributed on the radial direction from source point p to field point q.
For e-th element, the element integration can be expressed as follows: Ie1 ðxp Þ ¼
N node X
uαi
Z
α¼1
Γe
T ijk ðx; xp ÞN α ðξÞdΓðxÞ ¼
N node X
uαi
Z
α¼1
Γξ
T ijk ðx; xp ÞN α ðξÞj J ε j dΓðξÞ;
ð15aÞ Ie2 ðxp Þ ¼
N node X α¼1
t αi
Z Γe
U ijk ðx; xp ÞN α ðξÞdΓðxÞ ¼
N node X α¼1
t αi
Z Γξ
U ijk ðx; xp ÞN α ðξÞj J ε j dΓðξÞ:
ρðx;xp Þ ρξ ðεÞ
T ijk ðx; xp ÞN α ðξÞ J e ρ dρ; r β ðx; xp Þ
ð20Þ
In Eq. (19), x represents the field point which takes value on the boundary of the square as shown in Fig. 3; n0 is the outward normal to the boundary of the integration square in the intrinsic coordinate system (ξ; η); L is the contour of the intrinsic square of the element; and ρξ ðεÞ is an infinitesimal radius of a local circle isolating the source point (ξp ; ηp ); and ρ is the local distance defined as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ρ ¼ ðξ ξP Þ2 þ ðη ηP Þ2 ; ð21Þ where ðξp ; ηp Þ is intrinsic coordinate of the source point. From Fig. 3, it can be seen that ∂ρðx; xp Þ ¼ ρξ n0ξ þ ρη n0η ; ∂n0
ð22Þ
where x ξxp ρ;ξ ¼ ξρðx;x p Þ ¼ cos θ
ρη ¼
ηx ηxp ¼ sin θ: ρðx; xp Þ
ð23Þ
ð15bÞ In Eq. (15) Γ e is the global integration region of element e; Γ ξ is a straight line from 1 to þ 1 for 2D and a square shown in Fig. 2 for 3D problems, j J ε j is the transformation Jacobian from the global coordinates to intrinsic coordinates: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi j J ε j ¼ ð∂x=∂ξÞ2 þ ð∂y=∂ξÞ2 for 2D : ð16Þ ∂x for 3D j J ε j ¼ ∂x ∂ξ ∂η where x is the coordinate vector, and ð∂x=∂ξÞ can be determined by differentiating Eq. (12). For regular integrals, i.e., when xp is not located on the element under integration, Eq. (15) can be evaluated using standard Gaussian quadrature [1,2]. For singular integrals, i.e., xp is located on the element under integration, Gauss quadrature will generate unacceptable results and particular integral method is required.
4.1. Expanding global distance r as series of local distance ρ In order to evaluate the singular radial integral, i.e. Eq. (20), the global distance r need to be expressed as a power series in the local distance ρ. To derive the series, the intrinsic coordinates ξ and η need to be expressed in terms of ρ first. Referring to Fig. 3, the following relationships are written: ( ξ ¼ ξp þ ρρ;ξ ¼ ξp þ ρ cos θ : ð24Þ η ¼ ηp þ ρρ;η ¼ ηp þ ρ sin θ The expansion for the following power series is proposed vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u M N uX X r n ¼ρ¼ Cnρ ¼ t Gm ρm ; ð25Þ ρ n¼0 m¼0
W.-Z. Feng et al. / Engineering Analysis with Boundary Elements 61 (2015) 274–281
where N and M are the orders of power series, C n and Gm are coefficients to be determined. To determine the coefficients C n and Gm we rewrite the intrinsic coordinate as follows: ξ ¼ ξp þξs ρ:
ð26Þ
The global distance r can be determined as follows vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi )2 u β (N u M uX X node uX t p ¼ ρt N l ðξðρÞÞxli xi Gm ρm ;
ð27Þ
i¼1
m¼0
l¼1
where β is the dimension of the problem and N node is the number of element nodes. Substituting shape function N l into Eq. (27), and expanding the left-hand side of Eq.(27) into a power series in ρ, the coefficients Gm and order M can be determined by comparing equal order terms of ρ on the left- and right-hand sides of Eq.(27). The expressions for coefficients Gm can be seen in [16]. After obtaining the coefficients Gm , we rewrite Eq.(25) in the following form: !2 N M X X n Cnρ ¼ Gm ρm : ð28Þ n¼0
n¼0
Expanding the left-hand side of the above equation and comparing the coefficients for the same powers of ρ on the two sides, the following recursive relationship can be obtained: Gn ¼
n X
CiCn i:
ð29Þ
i¼0
Solving this equation set for C n , we can obtain the following recursive relationship: pffiffiffiffiffiffi C 0 ¼ G0 ; ð30aÞ !
Cn ¼
n 1 X 1 Gn CiCn i : 2C 0 i¼1
ð30bÞ
4.2. Evaluation of singular radial integrals Substituting Eq. (25) into the radial integrals Eq. (20) yields Z ρðx;xp Þ FðρÞ dρ; ð31Þ Fðx; xp Þ ¼ lim β1 ε-0 ρ ðεÞ ρ ξ where FðρÞ is the regular part, i.e. T ijk ðx; xp Þ J e Nα ðξÞ FðρÞ ¼ : β ρ ðρÞ
ð32Þ
In order to evaluate the singular integral shown in Eq. (31), the non-singular part of the integrand is expressed as a power series in the local distanceρ, i.e. FðρÞ ¼
N X
Bn ρn ;
ð33Þ
n¼0
where N is the order of the series, and Bn the coefficients which are determined by collocating N þ 1 points over the integration region ρðξρ ; ξs Þ, see Fig. 3. Eqs. (24) and (12) are used to determine the intrinsic and global coordinates needed for calculating the quantities appearing in the right-hand side of Eq. (32). The coefficient for the first point (n ¼0) is computed using T ijk ðxp ; xp Þj J ε j Nα ðξÞ ðC 0 Þβ
¼ B0 :
ð34Þ
Other coefficients can be solved using ½RðNÞfBg ¼ fYg;
ð35Þ
277
where [R(N)] is a square matrix of order N, i.e. 2 3 1 1; ρ1 ; ⋯ ρN 1 6 1 7 6 1; ρ2 ; ⋯ ρN 7 2 7: ½R ¼ 6 6⋮ ⋮ ⋯ ⋮ 7 4 5 1 1; ρN ; ⋯ ρN N fBg and fYg are vectors as follows: 9 8h i 9 8 > > > Fðρ1 Þ Bð0Þ =ρ1 > 1 > > > > > > B > > > > h i > > > > > > > = < Fðρ Þ Bð0Þ =ρ > < B2 = 2 2 : ; fYg ¼ fBg ¼ > > > > > > > > > >⋮ >⋮ > h i > > > > > > : N; > > ð0Þ B > ; : FðρN Þ B =ρN >
ð36Þ
ð37Þ
Once the coefficients Bn are obtained from above equations, substituting Eq. (33) into (31) yields Z ρðξp ;ξs Þ N N X X Bn lim ρn β dρ ¼ Bn E n ; ð38Þ Fðx; xp Þ ¼ n¼0
ε-0
ρξ ðεÞ
n¼0
where 8 > > < n 1β þ 1 ρβ n 11 ðξ ;ξ Þ limρ β n1 1 ðεÞ p s ε-0 ξ En ¼ > > ln ρðξ ; ξ Þ lim ln ρ p s ξ ðεÞ : ε-0
ðn aβ 1Þ ðn ¼ β 1Þ
ð39Þ
To separate the finite parts from the limiting terms in Eq. (39), we need to express the local infinitesimal distance ρξ ðεÞ in terms of the global one, say ε, and infinite terms [8], i.e. lim ln ρξ ðεÞ ¼ ln H 0 þ Inf inite terms;
ð40aÞ
ε-0
lim ε-0
1 ρβξ n 1 ðεÞ
¼ H β n 1 þ Inf inite terms
where 8 H 0 ¼ 1=C 0 > > > > > H1 ¼ C 1 > > > < 2 H 2 ¼ 2C 2 C 1 ; > 3 > > H ¼ 3C 3C C þ C > 3 3 1 2 1 > > > 2 2 4 > : H 4 ¼ 4C 4 þ 4C 1 C 2 4C 1 C 3 2C 2 C 1 C i ¼ C i =C 0 :
ð0 rn rβ 2Þ:
ð40bÞ
ð41Þ
ð42Þ
For a physical problem, the boundary integrals should exist. This means that the infinite terms in Eq. (40) should be eliminated after considering contributions of all adjacent elements around the source point. The final results can be obtained as
8 1 1 > 0 r n rβ 2 > β n 1 Hβ n 1 ; n β þ 1 > ρ < : ð43Þ En ¼ ln ρ ln H 0 ; n ¼ β1 > > > ρn β þ 1 ; : n4β1 nβþ1
Now, Eqs. (19) and (38) are regularized expressions for evaluating 3D singular boundary integrals, and can be applied to evaluate the singular integrals included in Eq. (17). When expanding non-singular part of the integrand (a function of r ;i , and r ;i ¼ ðxi xpi Þ=r) as a power series in the radial finite part integral on the intrinsic plane, see Eqs. (33)–(37), the first collocation point is taken as the source point, see Eq. (34) and Fig. 3, the distance variable r exactly equals to 0, therefore program bug arises in evaluating r ;i by using the function r ;i ¼ ðxi xpi Þ=r. In Ref. [16], an approximate way is taken to bypass this program bug, the first collocation point (p1 in Fig. 4) shrinks in a little distance ρο in the radial direction on the intrinsic plane, keeps a little distance away from the source point to evaluate the variable r ;i ¼ ðxqi xpi Þ=r
278
W.-Z. Feng et al. / Engineering Analysis with Boundary Elements 61 (2015) 274–281
Substituting Eqs. (47) and (49) into Eq. (44) yields r ;i ¼
Fig. 4. Collocation points distributed on the intrinsic radial direction, and pðξp ; ηp Þ is the source point, qðξq ; ηq Þ is the field point, p1 p5 are the collocation points for power series expansion, where p1 is the first collocation point with a little distance shrinking in from pðξp ; ηp Þ.
successfully. Because the value of ρο is very small, this approximation is acceptable for general case. But when the source point is very close to the counter of the element (e.g. when employing discontinuous elements), errors and unstabilities brought in by this approximation are obvious. In the next section, the exact analytical formulation for computing the spatial derivative of the global distance r ;i is derived, so that the collocation point can be located on the source point, and the expanding process (Eqs. (33)–(37)) is accurately satisfied.
5. Formulation for computing the spatial derivative of global distance r ;i In boundary element method, most kernel functions include the spatial derivative r ;i , which can be expressed as p
r ;i ¼
∂r xi xi ; ¼ ∂xi r
ð44Þ
where xpi and xi are the coordinates of source point and field point, and r is the distance between source point xpi and field xi , with the qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi p 2 expression r ¼ ðxi xi Þ . When xi is approaching to xpi , the expression of Eq. (44) is regular and can be directly computed. But when xi is coincide with xpi , r ¼ 0, Eq. (44) will lead to indeterminacy. In this case, a particular treatment needs to be performed in following explanation. Appling Taylor expansion, it follows that: xi xpi ¼
∂xi 1 ∂2 xi ðξK ξpK Þ þ ðξ ξpK ÞðξJ ξpJ Þ þ ⋯: 2 ∂ξK ∂ξJ K ∂ξK
ð45Þ
where K and J take the value of 1 for 2D problems, and 1 to 2 for 3D problems, and the repeated subscripts imply summation. From Eq. (24), we can write ξK ¼ ξpK þ ρρ;K ;
ð46Þ
in which, ρ; K is determined by Eqs. (23) for 3D problems, and ρ; 1 is 1 or þ 1 for 2D problems. Substituting Eq. (46) into Eq. (45) yields xi xpi ¼ Ai ρ þ Bi ρ2 þ Oðρ3 Þ;
Ai þBi ρ : MðρÞ
ð51Þ
When the collocation point is located on the source point,ρ ¼ 0, so that the expression of r ;i is as follows (repeated subscripts imply summation): Ai ffi: r ;i ¼ pffiffiffiffiffiffiffiffi Aj Aj
ð52Þ
Eq. (52) is the exact formulation for computing r ;i for the first collocation point, and this is also suitable for 2D case, in 2D i problems Ai ðφÞ ¼ ∂x ∂ξ ξs (ξs takes the value of 1 or þ1). In real application, when evaluating the radial integral using power series expansion and analytical approach, r ;i involved in the integrand for the first collocation point is evaluated using Eq. (52), for other collocation points, using Eq.(44).
6. Numerical examples 6.1. Singular integral over a 11-degree curved surface This example is to demonstrate the effectiveness of the modification in evaluating integrands in radial integrals. The superR singular integral Iðxp Þ ¼ Γ ½3r ;3 ∂r=∂n n3 =ð4πr β ÞdΓ over a cylindrical surface with radius 10, formed by rotating a line with length of 2 for 111 as shown in Fig. 5. The source point b, with intrinsic coordinates of (ξ ¼ 0:7η ¼ 0:7) are computed. Table 1 gives the results for singular order β ¼ 5, with varying power series order N in Eq. (33). For comparison, results using the method in Ref. [19] are also listed, using a same power series order N. From Table 1 we can see that, results of the improved method is very stable compared with the original method in Ref. [16], and coincide well with the results in Ref. [19], it is obvious that the improvement of the original method is effective. 6.2. Singular integral over a 90-degree cylindrical surface Comparing with the well-established intrinsic Laurent series expansion approach, this example is adopted to examine the correctness and stabilities of the proposed method. The following singular integral is operated on a 90-degree cylinder surface element (radius¼ 1, length ¼2) as shown in Fig. 6
Z 1 ∂r Iðxp Þ ¼ 3r ;3 n3 dΓ: ð53Þ β 4πr ∂n Γ
ð47Þ
where Ai ¼
∂xi ρ ; ∂ξK ; K
ð48aÞ
Bi ¼
1 ∂2 xi ρ ρ : 2 ∂ξK ∂ξJ ; K ; J
ð48bÞ
The global distance r can be evaluated using the following expression (repeated subscripts imply summation): qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r ¼ ðxi xpi Þðxi xpi Þ ¼ MðρÞρ; ð49Þ where MðρÞ ¼
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Ai Ai þ 2Ai Bi ρ þBi Bi ρ2 :
ð50Þ
Fig. 5. 11-degree cylindrical surface element.
W.-Z. Feng et al. / Engineering Analysis with Boundary Elements 61 (2015) 274–281
This example was analyzed by Guiggiani et al. [8], in which analytical Laurent series expansion in intrinsic space is adopted to eliminate singularities for the case β ¼ 3. Their results are used here as comparison. Table 2 gives the results for source points a, b and c with the intrinsic coordinates (0, 0), (0.66, 0) and (0.66, 0.66), respectively. From Table 2, it can be seen that the current results are in good agreement with the results given by Guiggiani el al. in Ref. [8], demonstrating the correctness of the formulations in this paper. To examine the influence of N, the order of power series in Eq. (38) determined by the number of collocation points distributed in the local distance ρ, on the computed results, Table 3 lists the computed results with varying power series order N. And relative errors are calculated taking the numerical results of Ref. [8] as the standard values. It can be seen that the results is stable, and when
Table 1 Results of supersingular integral on curved boundary element for source point b (β ¼ 5) N
Improved method
Original method
Ref. [19]
4 5 6 7 8
2.4055649 2.4055663 2.4055629 2.4055633 2.4055634
2.4068766 2.4113997 2.4210097 2.4377839 2.4638292
2.4055632 2.4055646 2.4055612 2.4055616 2.4055617
z 2
c b a
1
y
x Fig. 6. A 90-degree curved element.
the order of power series N is equal or more than 6, accurate results can be achieved. The number of Gauss points Ng used to evaluate the line integral shown in Eq. (19) is automatically determined in code based on a criterion described in Ref. [2]. To investigate the dependence of the computational results on the number of Gauss points, computations using different Gauss points were performed and the results are listed in Table 4, which shows that remarkably stable results can be achieved. 6.3. 3D high-order singular integral over a rectangular region This example is concerned with the evaluation of the integral R4 R1 r Iðxp Þ ¼ 0 1 r;1β dxdy over a rectangular region with x ¼ f 1; 1g and y ¼ f0; 4g, see Fig. 7. Source point a with the intrinsic coordinatesð0:99; 0:5Þ is computed. Table 5 lists the computed results for singularity orders of β from 3 to 6. Analytical results are easy to obtain for this example and are also listed in Table 5. From Table 5 we can see that the numerical results remarkably coincide with the analytical results. 6.4. 3D thick cylinder under internal pressure This example is about a thick cylinder under internal press. The internal diameter of the model is 200, and external diameter 400. The value of internal pressure is P ¼ 10, material constants: shear modulus μ ¼ 400, Poisson's ratio ν ¼ 0:25. Due to symmetry, a quarter of the model is selected for computing. The height of the cylinder is Z¼400, see Fig. 8. Two end sections are constrained in Z direction. Fig. 8 shows the BEM model for computation, in which the internal and external surfaces are discretized into 10 continuous equally-spaced quadratic elements in hoop direction, each straight line over sections is discretized into 10 continuous equally-spaced elements, and along thickness 10 elements are also used. Totally, there are 1082 nodes and 600 boundary elements. Because the mesh is fine, “Traction Recovery Method” (TRM) has corresponding high accuracy, so we uses the results of “Traction Recovery Method” as comparison. Table 6 shows the results of shear stress τxy over internal hoop line of the middle cross-section. Fig. 9 shows the distribution of normal stress σ xx over the internal line of the end cross-section. Table 4 Computational results for different numbers of Gauss points (β¼ 3).
Table 2 Computed results for three points a, b and c (β¼ 3).
Ref. [8] Current
279
Point a
Point b
Point c
0.343807 0.343808
0.497099 0.497095
0.877214 0.877370
Ng
Point a
Point b
Point c
4 8 10 14 20
0.343802427838 0.343802427838 0.343802427832 0.343802427838 0.343802427838
0.49709467538 0.49709467538 0.49709467532 0.49709467538 0.49709467538
0.87737047821 0.87737047821 0.87737047823 0.87737047821 0.87737047821
Table 3 Computed results with varying order of power series (β ¼ 3). Order N
4 5 6 7 8
Point a
Point b
Point c
Value
Error (%)
Value
Error (%)
Value
Error (%)
0.340241 0.343924 0.343986 0.343817 0.343802
1.037E0 3.403E 2 5.206E 2 2.909E 3 1.454E 3
0.493786 0.497201 0.497267 0.497109 0.497095
6.665E 1 2.052E 2 3.380E 2 2.012E 3 8.050E 4
0.881048 0.881320 0.876125 0.876885 0.877370
4.370E 1 4.680E 1 1.241E 1 3.751E 2 1.778E 2
280
W.-Z. Feng et al. / Engineering Analysis with Boundary Elements 61 (2015) 274–281
4
y
xp
Fig. 9. Normal stress σ xx over internal hoop line.
x
2 Fig. 7. 8-noded rectangular element.
Fig. 10. Rectangular plate with a hole. Table 5 Computed results for varying order of singularity. β
Analytical
Current
3 4 5 6
156.048 6666.37 392699 2.66667E7
156.048 6666.37 392699.0 2.66667E7
Fig. 11. Quarter of the model to be analyzed.
Fig. 8. 3D BEM model of the thick-wall cylinder. Fig. 12. BEM model.
Table 6 Shear stresses τxy over internal hoop line. Angle (1)
TRM
SIEPPEM
Error (%)
90 81 72 63 54 45 36 27 18 9 0
1.6742 4.1202 7.8371 10.7869 12.6807 13.3333 12.6807 10.7869 7.8371 4.1202 1.6742
1.6669 4.1205 7.8376 10.7875 12.6815 13.3341 12.6815 10.7875 7.8376 4.1205 1.6669
0.4391 0.0054 0.0058 0.0058 0.0059 0.0058 0.0059 0.0058 0.0058 0.0054 0.4391
From Table 6 and Fig. 9 we can see that the approach proposed in this paper coincide well with the classic method TRM, that demonstrates the accuracy of this method. 6.5. A plate with a circular hole under an extension load This example is about a tension P applied to a rectangle plate including a hole. The thickness of the plate is 0.4, and other geometrical dimensions have been shown in Fig. 10. The applied tension load is P ¼ 10. The shear modulus is μ ¼ 400, and Poisson's ratio is ν ¼ 0:25. Due to symmetry, a quarter of the model is selected for computing, see Fig. 11. Fig. 12 shows the BEM model, in which surface of the hole is dispersed into 6 quadratic elements, and 2 in
W.-Z. Feng et al. / Engineering Analysis with Boundary Elements 61 (2015) 274–281
281
results even using a very coarse mesh, demonstrating the accuracy of the presented method.
7. Conclusions
Fig. 13. FEM model.
A new technique for computing the spatial derivative of the global distance r is derived in the paper, which is inevitably used in the power series expansion. And by applying the improved singular boundary integral evaluation method to the stress integral equation, boundary stresses can be directly evaluated from the hypersingular stress boundary integral equations. Numerical examples show that accurate results can be obtained using the proposed method.
Acknowledgments The authors gratefully acknowledge the National Natural Science Foundation of China (11172055, 11202045). References
Fig. 14. Normal stress σ yy over symmetric arc line (z ¼ 0.2).
Fig. 15. Normal stress σ xx over arc line (z¼ 0.4).
the thickness direction. There are total of 338 nodes and 112 boundary elements. To test the accuracy, the FEM results computed by commercial software ANSYS is also given based on a very fine mesh shown in Fig. 13. And another hypersingular integral method proposed in Ref. [19] is also given for comparison. Fig. 14 shows the normal stress σ yy in y direction over the middle arc of the curved surface and Fig. 15 shows the normal stress σ xx over the end arc of the curved surface. Form Figs. 14 and 15, it can be seen that the current results coincide well with FEM
[1] Brebbia CA, Dominguez J. Boundary Elements: An introductory course. London: Mc Graw-Hill Book Co.; 1992. [2] Gao XW, Davies TG. Boundary Element Programming in Mechanics. Cambridge,UK: Cambridge University Press; 2002. [3] Banerjee PK, Raveendra ST. Advanced boundary blement analysis of twodimensional and 3-dimensional problems of elastoplasticity. Int J Numer Methods Eng 1986;23(6):985–1002. [4] Kane JH. Boundary Element Analysis in Engineering Continuum Mechanics. Englewood Cliffs, NJ: Prentice Hall; 1994. [5] Rudolphi TJ. The use of simple solutions in the regularization of hypersingular boundary integral equations. Math Comput Model 1991;15(3):269–78. [6] Chati MK, Mukherjee S. Evaluation of gradients on the boundary using fully regularized hypersingular boundary integral equations. Acta Mech 1999;135 (1–2):41–55. [7] Liu YJ, Rudolphi TJ. Some identities for fundamental solutions and their applications to weakly-singular boundary element formulations. Eng Anal Bound Elem 1991;8(6):301–11. [8] Guiggiani M, Krishnasamy G, Rudolphi T, Rizzo F. A general algorithm for the numerical solution of hypersingular boundary integral equations. ASME J Appl Mech 1992;59:604. [9] Guiggiani M. Hypersingular formulation for boundary stress evaluation. Eng Anal Bound Elem 1994;13(2):169–79. [10] Huber O, Lang A, Kuhn G. Evaluation of the stress tensor in 3D elastostatics by direct solving of hypersingular integrals. Comput Mech 1993;12(1-2):39–50. [11] Frangi A, Guiggiani M. Boundary element analysis of Kirchhoff plates with direct evaluation of hypersingular integrals. Int J Numer Methods Eng 1999;46 (11):1845–63. [12] Frangi A, Guiggiani M. A direct approach for boundary integral equations with high-order singularities. Int J Numer Methods Eng 2000;49(7):871–98. [13] Frangi A, Guiggiani M. Free terms and compatibility conditions for 3D hypersingular boundary integral equations. Z Angew Math Mech 2001;81 (10):651–64. [14] Chen JT, Hong HK. Review of dual boundary element methods with emphasis on hypersingular integrals and divergent series. ASME Appl Mech Rev 1999;52 (1):17–33. [15] Chen JT, Huang WS, Lee JW, Hong HK. On the free terms of the dual BIE for Ndimensional Laplace problems. Eng Anal Bound Elem 2015;59:123–8. [16] Gao XW. An effective method for numerical evaluation of general 2D and 3D high order singular boundary integrals. Comput Methods Appl Mech Eng 2010;199(45):2856–64. [17] Lachat JC. A Further Development of the Boundary Integral Technique for Elastostatics. Ph.D. thesis. Southampton: University of Southampton; 1975. [18] Gao XW. The radial integration method for evaluation of domain integrals with boundary-only discretization. Eng Anal Bound Elem 2002;26(10):905–16. [19] Gao XW, Feng WZ, Yang K, Cui M. Projection plane method for evaluation of arbitrary high order singular boundary integrals. Eng Anal Bound Elem 2015;50:265–74.