Singularity cancellation method for time-domain boundary element formulation of elastodynamics: A direct approach

Journal Pre-proof

Singularity cancellation method for time-domain boundary element formulation of elastodynamics: a direct approach Guizhong Xie , Yudong Zhong , Fenglin Zhou , Wenliao Du , Hao Li , Dehai Zhang PII: DOI: Reference:

S0307-904X(19)30736-X https://doi.org/10.1016/j.apm.2019.11.053 APM 13183

To appear in:

Applied Mathematical Modelling

Received date: Revised date: Accepted date:

14 April 2019 10 November 2019 28 November 2019

Please cite this article as: Guizhong Xie , Yudong Zhong , Fenglin Zhou , Wenliao Du , Hao Li , Dehai Zhang , Singularity cancellation method for time-domain boundary element formulation of elastodynamics: a direct approach, Applied Mathematical Modelling (2019), doi: https://doi.org/10.1016/j.apm.2019.11.053

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Inc.

Highlights  A kernel separation method is proposed to compute the singular integral accurately.  Four singular integrals with different ordered kernel are computed directly with the help of the kernel separation method.  A patch subdivision scheme is developed based on the relation between the product of wave velocity, time and the distance.

1

Singularity cancellation method for time-domain boundary element formulation of elastodynamics: a direct approach Guizhong Xiea, Yudong Zhonga, Fenglin Zhou b,*, Wenliao Dua, Hao Lia, Dehai Zhanga a

Henan Key Laboratory of Intelligent Manufacturing of Mechanical Equipment,

Mechanical and Electrical Engineering Institute, Zhengzhou University of Light Industry, Zhengzhou, 450002, Henan, China b

College of mechanical engineering, Hunan university of technology, Zhuzhou,

412007, Hunan, China * Corresponding Author: [email protected] Abstract In the implementation of time-domain boundary element method for elasto-dynamic problems, there are two types of singularities: the wave front singularity arising when the product of wave velocity and time is equal to the distance between the source point and the field point, and the spatial singularity arising when the source point coincides with the field point. In this paper, the singularity of the first type in the integrand is eliminated by an analytical integration over time, Cauchy principal value and Hadamard finite part integral. Four types of singularities with different orders appear in the integrand after analytical time integration. In order to accurately calculate the integral, in which the integrand is piecewise continuous, the integral domain is subdivided into several patches based on the relation between the product of wave velocity and time and the distance. In singular patches, the integrands are separated into a regular part and a singular part. The regular part can be computed by traditional numerical integration method such as Gaussian integration, while the singular part can be analytically integrated. Using the proposed method, the spatial singular integrals can be calculated directly. Numerical examples using various kinds of elements are presented to verify the proposed method. Keywords: time-domain boundary element method; elasto-dynamics; singularity; singular integrals. Introduction

In the structural design, the transient response of a structure should be considered carefully. Structural dynamic problems are usually governed by several second ordered partial differential equations of the hyperbolic type. In these problems, transient responses are driven by the boundary load and the initial condition [1]. Analytical solution to the structural dynamic problem is available only in the case of a very simple geometry and a very simple load. Thus, three types of numerical methods including the finite element method (FEM) [2-4], the boundary element method (BEM) [5-10] and meshless method [11-12] have been widely applied to solve these structural dynamic problems. In those numerical methods, the BEM exhibits many advantages such as dimension reduction, high accuracy and less restriction on elements continuity. Lots of applications of BEM in solution of engineering problems could be found. In the applications of BEM for transient problems, four types of approaches, which classified by J.A.M. Carrer, W.L.A. Pereira and W.J. Mansur in [13], are mostly adopted. These approaches are time-domain boundary element method (TD-BEM) in which the time-dependent fundamental solutions are employed [5-8], formulations in which the static fundamental solutions are applied [14], the BEM formulations based on Laplace and frequency domains [15, 16], and the Convolution Quadrature Method [17, 18], respectively. Among these approaches, the TD-BEM formulations are considered to be the most natural choice in the BEM application. In series works of Niwa et al [19], W.J. Mansur and Brebbia [5, 20], A. S. M. Israil and P.K. Banerjee [6, 21], a time-stepping scheme had been developed in the TDBEM formulation. In those works, two types of inherent singularities were found in the TD-BEM formulations, namely, the wave front singularity and the spatial singularity. The wave front singularity can be coped with through carrying out the time integration analytically. It should be noted that the Cauchy principal value and Hadamard finite part integral and the nature of Heaviside function are employed in the derivative process. To deal with the spatial singularity, two types of approaches are employed, namely, the indirect method [5, 20, 21, 22] and the direct method [22]. In the implementations of both indirect method and direct method, the time integrals are calculated analytically at first. In the indirect method, the spatial singularity is cancelled by subtracting elasto-static kernels which have the same singularity as the convoluted kernels. Then the integrals that is regarding to elasto-static kernels can be obtained by the rigid body displacement. This method is valid only in the case of the

large time steps. Moreover, for crack bodies or elastic-plastic problems, the rigid body displacement method may fail. In the direct method, the spatial singularity of the integrals with the convoluted kernels can be computed analytically. The computation of the singular integral, however, may be not an easy task. Until now, the available implementations could only be found in the cases of constant element or linear element. In the case of curved element such as the quadratic isoparametric element, the procedure of the integral calculation is very complicated. In this paper, a kernel separation method will be presented for computation of the singular integral that is involved in the direct method. In the kernel separation method, the integrals with the convoluted kernels are separated into two parts, which include the singular part and the regular part, through the Taylor expansion around the singular point in the isoparametric coordinate system. The singular part can be easily computed. It is worth noting that the convoluted kernel may be piecewise continuous. The consistency of the integrands in the singular integration interval should be guaranteed. In order to circumvent this problem, a subdivision method is proposed based on the relations between the product of wave velocity and time and the distance. Two benchmark numerical examples will be presented to illustrate the validity and accuracy of the proposed method. 2.Boundary integral formulations for TD-BEM The boundary integral equation for transient elasto-dynamics can be obtained by Graffi reciprocity theorem as follows, without considering initial conditions and body forces [5], clk uk (S, t )    ulk*  S, , F, t  pk ( , F)d d     plk*  S, , F, t  uk ( , F)d d  

t

  v 

* 0 lk

t

t

0

0



(1)

 S, , F, t vk ( , F)d d 

In equation (1), clk is a free term computed from Cauchy principal value integral, and clk  0.5 lk ,when  is smooth. S is the source point and F is the field point at the

boundary. uk , pk and vk , k  1..2 represent the displacements, tractions and velocities on the boundary  at time  . ulk*  S, , F, t  , plk*  S, , F, t  , and vlk*  S, , F, t  are the 2D fundamental solutions, which are

 1  r2 2 R22  r 2  H (c2 (t   )  r )  (  R )   r r    2 lk ,l , k  R2  r2 1  c2  R2  ulk*  S, , F, t     2  1  2 R12  r 2  H (c1 (t   )  r )   2  c   ik R1  r,l r,k  R r 1   1 

(2)

 3c1 (t   )r 2  2c13 (t   )3  r 3 2 R12  r 2  D  B lk  lk  H (c1 (t   )  r )  2 c1  R13 R1  (3) 2 2 2 3 3 3 c2 (t   )  r 2R  r 3c (t   )r  2c2 (t   )  r  1   Blk 2  Dlk 2   Alk  H (c2 (t   )  r ) 3 2 c2  R2 R2 R23  1

plk*  S, , F, t  

 r2  Dlk  2 R1   H (c1 (t   )  r ) 2 c R1   2  1  1 r    Alk  Dlk  2 R2    H (c2 (t   )  r ) 2  2 c2  R2 R2   

vlk*  S, , F, t  

1

2 1

(4)

