Next Article in Journal
Estimation of the Instantaneous Reproduction Number and Its Confidence Interval for Modeling the COVID-19 Pandemic
Next Article in Special Issue
Co-Circular Polarization Reflector Revisited: Reflection Properties, Polarization Transformations, and Matched Waves
Previous Article in Journal
Voting-Based Ensemble Learning Algorithm for Fault Detection in Photovoltaic Systems under Different Weather Conditions
Previous Article in Special Issue
Analysis of Eigenfrequencies of a Circular Interface Delamination in Elastic Media Based on the Boundary Integral Equation Method
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Direct Method for Solving Singular Integrals in Three-Dimensional Time-Domain Boundary Element Method for Elastodynamics

1
School of Science, Harbin Institute of Technology, Shenzhen 518055, China
2
Urban and Rural Construction Institute, Hebei Agricultural University, Baoding 071001, China
3
School of Civil and Environmental Engineering, Harbin Institute of Technology, Shenzhen 518055, China
*
Author to whom correspondence should be addressed.
Submission received: 20 December 2021 / Revised: 12 January 2022 / Accepted: 15 January 2022 / Published: 17 January 2022
(This article belongs to the Special Issue Analytical Methods in Wave Scattering and Diffraction)

Abstract

:
The analytically time integrable time-space domain (ATI-TSD) is discovered based on which the minimum time-space domain is identified for treatment on singularities in the three-dimensional time-domain boundary element method (3D TD-BEM) formulation. A direct method to solve singular integrals in the 3D TD-BEM formulation for elastodynamic problems is proposed. The wavefront singularity can be analytically eliminated in ATI-TSD, while the dual singularity can be treated by the direct method using Kutt’s quadrature in the identified minimum time-space domain. Three benchmark examples are presented to verify the correctness and the applicability of the direct method for solving the singular integrals in 3D TD-BEM.

1. Introduction

In general, practical dynamic problems can be approximated by numerical solutions [1,2,3]. Among numerical techniques, such as the finite element method (FEM) [4,5,6], the meshless method [7,8,9] and the boundary element method (BEM) [10,11,12,13,14,15], BEM has the advantages of high modelling accuracy, dimension reduction, and automatic satisfaction of boundary conditions at infinity. Therefore, BEM has been successfully applied to many engineering problems [16,17,18,19,20]. In the boundary element methods, the time domain boundary element method (TD-BEM, with TD for time domain) is considered to be the most natural choice for transient dynamic problems [12,13].
In the early historic stages of the development of TD-BEM, Niwa et al. [21], Mansur and Brebbia [22] made important contributions to the complete TD-BEM formulation based on time domain fundamental solutions for 2D scalar wave propagation problems. After that, Karabalis and Beskos [23,24] proposed the first simplification to the transient time domain boundary integral formulations for 3D elastodynamics in the interaction of soil and structure. Subsequently, Banerjee et al. [25] presented a systematic numerical implementation of the 3D TD-BEM formulation which is currently the most widely used and classic formulation. A review of research on accuracy and efficiency was undertaken by Beskos [26,27]. In the process of numerical implementation of the boundary integral equation in TD-BEM, it has been pointed out in the literature (e.g., [10] and [13]) that there are two types of inherent singularity in TD-BEM formulation, the wavefront singularity and the spatial singularity, respectively, arising in cases where the wave triggered from the source point approaches the field point, and the source point coincides with the field point. Mathematically, in 3D time domain fundamental solutions, the coefficient terms consisting of the terms 1 / r or 1 / r 2 cause the spatial singularity in space domain, while the terms consisting of the function terms (e.g., the Heaviside function term, the Dirac function term or the derivative term of the Dirac function) cause the wavefront singularity in the time domain. The spatial singularity exists in the form of a dual singularity, where the wavefront singularity and spatial singularity occur at the same time in both time and space domains. It is impossible for the sole spatial singularity to exist independently in the 3D TD-BEM formulation. However, the wavefront singularity exists in two forms, the dual singularity in the case of r → 0, and the independent wavefront singularity in the case of r 0 . Therefore, in the 3D TD-BEM formulation, there are only two types of singularity, the independent wavefront singularity and the dual singularity. All the coefficients in the coefficient matrixes are double-layered integrals. The independent wavefront singularity is incorporated in the singular integral in the layer in time domain, while the other layer of integral in the space domain is non-singular. In contrast, the dual singularity is incorporated in the singular integral in the double-layered integral in both time and space domains.
In order to solve the singular integrals, Aliabadi [28] further subdivided the integration domain, and then implemented analytical integration in the time domain to eliminate wavefront singularity. Coda and Venturini [29] adopted an improved time domain fundamental solution for constant approximation for both displacement and traction variables, where the piecewise expression of the singular integrals were given. In subsequent articles [30,31], the singular integrals were also solved analytically in the time domain. In contrast, for the dual singularity in both the time and space domains, after the wavefront singularity is analytically solved in time domain, only the weak spatial singularity 1/r and the hyper-spatial singularity 1/r2 are left. The degenerating element method was adopted to generate a Jacobian matrix, approaching zero at the field point, to help deal with the weak spatial singularity [32,33,34,35]. Three methods were proposed to solve the hyper-spatial singularity. The first one is the rigid-body displacement method [25,30,31,36], where the singular integrals are indirectly solved by subtracting the non-singular integral from the whole non-singular integral. Therefore, the rigid body displacement method is an indirect singularity treatment, with the numerical accuracy dependent on all elements involved, causing the accumulation of numerical errors. Additionally, there are the following two direct treatments on the hyper-singularity: In the second method of the treatment on the hyper-spatial singular integral, the hyper-spatial singular integral was processed to reduce the higher order to weak singular integral by special techniques, such as the polar coordinate transformation method [37,38,39] and the kernel function separation method [34,35]. The third method is called the method of Kutt’s quadrature [29,40,41]. Hildenbrand and Kuhn [40] applied Kutt’s quadrature to deal with hyper-singular integrals in the stress boundary integral equations for 2D linear elastostatics. The algorithm was further applied to 3D elastostatic cracks in anisotropic solids by BEM in [41]. For 3D elastodynamic problems, Coda and Venturini [29] adopted Kutt’s quadrature to treat the spatial singularity in the TD-BEM formulation for a small region r [ 0 , 1 ] near the singular point. In the second method, the reduction of the order of the hyper-spatial singularity did not finalize the singularity treatment on the hyper-spatial singularity. The processed singularity was suspended to be further treated by the degenerating element treatment on weak singularity. The third method is limited to a small integration range, and is obviously not applicable to the more general integration range in 3D TD-BEM. From the statement of the historic process of the treatments on the singularity in 3D TD-BEM, the treatment on singular integrals in the time domain, in the sense of the single-layered integration, is complete, direct and analytical; nevertheless, the dual singularity should be fully and analytically solved, in the sense of double-layered integral, in both time and space domains. Therefore, the direct and analytical treatment on singularities in 3D TD-BEM is worth investigating.
The current study is focused on a direct treatment on the weak and hyper singularities in 3D TD-BEM. Based on the characteristics of P- and S-wave propagation and the properties of the Heaviside function, the Dirac function and the derivative of the Dirac function, the time-space domain, where the dual singular integrals are analytically integrable in the time domain, is clearly discovered. Subsequently, the minimum time-space domain where the spatial singular integrals need to be analytically treated is identified. Both the weak and hyper-spatial singularities in the overall integration range are independently and explicitly expressed by Kutt’s quadrature. Finally, the correctness of the proposed treatment on the singularity in 3D TD-BEM is verified by three dynamic examples: a 1D rod and a restrained 2D beam in finite media and a spherical cavity in an infinite medium.
In this paper, as the important part of the 3D TD-BEM formulation, the boundary integral equation and the numerical discretization are briefly described in Section 2, while the main contribution of the current study, the treatment on singularity in boundary integral, is presented in Section 3. Three examples are given to validate the effectiveness of the proposed method for solving the singularity in the 3D TD-BEM formulation in Section 4, and conclusions are drawn in Section 5.

