Evaluation of singular integrals in the symmetric Galerkin boundary element method

Advances in Engineering Software 35 (2004) 781–789 www.elsevier.com/locate/advengsoft Evaluation of singular integrals in the symmetric Galerkin boundary element method Zhiye Zhao*, Weifeng Yuan School of Civil and Environmental Engineering, Nanyang Technological University, Nanyang Avenue, Singapore, Singapore 639798 Received 16 July 2003; received in revised form 3 June 2004; accepted 2 July 2004 Available online 25 August 2004 Abstract The implementation of the symmetric Galerkin boundary element method (SGBEM) involves extensive work on the evaluation of various integrals, ranging from regular integrals to hypersingular integrals. In this paper, the treatments of weak singular integrals in the time domain are reviewed, and analytical evaluations for the spatial double integrals which contain weak singular terms are derived. A special scheme on the allocation of Gaussian integration points for regular double integrals in the SGBEM is developed to improve the efficiency of the Gauss– Legendre rule. The proposed approach is implemented for the two-dimensional elastodynamic problems, and two numerical examples are presented to verify the accuracy of the numerical implementation. q 2004 Elsevier Ltd. All rights reserved. Keywords: Symmetric Galerkin boundary element; Elastodynamics; Singular integrals 1. Introduction One drawback of the traditional collocation based boundary element method (BEM) is that the system matrices are non-symmetric. Special matrix manipulation is needed to convert the system matrices into symmetric matrices if the symmetry is desirable such as in the coupling of the FEM/BEM approach. The symmetric Galerkin boundary element method (SGBEM), by employing both the displacement integral equation and the traction integral equation, produces a system of symmetric matrices. Due to its symmetric nature, the SGBEM has drawn much attention over the last decade, and it has been studied for a wide range of engineering fields, such as in elastostatics [1–4], coupling of the FEM and the BEM [5], and in elastodynamics [6–8]. One difficulty associated with the SGBEM is the higher order singularity of the kernel functions both in the time domain and in the space domain. For the two-dimensional elastodynamic SGBEM, the fundamental solution associated with the traction due to unit displacement discontinuity (Gpp) contains strong singularities in the time domain. In the space * Corresponding author. Tel.: C65-6790-5255; fax: C65-6791-0676. E-mail address: czzhao@ntu.edu.sg (Z. Zhao). 0965-9978/$ - see front matter q 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.advengsoft.2004.07.004 domain, Guu, Gup and Gpp, which are the fundamental solutions in the two-dimensional elastodynamics, contain singularities in the order of ln r, rK1 and rK2, respectively. In general, analytical approaches are required before the strong and hypersingular integrals can be evaluated numerically [9–12]. Recently, an alternative approach is proposed to deal with the singular double integrals in SGBEM [13], where the strong singular terms are indirectly evaluated through an artificial body force method and the hypersingular integrals are expressed in terms of the strong and weak singular terms. A direct evaluation approach has been just published [14] which deals with the singular integrals as a whole so the symmetry is preserved strictly. In this paper, the weak singular double integrals in the space domain are evaluated based on the locations of the source point and the field point. The double integrals are split into a few terms where the singular terms can be evaluated analytically. For regular double integrals, a new concept of valid segment is proposed to improve the efficiency of the Gauss–Legendre rule, in which the Gaussian integration points are placed on the valid segment only, instead of on the whole element. Numerical implementation based on the proposed evaluation of the singular integrals and the proposed Gaussian integration 782 Z. Zhao, W. Yuan / Advances in Engineering Software 35 (2004) 781–789 scheme for the regular double integrals is carried for the two-dimensional elastodynamic problems. Two examples are used to verify the correctness of the numerical implementation. where Amn ðtÞ Z ð ð ðt 0 JTa ðqÞGmn ðq; s; t K tÞJb ðsÞdtdGn dGm Gm Gn (5) 2. Formulation of the symmetric Galerkin boundary element method Bm ðtÞ Z ð JTa ðqÞg
m ðq; tÞdGm (6) Gm For two-dimensional elastodynamic SGBEM with zero initial conditions, the boundary integral equations can be written in the following form [6,15] ð ðt ð ðt Guu pðs; tÞdtdG K Gup uðs; tÞdtdG 0 0 G G Z uðq; tÞ K ð ðt
tÞdtdU Guu bðs; 0 ð ðt Gpp uðs; tÞdtdG 0 G G Z pðq; tÞ K ð ðt
tÞdtdU Gpu bðs; (2) 0 U mn where G Z Gmn ðq; s; tK tÞ; (m,nZu,p) are the fundamental solutions, u and p denote displacement and traction, up respectively. Guu ij and Gij represent the jth displacement component and traction component at q due to a concentrated unit force acting at s and time t in the ith pp direction; Gpu ij and Gij represent the jth displacement component and traction component at q due to a concentrated unit displacement discontinuity acting at s in pu the ith direction. It should be mentioned that Gup ij and Gij are closely linked by the symmetric relation up Gpu ij ðq; s; n; t; tÞ Z Gji ðs; q; n; t; tÞ (7) with the symmetric condition ATZA. 0 Gpu pðs; tÞdtdG K AX Z B (1) U and ð ðt where the subscripts and superscripts m, n, a, b can take the value of u and p; the term g
m accounts for the given load history, the initial conditions and the boundary conditions [15]. Eq. (4) can be simplified as (3) where n denotes the outward normal of the boundary at point s. In order to implement a numerical scheme to solve Eqs. (1) and (2), it is necessary to consider a set of discrete elements on the boundary G. In the space domain, the displacement u and traction p can be approximated using a set of shape functions Ju and Jp, respectively. In the time domain, the continuous time is divided into a set of discrete values tn, nZ1,2,.,N. Within each time step, u and p are modelled by functions Fu and Fp. As in the classical Galerkin’s weighted approach, shape functions Ju and Jp are selected as weight functions. The discretized state equations are finally obtained by enforcing Eqs. (1) and (2) in a Galerkin weighted-residual sense [6,15] " #( ) ( ) Auu KAup Bu Xp Z (4) KBp KApu App Xu The spatial double integral of Eq. (5) contains weak singular, strong singular and hyper-singular integrals for fundamental solutions Guu, Gup and Gpp, respectively. An efficient artificial body force method has been proposed to evaluate those singular double integrals [13], in which the strong singular and hyper-singular integrals are obtained indirectly through certain identity and the introduction of the artificial body forces. So in this paper, only the singular integrals in the time domain and the weak singular double integrals in the space domain will be discussed in the following sections. 3. Evaluation of temporal integrals in the SGBEM In this section, the temporal integrals derived analytically in Ref. [16] is presented. The results from those temporal integrals will be employed in Section 4 to derive the spatial double integrals encountered in the SGBEM. For the twodimensional elastodynamic SGBEM, all temporal integrals associated with fundamental solutions Guu and Gpu can be expressed in terms of the following four integrals [16] ðt 1 I1 Z Lw Hw dt (8) 0 cw ðt r;i r;j L N H dt (9) I2 Z 2 w w w 0 cw r ðt 1 K1 (10) L Hw dt I3 Z 2 w c 0 wr ðt r 3 I4 Z (11) Lw Hw dt 0 cw where 1 Lw Z pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 cw ðt K tÞ2 K r 2 Nw Z 2c2w ðt K tÞ2 K r 2 (12) (13) Z. Zhao, W. Yuan / Advances in Engineering Software 35 (2004) 781–789 Hw Z Hw ½cw ðt K tÞ K r (14) In the above expressions, cw denotes the wave velocity where the subscript w can take the value of either s or d corresponding to the rotational wave and dilatational wave, respectively; and r is the distance between the source point and the field point. To shorten the expressions, the following notations are defined: pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi T1 Z cw ðt K tnK1 Þ K r (15) T2 Z T3 Z T4 Z pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi cw ðt K tnK1 Þ C r pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi cw ðt K tn Þ K r pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi cw ðt K tn Þ C r (17) (18) (19) T6 Z cw ðt K tn Þ (20) In a time step [tnK1,tn], because the temporal integral kernels are discontinuous due to the characteristics of Heaviside function, the final expressions of the integrals I1 to I4 depend on three cases [16]. r Case 1 : t K ! tnK1 cw In this case, because the value of Heaviside function is zero, all temporal integrals I1 to I4 are zero. r Case 2 : tnK1 ! t K ! tn cw I1 Z I2 Z ð tn tnK1 tnK1 1 1 L H dt Z 2 ½lnðT5 C T1 T2 Þ K ln r cw w w cw r;i r;j r;i r;j Lw Nw Hw dt Z 2 2 T1 T2 T5 cw r 2 cw r Case 3 : tn ! t K (21) ð tn 1 1 T T C T5 L H dt Z 2 ln 1 2 cw w w cw T3 T4 C T4 I2 Z ð tn r;i r;j r;i r;j Lw Nw Hw dt Z 2 2 ðT1 T2 T5 K T3 T4 T6 Þ 2 cw r cw r tnK1 tnK1 4.1. Analytical evaluation of the spatial integrals The typical integral in the SGBEM formulation has the form of ðð Iw Z Kðr; t; tÞdGs dGf (29) Gf Gs where the integration kernel K(r,t,t) could be one of I1 to I4 defined in the previous section. By using the artificial body forces approach [13], the hypersingular double integrals associated with Gpp can be evaluated using an indirect method. So the discussion focuses on the evaluation of weak and strong singular integrals. The integration kernel of Iw is not continuous because of the characteristic of the Heaviside function. For easy evaluation, the integral is divided into three parts (30) where Iw1 Z (22) ð ð Kðr; t; tÞdGs dGf (31) ð ð Kðr; t; tÞdGs dGf (32) ð ð Kðr; t; tÞdGs dGf (33) Gf Gs1 Iw1 Z Gf Gs1 (23) Iw3 Z (24) r cw I1 Z (28) 4. Evaluation of spatial double integrals Iw Z Iw1 C Iw2 C Iw3 ð tn 1 K1 L Hw dt I3 Z 2 w tnK1 cw r  1 1 Z 2 2 T1 T2 T5 K lnðT5 C T1 T2 Þ C ln r 2cw r ð tn r 3 L 1 T1 Lw Hw dt C 2w Z K 2 I4 Z cw cw r T2 tnK1 cw (27) (16) T5 Z cw ðt K tnK1 Þ ð tn ð tn 1 K1 L Hw dt 2 w c tnK1 w r  1 1 T C T1 T2 Z 2 2 ðT1 T2 T5 K T3 T4 T6 Þ K ln 5 T6 C T3 T4 2cw r ð tn r 3 1 T6 T I4 Z Lw Hw dt Z 2 K 5 c T T T c r tnK1 w 3 4 1 T2 w I3 Z 783 (25) (26) Gf Gs3 It should be noted that cw(tKtnK1)!r is on the boundary pair GfwGs1, while on the boundary pairs GfwGs2 and GfwGs3 we have cw ðtK tn Þ! r! cw ðtK tnK1 Þ and r!cw(tKtn), respectively. The expressions of integration kernels change for different forms due to different integration domains. Two typical double integrals in SGBEM are discussed in detail. ðð 1 Type 1 : IT1 Z ln r dLs dLf (34) c2w Lf Ls 784 Z. Zhao, W. Yuan / Advances in Engineering Software 35 (2004) 781–789 Type 2 : IT2 Z ðð Lf Ls 1 T1 T2 T5 dLs dLf c2w r 2 (35) where T1, T2 and T5 are defined in Eqs. (15), (16) and (19), respectively.  l K q cos q lim q cos q ln q C q sin q arctan s q/0 q sin q  cos q Z0 Carctan sin q Therefore, the outer integral of IT1 can be computed by Gauss–Legendre rule. It must be noted that the special case for qZp is also covered in the derivation. 4.2. Evaluation of Type 1 integrals To evaluate IT1, three cases have to be dealt with depending on the reciprocal position of the two element segments: (a) distinct element segments; (b) adjacent element segments; (c) coincident element segments. 4.2.1. Type 1 integral for the distinct element segments The source segment of length ls and the field segment of length lf are separated. Because r is always greater than zero, the singularity is not activated thus the standard Gauss–Legendre rule may be used. 4.2.3. Type 1 integral for the coincident element segments The source element segment and the field element segment are the same one, as shown in Fig. 2. The distance between the source point and the field point can be expressed as r Z js K qj 4.2.2. Type 1 integral for the adjacent element segments The source segment and the field segment share one common end where r becomes zero, as shown in Fig. 1. The distance between the source point S and the field point F can be expressed as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r Z ðs K q cos qÞ2 C q2 sin2 q (36) (37) 0 1 Z ½ðls K q cos qÞlnðl2s K 2ls q cos q C q2 Þ K ls  2 C q cos q ln q  l K q cos q cos q C arctan C q sin q arctan s q sin q sin q (38) It should be noticed that when q/0, Eq. (38) does not contain any singularity, because (40) and IT1 Z then we have ð lf ð ls qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 ln ðs K q cos qÞ2 C q2 sin2 q dsdq IT1 Z 2 0 0 cw The analytical result of the inner integral is ð ls qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ln ðs K q cos qÞ2 C q2 sin2 q ds ð39Þ ðl ðl 0 1 lnjs K qjdsdq 2 c 0 w (41) The inner integral in IT1 is given ðl lnjs K qjds Z Kl C ðl K qÞlnjl K qj C q ln q (42) 0 Because limq/0 ðq ln qÞZ 0 and limq/l ½ðlK qÞlnjlK qj Z0; it is clear that no singularity is contained in the outer integral of IT1, thus the outer one can be calculated by means of the Gauss–Legendre rule. 4.3. Evaluation of Type 2 integrals The integral IT2 can be rewritten as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðð ðtKtnK1 Þ c2w ðtKtnK1 Þ2Kr 2Kcw ðtKtnK1 Þ dLs dLf IT2Z cw r 2 Lf Ls C ðð ðtKtnK1 Þ2 dLs dLf r2 ð43Þ Lf Ls On the right-hand side of Eq. (43), the first term is not singular, so only the second term will be discussed in the following. 4.3.1. Type 2 integral for the distinct element segment Since the distance between the source point and the field point rO0 always holds for this case, no singularities appear in the integrals. So IT2 can be evaluated numerically. Fig. 1. Reciprocal locations of two adjacent segments. Fig. 2. Reciprocal location of two coincident segments. 785 Z. Zhao, W. Yuan / Advances in Engineering Software 35 (2004) 781–789 4.3.2. Type 2 integral for the adjacent element segment The distance between the source point and the field point can be expressed as rZ and qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðs K q cos qÞ2 C q2 sin2 q Ihyp Z Z ð lf ð ls 0 ð lf 0 Z 0 (44) ð lf ls K k cos q p Z k sin q 2 cos q z2 Z arctan sin q (45) (46) 3/0 ðz1 C z2 Þ ln 3 sin q ðl ðl 1 Ihyp Z dsdq Z ðs K qÞ2 0 0 l  Z ðlnjl K qj K ln qÞ ðl l dq qðq K lÞ 0 (48) (53) o/0 So the finite-part integral is defined by: ðl ðl 1 Ifin Z dsdq 2 0 0 ðs K qÞ Z Ihyp K limC ln 3 K limC ln o Z K2 ln l 3/0 The term lim3/0C ðz1 C z2 Þ=sin q ln 3 is infinite, but it should be deleted from the formula if two kinds of solid waves are considered together. For convenience, a finite part integral is defined ð lf ð ls 1 Ifin Z dsdq 2 2 2 0 0 ðs K q cos qÞ C q sin q Z Ihyp C limC (52) 3/0 (47) ðz1 C z2 Þ ðz C z2 Þ ln lf K limC 1 ln 3 sin q sin q 3/0 4.3.3. Type 2 integral for the coincident element segments The distance between the source point and the field point is Z limC ln 3 C limC ln o K 2 ln l On the right-hand side of Eq. (45), the first term is regular thus it can be evaluated numerically. The second term is simple enough to be evaluated analytically as lf ð lf  1 1 ðz1 C z2 Þdq Z ðz1 C z2 Þln q sin q 0 q sin q 0 Z 0 0 where k/0 0 1 ll dsdq Z Ihyp C limC 3 Z ln s f 2 ðls C lf Þ ðs C qÞ 3/0 (51) and 1 l K q cos q cos q arctan s C arctan dq q sin q q sin q sin q z1 Z lim arctan ð lf ð ls r Z js K qj 1 dsdq ðs K q cos qÞ2 C q2 sin2 q 1 l K q cos q arctan s K z1 dq q sin q 0 q sin q ð lf 1 ðz1 C z2 Þdq C 0 q sin q Ifin Z (54) o/0 It must be noticed that the definition of finite part integral has a clear physical meaning. To define the finite part integral, some infinite terms, for instance, lim3/0C 3 in Eq. (50), are directly deleted. This operation actually indicates the superposition of two solid waves. As limr/0 H½cd ðtK tÞK rZ limr/0 H½cs ðtK tÞK r; the influence due to the dilatational wave and the influence due to rotational wave will overlap at the field point in the elastic domain. Each kind of solid wave results in an infinite integral. The two infinite terms have the same value but with different signs, thus they can be eliminated. 4.4. A special scheme for Gaussian integration If the integration kernel is not singular or it is only weak singular, the outer integrals of Eq. (29) can be carried out using the Gauss–Legendre rule. In order to improve (49) The special case for two adjacent elements is rZsCq, when qZp. The related Ihyp and Ifin are given as: ð lf ð lf ð ls 1 ls dq Ihyp Z dsdq Z 2 0 tðls C qÞ 0 0 ðs C qÞ Z ln ls lf K lim 3 ðls C lf Þ 3/0C (50) Fig. 3. An example on the determination of valid field segment. 786 Z. Zhao, W. Yuan / Advances in Engineering Software 35 (2004) 781–789 Table 1 Comparison between two integration schemes Scheme 2!2 Nodes 3!3 Nodes 4!4 Nodes 5!5 Nodes Standard Gaussian scheme Proposed Gaussian scheme 0.8452995 0.1391379 0.2149339 0.3479466 0.3333333 0.3333333 0.3333333 0.3333333 the determination of valid segments for both the source element and the field element can be found in Ref. [17]. For testing purpose, one numerical example is presented to check the validity and accuracy of the proposed scheme. Consider an integral ð2 ð2 r$Hð1 K rÞdsdq (55) Itest Z 0 0 Suppose the source element and the field element are adjacent to each other, and the angle between them is p, then the test integral is simplified to be: ð2 ð2 Itest Z ðs C qÞ$Hð1 K s K qÞdsdq (56) 0 Fig. 4. Loading history of example 1. the integration accuracy, a scheme is proposed to place the Gaussian points on certain part (called valid segment) of the element instead of the whole element. In the proposed integration scheme, the field element is divided into several segments among which one is named valid segment, the rests are called invalid segments. Gaussian integration points are only generated on the valid segment. One example on the delimitation of the valid segment is depicted in Fig. 3, where Es and Ee are the start point and the end point of the source element, respectively, while GsGe is the field element. The radius of the two circles are RwZcwt 0 , where t 0 Z(tKtnK1) and t 0 Z(tKtn), respectively. The field element and the boundary of the shadow domain have two intersections js and je, so the valid element is jsKje. As the corresponding integration kernels in the double integrals have non-zero values only on the valid segment of the element, so the integration points placed outside the valid segment of the element will have no effect on the final value of the double integral. Therefore, in this study, the Gaussian integration points are placed on the valid segment only to improve the efficiency of the Gauss–Legendre rule. Full details on 0 This double integral is simple enough to be evaluated analytically, and the exact final value is 1/3. Table 1 shows the comparison between the standard Gaussian integration scheme (i.e. placing the integration points on the whole element) and the proposed integration scheme (i.e. placing the integration points on the valid segment only). Based on the Guassian–Legendre rule, the more Gaussian integration nodes are hired, the better the numerical result is. For the standard Gaussian integration scheme, the results have been improved from very poor when 2!2 nodes are used to reasonable when 5!5 nodes are used. For the proposed scheme, the exact value can be achieved even for the case of 2!2 nodes. So the proposed integration scheme is highly efficient as compared with the standard integration scheme. 5. Numerical examples A SGBEM computer code is developed for the twodimensional elastodynamic problems. In this code, the singular integrals are evaluated by different approaches: (a) the temporal integrals are evaluated analytically as in Ref. [16]; (b) the strong singular double integrals and the hypersingular double integrals are carried out by an indirect method as presented in Ref. [13]; (c) the weak singular double integrals are evaluated based on the method presented in the previous section. Two numerical examples are presented to verify the numerical implementation. Fig. 5. Boundary discretization and distribution of artificial nodes: (a) 40 artificial nodes; (b) 60 artificial nodes. Z. Zhao, W. Yuan / Advances in Engineering Software 35 (2004) 781–789 787 the strong singular and hypersingular double integrals [13]. In order to show that the numerical results are not sensitive to the location and the number of the artificial nodes, different d values are used and two sets of artificial nodes (40 nodes and 60 nodes) are employed. For comparison purpose, two FEM models for the same problem is generated using software ABAQUS. One FEM model, one-quarter of the region with 300 linear plain strain elements (CPE4), and the other model uses 120 CAX4 type elements under the axis-symmetric problem. In this example, the time steps used for the SGBEM and the FEM are 25 and 4 s, respectively. The radial displacement history of the cavity’s surface is plotted in Fig. 6. Considering the wave reflection, the FEM analysis stops when the wave reaches the outer edge of the region. The results show a very good agreement between the FEM results and the SGBEM results. The differences by using different number of the artificial nodes and different locations of the artificial loads are negligible. After the internal pressure vanishes, the displacement finally tends to zero because the input energy spreads to remote area as time increases. Fig. 6. Nodal displacements. 5.2. A plane strip under a Heaviside type loading Fig. 7. Illustration of plane strip with distributed Heaviside traction. 5.1. Cylindrical cavity under uniform pressure The example simulates the propagation of the transient compressive wave emanating from a cylindrical cavity of radius RZ100 m to the surrounding unbounded isotropic elastic solid. The material constants are set as: Young’s modulus of EZ1.0 Pa, Poisson ratio nZ0.0, density rZ1.0 kg/m3. This problem is regarded as a plain strain problem. The pressure p(t), applied on the boundary surface of the cavity is supposed to be uniform and modeled in time by a piecewise-defined function, as plotted in Fig. 4 where p0Z1.0!10K3 Pa and t2Z2, t1Z200 s. As shown in Fig. 5, the boundary of the cavity is discretized into 20 constant elements, thus the length of each element can be approximately expressed as lzRp/10. Around the boundary, the artificial nodes are distributed on a circle of the radius RaZRCd for the purpose of evaluating The problem considered is a rectangular domain with the side lengths a and b (bZ4a) as shown in Fig. 7. On the boundary yZ0, the displacement u is fixed in the x direction. At the boundary yZa, a loading p0 is suddenly applied at tZ0 and kept constant until the analysis stops. On the other two boundaries, where xZKb/2 and b/2, displacement is free and traction is taken as zero. The Poisson ratio is taken as nZ0.0. The boundary is discretized into 80 constant elements with element’s length lZa/8, as shown in Fig. 8. The constant time step Dt is chosen by the following formula Dt Z b l cd (57) where b is a constant factor. The displacement responses of four elements e1, e2, e3 and e4 for bZ1.0 are shown in Fig. 9, and are compared with the analytical results for a plane strip under a Heaviside type loading. It can be observed that the results by the SGBEM forms excellent periodical patterns, and are close to the analytical results. The traction time-history of element Fig. 8. Boundary discretization of the plane strip. 788 Z. Zhao, W. Yuan / Advances in Engineering Software 35 (2004) 781–789 Fig. 9. Displacement time-history at nodes e1, e2, e3 and e4. e5 are plotted in Fig. 10 which shows very good accuracy in the first cycle. The accuracy drops slightly after the wave reflection reaches the element again, but still follows the cyclic pattern well. 6. Conclusions This paper concerns with the evaluation of the singular integrals in the two-dimensional elastodynamic SGBEM. Though various approaches for the treatments of the singular integrals have been developed, very few numerical examples have been published, mainly due to the complexity of the various integrals involved in the SGBEM. The strong singular and the hypersingular singular double integrals have been successfully evaluated in Ref. [13]. In this paper, the weak singular double integrals are split into two parts: the singular part is evaluated by analytical Fig. 10. Traction time-history at node e5. approach and the regular part is calculated by direct numerical evaluation. The numerical examples show that good accuracy can be obtained by the proposed SGBEM implementation. The authors are currently working on further improvement of the SGBEM implementation in the following two areas: (a) to extend the current constant element into linear elements in both the space domain and the time domain; (2) to use a more stable time integration scheme. References [1] Hartmann F, Katz C, Protopsaltis B. Boundary elements and symmetry. Ing Arch 1985;55:440–9. [2] Polizzotto C. An energy approach to the boundary element method. Part I. Elastic solids. Comput Meth Appl Mech Eng 1988;69: 167–84. [3] Sirtori S, Maier G, Novati G, Miccoli S. A Galerkin symmetric boundary element method in elasticity: formulation and implementation. Int J Numer Meth Eng 1992;35:255–82. [4] Bonnet M. Regularized direct and indirect symmetric variational BIE formulation for three dimensional elasticity. Eng Anal Bound Elem 1995;15:93–102. [5] Sirtori S, Miccoli S, Korach E. Symmetric coupling of finite elements and boundary elements Advances in boundary element techniques. Berlin: Springer; 1993, p. 407–27. [6] Maier G, Diligenti M, Carini A. A variational approach to boundary element elastodynamic analysis and extension to multidomain problems. Computer Meth Appl Mech Eng 1991;92: 193–213. [7] Yu G, Mansur W, Carrer J, Gong L. Time weighting in time domain BEM. Eng Anal Bound Elem 1998;22:175–81. [8] Zhao Z, Yuan W, Lie S, Yu G. Symmetric Galerkin BEM for dynamic analysis Proceedings of the First Asian-pacific Congress on Computational Mechanics. Sydney: Elsevier; 2001, p. 1547–52. [9] Carini A, Diligenti M, Maranesi P, Zanella M. Analytical integrations for two- dimensional elastic analysis by the symmetric Galerkin boundary element method. Comput Mech 1999;23:308–23. Z. Zhao, W. Yuan / Advances in Engineering Software 35 (2004) 781–789 [10] Frangi A, Novati G. Symmetric BE method in two-dimensional elasticity: evaluation of double integrals for curved elements. Comput Mech 1996;19:58–68. [11] Toh K, Mukherjee S. Hypersingular and finite part integrals in the boundary element method. Int J Solids Struct 1994;31:2299–312. [12] Aimi A, Diligenti M, Monegato G. Numerical integration schemes for the BEM solution of hypersingular integral equations. Int J Numer Meth Eng 1999;45:1807–30. [13] Yuan W, Zhao Z, Lie S, Yu G. Numerical implementation of the symmetric Galerkin boundary element method in 2D elastodynamics. Int J Numer Meth Eng 2003;58(7):1049–69. 789 [14] Bonnet M, Guiggiani M. Direct evaluation of double singular integrals and new free terms in 2D (symmetric) Galerkin BEM. Computer Meth Appl Mech Eng 2003;192:2565–96. [15] Bonnet M, Maier G, Polizzotto C. Symmetric Galerkin boundary element methods. Appl Mech Rev 1998;51:669–704. [16] Carrer J, Mansur W. Stress and velocity in 2D transient elastodynamic analysis by the boundary element method. Eng Anal Bound Elem 1999;23:233–45. [17] Yuan W. Symmetric Galerkin boundary element method in elastodynamic analysis. PhD Thesis, Nanyang Technological University, Singapore; 2003.