In equations (2), (3) and (4), c1 and c2 represent the dilatational and equivoluminal body wave velocities. R1  c12 (t   )2  r 2 , R2  c22 (t   )2  r 2 , H (c1 (t   )  r ) and

H (c2 (t   )  r ) are the Heaviside functions.   Blk 

2G  r r    4 r,l r,k  nk r,l  nl r,k  , 3  lk r  n n 

Dlk 



r  , Alk  G  2 nk r,l   lk  nl r,k  , 2G n  

2G r ( nk r,l  r,l r,k ) . 2 r n

Note that unlike many other applications of the TD-BEM for elasto-dynamic prob lems, that fundamental solution for the velocity is not involved. From references [5, 7], the fundamental solution for the velocity comes into being through removing the derivatives of Heaviside function in the integrals. It should be noted that the time t is divided into n steps, assume that the displacements vary linearly with respect to time and the tractions stay steady in each time scan. The following integral with respect to velocity fundamental solutions is  r2  Dlk  2 R1   H (c1 (t   )  r )vk d tm 2 c tm R1   2 tm1   1 1 r    Dlk  2 R2    H (c2 (t   )  r )vk d   Alk 2 t 2 c2 m  R2 R2    I 

tm1

vlk*  S, , F, t  vk d 

1

2 1



tm1

(5)

where vk ( , F) can be obtained by vk    (m  1)t   ukm    mt ukm1     ukm1  ukm  

t

t



t

,

in which ukm 1 and ukm are the displacements at the time scan m  1 and m , and t is the length of the time scan.. It should be noted that from the work of Cole and his cooperators and Dominguez J [7, 24], they consider that if the piecewise linear time interpolation functions are used

for tractions, as for displacements, the numerical solution process becomes unstable in certain cases. Thus, we have given up both the linear time interpolation functions for the tractions and displacements. From the work of Peirce A and Siebrits E, they consider that there are restrictions governing the choice of functional variation in space and time. A piecewise constant functional variation in time (in two and three dimensions) is not possible because it leads to infinite stresses at the wave fronts [25]. Moreover, from the work of Frangi A and Novati G, they have found that when a longitudinal wave keeps travelling (and thus reflecting) between two sides of a strip accelerations are concentrated along the moving line of the wave front (where they assume infinite values) and are zero elsewhere, which could accounts for the piecewise linear and constant behaviour of displacements and tractions, respectively [26]. Thus we have chosen two different approximations for the boundary displacement and boundary traction. After the time discretization for Eq. (1) and using Eq. (5), the following integrals of time considering the 2D fundamental solution can be obtained: tm1

I1  

tm

I2  

tm1

I3  

tm1

I4  

tm1

tm

tm

tm

ulk*  S, , F, t  d

(6)

plk*  S, , F, t 

  mt

plk*  S, , F, t 

(m  1)t   d t

t

d

(7) (8)

vlk*  S, , F, t  d t

(9)

It should be noted that another integral form still existed, which are the linear combinations of the four integral forms above. In order to simplify the indication, four symbols are used below,

Sc1tm  c1 (t  tm ) Sc1tm1  c1 (t  tm1 )

Sc 2tm  c2 (t  tm ) Sc 2tm1  c2 (t  tm1 )

S1  c1 (t   ) S2  c2 (t   )

(10) (11) (12) (13) (14) (15)

Exact expressions of I1 , I 2 , I 3 and I 4 are listed in Appendix A, B, C, D.

3. Singularity cancellation method for TD-BEM It should be noted that the formulation in appendix A, B, C and D are deduced in the sense of Cauchy principal value and Hadamard finite part integrals. The divergent terms related to S1  r  0 or S2  r  0 are removed directly, thus the wave front singularities disappear. When carrying out the analytical integrals fully, the space discretization can be performed by employing constant, linear or quadratic isoparametric elements. Coordinates xi , i  1..2or 3 at arbitrary point inside those elements can be computed through

xi  Nk ( ) xik

(40)

where k  1..2 for constant and linear elements, while k  1..3 for quadratic elements.

 is the parametric coordinate.  xk , yk  is the nodal global coordinate. N k ( ) is the shape functions. After the spatial discretization, the spatial integrals involved in I1 , I 2 , I 3 and I 4 can be expressed as I1_ s   I1   N k   Jbd

(41)

I 2_ s   I 2   N k   Jbd

(42)

I3_ s   I3   N k   Jbd

(43)

1

1 1

1

1

1

I 4_ s   I 4   N k   Jbd 1

(44)

1

Where Jb stands for the Jacobian of the parametric transformation. In the following subsections, we will apply a singularity separation scheme to separate the singular terms from the kernel of the integral. During the separation, the following approximation should be applied

lim(1  x)q  1  qx x 0

(45)

3.1 The singular terms in I1_ s In the expressions of I1_ s in Eq. (46)-Eq. (51), the singular kernel of the integrand can be described in appendix E. It can be found that log(r) or 1/r2 type singular terms appear in the I1 . In the integral calculation, both types of singularity should be considered. 3.2 The singular terms in I 2 _ s

The singular kernel in I 2 _ s is listed in appendix F. It is much more complicated than that in I1_ s . The influence factors Blk and Dlk exhibit singularities with orders of 1/r3 and 1/r2, respectively. As illustrated in Eq. (52)-Eq. (57), four different orders of singularity as log(r), 1/r, 1/r2 and 1/r3 may be involved in the kernel of I1_ s . Thus, it is not an easy task to compute I 2 _ s accurately. 3.3 The singular terms in I 3_ s The singular kernel in I 3_ s is listed in Eq. (58)-Eq. (63). The order of singular kernel in I 3_ s is the same as that I 2 _ s . 3.4 The singular terms in I 4 _ s The singular kernel in I 4 _ s is listed in Eq. (64)-Eq. (69) in Appendix H. The influence factors Dlk exhibits singular order of 1/r2. In the kernel of I 4 _ s , singularities of orders of log(r) and 1/r2 are involved. 3.5 Taylor expansion of the singular kernels at the source point in the intrinsic coordinate system A source point is denoted by y with intrinsic coordinate  0 . One arbitrary point in the element is denoted by x . xi , yi are components in coordinates of x and y . The difference of these components can be expanded into Taylor series [27-30]. xi  yi 

dx d

 0

  0  

where   0   , Ai 

dx d

1 d 2x 2 d 2

  0

 0

, Bi 

  0 

1 d 2x 2 d 2

2

 0

 O   0   Ai  Bi 2  O( 3 ) (70) 3

.

Thus, the square of the distance between x and y can be described in Eq. (71)

r 2   xi  yi   A2 2 +2C 3 +O( 4 ) 2

(71)

1

where A   Ai Ai  2  J (0 ) , C  Ai Bi . The distance r and its n powers can be further expressed as r  A2 2 +2C 3 +O( 4 )  A 

1+

2C C  +O( 2 )  A  (1  2 )  O( 3 ) 2 A A

C C n   r   A  (1  2 )  O( 3 )   An  (1  n 2  )  O( n  2 ) A A   n

n

The component of r in the global coordinate system can be written as

(73)

(72)

r,i 

xi  yi A   Bi 2  O( 3 ) A B AC  i  i sgn( )  ( i  i 3 )   O( 2 ) C r A A A  (1  2 )  O( 3 ) A A

(74)

The shape function can also be expanded in the center of  0

N ( )  N (0 ) 

dN d

 0

(  0 )  O[(  0 ) 2 ]  N0 (0 )  N1 (0 )  O( 2 ) (75)

Considering the geometry relations between the normal and Jacobian, the following equations can be obtained J1 ( )  J ( )n1 (x) 

dx2 d

J 2 ( )  J ( )n2 (x)  

(76a)

dx1 d

(76b)

Using the above basic formulas from (69) to (75), the integrands of I1_ s , I 2 _ s , I 3_ s and I 4 _ s can be reorganized as ExI1_ s  F1ij (0 ) ln   0  G1ij (0 )