2. Boundary Integral Equation and Numerical Discretization

Neglecting the initial conditions and the body force, the displacement boundary integral equation for 3D elastodynamic problems can be obtained by Graffi’s reciprocity theorem, as expressed in [25,28]:
c k i ( P ) u i ( P , t ) = S 0 t u k i s ( P , τ ; Q , t ) p i ( Q , τ ) d τ d S ( Q ) S 0 t p k i s ( P , τ ; Q , t ) u i ( Q , τ ) d τ d S ( Q )
In Equation (1), c k i is the free term computed by Hartmann [42]. Q and P are, respectively, the field point and the source point. p i and u i (i = 1, 2, 3) represent the traction and displacement at instant τ. u k i s and p k i s are fundamental solutions in the 3D time domain, which can be found in the well known papers [28] and [43], as follows:
u k i s ( P , τ ; Q , t ) = 1 4 π μ ( ψ δ k i χ r , k r , i ) ,
p k i s ( P , τ ; Q , t ) = 1 4 π ψ r χ r r n δ k i + r , i n k 2   χ r r , k n i 2 r , k r , i r n 2 χ r r , k r , i r n + c 1 2 c 2 2 2 ψ r χ r 2   χ r r , k n i .
In Equations (2) and (3), the related parameters are expressed, as:
ψ = c 2 2 r 3 t   H t r c 2 H t r c 1 + 1 r δ t r c 2 , χ = 3 ψ 2 r δ t r c 2 c 2 2 c 1 2 1 r δ t r c 1 , ψ r = χ r 1 r 2 δ t r c 2 + r c 2 δ · t r c 2 , χ r = 3 χ r 1 r 2 δ t r c 2 + r c 2 δ · t r c 2 + c 2 2 c 1 2 1 r 2 δ t r c 1 + r c 1 δ · t r c 1 ,
r i = r i Q r i P , r , i = r x i Q = r x i P = r i r , n i = x i n , r n = r i n i .
In Equations (2)–(5), t = t τ . The variable r represents the distance from the source point to the field point. The notations H t r / c w , δ t r / c w and δ · t r / c w , respectively, stand for the Heaviside function, the Delta function and the derivative of the Delta function, where w = 1, 2. The constants c1 and c2 represent the P-wave and S-wave velocities, respectively, expressed by:
c 1 = λ + 2 μ ρ , c 2 = μ ρ .
In Equation (6), ρ is the material density. The constants μ and λ are Lame constants, expressed as:
μ = E 2 ( 1 + ν ) , λ = ν E ( 1 + ν ) ( 1 2 ν ) .
In Equation (7), E and ν are Young’s modulus and Poisson’s ratio, respectively.
The established boundary integral equation Equation (1) needs to be discretized and after that assembled into the solvable form. The calculated time interval [0, t] is divided into M equal segments with time step Δ t = t / M , such that tm = mΔt. The boundary S is discretized into Ne linear triangular elements. After that, the terms of displacement and traction in Equation (1) can be expressed, as:
S 0 t u k i s p i ( Q , τ ) d τ d S ( Q ) = m = 1 M e = 1 N e a = 1 N q G k i m ; e , a p i m ; e , a , S 0 t p k i s u i ( Q , τ ) d τ d S ( Q ) = m = 1 M e = 1 N e a = 1 N q H k i m ; e , a u i m ; e , a .
In Equation (8), p i m ; e , a and u i m ; e , a , respectively, represent the traction and displacement at point a on element e at instant tm. The notation N q is the amount of nodes of element e. The terms G k i m ; e , a and H k i m ; e , a are the influence coefficient elements of the traction and displacement, expressed by Equations (9) and (10), respectively, as:
G k i m ; e , a =   S t m 1 t m u k i s τ t m 1 Δ t M a d τ d S + S t m t m + 1 u k i s t m + 1 τ Δ t M a d τ d S , when   m   M , S t m 1 t m u k i s τ t m 1 Δ t M a d τ d S ,   when   m = M ,
H k i m ; e , a =   S t m 1 t m p k i s τ t m 1 Δ t M a d τ d S + S t m t m + 1 p k i s t m + 1 τ Δ t M a d τ d S , when   m   M , S t m 1 t m p k i s τ t m 1 Δ t M a d τ d S ,   when   m = M .
In Equations (9) and (10), M a Q is the shape function of the spatial discrete linear triangular element.
After discretization and assemblage in time and space domains, the boundary integral equation Equation (1) is transformed into the equation set in matrix form, expressed as follows:
( C + H M M ) u M = G M M p M + m = 1 M 1 H M m u m + G M m p m .
In Equation (11), u M and p M are the 3 N × 1 overall displacement and traction vector, with N for the number of boundary nodes. C ,   G M M and H M M are the overall matrices with the 3 N × 3 N elements at instant t, formed in the process of assembling the coefficient c k i , the influence coefficient elements G k i m ; e , a and H k i m ; e , a . The matrices u m and p m ( m = 1 ,   2   , , M 1 ) represent the displacement and traction at instant tm, while H M m and G M m are the overall matrices at instant tm, respectively.
Both matrices H M m and G M m are solely dependent on M m . Therefore, once the differential value M m is determined, the matrices H M m and G M m are accordingly determined.

3. Analytical Treatment on Singularity

In the expressions of coefficients of the traction and displacement in Equation (4), if the wave triggered from the source point is approaching the field point, in the case of c w t r , the first type of singularity of wavefront singularity appears. In other words, in the case of t r / c w 0 , the Heaviside function term H w = H t r / c 2 H t r / c 1 , the Dirac function term δ ( t r / c w ) and the derivative of the Dirac function  δ · ( t r / c w ) containing the term ( t r / c w ) are of the wavefront singularity. Therefore, the wavefront singularity is the singularity related only to time. On the other hand, in the case of τ t r / c w and r 0 , the second category singularity of the dual singularity appears. For example, the integral terms, containing one of terms of c 2 2 r 3 t H w ,   1 r 2 δ t r / c 2 + r c 2 δ · t r / c 2 and 1 r δ t r / c 2 , are singular integrals with the dual singularity. So, the dual singularity is the singularity related to both time and space.
In this section, different approaches are employed to handle the two categories of singularities. The treatments on singularities for the two scenarios are independently presented.