1 1  H1ij (0 ) 2   0   0 

1 1 1  H 2ij (0 )  K 2ij (0 ) 2 3   0   0     0 

(78)

1 1 1  H 3ij (0 )  K 3ij (0 ) 2 3   0    0     0 

(79)

ExI 2_ s  F 2ij (0 ) ln   0  G 2ij (0 ) ExI 3_ s  F 3ij (0 ) ln   0  G3ij (0 )

(77)

ExI 4_ s  F 4ij (0 ) ln   0  G 4ij (0 )

1

  0

 H 4ij (0 )

1

  0 

(80)

2

In the sense of Cauchy principal value and Hadamard finite part integrals, singular terms in Eqs. (77)-(80) can be integrated analytically when 1  0  1 .



1

1

ln   0 d   (1  0 ) ln(1  0 )  (1  0 ) ln(1  0 )  2 

1  0 1 d  ln 1    1  0 0



1

1

1

     1

2

d  

0

1

1

     1

0

3

d 

(82a)

1 1  1  0 1  0 1

2 1  0 

(81a)

2



(83a)

1 2 1  0 

2

(84a)

In the sense of Cauchy principal value and Hadamard finite part integrals, singular terms in Eqs. (77)-(80) can be integrated analytically when 0  1 .



1



1

1

1

1

ln   1d  2ln  2   2 1 d   ln  2    0

1

     1

2

d  

0

1

1

     1

3

0

d 

1 2

1 8

(81b) (82b)

(83b)

(84b)

In the sense of Cauchy principal value and Hadamard finite part integrals, singular terms in Eqs. (77)-(80) can be integrated analytically when 0  1 .



1



1

1

1

1

ln   0 d  2ln  2   2 1 d  ln  2    0

1

     1

2

1

1

1

0

3

(82c)

d  

1 2

(83c)

d  

1 8

(84c)

0

    

(81c)

More details about the truncation scheme can be found in [27-30]. In order to numerically calculate the singular integrals, we need further to subdivide the intervals for singular integrals considering Sc1tm , Sc1tm1 , Sc 2tm , Sc 2tm1 and the length of the integration element. The ratio of minimum non-zero absolute value of the five parameters over the length of the integration element was usually taken as the subdivision criterion. This subdivision scheme can keep the consistency of I1 , I 2 , I 3 and I 4 in the intervals for singular integrals. Note that The expressions of I1 , I 2 , I 3 and

I 4 are divided into 6 segments considering the relations of Sc1tm , Sc1tm1 , Sc 2tm , Sc 2tm1 and r, respectively. When a space integral element is cut by the wave front, the integral element inside and outside the front edge should employ the fundamental solutions over different moments of time, thus the integration accuracy can be maintained. Using Eqs. (77)-(80), now I1_ s , I 2 _ s , I 3_ s and I 4 _ s can be regularized into

  I   N   Jbd   ExI d    ExI d    I   N   Jbd   ExI d    ExI d    I   N   Jbd   ExI d    ExI d    I   N   Jbd   ExI d    ExI d 1

I1_ s 

1 1

I 2_ s

1 2

I 3_ s

I 4_ s

1

k

1

1

k

1

1 3

1

1

1

1_ s

1

2_ s

1

3_ s

1

k

1

1

2_ s

1

k

1

1 4

1

1

1_ s

1

3_ s

1

4_ s

1

4_ s

(85) (86) (87) (88)

From Eqs. (85)-(88), the first terms of these four equations which are regularized can be directly computed by the traditional Gauss integration method. The second terms can be computed analytically through Eqs. (81)-(84). In our method, five types of elements are employed as illustrated in Figure 1.  is the position parameter of the collation nodes in the intrinsic coordinates system. The range of its values varies from 0.125 to 0.5 in our method. In numerical examples, we will compare results obtained by using different value of  . The comparisons can be found in the next section.

Figure 1. Five types of elements (a) Constant element; (b) Linear continuous element; (c) Linear semi-discontinuous element; (d) Quadratic continuous element; (e) Quadratic semi-discontinuous element 4. Numerical examples Two numerical examples considering a rectangular shaped structure and a plate with a hole subject to an impact loading are presented in this section. Comparisons with available results are made to demonstrate the validity and accuracy of the proposed method. 4.1 Comparison with the existing results As shown in Figure 2, a rod with a traction load of the Heaviside type at the right end is considered in this numerical example. The geometry and boundary conditions are

also illustrated in Figure 2. a and b are the length and the width of the rod, respectively. Firstly, the results which are obtained by using the present method are compared with these in Li and Zhang’s work [31]. In their paper, a = 8m and b = 2m. The Young’s modulus E = 1.1e5N/m2, the Poisson’s ratio v  0 and the mass density ρ = 2.0kg/m3. The load parameter P(t) = 1000N. In the analysis, we employ three different types of element including constant element, linear element and quadratic element. Furthermore, the influence of the position parameter of  is also studied. The sample point is located in the middle of the right edge. In the application of constant element, totally 22 constant elements including 3 elements in the left edge, 3 elements in the right edge, 8 elements in the top and 8 elements in the bottom are employed. In the application of linear element, 20 linear elements including 2 in the left edge, 2 in the right edge, 8 in the top and 8 in the bottom are employed. It is worth noting that linear semi-discontinuous elements are employed in the neighborhood of the corners to avoid the corner problem as mentioned in [27, 29]. Thus, the position parameter  should be considered. In the application of quadratic element, the same mesh as that in the application of linear element is applied for discretization. The time scan size that we adopted in all applications is t  0.00356 . The displacements on the sample point are computed by using three types of elements. The variation history of the computed displacements is compared with the analytical solution and that in Li and Zhang’s paper. The comparison is illustrated in Fig. 3. It can be seen from Fig. 3 that the variation history that is obtained by the present method with different elements coincides nicely with the analytical solutions. The error of the results which are obtained by the present method is much lower than that in the referred work. In three applications of different elements, the best result is obtained by using quadratic element. And the constant element gives results with a lower accuracy. It is interesting to find in Fig. 3 that the error of variation frequency seems to be much less than that in the referred work even in the application of constant element. This is very important in some engineering cases since much information is involved in the frequency. However, the amplitude attenuation of the variation history in the present method seems to be larger than that in the referred work. In applications of linear element and quadratic element, three evaluations of the position parameter  are considered and compared in Fig. 4 and Fig. 5. It should be

noted that the parameter  is given to determine the location of the source point in an element, specially for the element which is located at a corner, which is used to keep the source point in the element and avoid the corner problem in BEM. We have studied the influence of the parameter  through numerical experiments and found that the values of  have little influence on the final results. From Fig. 4 and Fig. 5, the numerical result is not sensitive to the position parameter  which describes the location of the source point. Thus, no optimizing method has been developed for evaluating this parameter.

Figure 2. Geometry and boundary conditions for a rod under Heaviside type traction

Figure 3. Displacements u1 at the middle of right end compared with Li and Zhang's method

Figure 4. Displacements u1 at the middle of right end using linear elements with different  values

Figure 5. Displacements u1 at the middle of right end using quadratic elements with different  values Note that the example of Fig. 2 is a classical benchmark [25, 26]. It is known to have severe instability problems. Take using the linear elements as an example, smaller DTs are also employed with the same number of mesh, in which t  0.00356 , t  0.00178 and t  0.00089 , respectively. And more cycles are plotted. The time cost of the presented method is 71841ms for 300 steps when t  0.00356 , and the time cost of the presented method is 152270ms for 600 steps when t  0.00178 , while the time cost of the presented method is 150020ms for 600 steps when t  0.00089 . It can be found that from Fig. 6 when t  0.00089 , the instability arises in the second cycles,

while for the two cases t  0.00356 and t  0.00178 , the results agree well with the analytical solutions, in which   0.8 and   0.4 . The results when t  0.00356 show a higher accuracy than those when t  0.00178 . In the future work, we will study the method for improving stability.

Figure 6. Displacements u1 at the middle of right end with different time scan sizes Then, we compare with the results of method of Lei, Li and Qin et al [23]. Assume

a  2m and b  1m . The Young’s modulus E  2.11011 N m2 , the Poisson’s ratio v  0 and the mass density   7900kg / m3 . p=100MN. Using about the same number of elements as their method, also different kinds of elements with different position parameters are considered. In their method,   c1t / l =1.0317. In Fig. 7 and Fig. 8, it can be found that the results also agree well with the analytical solutions and these of Lei, Li and Qin et al. The time cost of the presented method is 229086ms for 300 steps when using constant elements. The time cost of the presented method is 443777ms for 300 steps when using linear elements, while the time cost of the presented method is 478734ms for 300 steps when using quadratic elements with the same nodes with these using linear elements. The accuracy and efficiency of our method are also can be found in Fig. 8.

Figure 7. Displacements u1 at the middle of right end compared with method of Lei, Li and Qin et al

Figure 8. Displacements Eu1 /  p1a  at the middle of right end compared with method of Lei, Li and Qin et al 4.2 A plate with a hole which is subject to an impact loading The performance of the proposed method is further demonstrated through a more realistic problem. A plate with a hole is subject to a force as illustrated in Fig. 9. The material properties of the plate are: E = 69Gpa, ρ = 2700kg/m3, ν = 0.3 and P(t) = 1000N, respectively.

Figure 9. A plate with a hole (unit: m) The mechanical response of the plate is simulated using both the proposed method and the FEM. In both cases, the impact force is assumed to be applied uniformly on the top of the plate and the bottom of the plate is fixed. The time scan size is set to be 0.0001s. The results obtained by FEM with 12645 linear quadrilateral elements and 12776 nodes are denoted as reference solutions. Fig. 10 shows the displacement response in the load direction at point A. Fig. 11 and Fig. 12 show the displacement responses in the x and y directions at point B. These numerical results are obtained by proposed method with 100 nodes, in which 17 nodes and 8 quadratic elements are using in each straight edge, while 32 nodes and 16 quadratic elements are using in circular edge. The time cost of the proposed method is 758859ms for 300 steps using our method. From Fig.10 to Fig.12, the numerical results of the proposed method agree very well with the FEM results.

Figure 10. The displacement response in the load direction at point A

Figure 11. The displacement response in the x direction at point B

Figure 12 The displacement response in the y direction at point B 5. Conclusions This paper has developed a direct singularity cancellation method for time-domain boundary element formulation for elastodynamic problems. The singularities in the integrals of the time-dependent fundamental solutions are resolved by the Taylor expansion. Thus, the integrals with the similar singularities which can be easily treated are obtained. Using the integrals with simple forms, the original integrals with singularities can be separated into a regular part and a singular part. The singular part can be treated with standard Gaussian Integration, while the singular part can be obtained by analytical integrals. In our numerical implementation, the singular integration elements are subdivided into several parts to assume the consistency of the integrands. Numerical examples are proposed to compared with analytical solutions and some pioneering works. The numerical results of our method are in good agreement with these published results. Acknowledgments This work was supported by the National Natural Science Foundation of China (11602229 and 11602082) and the Hunan Provincial Natural Science Foundation of China (2017JJ3061).

References [1] Eringer A C. Elastodynamics[M]. Рипол Классик, 1974.

[2] Zienkiewicz O C, Taylor R L, Zienkiewicz O C, et al. The finite element method[M]. London: McGraw-hill, 1977. [3] Feng S Z, Han X. A novel multi-grid based reanalysis approach for efficient prediction of fatigue crack propagation[J]. Computer Methods in Applied Mechanics and Engineering, 2019, 353: 107-122. [4] Lv J H, Jiao Y Y, Wriggers P, et al. Efficient integration of crack singularities in the extended finite element method: Duffy-distance transformation and conformal preconditioning strategy[J]. Computer Methods in Applied Mechanics and Engineering, 2018, 340: 559-576. [5] Mansur W J. A time-stepping technique to solve wave propagation problems using the boundary element method[D]. University of Southampton, 1983. [6] Banerjee P K, Ahmad S, Manolis G D. Transient elastodynamic analysis of three‐dimensional problems by boundary element method[J]. Earthquake engineering & structural dynamics, 1986, 14(6): 933-949. [7] Dominguez J. Boundary elements in dynamics[M]. Wit Press, 1993. [8] Beskos D E. Boundary element methods in dynamic analysis[J]. Applied Mechanics Reviews, 1987, 40(1): 1-23. [9] Wen P H, Aliabadi M H, Rooke D P. Cracks in three dimensions: A dynamic dual boundary element analysis[J]. Computer methods in applied mechanics and engineering, 1998, 167(1-2): 139-151. [10] Carrer J A M, Mansur W J. Stress and velocity in 2D transient elastodynamic analysis by the boundary element method[J]. Engineering analysis with boundary elements, 1999, 23(3): 233245. [11] Liu Y H, Chen S S, Li J, et al. A meshless local natural neighbour interpolation method applied to structural dynamic analysis[J]. CMES: Computer Modeling in Engineering & Sciences, 2008, 31(3): 145-156. [12]Chen S S, Wang W, Zhao X. An interpolating element-free Galerkin scaled boundary method applied to structural dynamic analysis[J]. Applied Mathematical Modelling, 2019. [13] Carrer J A M, Pereira W L A, Mansur W J. Two-dimensional elastodynamics by the timedomain boundary element method: Lagrange interpolation strategy in time integration[J]. Engineering Analysis with Boundary Elements, 2012, 36(7): 1164-1172. [14] Hatzigeorgiou G D, Beskos D E. Transient dynamic response of 3-D elastoplastic structures by the D/BEM[J]. WIT Transactions on Modelling and Simulation, 2001, 28. [15] Manolis GD. A comparative study on three boundary element method approaches to problems in elastodynamics. International Journal for Numerical Methods in Engineering, 1983, 19:73–91. [16] De Lacerda LA, Wrobel LC, Mansur WJ. A boundary integral formulation for twodimensional acoustic radiation in a subsonic uniform flow. J Acoustical Soc Am 1996, 100:98– 107.

[17] Lubich C. Convolution quadrature and discretized operational calculus I. Numer Math 1988, 52:129–45. [18] Lubich C. Convolution quadrature and discretized operational calculus II. Numer Math 1988, 52:413–25. [19] Niwa Y. An application of the integral equation method to two-dimensional elastodynamics[J]. Theoretical and applied mechanics, 1980, 28: 281-290. [20] Mansur W J, Brebbia C A. Formulation of the boundary element method for transient problems governed by the scalar wave equation[J]. Applied Mathematical Modelling, 1982, 6(4): 307-311. [21] Israil A S M, Banerjee P K. Advanced time‐domain formulation of BEM for two‐ dimensional transient elastodynamics[J]. International Journal for Numerical Methods in Engineering, 1990, 29(7): 1421-1440. [22] Lei W, Ji D, Li H, et al. On an analytical method to solve singular integrals both in space and time for 2-D elastodynamics by TD-BEM[J]. Applied Mathematical Modelling, 2015, 39(20): 6307-6318. [23] Lei W, Li H, Qin X, et al. Dynamics-based analytical solutions to singular integrals for elastodynamics by time domain boundary element method[J]. Applied Mathematical Modelling, 2018, 56: 612-625. [24]Cole D M, Kosloff D D, Minster J B. A numerical boundary integral equation method for elastodynamics. I[J]. Bulletin of the Seismological Society of America, 1978, 68(5): 1331-1357. [25]Peirce A, Siebrits E. Stability analysis and design of time‐stepping schemes for general elastodynamic boundary element models[J]. International Journal for Numerical Methods in Engineering, 1997, 40(2): 319-342. [26]Frangi A, Novati G. On the numerical stability of time-domain elastodynamic analyses by BEM[J]. Computer Methods in Applied Mechanics and Engineering, 1999, 173(3-4): 403-417. [27] Aliabadi M H, Hall W S, Phemister T G. Taylor expansions for singular kernels in the boundary element method[J]. International Journal for Numerical Methods in Engineering, 1985, 21(12): 2221-2236. [28] Guiggiani M. Formulation and numerical treatment of boundary integral equations with hypersingular kernels[J]. Singular integrals in boundary element methods, 1998: 85-124. [29] Xie G Z, Zhang J, Huang C, et al. A direct traction boundary integral equation method for three-dimension crack problems in infinite and finite domains[J]. Computational Mechanics, 2014, 53(4): 575-586. [30] Xie G, Zhang D, Meng F, et al. Calculation of stress intensity factor along the 3D crack front by dual BIE with new crack front elements[J]. Acta Mechanica, 2017, 228(9): 3135-3153. [31] Li Y, Zhang J. A comparative study of time domain BEM for 3D elastodynamic analysis[J]. WIT Transactions on Modelling and Simulation, 2013, 56: 503-514.