3.1. Analytical Treatment on Wavefront Singular Integrals

In Equation (4), the coefficients ( ψ ,   χ ,   ψ r and χ r ) contain the wavefront singularity relevant functions (the Heaviside function term H w = H t r / c 2 H t r / c 1 , the Dirac function term δ ( t r / c w ) and the derivative term of Dirac function δ · ( t r / c w ) , so, the coefficients ( ψ ,   χ ,   ψ r and χ r ) are wavefront singular coefficients. For the case where the wave has not yet arrived at the field point, the influence value should be mathematically assigned 0. In the process of the integration, the expressions are different in different integration domains, therefore, the entire integration domain needs to be divided. Looking into the composition of the function, the three wavefront singularity relevant functions ( H w ,   δ ( t r / c w )   and δ · ( t r / c w ) ) are described by P-wave velocity c1, S-wave velocity c2, the distance r between the source point and the field point, and the time duration t of wave propogation. Therefore, division of the integration domain should be dependent with the variables, r, c1, c2 and t . By introducing the time-space coordinate system, the division of the integration domain, relevant with the P-wave velocity c1 and S-wave velocity c2, is shown in Figure 1.
In Figure 1, the horizontal axis τ is the travelling time of the impulse originated from the source point P, and the vertical axis r is the distance from the field point to the source point. The origin  o is the source point and instant zero, as well. t is the total analysis time, and the coordinates r 1 = c 1 Δ t and r 2 = c 2 Δ t are, respectively, the limits of the distance, that the P-wave and S-wave can propagate from the source point P from instant 0 to instant t. The oblique segments r = c 1 t τ and r = c 2 t τ are the farthest distances that the P-wave and S-wave pulses triggered at the source point P at instant τ can travel during period t τ . Therefore, on the one hand, the upper part r > c 1 ( t τ ) of the oblique segment r = c 1 ( t τ ) is the domain, where the P-wave has not yet reached, so that, H ( c 1 ( t τ ) r ) ) = 0 , δ ( t r / c 1 ) = 0 and δ · ( t r / c 1 ) = 0 . On the other hand, the lower part r < c 1 ( t τ )  of the oblique line segment r = c 1 t τ indicates the domain, where the P-wave has passed over, so that, H ( c 1 ( t τ ) r ) ) = 1 . Similarly, in the upper domain of r = c 2 t τ ,  one has, H ( c 2 ( t τ ) r ) ) = 0 ,   δ ( t r / c 2 ) = 0 and δ · ( t r / c 2 ) = 0 , while in the lower domain of r = c 2 t τ , one has, H ( c 2 ( t τ ) r ) ) = 1 . The prerequisite of the existence and being non zero for Equations (9) and (10) is that the equalities H w = H t r / c 2 H t r / c 1 = 0 ,   δ ( t r / c w ) = 0 and δ · ( t r / c w ) = 0 are not tenable at the same time. This prerequisite is satisfied, when, and only when, the time-space integration domain of the wavefront singular coefficients is located in the shaded part shown in Figure 1. Therefore, in the shaded triangle in Figure 1, the coefficients in the time-space domain are of wavefront singularity, but the singular coefficients in the time domain can be analytically expressed. In this paper the shaded triangle in Figure 1 is defined as the analytically time integrable time-space domain, abbreviated to ATI-TSD, with ATI meaning analytically time integrable, and TSD meaning time-space domain.
In the process of the discretization of the boundary integral equation, the ATI-TSD must be correspondingly discretized into every time element. The ATI-TSD, limited by time element [ t m 1 , t m ], is shown in Figure 2, which is mathematically expressed, as:
D τ r =   τ , r | τ t m 1 , t m , r c 2 t τ , c 1 t τ
In Figure 2, the following relationships are tenable, r 3 = c 1 ( t t m 1 ) ,   r 4 = c 1 ( t t m ) , r 5 = c 2 ( t t m 1 ) , r 6 = c 2 ( t t m ) . Due to the different performance of the wavefront singularity relevant functions (the Heaviside function H w , the Dirac function δ ( t r / c w ) , and the derivative of Dirac function  δ · ( t r / c w ) ) in different integration domains, according to the amount relationship between r 4 and r 5 , the ATI-TSD in time element [ t m 1 , t m ] in Figure 2 is divided into two cases, the case of r 4 > r 5 in Figure 3a and the case of r 4 < r 5 in Figure 3b. Considering the position of r in the shaded area, the ATI-TSD in time element [ t m 1 , t m ] in Figure 3a,b are divided into three domains, respectively denoted by , and in Figure 3a,b.
Taking the integration process of the wavefront singular coefficient ψ as an example, the division mechanism in the sense of mathematics of the ATI-TSD is illustrated. The integration of ψ in time element [ t m 1 , t m ] in Figure 3a,b can be mathematically expressed, respectively, as:
  S t m 1 t m ψ τ t m 1 Δ t d τ d S = r 4 r 3 t m 1 t r c 1 ψ τ t m 1 Δ t d τ d S + r 5 r 4 t m 1 t m ψ τ t m 1 Δ t d τ d S + r 6 r 5 t r c 2 t m ψ τ t m 1 Δ t d τ d S ,
  S t m 1 t m ψ τ t m 1 Δ t d τ d S = r 5 r 3 t m 1 t r c 1 ψ τ t m 1 Δ t d τ d S + r 4 r 5 t r c 2 t r c 1 ψ τ t m 1 Δ t d τ d S + r 6 r 4 t r c 2 t m ψ τ t m 1 Δ t d τ d S .