Appendix A The time integrations for I1 may be piecewise functions, which are as follows 1. Sc1tm1  r , Sc 2tm1  r



 1  1  lk  2 2 2 2 2  2  2  S1 S1  r  r ln S1  S1  r c 2 r 1  1  I1   2  1  1  lk rr   S S 2  r 2  ,l 2,k S2 S22  r 2  c 2  2 r 2 2 2 r  2



  rrr  S ,l , k 2



1

  



 S12  r 2    tm 1   tm  

(16)

2. Sc1tm1  r , Sc 2tm  r & &Sc 2tm1  r I1    



1  1 1  lk ln Sc1tm 1   2  c12  2

1  1 1  lk ln( Sc 2tm   2  c22  2 1 4 c

2 2

 lk ln(r ) 

 Sc1tm1 

 Sc 2tm 

2



 Sc1tm 

2



  r 2  

(17)

  r 2 )  

1 1  1  lk  2r,l r,k  Sc1tm1  2  c12  2r 2

1 1  1  lk  2r,l r,k  Sc1tm  2  c12  2r 2



1  1 1   r2    2  lk ln Sc1tm    2  c1  2

2

 Sc1tm 

2

 Sc1tm1 

2

  r2  



 1  1  1  r 2     2r,l r,k  Sc 2tm  2 2  lk   2  c2  2r

 S c 2t m 



  r 2  

2

3. Sc1tm1  r , Sc 2tm  r I1  



1  1 1  lk ln Sc1tm1   2  c12  2

 Sc1tm1 

1 1  1  lk  2r,l r,k  Sc1tm1  2  c12  2r 2

2





1  1 1   r2    2  lk ln Sc1tm    2  c1  2

 Sc1tm1 

2

 Sc1tm 

 1  1  1  r2     2r,l r,k  Sc1tm  2 2  lk   2  c1  2r

2



  r 2  

 Sc1tm 

2

(18)

  r 2  

4. Sc1tm  r & &Sc1tm1  r Sc 2tm  r & &Sc 2tm1  r I1  

1 1  lk ln(r )   lk ln(r ) 4 c12 4 c22

1 1   2  c12 1 1   2  c12



1  2  lk ln Sc1tm 

 Sc1tm 

 1  2r 2  lk  2r,l r,k  Sc1tm

2



 1  1  r    2   2  c2 2

 Sc1tm 

2

1  2  lk ln( Sc 2tm 

 S c 2t m 



2

 1  1  1  r 2     2r,l r,k  Sc 2tm  2 2  lk   2  c2  2r

 S c 2t m 

5. Sc1tm  r & &Sc1tm1  r Sc 2tm  0 I1   

1 4 c

1 1  2  c12

2 1

 lk ln(r ) 



1  1 1  lk ln Sc1tm   2  c12  2

 1  2r 2  lk  2r,l r,k  Sc1tm

 Sc1tm 

2

 Sc1tm 

2

  r2  

6. Sc1tm  r Sc 2tm  r

I1  0



  r2  

(21)

(19)

  r )   2

(20)

2



  r2  

Appendix B The time integrations for I2 are also piecewise functions, which are as, 1. Sc1tm1  r , Sc 2tm1  r  2  2  t  mt    r  t  mt    1 2  r  t  mt   S1   r 2      2  c1t  c1t   S1   r 2  c1t  1  I2  Dlk  2 c1    3  r2  S1 1    r     1  S  S 2  r 2  2   2   2  1 1 2 2   2 2  c1 t  S1   r  c1 t   S1   r   c1 t      





 2  2  t  mt    r  t  mt    1 2  r  t  mt   S2   r 2      2  c2 t  c2 t    S2   r 2  c2 t  1   Dlk  2 c2    3  r2  S2 1    r     1  S  S 2  r 2  2     2  2 2 2 2 2  c2 t   S2   r 2   c2 t    S 2   r 2   c2 t      





 t  mt Blk  S1 2 c1  c1t 1

 t  mt 1  Blk  S2 2 c2 ct 

 S1 

2

 S2 

   crt   S   r  tt  t  r  2  S   r     S   r  3c t  c  t    t

 r2  2

 r2

2 3c 2 t





 S1 

2

2

3

 r2

2

  2   S1   r 2  tm 1   tm     S2  2   S2   r 2  tm 1   tm    S1

m 1

2

1

2

m

2

3

2

2

2

2

2 1

    t  mt     1  t  mt        2  c2 t   2  c  t 2  r 1    S2   r    Alk  2 c2  1    S2 2  2    ln S2   S 2   r 2  c2 t    S 2 2  r 2  



2

m 1

2

2 1

 S2  2  S2   r 2        r    c22 t    



m

    tm 1    tm 1  2  S2   r 2  

(22)

2. Sc1tm1  r , Sc 2tm  r & &Sc 2tm1  r  2  2  t  mt    r  t  mt    1 2  r  t  mt   S1   r 2      2  c1t  c1t   S1   r 2  c1t  1  I2  Dlk  2 c1     3  r2 S1 1    r     1  S  S 2  r 2  2   2   2  1 1  c1 t   S1 2  r 2   c1 t    S1 2  r 2   c1 t      





  2   S1   r 2  tm 1   tm    S1

 2  2  t  mt    r  t  mt    Sc 2tm 1 2 2  r  t  mt     Sc 2tm   r    2  2 c2 t c2 t c2 t    S t 2  r 2    S t  r   c2 m c2 m 1    Dlk  2 c2    3  r2  S c 2tm 1    r     1  S t  S t 2  r 2  2   2   2  c2 m c2 m  c2 t   Sc 2tm 2  r 2   c2 t    Sc 2tm 2  r 2   c2 t      









3  t  mt t  r2  1 2 2 2 2 Blk  S1  S1   r 2  2  S1   r 2   2   S1   r 2  m1 2 c1 3c1 t c  t  c1t  tm  1  2 3   r  1 2 2 2 2  t  mt  Blk  S c 2t m  S c 2t m   r 2  2  Sc 2tm   r 2   2   Sc 2tm   r 2  2 c2 3c2 t  c2 t   c2 t 









1 Alk ln  r  2 c23t

     Sc 2tm 1  t  mt      t  mt      2 2  c2 t   2  2  c2 t  r  Sc 2tm   r  1    Sc 2tm   r     Alk  2 c2    1     S c 2tm r  2  2    ln Sc 2tm   Sc 2tm   r 2     2   2   2  c2 t   Sc 2tm   r  c2 t      



3. Sc1tm1  r , Sc 2tm  r



      1  2  Sc 2tm   r 2  

        

(23)

 2  2  t  mt    r  t  mt    1 2  r  t  mt   S1   r 2      2  c1t c1t    S1   r 2  c1t  1  I2  Dlk  2 c1    3  r2  S1 1    r     1  S  S 2  r 2  2     2  1 1 2  c1 t   S1 2  r 2   c1 t    S1 2  r 2   c1 t      





 t  mt  Blk  S1 2 c1   c1t 1

 S1 

2

 r2 

2 3c12 t



 S1 

2

 r2



   crt  2

3

2 1

 S1 

2

  2   S1   r 2  tm 1   tm    S1

(24)

t   r 2  m 1   tm

4. Sc1tm  r & &Sc1tm1  r Sc 2tm  r & &Sc 2tm1  r  2  2  t  mt    r  t  mt    Sc1tm 1 2  r  t  mt   Sc1tm   r 2      2 2  c1t  c1t    Sc1tm   r 2  Sc1tm   r 2  c1t  1  I2   Dlk  2 c1    3  r2  Sc1tm 1    r     1  S t  S t 2  r 2  2   2   2  c1 m c1 m 2 2   2 2  c1 t  Sc1tm   r  c1 t   Sc1tm   r   c1 t      





 2  2  t  mt    r  t  mt    Sc 2tm 1 2 2  r  t  mt     S c 2tm   r    2  2 c2 t c2 t c2 t    S t 2  r 2    S t  r   c2 m c2 m 1    Dlk  2 c2    3  r2      Sc 2tm 1  r    1 S t  S t 2  r 2  2   2   2  c2 m c2 m  c2 t   Sc 2tm 2  r 2   c2 t    Sc 2tm 2  r 2   c2 t      





 t  mt 1 Blk  Sc1tm 2 c1  c1t

 t  mt 1 Blk  Sc 2tm 2 c2  c2 t 1  Alk ln  r  2 c23 t 

 Sc1tm 

2

 S c 2tm 

  S t   r    crt   S t   r    r  2  S t   r     S t   r  3c t  c  t   

 r2  2

 r2



2 3c12 t

2

2

2

3

2

c2 m

2 2

     Sc 2tm 1  t  mt      t  mt      2 2  c2 t   2  2  c  t  2  S t  r r S t  r     c 2 m c 2 m 1       Alk  2 c2    1     S t r  2 c2 m  2    ln Sc 2tm   Sc 2tm   r 2     2     c2 t    c2 t    Sc 2tm 2  r 2    



2

2

2

c2 m

        

c1 m

2 1

2

2 2

2

3

2

c1 m

        



(25)

      1  2  Sc 2tm   r 2  

5. Sc1tm  r & &Sc1tm1  r Sc 2tm  0  2  2  t  mt    r  t  mt    Sc1tm 1 2  r  t  mt   Sc1tm   r 2      2 2  2 2 c1t c  t c  t 1 1     S t  r S  c1 m  1    c1tm   r I2   Dlk  2 c1    3  r2  Sc1tm 1    r     1  S t  S t 2  r 2  2     2  c1 m c1 m 2 2 2    2 2 c  t c  t  1  1   Sc1tm   r   c1 t    Sc1tm   r    





 t  mt  Blk  Sc1tm 2 c1  c1t  1

 Sc1tm 

2

 r2 

2 3c12 t



 Sc1tm 

6. Sc1tm  r Sc 2tm  r

I2  0 Appendix C The time integrations for I3 are as, 1. Sc1tm1  r , Sc 2tm1  r

(27)

2

 r2



   crt  2

3

2 1

 Sc1tm 

2

   r2   

        

(26)

 2   r   m  1 t  t       c12 t    1  I3  Dlk  2 c1    2 S1   r  2   c12 t  2   S  1  r     2   r   m  1 t  t          c2 t    1   Dlk  2 c2    2 S2   r  2   c22 t  2   S  2  r        m  1 t  t  1 Blk   S1 2 c1 c1t       m  1 t  t  1  Blk   S2 2 c2 c2 t   

  2   m  1 t  t     2 c1t   S1   r 2    3   r     2  c t    1   

   1     2  S1 2  S1   r 2   c1 t 

  2   m  1 t  t     2 c2 t   S2   r 2    3   r     2  c t    2   

S 

2

1

S 

2

S 

2

2

   1     2  S2 2  S2   r 2   c2 t 

S 

S 

  2   r2    2    3c2 t 

S 

2

2

   r2   

2

2

S 

2   r   r2    2    c2 t 

S 

2

1

3

S2

   r2   

2   r   r2    2    c1 t 

2

2

    m  1 t  t  S2 1    m  1 t  t        2 2  c  t c  t 2 2 2 2    S2   r r  S2   r    1   Alk  2 c2     2  S2 1    r    1    ln  S2   S 2   r 2     2  2   c22 t   2     c2 t  S 2  r 2  r S  r  2     2 

      2  S1   r 2  tm1   tm    S1

      2  S2   r 2  tm1   tm   

  r   m  1 t  t    r2    c2 t   

3

2

1

2

1

S 

1

  2   r2    2    3c1 t 

  r   m  1 t  t    r2    c1t   

2

1

1

1

2

S 

1

 t   r 2  m 1  tm  t   r 2  m 1  tm

   t  m 1    tm    

(28)

2. Sc1tm1  r , Sc 2tm  r & &Sc 2tm1  r  2   r   m  1 t  t       c12 t    1  I3  Dlk  2 c1    2 S1   r  2   c12 t    S  r2     1

  2   m  1 t  t     2 c1t   S1   r 2  1

  3   r     c 2 t     1   

1

S  1

2

S  1

   1     2  S1  r 2   c1 t  

2

      2  S1   r 2  tm1   tm      Sc 2tm   2 S  c 2tm   r 2       

  r   m  1 t  t    r2    c1t   

S  1

2

S1

   r2   

 2   2   m  1 t  t  r   m  1 t  t   1 2  r   m  1 t  t      S c 2tm   r 2   2 2   S t  r   c2 t c2 t c2 t   c2 m    1  Dlk  2 c2       r2  3 S c 2tm 1    r     1  S t  S t 2  r 2    2   2   2  c2 m c2 m   c2 t   S t 2  r 2   c2 t    S t 2  r 2   c2 t   c 2 m c 2 m      





3 2     m  1 t  t   t 2 2 2 1   2     r  Blk   S1  S1   r 2    2   S1   r 2    2   S1   r 2  m 1 2 c1 c  t 3 c  t c  t   1  1  tm   1    3   r2  1 2 2 2     m  1 t  t    2  Blk   S c 2tm  S c 2t m   r 2    2  Sc 2tm   r 2   2   Sc 2tm   r 2   2 c2 c2 t  c2 t       3c2 t 1  ln  r  Alk 2 c23 t





(29)



 m  1 t  t       m  1 t  t  Sc 2tm 1       2 2    2 2 c2 t c2 t   S t  r r S t  r       c2 m c2 m 1   Alk  2 c2    1    r   S c 2tm 1 2   2    ln Sc 2tm   Sc 2tm   r 2    2    c2 t  S t 2  r 2   c2 t   r  Sc 2tm 2  r 2  c2 m     





         

3. Sc1tm1  r , Sc 2tm  r  2   r   m  1 t  t      2  c1 t    1  I3  Dlk  2 c1    2 S1   r  2   c12 t  2   S  1  r    

    m  1 t  t  Blk   S1 2 c1 c1t   1

  2   m  1 t  t     2 c1t 2   S1   r  1

  3   r     c 2 t     1   

S  1

2

S  1

   1     2  S1 2 2  S1   r   c1 t  1

  2   r2    2    3c1 t 

S  1

2

2

  r   m  1 t  t    r2    c1t   

S  1