The integration domains of the first term on the right-hand side of Equations (13) and (14) are mathematically expressed, respectively as, D τ r 1 a =   τ , r | τ t m 1 , t r / c 1 ,   r   r 4 , r 3 and D τ r 1 b =   τ , r | τ t m 1 , t r / c 1 ,   r   r 5 , r 3 , corresponding to ATI-TSD ① in Figure 3a and ATI-TSD ① in Figure 3b. Similarly, the integration domains of the second and the third terms on the right-hand side of Equations (13) and (14) are expressed, respectively, as, D τ r 2 a =   τ , r | τ   t m 1 , t m ,   r   r 5 , r 4 ,   D τ r 2 b =   τ , r | τ   t r / c 2 , t r / c 1 ,   r   r 4 , r 5 , and D τ r 3 a =   τ , r | τ   t r / c 2 , t m ,   r   r 6 , r 5 ,   D τ r 3 b = τ , r | τ   t r / c 2 , t m ,   r   r 6 , r 4 , respectively corresponding to ATI-TSDs ② and ③ in Figure 3a,b.
From the above analysis of the analytical integration of the integral term in time domain for time interval [ t m 1 , t m ], the integral terms containing the wavefront singularity relevant terms, H w ,   δ ( t r / c w ) and δ · ( t r / c w ) , can be analytically expressed in different time integrable domains, tabulated in Table 1 and Table 2. In the tables, the upper and lower limits for every integration domain are given in the first rows, f ( τ ) is a linear time interpolation function. To reduce the size, the following notations are adopted:
δ w 1 = δ ( t r / c 1 ) , δ w 2 = δ ( t r / c 2 ) , δ · w 1 = δ · ( t r / c 1 ) , δ · w 2 = δ · ( t r / c 2 ) .
The integration results of all the wavefront singularity relevant functions ( H w , δ w 1 ,   δ w 2 , δ ˙ w 1 and δ ˙ w 2 ) are determined in Table 1 and Table 2. The analytical integration values in the time domain of coefficient  ψ in different ATI-TSDs can be obtained by recalling the expression of ψ in Equation (4), and the corresponding analytical integration results shown in Table 1 and Table 2. Similarly, the analytical integration results in the time domain of the other coefficients ( χ ,   ψ r and χ r ) can also be obtained. The independent wavefront singularity is only the single-layered singularity in time domain. Therefore, by analytical integration in the time domain, the wavefront singularity in the coefficients ( ψ ,   χ ,   ψ r and χ r ) is analytically eliminated.

3.2. Treatment on Dual Singular Integrals in Both Time and Space Domains

For the dual singularity, in the sense of double-layered integration both in time and space domains in the coefficient matrices H and G, the analytical integration in the time domain is carried out to eliminate the wavefront singularity in advance, and the spatial singularity is treated by the direct method of Kutt’s quadrature.
When τ [tm−1,tm] and m M are satisfied, meaning that  r c 2 ( t t m 1 ) , c 1 ( t t m ) 0 is tenable, therefore, no spatial singularity exists. In the case of m = M, the calculation time interval [tM−1,t] is the last time period in the total analysis duration t. The corresponding ATI-TSD in time element [tM−1,t] is shown in the shaded triangle in Figure 4.
The ATI-TSD in Figure 4 is mathematically expressed, as:
D τ r =   τ , r | τ t M 1 , t ,   r c 2 t τ , c 1 t τ
In Figure 4, the following relationships exist, r 3 = c 1 ( t t M 1 ) = c 1 Δ t ,   r 4 = c 2 ( t t M 1 ) = c 2 Δ t . The ATI-TSD in time element [ t M 1 , t ] in Figure 4 needs to be further subdivided in the space domain to analytically express the integral terms in different integration domains. According to the position of variable r in the shaded triangle in Figure 4, the ATI-TSD is subdivided into two parts, denoted by ① and ②.
Through the above subdivision of the ATI-TSD, the integrals of the dual singular coefficients in Equation (4) ( ψ ,   χ ,   ψ r and χ r ) in different ATI-TSDs in time interval [ t M 1 , t ] are analytically expressed and tabulated in Table 3.
It can be seen from Table 3 that, after the analytical integration in the time domain, the spatial singular integrals in the space domain can be analytically expressed. In the sense of the single-layered integration over space, for the second column in Table 3, the coefficients are not spatially singular, due to r 0 . Nevertheless, for the third column in Table 3, the coefficients are spatially singular due to r 0 . Therefore, only the single-layered singular integral term over space in the third column in the range of D τ r 2 needs to be dealt with.
The treatment on the spatial singularity of the single-layered integrals in the third column in Table 3 is taken as an example to present the treatment method.
When r indefinitely approaches 0, the source point P coincides with the node of the analyzed element. In this senario, the element is called the singular element, and the node, which is coincident with the source point P, is called the singular point. Setting the singular point as the origin, the local oblique coordinate system for the singular element is established on the plane, defined by the singular element, shown in Figure 5a. In the coordinate system, the direction of the z axis is the outer normal direction of the plane defined by the singular element, while the direction of the x axis is the direction from Q1 to Q2, and the direction of the y axis is derived from the cross-product of the direction of z and the direction of x. To facilitate the processing of the spatial singularity, the local oblique coordinate system is transformed into the polar coordinate system, shown in Figure 5b.
In the regularized solution of the coefficient integrals on the boundary elements in the local coordinate system, the dimensionless coordinates of the nodes on the boundary elements are adopted [29]. The dimensionless coordinate of any point Q in the local polar coordinate system in Figure 5b can be expressed as:
η i Q = 1 2 A ( 2 A i 0 + b i x 1 P + a i x 2 P ) + 1 2 A ( b i r cos θ + a i r sin θ ) ,
where,
a i = x 1 k x 1 j ,
b i = x 2 j x 2 k ,
2 A i 0 = x 1 j x 2 k x 1 k x 2 j ,
A = 1 2 ( b 1 a 2 b 2 a 1 ) ,
for i = 1, 2, 3, j = 2, 3, 1, and k = 3, 1, 2.
After relations between the coordinate systems are determined, the single-layered singular integral over space can be performed over two dimensions of r and θ .
By letting the generic function g ( r ) = f ( r ) / r n represent the expressions in terms of r in the third column in Table 3 for time interval [ t M 1 , t ], the integration of g ( r ) in the integration range of D τ r 2 in the polar coordinate system is expressed by the following equation:
S g ( r ) d S = θ 1 θ 2 0 R g ( r ) J d r d θ .
In Equation (19), J is the Jacoby matrix regarding the transformation from the local oblique coordinate system to the element’s polar system, and the upper limit of the integration R is expressed by R = m i n ( L , c 2 Δ t ) , with L = r ( θ ) meaning the radial boundary of the singular element in the polar coordinate system along the wave propagation direction in Figure 5b.
After the coefficients are converted to the local polar coordinate system, Kutt’s quadrature is adopted in this paper to directly solve the singular integral in Equation (19). Kutt’s quadrature can be directly adopted to evaluate any order singular integrals. Interested readers can refer to [44] for more details for Kutt’s quadrature.
Applying Kutt’s N-point equispace quadrature, one has:
0 R f ( r ) r n d r R 1 n i = 1 N ( ω i + c i ln R ( n 1 ) ! ) f ( i 1 N R ) ,
where, ω i and c i are the weights and the coefficients given by Kutt [44].
It can be seen that Kutt’s quadrature retains the finite integrable parts of the Cauchy principal integral. As the optimal formula, the high accuracy of Kutt’s quadrature was verified in several papers [29,40,41]. The single-layered singular integrals over space in the third column of Table 3 are expressed in the form of the integrands in Equation (20). Therefore, by adopting Kutt’s quadrature, the single-layered singular integral terms over space in the third column of Table 3 can be solved with high accuracy.

4. Verification

To validate the proposed direct method to solve singular integrals in the 3D TD-BEM formulation, three benchmark problems are conducted. The first and the second are a 1D rod and a 2D restrained beam under the same Heaviside-type load in finite elastic media, and the third is a spherical cavity under an explosive load in an infinite elastic medium.