2   r   r2    2    c1 t  3

2

   r2   

S  1

2

 t   r 2  m 1  tm

      2  S1   r 2  tm1   tm    S1

(30)

4. Sc1tm  r & &Sc1tm1  r Sc 2tm  r & &Sc 2tm1  r  2   2  m  1 t  t  S t 2  r 2  r   m  1 t  t   1  r   m  1 t  t       c1 m  2 2      2 c1 t c1t c1t    Sc1tm   r    1 I3   Dlk  2 c1       r2  3     S t r 1 1 2   2 c1 m      2   2   2  Sc1tm  Sc1tm   r    c1 t   S t 2  r 2   c1 t    S t 2  r 2   c1 t   c1 m c1 m       





 2   2  m  1 t  t  S t 2  r 2  r   m  1 t  t   1  r   m  1 t  t       c2 m  2    2  c  t c  t c2 t 2 2  S t  r    c2 m     1  Dlk  2 c2       r2  3 S c 2tm 1    r     1  S t  S t 2  r 2    2   2   2  c2 m c2 m   c2 t   S t 2  r 2   c2 t    S t 2  r 2   c2 t   c2 m c2 m      



1     m  1 t  t  Blk   Sc1tm 2 c1 c1t   1     m  1 t  t   Blk   S c 2tm 2 c2 c2 t   



 Sc1tm 

2

 S c 2t m 

  2  r2    2   3c1 t 2



 Sc1tm 



 S c 2t m 

  2  r2    2   3c2 t

2



 r2

2

   crt 

 r2

2

3

2 1

   crt  2

3

2 2

 Sc1tm 

2

2

   r 2  

2

 S c 2t m 

     Sc1tm   r          S c 2tm   2  Sc 2tm   r 2        Sc1tm

2

(31)

   r 2  

1 ln  r  Alk 2 c23 t

 m  1 t  t       m  1 t  t  S c 2tm 1       2 2   2 2 c2 t c2 t  S t  r r S t  r       c2 m c2 m 1   Alk  2 c2      1   r   S c 2tm 1 2   2    ln Sc 2tm   Sc 2tm   r 2    2    c2 t  S t 2  r 2   c2 t   r  Sc 2tm 2  r 2  c2 m     





         

5. Sc1tm  r & &Sc1tm1  r Sc 2tm  0  2   2  m  1 t  t  S t 2  r 2  r   m  1 t  t   1  r   m  1 t  t       c1 m  2 2     2  c  t c  t c1t 1 1    Sc1tm   r    1 I3   Dlk  2 c1       r 2  3     S t 1 c1 m  r    1 S t  S t 2  r 2    2   2   2  c1 m c1 m   c1 t   S t 2  r 2   c1 t    S t 2  r 2   c1 t  c1 m c1 m        2 3    m  1  t  t     r   S t S t 2  r2    2 1 2 2     Blk    Sc1tm   r 2   2   Sc1tm   r 2    2  c1 m  c1 m  2 c1 c  t 3 c  t c  t   1  1    1  







Sc1tm

 Sc1tm 

2

      r       

(32)

2



6. Sc1tm  r Sc 2tm  r

I3  0

(33)

It should be noted that in appendix B and C the concept of Cauchy principal value and Hadamard finite part integral are employed in this process. The divergent terms related to S1  r  0 or S2  r  0 are removed directly, thus the wave front singularities disappear. Appendix D The time integrations for I4 are as, 1. Sc1tm1  r , Sc 2tm1  r   1   1  1 1 Dlk   S1 S12  r 2   Alk   ln( S2  S22  r 2 )    2 2  2 c1 t  c1  2 c2 t  c2   tm1 I4     1 1   tm 2 2  D  S S  r   lk 2 2 2  2 c t  2  c2   

2. Sc1tm1  r , Sc 2tm  r & &Sc 2tm1  r

(34)

 1   t 1 1  I4   Dlk   S1 S12  r 2   m1  Alk ln(r ) 2 t 2  c  t c 2  c23t 1  1   m  1 1 2 2  Alk ln( Sc 2tm   Sc 2tm   r 2 )  Dlk Sc 2tm  Sc 2tm   r 2 2 c23t 2 c23 t



(35)



3. Sc1tm1  r , Sc 2tm  r  1   t  1 I4   Dlk   S1 S12  r 2   m1 2  2 c1 t  c1   tm

(36)

4. Sc1tm  r & &Sc1tm1  r Sc 2tm  r & &Sc 2tm1  r I4 



1 Dlk Sc1tm 2 c13t

 Sc1tm 



1  Alk ln( Sc 2tm  2 c23t

2



 r2 

 Sc 2tm 

2

1 Alk ln(r ) 2 c23t



1  r2 )  Dlk Sc 2tm 2 c23 t

(37)

 S c 2t m 

2

 r2

5. Sc1tm  r & &Sc1tm1  r Sc 2tm  0 I4 



1 Dlk Sc1tm 2 c13t

 Sc1tm 

2

 r2



(38)

6. Sc1tm  r Sc 2tm  r

I4  0

(39)

It should be noted that in Appendix A and D, the concept of Cauchy principal value are employed in this process. The divergent terms related to S1  r  0 or S2  r  0 is removed directly, thus the wave front singularities disappear.

Appendix E The singularities existing in the time integrations I1 can be obtained by using approximate expansion, which are also piecewise functions. The expressions are as, 1. Sc1tm1  r , Sc 2tm1  r I1_ st 

1 1   lk   1 2 1 2  tm1   r,l r,k   2 S1  2 S2  t  0 2 r 2  2 c2  m   c1

(46)

2. Sc1tm1  r , Sc 2tm  r & &Sc 2tm1  r I1_ st  

1 4 c

2 2

 lk ln(r ) 

1 1  1 2   2  2  lk  2r,l r,k   Sc1tm1    2  c1  2r 

3. Sc1tm1  r , Sc 2tm  r

(47)

I1_ st 

1 1  2  c12

1 1 2   1  2r 2  lk  2r,l r,k   Sc1tm1     2  c 2   1

2   1  2r 2  lk  2r,l r,k   Sc1tm    

4. Sc1tm  r & &Sc1tm1  r Sc 2tm  r & &Sc 2tm1  r I1_ st  

1 4 c

2 1

 lk ln(r ) 

1 4 c22

(49)

 lk ln(r )

5. Sc1tm  r & &Sc1tm1  r Sc 2tm  0 I1_ st  

1 4 c

2 1

 lk ln(r ) 

1 1  1 2   2  2  lk  2r,l r,k   Sc1tm    2  c1  2r 

6. Sc1tm  r Sc 2tm  r

I1  0

(51)

(50)

(48)

Appendix F The singularities existing in the time integrations I2 can be obtained by using approximate expansion, which are also piecewise functions. The expressions are as, 1. Sc1tm1  r , Sc 2tm1  r I 2 _ st  

t r  t  mt   1  2  2  t  mt   Dlk    2   S1   m 1  S1  2 c1 c1t c1t  c1 t     tm 1

 2  t  mt    r  t  mt    1  2  tm 1 Dlk   S2      2   S2   2 c2  c2 t  tm   c2 t   c2 t  1

2  3 3S1r 2  r 2  t  mt  2 r 2   tm 1  Blk    2 S1   S1 -   2  S1 2 c1 c  t 2 3 c  t 2 c  t     1   tm

(52)

1



2  3 3S 2 r 2  r 2  t  mt  2 r 2   tm 1 Blk    2 S2   S2 -   2  S2 2 c2 ct  2  3c1 t  2  c1 t  tm 



  t  mt   tm 1 Alk   2 c2  rc2 t  tm

1

1

2. Sc1tm1  r , Sc 2tm  r & &Sc 2tm1  r I 2 _ st  

t r  t  mt   1  2  2  t  mt   Dlk    2   S1   m 1  S1  2 c1 c1t c1t   tm  c1 t   1

 2  t  mt    r  t  mt    1  2 2  Dlk    S c 2t m       2   S c 2t m   2 c2 c  t c  t c  t 2 2     2    1

2  3 3S1r 2  r 2  t  mt  2 r 2   tm 1  Blk   S1    2  S1    2 S1  2 c1 2S1  3c t  2  c t  tm  c1t 

(53)

1



 t  mt   3Sc 2tm r 2   r 2  r2  2  2 3 Blk  S t   S t       c2 m   c2 m    2  S c 2t m  2 2 c2 2Sc 2tm  3c2 t  2   c2 t   c2 t  



 t  mt  1 1 1 Alk ln  r   Alk   2 c23 t 2 c2  c2 t  r

1

3. Sc1tm1  r , Sc 2tm  r  2  t  mt   t r  t  mt   1  2  Dlk    2   S1   m 1  S1  2 c1 c1t c1t    c1 t     tm  t  mt   tm1 1 r2  2  3 3S1r 2  r 2    Blk   S1    2  S1    2 S1  2 c1 2S1  3c t  2  c t    c1t   tm I 2 _ st 

1

(54)

4. Sc1tm  r & &Sc1tm1  r Sc 2tm  r & &Sc 2tm1  r  r  t  mt   1  2  2  t  mt   Dlk    2   Sc1tm    Sc1tm  2 c1 c1t c1t    c1 t    2  t  mt    r  t  mt    1  1 2 2   Dlk    S c 2t m       2   S c 2t m   2 c2  c2 t    c2 t   c2 t  2 2  t  mt   3Sc1tm r  r 2 1 r  2  2 3 + Blk    Sc1tm     2   Sc1tm     2 Sc1tm  2 c1 2 Sc1tm  3c t  2  c1t    c t 2 2 2  t  mt   3Sc 2tm r   r  1 r  2  2 3  Blk    S c 2t m     2   S c 2t m      2  S c 2t m  2 c2 2Sc 2tm  3c2 t  2  c2 t     c2 t  1

I 2 _ st  



(55)

 t  mt  1 1 1 Alk ln  r   Alk   3 2 c2 t 2 c2  c2 t  r

5. Sc1tm  r & &Sc1tm1  r Sc 2tm  0 I 2 _ st   +

 2  t  mt    r  t  mt   1  2  Dlk    2   Sc1tm    Sc1tm  2 c1 c1t c1t    c1 t     1

(56)

 t  mt   3Sc1tm r  r r  2  2 3   Blk    Sc1tm     2   Sc1tm     2 Sc1tm  2 c1 2Sc1tm  3c1 t  2    c1 t  c1t   2

2

1

2

6. Sc1tm  r Sc 2tm  r I 2 _ st  0

(57)

Appendix G The singularities existing in the time integrations I3 can be obtained by using approximate expansion, which are also piecewise functions. The expressions are as, 1. Sc1tm1  r , Sc 2tm1  r I 3_ st   

 2   m  1 t  t  r   m  1 t  t   1  2 t Dlk  S1  +  2   S1   m 1 2 c1 c  t c  t c  t 1 1  1    tm 1

 2   m  1 t  t  r   m  1 t  t   1  2 t Dlk  S2  +  2   S 2   m 1 2 c2 c  t c  t c  t 2 2  2    tm 1

   m  1 t  t   2 r 2  2  3 3S1r 2   tm 1 Blk    S1    2  S1  2 c1 c1t 2  3c1 t  2   tm   1

   m  1 t  t   2 r 2   tm 1 2  3 3S 2 r 2    Blk     S2    2  S2  2 c2 c  t 2 3 c  t 2      tm 2 2     m  1 t  t   tm 1 1  Alk   2 c2 rc2 t   tm 

1

2. Sc1tm1  r , Sc 2tm  r & &Sc 2tm1  r

(58)

I 3_ st   

 2   m  1 t  t  r   m  1 t  t   1  2 t Dlk  S1  +  2   S1   m 1 2 c1 c1t c1t  c1 t    tm 1

 2   m  1 t  t  r   m  1 t  t   1  2  Dlk  Sc 2tm  +  2   Sc 2tm   2 c2 c  t c  t c  t 2 2  2    1

   m  1 t  t   2 r 2  2  3 3S1r 2   tm 1 Blk    S1    2  S1  2 c1 c1t 2  3c1 t  2   tm   1

(59)

   m  1 t  t   2 3 3Sc 2tm r 2   r2  2  Blk     S c 2tm     2   S c 2tm   2 c2 c2 t 2  3c2 t  2     1  ln  r  Alk 2 c23 t 

1

3. Sc1tm1  r , Sc 2tm  r  2   m  1 t  t  t r   m  1 t  t   1  2  Dlk  S1  +  2   S1   m1 2 c1 c1t c1t  c1 t      tm    m  1 t  t   2 r 2   tm1 1 2  3 3S1r 2     Blk    S1    2  S1  2 c1 c1t 2  3c1 t  2      tm I 3_ st 

1

(60)

4. Sc1tm  r & &Sc1tm1  r Sc 2tm  r & &Sc 2tm1  r  2   m  1 t  t  r   m  1 t  t   1  2  Dlk  S1 tm  +  2   S1 tm   2 c1 c1t c1t  c1 t     r   m  1 t  t   1  2 1  2   m  1 t  t   Dlk  S c 2tm  +  2   S c 2tm   2 c2 c2 t c2 t  c2 t    (61) 2 3 3S1 tm r 2   1 r2  2     m  1 t  t    Blk     S1 tm     2   S1 tm   2 c1 c1t 2  3c1 t  2     1

I 3_ st  

2 3Sc 2tm r 2   r2  2     m  1 t  t   Blk  S t   S t       c2 m   c2 m  2 c2 c2 t 2  3c22 t  2     1  ln  r  Alk 2 c23 t



1

5. Sc1tm  r & &Sc1tm1  r Sc 2tm  0 I 3_ st   

 r   m  1 t  t   1  2  2   m  1 t  t   Dlk  S1 tm  +  2   S1 tm   2 c1 c  t c  t c  t 1 1  1      1

 2 3 3S1 tm r 2   r2  2     m  1 t  t    Blk     S1 tm     2   S1 tm   2 c1 c  t 2 3 c  t 2     1 1    1

6. Sc1tm  r Sc 2tm  r I 3_ st  0

(63)

(62)

Appendix H The singularities existing in the time integrations I4 can be obtained by using approximate expansion, which are also piecewise functions. The expressions are as, 1. Sc1tm1  r , Sc 2tm1  r

 1  1 2  1 2   tm1 1 I 4 _ st   D  S  D     S2   lk 1 lk 2 2  c1  2 c2 t  c2   tm  2 c1 t

(64)

2. Sc1tm1  r , Sc 2tm  r & &Sc 2tm1  r  t  1  1 1 1 2  I 4 _ st   Dlk   S12   m1  Alk ln(r )  Dlk  Sc 2tm  2 3 3 2 c2 t   c1    2 c1 t  tm 2 c2 t

(65)

3. Sc1tm1  r , Sc 2tm  r   1  1  t I 4 _ st   Dlk   S12   m1 2  2 c1 t  tm  c1   

(66)

4. Sc1tm  r & &Sc1tm1  r Sc 2tm  r & &Sc 2tm1  r I 4 _ st 

1 1 1 2 2 Dlk  Sc1tm   Alk ln(r )  Dlk  Sc 2tm  3 3 2 c1 t 2 c2 t 2 c23t

5. Sc1tm  r & &Sc1tm1  r Sc 2tm  0 I 4 _ st 

1 2 Dlk  Sc1tm  3 2 c1 t

(68)

6. Sc1tm  r Sc 2tm  r I 4 _ st  0

(69)

(67)