4.1. 1D Rod under Heaviside-Type Load in Finite Elastic Medium

The 1D rod under the Heaviside-type load along axis x is chosen as the first benchmark test in the finite elastic medium, shown in Figure 6a. The 1D benchmark example is very classical and was used in several TD-BEM verifications [29,30,35,38]. The length of the rod (l1 = 4 m) is four times the width (l2 = 1 m) and the height (l3 = 1 m). The left end of the rod (x = 0) is fully fixed, and the right end (x= l1) is subjected to the Heaviside-type load p = p 0 H ( t 0 ) , with p0 = 200 Mpa, shown in Figure 6b. The rod is characterized by Young’s modulus E = 2.1 × 10 11 Pa, Poisson’s ratio ν = 0 and mass density ρ = 7.9 × 10 3 kg/m3.
The surface of the rod is uniformly divided into 144 isosceles right triangular linear boundary elements with the right-angle side of 500 mm, resulting in 74 boundary nodes, shown in Figure 7. According to the dimensionless time step parameter β = c1Δt/le [10,11,12,13,14,15], where le represent the length of the elements in the direction of wave propagation, in this example and the following ones, the arrangement β = 0.4~1.0 is adopted to guarantee numerical accuracy.
The displacements in direction x in three periods at the center of the loaded end A (l1, l2/2, l3/2) and at the center of the rod B (l1/2, l2/2, l3/2) are calculated from current 3D TD-BEM and traditional 3D TD-BEM (in which singular integrals are analyzed by the indirect numerical method of rigid body displacement [25,30]), and shown in Figure 8, together with the analytical solutions from [45]. The time history of the normal traction at the center of the fixed end of the rod C (0, l2/2, l3/2) from current 3D TD-BEM, traditional 3D TD-BEM and the analytical solutions from [45] are shown in Figure 9. In the figures, the vertical axis is normalized by E / p 0 l 1 and 1 / p 0 , respectively, while the horizontal axis is normalized by c 1 / l 1 .
It can be found from the figures that the numerical results from the proposed method of solving the singularities in 3D TD-BEM agree well with the analytical solutions. The error of the results which are obtained by the current method is smaller than the traditional 3D TD-BEM in terms of both displacement and traction.

4.2. A 2D Restrained Beam under the Heaviside Load in Finite Elastic Medium

As shown in Figure 10, the 2D restrained beam subjected to the same Heaviside uniform load in the 1D rod example on the upper boundary z = l3 is chosen to verify the proposed 3D TD-BEM. The Poisson’s ratio of the beam is ν = 0.3, while the Young’s modulus and the density of the beam and meshes of the TD-BEM model are the same as those in the 1D rod example.
Because the analytical solution for the problem of the 2D example does not exist, the results from the current TD-BEM formulation and traditional 3D TD-BEM are compared with the FEM solutions. In the FEM model, a much finer mesh size 0.1 m × 0.1 m × 0.1 m and the much shorter time step 1 × 10−5 s are used. The transient displacements uz in three periods at boundary point D (l1/2, l2/2, 0) in Figure 10 from current 3D TD-BEM, traditional 3D TD-BEM, and FEM formulations are shown in Figure 11.
It can be seen that a good agreement between the current 3D TD-BEM and FEM has also been achieved, with a much more favorable mesh size and time step for FEM. Comparison of results obtained by the traditional 3D TD-BEM is performed with FEM which shows that the results from the current 3D TD-BEM are more consistent and accurate.

4.3. A Spherical Cavity under an Explosive Load in an Infinite Elastic Medium

To further explore the applicability to the infinite domain problem of the proposed method for solving singularity in the 3D TD-BEM formulation, the spherical cavity of inner radius r 0 = 1 m in the 3D infinite elastic medium is considered, where the explosive load is applied on the boundary of the cavity. The explosive load is p t = k p 1 e a t e b t , with k = 1.435, a = 1279 , b = 12,792 and p1 = 100 Mpa, shown in Figure 12c. The Young’s modulus and mass density in the example are the same as those of the 1D rod example, nevertheless with the Poisson’s ratio of ν = 0.3 .
The geometry and boundary discretization of the spherical cavity are shown in Figure 12a,b. The inner boundary is divided into 134 triangular linear elements with the mesh size 250   mm × 250   mm × 250   mm , resulting in 120 boundary nodes.
The displacement responses of two monitoring points with the radial distances from the center of the cavity, r = 2 m and r = 8 m, are calculated by the 3D TD-BEM formulation based on the proposed treatment on singularities, the traditional 3D TD-BEM based on the numerical method of rigid body displacement and the method of characteristics [46]. Figure 13 and Figure 14, respectively, show the comparisons between the calculated radial displacement from the current study, traditional 3D TD-BEM and the method of characteristics, where the time is normalized by c 1 / r 0 .
From the comparisons, it can be seen that the results from the proposed method for solving singular integrals in the 3D TD-BEM formulation coincide nicely with the analytical solutions, for both the internal points in the near field and the far field in the surrounding medium. However, the traditional 3D TD-BEM gives results with a lower accuracy.
The results from the above three dynamic examples show that computational accuracy of the current 3D TD-BEM based on the proposed treatment on singularities in this paper is higher than that of the traditional 3D TD-BEM based on the method of rigid body displacement. Next, the calculation efficiency of the current 3D TD-BEM and the traditional 3D TD-BEM is further compared. The time cost of the current 3D TD-BEM and the traditional 3D TD-BEM in the three examples is calculated as shown in Table 4.
It can be seen from Table 4 that the time cost of the current 3D TD-BEM does not increase significantly compared to that of the traditional 3D TD-BEM. The main reason is that the traditional 3D TD-BEM is based on the indirect numerical method of rigid body displacement, in which the singular integrals are indirectly solved by subtracting the non-singular integral from the whole non-singular integral. Therefore, the solution of singular integrals in the traditional 3D TD-BEM depends on the numerical solution of all elements involved. The current 3D TD-BEM is based on the proposed treatment on singularities in this paper, which directly solves the singularity on the singular element, without being dependent on the calculation results of all elements involved.

5. Conclusions

The direct method using Kutt’s quadrature to solve the singular integrals in time and space domains in the 3D TD-BEM formulation is developed, with the following conclusions.
  • The analytically time integrable time-space domain (ATI-TSD) is discovered, and the minimum time-space domain is identified for the analytical treatment on the spatial singularity. Based on the analysis of the integration domain, in the sense of single-layered integration, the wavefront singularity can be analytically eliminated in ATI-TSD, while the dual singularity can be analytically eliminated, in the sense of double-layered integration, in the identified minimum time space domain.
  • Regardless of the wavefront singularity or the dual singularity, by analytically treating the wavefront singularity in the time domain in advance, the spatial singularity is then treated in the space domain in the last time step, realizing the full treatment on the singularities in double-layered integration in both time and space domains.
  • On the comparison with the traditional 3D TD-BEM based on the rigid body displacement method, the current TD-BEM based on the proposed treatment on singularities in this paper has higher calculation accuracy, without increasing the time cost of calculation.
  • The correctness of the proposed direct method using Kutt’s quadrature to solve the singularity in the 3D TD-BEM formulation is verified by three different examples with satisfactory results.

Author Contributions

Conceptualization, X.Q. and W.L.; methodology, X.Q., W.L., H.L. and Y.F.; software, X.Q.; project administration, Y.F. and W.L.; validation, X.Q., W.L. and H.L.; writing—original draft preparation, X.Q. and W.L.; writing—review and editing, X.Q., W.L., H.L. and Y.F. All authors have read and agreed to the published version of the manuscript.

Funding

This research received financial support from research grants, No. 2021YFC3001002 from the National Key R&D Program of China, No. JCYJ20210324121402008 from the Shenzhen Science and Technology Innovation Commission, and No. 51778193 from the National Natural Science Foundation of China.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data sharing not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Liu, X.; Ahsan, M.; Ahmad, M.; Nisar, M.; Liu, X.L.; Ahmad, I.; Ahmad, H. Applications of Haar wavelet-finite difference hybrid method and its convergence for hyperbolic nonlinear Schrödinger equation with energy and mass conversion. Energies 2021, 14, 7831. [Google Scholar] [CrossRef]
  2. Feng, S.Z.; Han, X.; Li, Z.X.; Incecik, A. Ensemble learning for remaining fatigue life prediction of structures with stochastic parameters: A data-driven approach. Appl. Math. Model. 2022, 101, 420–431. [Google Scholar] [CrossRef]
  3. Liu, X.; Ahsan, M.; Ahmad, M.; Hussian, I.; Alqarni, M.M.; Mahmoud, E.E. Haar wavelets multi-resolution collocation procedures for two-dimensional nonlinear Schrödinger equation. Alex. Eng. J. 2021, 60, 3057–3071. [Google Scholar] [CrossRef]
  4. Feng, S.Z.; Han, X. A novel multi-grid based reanalysis approach for efficient prediction of fatigue crack propagation. Comput. Methods Appl. Mech. Eng. 2019, 353, 107–122. [Google Scholar] [CrossRef]
  5. Zienkiewicz, O.C.; Taylor, R.L. The Finite Element Method, 3rd ed.; McGraw-Hill: London, UK, 1977. [Google Scholar]
  6. Lv, J.H.; Jiao, Y.Y.; Wriggers, P.; Rabczuk, T.; Feng, X.T.; Tan, F. Efficient integration of crack singularities in the extended finite element method: Duffy-distance transformation and conformal preconditioning strategy. Comput. Methods Appl. Mech. Eng. 2018, 340, 559–576. [Google Scholar] [CrossRef]
  7. Liu, Y.H.; Chen, S.S.; Li, J.; Cen, Z.Z. A meshless local natural neighbour interpolation method applied to structural dynamic analysis. Comput. Model. Eng. Sci. 2008, 31, 145–156. [Google Scholar]
  8. Zaheer-ud-Din; Ahsan, M.; Ahmad, M.; Khan, W.; Mahmoud, E.E.; Abdel-Aty, A.-H. Meshless analysis of nonlocal boundary value problems in anisotropic and inhomogeneous media. Mathematics 2020, 8, 2045. [Google Scholar] [CrossRef]
  9. Chen, S.S.; Wang, W.; Zhao, X.S. An interpolating element-free Galerkin scaled boundary method applied to structural dynamic analysis. Appl. Math. Model. 2019, 75, 494–505. [Google Scholar] [CrossRef]
  10. Mansur, W.J. A Time-Stepping Technique to Solve Wave Propagation Problems Using the Boundary Element Method. Ph.D. Thesis, University of Southampton, Southampton, UK, 1983. [Google Scholar]
  11. Lei, W.D.; Li, H.J.; Qin, X.F.; Chen, R.; Ji, D.F. Dynamics-based analytical solutions to singular integrals for elastodynamics by time domain boundary element method. Appl. Math. Model. 2018, 56, 612–625. [Google Scholar] [CrossRef]
  12. Dominguez, J. Boundary Elements in Dynamics; Wit Press: Southampton, UK, 1993. [Google Scholar]
  13. Xie, G.Z.; Zhong, Y.D.; Zhou, F.L.; Du, W.L.; Li, H.; Zhang, D.H. Singularity cancellation method for time-domain boundary element formulation of elastodynamics: A direct approach. Appl. Math. Model. 2020, 80, 647–667. [Google Scholar] [CrossRef]
  14. Carrer, J.A.M.; Mansur, W.J. Stress and velocity in 2D transient elastodynamic analysis by the boundary element method. Eng. Anal. Bound. Elem. 1999, 23, 233–245. [Google Scholar] [CrossRef]
  15. Lei, W.D.; Ji, D.F.; Li, H.J.; Li, Q.X. On an analytical method to solve singular integrals both in space and time for 2-D elastodynamics by TD-BEM. Appl. Math. Model. 2015, 39, 6307–6318. [Google Scholar] [CrossRef]
  16. Moro, F.; Codecasa, L. Coupling the Cell Method with the Boundary Element Method in static and quasi-static electromagnetic problems. Mathematics 2021, 9, 1426. [Google Scholar] [CrossRef]
  17. Peng, X.; Atroshchenko, E.; Kerfriden, P.; Bordas, S.P.A. Isogeometric boundary element methods for three dimensional static fracture and fatigue crack growth. Comput. Methods Appl. Mech. Eng. 2017, 316, 151–185. [Google Scholar] [CrossRef] [Green Version]
  18. Panji, M.; Ansari, B. Modeling pressure pipe embedded in two-layer soil by a half-plane BEM. Comput. Geotech. 2017, 81, 360–367. [Google Scholar] [CrossRef]
  19. Panji, M.; Mojtabazadeh-Hasanlouei, S.; Yasemi, F. A half-plane time-domain BEM for SH-wave scattering by a subsurface inclusion. Comput. Geosci. 2020, 134, 104342. [Google Scholar] [CrossRef]
  20. Li, Y.B.; Mulani, S.B.; Fei, Q.G.; Wu, S.Q.; Zhang, P. Vibro-acoustic analysis under stationary and non-stationary random excitations with KLE/FEM/BEM. Aerosp. Sci. Technol. 2017, 66, 203–215. [Google Scholar] [CrossRef]
  21. Niwa, Y.; Fukui, T.; Kato, S.; Fujiki, K. An application of the integral equation method to two-dimensional elastodynamics. Theor. Appl. Mech. 1980, 28, 281–290. [Google Scholar]
  22. Mansur, W.J.; Brebbia, C.A. Formulation of the boundary element method for transient problems governed by the scalar wave equation. Appl. Math. Model. 1982, 6, 307–311. [Google Scholar] [CrossRef]
  23. Karabalis, D.L.; Beskos, D.E. Dynamic response of 3-D rigid surface foundamentions by time domain element method. Earthq. Eng. Struct. Dyn. 1984, 12, 73–93. [Google Scholar] [CrossRef]
  24. Karabalis, D.L.; Beskos, D.E. Dynamic response of 3-D flexible foundamentions by time domain BEM and FEM. Int. J. Soil Dyn. Earthq. Eng. 1985, 22, 91–101. [Google Scholar]
  25. Banerjee, P.K.; Ahmad, S.; Manolis, G.D. Transient elastodynamic analysis of three-dimensional problems by boundary element method. Earthq. Eng. Struct. Dyn. 1986, 14, 933–949. [Google Scholar] [CrossRef]
  26. Beskos, D.E. Boundary element methods in dynamic analysis. Appl. Mech. Rev. 1987, 40, 1–23. [Google Scholar] [CrossRef]
  27. Beskos, D.E. Boundary element methods in dynamic analysis: Part II (1986–1996). Appl. Mech. Rev. 1997, 50, 149–197. [Google Scholar] [CrossRef]
  28. Aliabadi, M.H. The Boundary Element Method. Volume 2: Applications in Solids and Structures; Wiley: New York, NY, USA, 2002. [Google Scholar]
  29. Coda, H.B.; Venturini, W.S. Three-dimensional transient BEM analysis. Comput. Struct. 1995, 56, 751–768. [Google Scholar] [CrossRef]
  30. Marrero, M.; Dominguez, J. Numerical behavior of time domain BEM for three-dimensional transient elastodynamic problems. Eng. Anal. Bound. Elem. 2003, 27, 39–48. [Google Scholar] [CrossRef]
  31. Panagiotopoulos, C.G.; Manolis, G.D. Three-dimensional BEM for transient elastodynamics based on the velocity reciprocal theorem. Eng. Anal. Bound. Elem. 2011, 35, 507–516. [Google Scholar] [CrossRef]
  32. Ashlock, J.C. Computational and Experimental Modeling of Dynamic Foundation Interactions with Sand. Ph.D. Thesis, University of Colorado at Boulder, Boulder, CO, USA, 2006. [Google Scholar]
  33. Luchi, M.L.; Rizzuti, S. Boundary elements for three-dimensional elastic crack analysis. Int. J. Numer. Methods Eng. 1987, 24, 2253–2271. [Google Scholar] [CrossRef]
  34. Pak, R.Y.S.; Guzina, B.B. Seismic soil-structure interaction analysis by direct boundary element methods. Int. J. Solids. Struct. 1999, 36, 4743–4766. [Google Scholar] [CrossRef]
  35. Pak, R.Y.S.; Bai, X.Y. A regularized boundary element formulation with weighted-collocation and higher-order projection for 3D time-domain elastodynamics. Eng. Anal. Bound. Elem. 2018, 93, 135–142. [Google Scholar] [CrossRef]
  36. Manolis, G.D.; Beskos, D.E. Boundary Element Methods in Elastodynamics; Unwin Hyman: London, UK, 1988. [Google Scholar]
  37. Rizos, D.C.; Karabalis, D.L. An advanced direct time domain BEM formulation for general 3-D elastodynamic problems. Comput. Mech. 1994, 15, 249–269. [Google Scholar] [CrossRef]
  38. Li, H.B.; Han, G.M.; Mang, H.A. A new method for evaluating singular integrals in stress analysis of solids by the direct boundary element method. Int. J. Numer. Meth. Eng. 1985, 21, 2071–2098. [Google Scholar] [CrossRef]
  39. Brebbia, C.A.; Telles, J.C.F.; Wrobel, L.C. Boundary Element Techniques: Theory and Applications in Engineering; Springer: Berlin/Heidelberg, Germany, 1984. [Google Scholar]
  40. Hildenbrand, J.; Kuhn, G. Numerical computation of hypersingular integrals and application to the boundary integral equation for the stress tensor. Eng. Anal. Bound. Elem. 1992, 10, 209–217. [Google Scholar] [CrossRef]
  41. Pan, E.; Yuan, F.G. Boundary element analysis of three-dimensional cracks in anisotropic solids. Int. J. Numer. Meth. Eng. 2000, 48, 211–237. [Google Scholar] [CrossRef]
  42. Hartmann, F. Computing the C-matrix on non-smooth boundary points. In New Developments in Boundary Element Methods; Brebbia, C.A., Ed.; Butterworths: London, UK, 1980; pp. 367–379. [Google Scholar]
  43. Wen, P.H.; Aliabadi, M.H.; Rooke, D.P. The influence of elastic waves on dynamic stress intensity factors (three-dimensional problems). Arch. Appl. Mech. 1996, 66, 385–394. [Google Scholar] [CrossRef]
  44. Kutt, H.R. The numerical evaluation of principal-value integrals by finite-part integration. Numer. Math. 1975, 24, 204–210. [Google Scholar] [CrossRef]
  45. Eringen, A.C.; Suhubi, E.S. Elastodynamics. Volume 2: Linear Theory; Academic Press: New York, NY, USA, 1975. [Google Scholar]
  46. Chou, P.C.; Koenig, H.A. A unified approach to cylindrical and spherical elastic waves by method of characteristics. J. Appl. Mech. 1966, 33, 159–167. [Google Scholar] [CrossRef]
Figure 1. Sketch of time-space integration domain.
Figure 1. Sketch of time-space integration domain.
Mathematics 10 00286 g001
Figure 2. Division of ATI-TSD.
Figure 2. Division of ATI-TSD.
Mathematics 10 00286 g002
Figure 3. Subdivision of ATI-TSD in [ t m 1 , t m ]: (a) the case of r 4 > r 5 ; (b) the case of r 4 < r 5 .
Figure 3. Subdivision of ATI-TSD in [ t m 1 , t m ]: (a) the case of r 4 > r 5 ; (b) the case of r 4 < r 5 .
Mathematics 10 00286 g003
Figure 4. Subdivision of ATI-TSD in time element [ t M 1 , t ].
Figure 4. Subdivision of ATI-TSD in time element [ t M 1 , t ].
Mathematics 10 00286 g004
Figure 5. Local coordinate systems: (a) local oblique coordinate system; (b) local polar coordinate system.
Figure 5. Local coordinate systems: (a) local oblique coordinate system; (b) local polar coordinate system.
Mathematics 10 00286 g005
Figure 6. Sketch of the example of 1D rod: (a) geometry and distribution of the load; (b) Heaviside-type load.
Figure 6. Sketch of the example of 1D rod: (a) geometry and distribution of the load; (b) Heaviside-type load.
Mathematics 10 00286 g006
Figure 7. Boundary discretization for 1D rod.
Figure 7. Boundary discretization for 1D rod.
Mathematics 10 00286 g007
Figure 8. Comparison between the results of displacement from current 3D TD-BEM, traditional 3D TD-BEM, and the analytical solutions.
Figure 8. Comparison between the results of displacement from current 3D TD-BEM, traditional 3D TD-BEM, and the analytical solutions.
Mathematics 10 00286 g008
Figure 9. Comparison between the results of traction from current 3D TD-BEM, traditional 3D TD-BEM, and the analytical solutions.
Figure 9. Comparison between the results of traction from current 3D TD-BEM, traditional 3D TD-BEM, and the analytical solutions.
Mathematics 10 00286 g009
Figure 10. Sketch of the example of the restrained beam.
Figure 10. Sketch of the example of the restrained beam.
Mathematics 10 00286 g010
Figure 11. Comparison between the results from current 3D TD-BEM, traditional 3D TD-BEM and FEM in terms of the displacement uz at boundary point D.
Figure 11. Comparison between the results from current 3D TD-BEM, traditional 3D TD-BEM and FEM in terms of the displacement uz at boundary point D.
Mathematics 10 00286 g011
Figure 12. Sketch of the example of spherical cavity: (a) geometry; (b) boundary discretization; (c) explosive load.
Figure 12. Sketch of the example of spherical cavity: (a) geometry; (b) boundary discretization; (c) explosive load.
Mathematics 10 00286 g012
Figure 13. Calculated radial displacements at r = 2 m from the proposed 3D TD-BEM, the traditional 3D TD-BEM and the analytical solution.
Figure 13. Calculated radial displacements at r = 2 m from the proposed 3D TD-BEM, the traditional 3D TD-BEM and the analytical solution.
Mathematics 10 00286 g013
Figure 14. Calculated radial displacements at r = 8 m from the proposed 3D TD-BEM, the traditional 3D TD-BEM and the analytical solution.
Figure 14. Calculated radial displacements at r = 8 m from the proposed 3D TD-BEM, the traditional 3D TD-BEM and the analytical solution.
Mathematics 10 00286 g014
Table 1. Integration results of H w , δ w 1 , δ w 2 , δ · w 1 and δ · w 2 in time element [ t m 1 , t m ] in the case of r 4 > r 5 .
Table 1. Integration results of H w , δ w 1 , δ w 2 , δ · w 1 and δ · w 2 in time element [ t m 1 , t m ] in the case of r 4 > r 5 .
(1)
τ [ t m 1 , t r / c 1 ] r [ r 4 , r 3 ]
(2)
τ [ t m 1 , t m ] r [ r 5 , r 4 ]
(3)
τ [ t r / c 2 , t m ] r [ r 6 , r 5 ]
t m 1 t m H w d τ 1 1 1
t m 1 t m δ w 1 f ( τ ) d τ f ( t r / c 1 ) 00
t m 1 t m δ w 2 f ( τ ) d τ 00 f ( t r / c 2 )
t m 1 t m δ · w 1 f ( τ ) d τ f ( t r / c 1 ) 00
t m 1 t m δ · w 2 f ( τ ) d τ 00 f ( t r / c 2 )
Table 2. Integration results of H w ,   δ w 1 , δ w 2 , δ · w 1 and δ · w 2 in time element [ t m 1 , t m ] in the case of r 4 < r 5 .
Table 2. Integration results of H w ,   δ w 1 , δ w 2 , δ · w 1 and δ · w 2 in time element [ t m 1 , t m ] in the case of r 4 < r 5 .
(1)
τ [ t m 1 , t r / c 1 ] r [ r 5 , r 3 ]
(2)
τ [ t r / c 2 , t r / c 1 ] r [ r 4 , r 5 ]
(3)
τ [ t r / c 2 , t m ] r [ r 6 , r 4 ]
t m 1 t m H w d τ 1 1 1
t m 1 t m δ w 1 f ( τ ) d τ f ( t r / c 1 ) f ( t r / c 1 ) 0
t m 1 t m δ w 2 f ( τ ) d τ 0 f ( t r / c 2 ) f ( t r / c 2 )
t m 1 t m δ · w 1 f ( τ ) d τ f ( t r / c 1 ) f ( t r / c 1 ) 0
t m 1 t m δ · w 2 f ( τ ) d τ 0 f ( t r / c 2 ) f ( t r / c 2 )
Table 3. Integration results of the dual singular coefficients ( ψ , χ ,   ψ r and χ r ) in time interval [ t M 1 , t ].
Table 3. Integration results of the dual singular coefficients ( ψ , χ ,   ψ r and χ r ) in time interval [ t M 1 , t ].
D τ r 1 =   τ , r τ t M 1 , t r / c 1 r r 4 , r 3 D τ r 2 =   τ , r τ t r / c 2 , t r / c 1 r 0 , r 4
t M 1 t ψ f ( τ ) d τ ( c 2 Δ t ) 2 6 1 r 3 + c 2 2 2 c 1 2 1 r c 2 2 3 c 1 3 Δ t c 2 2 2 c 1 2 + 1 2 1 r 1 3 c 2 2 c 1 3 Δ t + 2 c 2 Δ t
t M 1 t χ f ( τ ) d τ ( c 2 Δ t ) 2 2 1 r 3 + c 2 2 2 c 1 2 1 r c 2 2 2 c 1 2 1 2 1 r
t M 1 t ψ r f ( τ ) d τ ( c 2 Δ t ) 2 2 1 r 4 c 2 2 2 c 1 2 1 r 2 c 2 2 2 c 1 2 + 1 2 1 r 2
t M 1 t χ r f ( τ ) d τ 3 ( c 2 Δ t ) 2 2 1 r 4 c 2 2 2 c 1 2 1 r 2 c 2 2 2 c 1 2 1 2 1 r 2
Table 4. Time cost of the current 3D TD-BEM and the traditional 3D TD-BEM in the three examples.
Table 4. Time cost of the current 3D TD-BEM and the traditional 3D TD-BEM in the three examples.
1-D Rod2D Restrained BeamSpherical Cavity
Current 3D TD-BEM123.461 s118.224 s71.137 s
Traditional 3D TD-BEM122.476 s114.137 s70.339 s
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Qin, X.; Fan, Y.; Li, H.; Lei, W. A Direct Method for Solving Singular Integrals in Three-Dimensional Time-Domain Boundary Element Method for Elastodynamics. Mathematics 2022, 10, 286. https://0-doi-org.brum.beds.ac.uk/10.3390/math10020286

AMA Style

Qin X, Fan Y, Li H, Lei W. A Direct Method for Solving Singular Integrals in Three-Dimensional Time-Domain Boundary Element Method for Elastodynamics. Mathematics. 2022; 10(2):286. https://0-doi-org.brum.beds.ac.uk/10.3390/math10020286

Chicago/Turabian Style

Qin, Xiaofei, Youhua Fan, Hongjun Li, and Weidong Lei. 2022. "A Direct Method for Solving Singular Integrals in Three-Dimensional Time-Domain Boundary Element Method for Elastodynamics" Mathematics 10, no. 2: 286. https://0-doi-org.brum.beds.ac.uk/10.3390/math10020286

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop