Three-Dimensional Elastodynamic Analysis Employing Partially Discontinuous Boundary Elements

: Compared with continuous elements, discontinuous elements advance in processing the discontinuity of physical variables at corner points and discretized models with complex boundaries. However, the computational accuracy of discontinuous elements is sensitive to the positions of element nodes. To reduce the side effect of the node position on the results, this paper proposes employing partially discontinuous elements to compute the time-domain boundary integral equation of 3D elastodynamics. Using the partially discontinuous element, the nodes located at the corner points will be shrunk into the element, whereas the nodes at the non-corner points remain unchanged. As such, a discrete model that is continuous on surfaces and discontinuous between adjacent surfaces can be generated. First, we present a numerical integration scheme of the partially discontinuous element. For the singular integral, an improved element subdivision method is proposed to reduce the side effect of the time step on the integral accuracy. Then, the effectiveness of the proposed method is veriﬁed by two numerical examples. Meanwhile, we study the inﬂuence of the positions of the nodes on the stability and accuracy of the computation results by cases. Finally, the recommended value range of the inward shrink ratio of the element nodes is provided.


Introduction
Due to the advantages of dimensionality reduction, semi-analysis, and being suitable for dealing with stress concentration and infinite domain problems, the boundary element method (BEM) [1] has been successfully adopted in different environments, such as fracture mechanics [2][3][4] and acoustics [5][6][7][8]. In contrast to the finite element method (FEM) [9,10], the BEM does not require the continuity of physical quantities between elements [11]. This point can be confirmed from the successful application of constant boundary elements. Therefore, the BEM is adopted to analyze the elastodynamic problem in this paper.
According to the positions of the nodes, the boundary elements are classified as continuous (conforming) or discontinuous (non-conforming) elements. Continuous elements share common nodes with their neighbors, whereas discontinuous elements do not [12]. As discontinuous boundary elements have an advantage in processing the geometric corner and boundary conditions, they are widely employed in the BEM [13][14][15]. Moreover, the discontinuous element has good adaptability to the boundary shape of the analysis domain so that it is easy to realize the mesh discretization of complex models. However, the positions of nodes in the discontinuous element significantly affect the computational accuracy. Xu and Brebbia [16] first discussed the problem through the computation of the two-dimensional elastostatic problem and made the conclusion that the optimal position of (1) The numerical integral scheme of the partially discontinuous elements is studied, in which the calculation of a singular integral [27][28][29] is carefully examined. An improved singular element subdivision method is proposed to reduce the side effect of the time step on the integral accuracy. (2) The influence of the node position of the partially discontinuous elements on the computational accuracy and result stability is studied through two classical numerical examples. The recommended value range of the inward shift ratio of the element node is provided.
The remainder of the paper is organized as follows: We introduce the boundary integral equation of 3D elastodynamics and its numerical implementation in Section 2. The numerical integral scheme of partially discontinuous elements together with the processing of the singular integral is described in Section 3. The experiments and analysis are shown in Section 4. Finally, we summarize the paper in the last section.

Time-Domain Boundary Integral Equation of Elastodynamics and Its Numerical Implementation
In this section, the time-domain boundary integral equation of elastodynamics and its numerical implementation are presented. As the basis of the research in this paper, the formulas shown in this section are well known and refer to the classic literature [23] in this field.
with boundary conditions u i (q, t) = u i (q, t), ∀q ∈ S ui p i (q, t) = p i (q, t), ∀q ∈ S pi (2) and initial conditions in which the symbols used are as follows: • i and j represent the three directions of the 3D space; • u i represents the displacement component; u i is the second-order derivative of displacement u i versus time t; • u i and p i represent the boundary displacement on S ui and the boundary traction on S pi , respectively, where the whole boundary S = S ui ∪ S pi ; • V represents the whole domain of the analysis object; • ρ is the medium density; • λ and G are the Lame constants given by where E and ν represent Young's modulus and Poisson's ratio, respectively.
The time-domain boundary integral equation (BIE) corresponding to the governing Equation (1) can be expressed as C ij (p)u j (p, t) = s t 0 u s ij (p, q; t − τ)p j (q, τ)dτdS(q) − s t 0 p s ij (p, q; t − τ)u j (q, τ)dτdS (q) +ρ V t 0 u s ij (p, q; t − τ)b j (q, τ)dτdV(q) + ρ V u s ij (p, q; t)v 0 i (q)dV(q) +ρ V . u s ij (p, q; t)u 0 i (q)dV(q) (5) in which the meanings of the symbols used are as follows: • C ij = δ ij /2 when the source point p is located on smooth surfaces, and C ij = δ ij when the source point p is located in the domain (δ ij = 1 when i = j, otherwise δ ij = 0); • u s ij (p, q; t − τ) and p s ij (p, q; t − τ) are the time-dependent displacement fundamental solution and the traction fundamental solution, respectively, whose expressions can be found in [23]; • . u s ij (p, q; t) represents the velocity fundamental solution, which is the first-order derivative of u s ij (p, q; t − τ) versus time. Assuming that the body force and initial conditions are both zero, Equation (5) can be simplified to C ij (p)u j (p, t) = s t 0 u s ij (p, q; t − τ)p j (q, τ)dτdS(q) − s t 0 p s ij (p, q; t − τ)u j (q, τ)dτdS(q) (6) (6) can be solved numerically by using the time step method to deal with the time integral and the collocation point method to process the space integral. If we discretize the time interval [0, t] into M time steps of duration ∆t, the discretized time nodal t M = M∆t, where M = 1, 2, . . . , N. Lagrange interpolation is used to approximate the dynamic response at any time in the discrete time interval [t M−1 , t M ]. In this paper, both displacement and traction adopt linear interpolation in the time interval, that is, for a specific time step M, the time interpolation formulas of displacement and traction are as follows: in which Substituting Equation (7) into Equation (6), the BIE of the time-discrete format can be obtained as follows: The time integral in Equation (10) can be calculated analytically. After time integration, the BIE without the time variable τ can be obtained, as shown below: (11) in which U Mm ij and P Mm ij represent the displacement and traction kernel functions after time integration, respectively, and u 0 j (q) and p 0 j (q) represent the displacement and traction at the initial moment, respectively. Under zero initial conditions, u 0 j (q) = p 0 j (q) = 0. After discretizing the boundary S into L elements, we use the K-th Lagrange interpolation function to formulate the BIE of the spatial discretization: in which u mlk j and p mlk j , respectively, represent the displacement and traction at the kth node in the l-th element at t = m∆t, and N k (q) represents the shape function at the integration point q. Equation (12) can use the Gauss-Legendre quadrature method to perform numerical integration on each element. Then, the physical variables at all element nodes can be obtained by solving the linear equations derived from Equation (12).
In this paper, the spatial discretization of Equation (12) adopts partially discontinuous elements. As shown in Figure 1, the blue lines represent the boundary of the surface, the black lines represent the edges of the mesh elements, the hollow dots represent the original element nodes, and the red solid dots represent the element collocation points which are moved inward. The nodes located on the surface but not on the edges (such as black solid dots in Figure 1) remain unchanged. It can be seen that in the partially discontinuous element, only the nodes located at the corner points are shrunk to the surface, and other nodes remain unchanged. Therefore, the partially discontinuous element can avoid the discontinuity problem of physical variables at the corner points without adding too many nodes. It is a compromise between the continuous element and the completely discontinuous element. gration point q. Equation (12) can use the Gauss-Legendre quadrature method to perform numerical integration on each element. Then, the physical variables at all element nodes can be obtained by solving the linear equations derived from Equation (12).
In this paper, the spatial discretization of Equation (12) adopts partially discontinuous elements. As shown in Figure 1, the blue lines represent the boundary of the surface, the black lines represent the edges of the mesh elements, the hollow dots represent the original element nodes, and the red solid dots represent the element collocation points which are moved inward. The nodes located on the surface but not on the edges (such as black solid dots in Figure 1) remain unchanged. It can be seen that in the partially discontinuous element, only the nodes located at the corner points are shrunk to the surface, and other nodes remain unchanged. Therefore, the partially discontinuous element can avoid the discontinuity problem of physical variables at the corner points without adding too many nodes. It is a compromise between the continuous element and the completely discontinuous element.

Numerical Integration Scheme of the Partially Discontinuous Element and the Treatment of the Singular Integral
In this section, the numerical integration scheme of the partially discontinuous element is presented according to [20][21][22]. To reduce the side effect of the time step on the singular integral accuracy, an improved singular element subdivision method is proposed in this section.

Numerical Integration Scheme
For convenience, we take the numerical integration of an element in Equation (12) as an example. The integral of the l-th element for the source point p is as follows: (13) in which * ( , )= Mm ij ij FU pq or Mm ij P , () k Q q represents the shape function of the partially discontinuous element. Equation (13) is usually numerically calculated by the Gauss-Legendre quadrature method. Before applying the Gauss-Legendre quadrature method to

Numerical Integration Scheme of the Partially Discontinuous Element and the Treatment of the Singular Integral
In this section, the numerical integration scheme of the partially discontinuous element is presented according to [20][21][22]. To reduce the side effect of the time step on the singular integral accuracy, an improved singular element subdivision method is proposed in this section.

Numerical Integration Scheme
For convenience, we take the numerical integration of an element in Equation (12) as an example. The integral of the l-th element for the source point p is as follows: (13) in which F * ij (p, q) = U Mm ij or P Mm ij , Q k (q) represents the shape function of the partially discontinuous element. Equation (13) is usually numerically calculated by the Gauss-Legendre quadrature method. Before applying the Gauss-Legendre quadrature method to Equation (13), the coordinate system should be transformed from the real 3D space xyz to the regular space ξηw, as shown in Figure 2.  The transformation formula is formulated as follows:  The transformation formula is formulated as follows: where |J(q)| local is the determinant of the Jacobian matrix of the transformation from the xyz space to the uvw space, and |J(q)| element is the determinant of the Jacobian matrix of the transformation from the uvw space to the ξηw space. The calculation method of |J(q)| local and |J(q)| element in the partially discontinuous element is the same as that in the continuous element [1].
Since the interpolation points in partially discontinuous elements are the element nodes after inward shifting, the shape function Q k (q) is different from the standard shape function N k (q). As shown in Figure 3, d k is the position of node k after moving inside the element.
The transformation formula is formulated as follows: Jq in the partially discontinuous element is the same continuous element [1].
Since the interpolation points in partially discontinuous elements are nodes after inward shifting, the shape function Figure 3, dk is the position of node k inside the element. The physical variable at any point in the element can be calculated by shape function as The u at the interpolation point dk can be calculated by The physical variable at any point in the element can be calculated by the standard shape function as The u at the interpolation point d k can be calculated by which can be rewritten in a matrix form as follows: We can reverse Equation (17) to obtain and substituting Equation (18) into Equation (15), we can obtain The calculation formula of the new shape function Q k (q) can be obtained as shown in Equation (21). For the interpolation point d k that does not move inward, just take ξ d k , η d k as the original node coordinates.

Singular Integral
In this paper, the coordinate transformation method combined [29] with an improved element subdivision method is adopted to deal with weakly singular integrals (f (1/r)), while the rigid body displacement method [23] is used to compute strongly singular integrals (f (1/r 2 )) indirectly. In the traditional element subdivision method, the integral element is divided into several triangular patches according to the position of the source point p, as shown in Figure 4. Then, the α-β transformation method [29] is used to eliminate the weak singularity in each patch. point p, as shown in Figure 4. Then, the α-β transformation method [29] is used to eliminate the weak singularity in each patch. However, the value of the fundamental solutions of elastodynamics is related to the time step. As shown in Figure 5, the curves of the displacement kernel function MM ij U are different under different time steps. In Figure 5, r represents the distance between the source point p and the field point q. It can be seen that the function value changes rapidly when r < c2Δt , in which c2 represents the shear wave velocity and Δt represents the time step. The smaller the time step Δt, the more drastically the function value changes. However, the value of the fundamental solutions of elastodynamics is related to the time step. As shown in Figure 5, the curves of the displacement kernel function U MM ij are different under different time steps. In Figure 5, r represents the distance between the source point p and the field point q. It can be seen that the function value changes rapidly when r < c 2 ∆t, in which c 2 represents the shear wave velocity and ∆t represents the time step. The smaller the time step ∆t, the more drastically the function value changes. Figure 4. The traditional element subdivision method: (a) p is located at the corner; (b) p is located on the edge; (c) p is located in the element.
However, the value of the fundamental solutions of elastodynamics is related to the time step. As shown in Figure 5, the curves of the displacement kernel function MM ij U are different under different time steps. In Figure 5, r represents the distance between the source point p and the field point q. It can be seen that the function value changes rapidly when r < c2Δt , in which c2 represents the shear wave velocity and Δt represents the time step. The smaller the time step Δt, the more drastically the function value changes. Therefore, this paper proposes an improved element subdivision scheme related to the time step as follows: (1) Create a square with the source point p as the center and 2 c2Δt as the side length. If the square exceeds the boundary of the element (as shown by the red dashed lines in Figure 6), the boundary of the element is taken as the boundary of the square; (2) According to the position of the source point p and the square area established, the element is subdivided into several triangular patches containing p and quadrilateral patches without p, as shown by the black dashed lines in Figure 6; (3) The α-β transformation method is used for integrals on the triangular patches, and the standard Gauss-Legendre quadrature method is used on the quadrilateral patches. Therefore, this paper proposes an improved element subdivision scheme related to the time step as follows: (1) Create a square with the source point p as the center and 2 c 2 ∆t as the side length. If the square exceeds the boundary of the element (as shown by the red dashed lines in Figure 6), the boundary of the element is taken as the boundary of the square; (2) According to the position of the source point p and the square area established, the element is subdivided into several triangular patches containing p and quadrilateral patches without p, as shown by the black dashed lines in Figure 6; (3) The α-β transformation method is used for integrals on the triangular patches, and the standard Gauss-Legendre quadrature method is used on the quadrilateral patches. Compared with the traditional element subdivision method, the improved element subdivision method related to the time step proposed in this paper can cause the Gauss points to be mainly distributed in the region of r < c2Δt (the region where the kernel function value changes drastically). Therefore, in the case of the same number of Gauss points, the improved subdivision method can significantly improve the accuracy of the singular integral. To verify the computational accuracy of the improved element subdivision method, a numerical example is given next to compare the computational accuracy before and after the improvement.
As shown in Figure 7, the coordinate of the source point p is (0, 0), and the coordinates of the quadrilateral element are shown in the figure. The dotted circles represent the boundaries of the wavefronts. The radius of the oscillation area changes with the time step Δt. Consider the following integral function: Figure 6. The improved element subdivision method: (a) p is located at the corner; (b) p is located on the edge; (c) p is located in the element.

MM
Compared with the traditional element subdivision method, the improved element subdivision method related to the time step proposed in this paper can cause the Gauss points to be mainly distributed in the region of r < c 2 ∆t (the region where the kernel function value changes drastically). Therefore, in the case of the same number of Gauss points, the improved subdivision method can significantly improve the accuracy of the singular integral. To verify the computational accuracy of the improved element subdivision method, a numerical example is given next to compare the computational accuracy before and after the improvement.
As shown in Figure 7, the coordinate of the source point p is (0, 0), and the coordinates of the quadrilateral element are shown in the figure. The dotted circles represent the boundaries of the wavefronts. The radius of the oscillation area changes with the time step ∆t. Consider the following integral function: where U MM ij (p, q) is the displacement kernel function after time integration. The material parameters in the function U MM ij (p, q) are defined as follows: Young's modulus E = 1.1 × 10 5 N/m 2 , density ρ = 2.0 kg/m 3 , and Poisson's ratio ν = 0.0. The average relative error of the integration result is defined as follows: in which I * ij represents the reference result. In this example, the integration results of 50 * 50 Gauss points for each triangle patch and 25 * 25 Gauss points for each quadrilateral patch are used as reference results.
points to be mainly distributed in the region of r < c2Δt (the region where the kernel function value changes drastically). Therefore, in the case of the same number of Gauss points, the improved subdivision method can significantly improve the accuracy of the singular integral. To verify the computational accuracy of the improved element subdivision method, a numerical example is given next to compare the computational accuracy before and after the improvement.
As shown in Figure 7, the coordinate of the source point p is (0, 0), and the coordinates of the quadrilateral element are shown in the figure. The dotted circles represent the boundaries of the wavefronts. The radius of the oscillation area changes with the time step Δt. Consider the following integral function:   Table 1. It can be seen that when the number of Gauss points used is the same, the computation error of the traditional element subdivision method becomes larger and larger as the time step ∆t decreases. However, the computation error of the improved method can be kept at the same order of magnitude. In the case of the same time step, whether it is improved or not, the computation error will decrease as the number of Gauss points increases. In general, the computation error of the improved method is an order of magnitude less than that of the traditional method. Therefore, in the following computation examples of elastodynamic problems, we will use the improved element subdivision method to deal with singular integrals.

Numerical Examples
In this section, two classical elastodynamic problems will be examined to verify the effectiveness of the proposed method. For the first example, the dynamic response of a cantilever beam under a longitudinal load is computed. The influence of the positions of element nodes on the computational stability and accuracy of the result is studied by cases. For the second example, the response of the plate with a square hole under a dynamic load is computed. The results under different inward shrink ratios are compared. According to the results of the two examples, the recommended value range of the inward shrink ratios of element nodes when using partially discontinuous elements is given.

Longitudinal Forced Vibration of a Cantilevered Beam
The geometric model used in this example is shown in Figure 8. The left end of the beam is fixed, and the right end is subjected to a Heaviside-type load p(t) = 1000 H(t) Pa. The material parameters of the beam are Young's modulus E = 1.1 × 10 5 N/m 2 , density ρ = 2.0 kg/m 3 , and Poisson's ratio ν = 0.0.

Numerical Examples
In this section, two classical elastodynamic problems will be examined to verify the effectiveness of the proposed method. For the first example, the dynamic response of a cantilever beam under a longitudinal load is computed. The influence of the positions of element nodes on the computational stability and accuracy of the result is studied by cases. For the second example, the response of the plate with a square hole under a dynamic load is computed. The results under different inward shrink ratios are compared. According to the results of the two examples, the recommended value range of the inward shrink ratios of element nodes when using partially discontinuous elements is given.

Longitudinal Forced Vibration of a Cantilevered Beam
The geometric model used in this example is shown in Figure 8. The left end of the beam is fixed, and the right end is subjected to a Heaviside-type load p(t) = 1000 H(t)Pa. The material parameters of the beam are Young's modulus and respectively, where c = E/ρ is the wave speed. The discretized mesh model is shown in Figure 9a, which includes 72 linear quadrilateral elements and 126 nodes. In this model, the nodes located at the corner points are shrunk to the surfaces in proportion λ = a/b(0 < λ < 0.5) (as shown in Figure 9b), where a represents the moving distance of the node and b represents the edge length of the element.
The computation results of the time-domain BEM are sensitive to the time step. An improper time step length will cause the results to diverge, which is unstable. Therefore, in many papers [24][25][26], the time step is selected according to the dimensionless parameters β = c 1 ∆t/d, where c 1 is the pressure wave velocity, ∆t is the time step, and d is the element edge length. The smaller the value of β is, the more unstable the result will be. The larger the value of β is, the lower the computational accuracy. In general, acceptable results can be obtained when β ≈ 1. In this experiment, we choose the time step that can obtain stable results for computation. On this basis, the influence of the inward shrink ratio of the element node on the computation result is investigated. respectively, where / cE   is the wave speed. The discretized mesh model is shown in Figure 9a, which includes 72 linear quadrilateral elements and 126 nodes. In this model, the nodes located at the corner points are shrunk to the surfaces in proportion   The computation results of the time-domain BEM are sensitive to the time step. An improper time step length will cause the results to diverge, which is unstable. Therefore, in many papers [24][25][26], the time step is selected according to the dimensionless parameters  When ∆t = 0.00445 s(β ≈ 1), the displacement response at the free end and the traction response at the fixed end under different λ values are shown in Figures 10 and 11, respectively. In these two figures, the left picture and the right picture show unstable results and stable results, respectively. It can be seen that the results are unstable when λ = 0.05 and λ = 0.45. The results are stable when 0.1 ≤ λ ≤ 0.4, and the result curves overlap each other in these cases. Therefore, the value of λ should avoid being too large or too small, and the computational accuracy is insensitive to the value of λ if the results converge. . Displacement response at the free end computed by using partially discontinuous elements in terms of (a) the cases of unstable results and (b) the cases of stable results.
(a) (b) Figure 11. Traction response at the fixed end computed by using partially discontinuous elements in terms of (a) the cases of unstable results and (b) the cases of stable results.
Under the same conditions, if completely discontinuous elements are used for discretization, the mesh model contains 72 linear quadrilateral elements and 288 nodes, which are greater values than those using partially discontinuous elements. The displacement response at the free end and the traction response at the fixed end under different λ values are shown in Figures 12 and 13, respectively. It can be seen that the results tend to be stable only when 0.1   , and the results are unstable in other cases. This demonstrates that the performance of the completely discontinuous element is very sensitive to the position of the element node. An inappropriate node inward shrink ratio can easily cause . Displacement response at the free end computed by using partially discontinuous elements in terms of (a) the cases of unstable results and (b) the cases of stable results.
(a) (b) Figure 11. Traction response at the fixed end computed by using partially discontinuous elements in terms of (a) the cases of unstable results and (b) the cases of stable results.
Under the same conditions, if completely discontinuous elements are used for discretization, the mesh model contains 72 linear quadrilateral elements and 288 nodes, which are greater values than those using partially discontinuous elements. The displacement response at the free end and the traction response at the fixed end under different λ values are shown in Figures 12 and 13, respectively. It can be seen that the results tend to be stable only when 0.1   , and the results are unstable in other cases. This demonstrates that the performance of the completely discontinuous element is very sensitive to the position of the element node. An inappropriate node inward shrink ratio can easily cause Figure 11. Traction response at the fixed end computed by using partially discontinuous elements in terms of (a) the cases of unstable results and (b) the cases of stable results.
Under the same conditions, if completely discontinuous elements are used for discretization, the mesh model contains 72 linear quadrilateral elements and 288 nodes, which are greater values than those using partially discontinuous elements. The displacement response at the free end and the traction response at the fixed end under different λ values are shown in Figures 12 and 13, respectively. It can be seen that the results tend to be stable only when λ = 0.1, and the results are unstable in other cases. This demonstrates that the performance of the completely discontinuous element is very sensitive to the position of the element node. An inappropriate node inward shrink ratio can easily cause unstable results, so it is not suitable for the time-domain BEM. Figure 11. Traction response at the fixed end computed by using partially discontinuous elements in terms of (a) the cases of unstable results and (b) the cases of stable results.
Under the same conditions, if completely discontinuous elements are used for discretization, the mesh model contains 72 linear quadrilateral elements and 288 nodes, which are greater values than those using partially discontinuous elements. The displacement response at the free end and the traction response at the fixed end under different λ values are shown in Figures 12 and 13, respectively. It can be seen that the results tend to be stable only when 0.1   , and the results are unstable in other cases. This demonstrates that the performance of the completely discontinuous element is very sensitive to the position of the element node. An inappropriate node inward shrink ratio can easily cause unstable results, so it is not suitable for the time-domain BEM.
(a) (b) Figure 12. Displacement response at the free end computed by using completely discontinuous elements in terms of (a) the cases of unstable results and (b) the cases of stable results. Figure 12. Displacement response at the free end computed by using completely discontinuous elements in terms of (a) the cases of unstable results and (b) the cases of stable results.
Algorithms 2021, 14, x FOR PEER REVIEW 13 of 19 (a) (b) Figure 13. Traction response at the fixed end computed by using completely discontinuous elements in terms of (a) the cases of unstable results and (b) the cases of stable results.
In order to verify the effectiveness of the improved singular element subdivision method proposed in Sub-Section 3.2, we compared the accuracy of the results before and after the improvement. It should be mentioned that the partially discontinuous element is adopted in the following study. Since the singular integral only appears in the first analysis step, we select the results of the first ten analysis steps to clearly observe the difference between the results before and after the improvement. The comparison results under the time steps , and shown in Figure 14a-c, respectively. In the figures, the "Exact" curve represents the exact solution, the "Traditional" curve represents the results before the improvement, and the "Improved" curve represents the results after the improvement. It can be seen that the improved results are in better agreement with the exact solutions. The results before the improvement become unstable as the time step decreases. This demonstrates that the improved singular element subdivision method can reduce the side effect of the time step on the stability of the result and effectively improve the accuracy of the singular integration. However, it should be noted that although the proposed subdivision method has improved the result accuracy in the first ten steps, the computational error becomes larger and larger with the accumulation of time steps, and instability results may still appear. The subdivision method proposed in this paper cannot completely solve the stability problem of the time-domain BEM but can improve the stability of the result. In order to verify the effectiveness of the improved singular element subdivision method proposed in Section 3.2, we compared the accuracy of the results before and after the improvement. It should be mentioned that the partially discontinuous element is adopted in the following study. Since the singular integral only appears in the first analysis step, we select the results of the first ten analysis steps to clearly observe the difference between the results before and after the improvement. The comparison results under the time steps ∆t = 0.00356 s(β ≈ 0.8), ∆t = 0.00178 s(β ≈ 0.4), and ∆t = 0.00089 s(β ≈ 0.2) are shown in Figure 14a-c, respectively. In the figures, the "Exact" curve represents the exact solution, the "Traditional" curve represents the results before the improvement, and the "Improved" curve represents the results after the improvement. It can be seen that the improved results are in better agreement with the exact solutions. The results before the improvement become unstable as the time step decreases. This demonstrates that the improved singular element subdivision method can reduce the side effect of the time step on the stability of the result and effectively improve the accuracy of the singular integration. However, it should be noted that although the proposed subdivision method has improved the result accuracy in the first ten steps, the computational error becomes larger and larger with the accumulation of time steps, and instability results may still appear. The subdivision method proposed in this paper cannot completely solve the stability problem of the time-domain BEM but can improve the stability of the result.
proved singular element subdivision method can reduce the side effect of the time step on the stability of the result and effectively improve the accuracy of the singular integration. However, it should be noted that although the proposed subdivision method has improved the result accuracy in the first ten steps, the computational error becomes larger and larger with the accumulation of time steps, and instability results may still appear. The subdivision method proposed in this paper cannot completely solve the stability problem of the time-domain BEM but can improve the stability of the result. Finally, the condition number of matrix A in Equation (21) with a different inward shrink ratio λ is studied. The condition number of matrix A is computed by Finally, the condition number of matrix A in Equation (21) with a different inward shrink ratio λ is studied. The condition number of matrix A is computed by (26) in which A 2 represents the two-norm of matrix A.
Three common types of discontinuous element are considered, as shown in Figure 15. In the figure, case 1 shows that four nodes are shrunk into the element, case 2 shows that three nodes are shrunk into the element, and case 3 shows that two nodes are shrunk into the element. Case 1 corresponds to the completely discontinuous element, and cases 2 and 3 correspond to the partially discontinuous element. Through computing the condition number of matrix A in Equation (21) with a different inward shrink ratio λ, Figure 16 is plotted for a results comparison. From Figure 16, the following can be seen: (1) When the λ = 0, the condition number is equal to 1. The condition number of A becomes larger as the value of λ increases. In case 1, the matrix A becomes singular when λ = 0.5. The result verifies the conclusion that the value of λ should not be too large.
(2) The condition numbers of A in case 1 are larger than those in the other two cases. This verifies the conclusion that the result stability of the completely discontinuous element (corresponding to case 1) is poorer than that of the partially discontinuous element (corresponding to case 2 and case 3).
Algorithms 2021, 14, x FOR PEER REVIEW 14 of 19 (26) in which 2 A represents the two-norm of matrix A.
Three common types of discontinuous element are considered, as shown in Figure  15. In the figure, case 1 shows that four nodes are shrunk into the element, case 2 shows that three nodes are shrunk into the element, and case 3 shows that two nodes are shrunk into the element. Case 1 corresponds to the completely discontinuous element, and cases 2 and 3 correspond to the partially discontinuous element. Through computing the condition number of matrix A in Equation (21) with a different inward shrink ratio λ, Figure  16 is plotted for a results comparison. From Figure 16, the following can be seen: (1) When the λ = 0, the condition number is equal to 1. The condition number of A becomes larger as the value of λ increases. In case 1, the matrix A becomes singular when λ = 0.5. The result verifies the conclusion that the value of λ should not be too large.
(2) The condition numbers of A in case 1 are larger than those in the other two cases. This verifies the conclusion that the result stability of the completely discontinuous element (corresponding to case 1) is poorer than that of the partially discontinuous element (corresponding to case 2 and case 3).     11 show that when the value of λ is too small, the result will also be unstable. However, from the point of view of the condition number of A, the smaller the value of λ, the smaller the condition number, which is inconsistent with the experimental results. The reason may be that a too small value of λ will reduce the accuracy of the nearly singular integration. In fact, many discrete parameters other than the value of λ, such as the element size, interpolation type, and time step, will affect the stability of the result. Further in-depth research is needed to solve the instability problem of the time-domain BEM.  Figures 10 and 11 show that when the value of λ is too small, the result will also be unstable. However, from the point of view of the condition number of A, the smaller the value of λ, the smaller the condition number, which is inconsistent with the experimental results. The reason may be that a too small value of λ will reduce the accuracy of the nearly singular integration. In fact, many discrete parameters other than the value of λ, such as the element size, interpolation type, and time step, will affect the stability of the result. Further in-depth research is needed to solve the instability problem of the time-domain BEM.

Response of a Plate with a Square Hole under Dynamic Load
In this example, a more realistic problem, which is the dynamic response of a plate with a square hole under an impact load, is computed. The geometric model used in this example is shown in Figure 17. The bottom of the plate is fixed, and the top is subjected to an impact loading p(t) = −8.4 × 10 10 t 2 + 6.3 × 10 6 t Pa. The material parameters of the plate are Young's modulus E = 6.9 × 10 4 N/mm 2 , density ρ = 2.7 × 10 −9 t/mm 3 , and Poisson's ratio ν = 0.3. The dynamic responses at points S1 and S2 in Figure 17a are computed in this example.

Response of a Plate with a Square Hole under Dynamic Load
In this example, a more realistic problem, which is the dynamic response of a plate with a square hole under an impact load, is computed. The geometric model used in this example is shown in Figure 17. The bottom of the plate is fixed, and the top is subjected to an impact loading The discretized mesh model with an element size of 5 mm is shown in Figure 18a. The exact solution of this example cannot be derived, so we adopt the results computed by FEM software to be the reference solution. The FEM mesh model is shown in Figure  18b, which uses a finer mesh model to ensure the reliability of the results. The discretized mesh model with an element size of 5 mm is shown in Figure 18a. The exact solution of this example cannot be derived, so we adopt the results computed by FEM software to be the reference solution. The FEM mesh model is shown in Figure 18b, which uses a finer mesh model to ensure the reliability of the results. The discretized mesh model with an element size of 5 mm is shown in Figure 18a. The exact solution of this example cannot be derived, so we adopt the results computed by FEM software to be the reference solution. The FEM mesh model is shown in Figure  18b, which uses a finer mesh model to ensure the reliability of the results. , the displacement response and the traction response at point S1, which is next to the hole and denoted in Figure 17a, under different λ values are shown in Figures 19 and 20, respectively. It can be seen from these figures that the results are stable when 0.15 0.35

 
, and the result curves overlap each other in these cases. In other cases, the results are unstable. Therefore, a conclusion similar to the previous When ∆t = 1.2 × 10 −6 s(β ≈ 1.4), the displacement response and the traction response at point S1, which is next to the hole and denoted in Figure 17a, under different λ values are shown in Figures 19 and 20, respectively. It can be seen from these figures that the results are stable when 0.15 ≤ λ ≤ 0.35, and the result curves overlap each other in these cases. In other cases, the results are unstable. Therefore, a conclusion similar to the previous example can be obtained: the value of λ should avoid being too large or too small, and the computational accuracy is insensitive to the value of λ if the results converge. In this example, the value range of λ for stable results is smaller than that in the previous example. The reason is that the boundary of the geometric model in this example is more complicated than that in the previous example. Therefore, instability is more likely to occur in this example, and it is more stringent for the selection of the value of λ.
Similarly, the displacement and traction responses at point S2, which is next to the corner of the square hole, are shown in Figures 21 and 22 example can be obtained: the value of λ should avoid being too large or too small, and the computational accuracy is insensitive to the value of λ if the results converge. In this example, the value range of λ for stable results is smaller than that in the previous example. The reason is that the boundary of the geometric model in this example is more complicated than that in the previous example. Therefore, instability is more likely to occur in this example, and it is more stringent for the selection of the value of λ.
(a) (b) Figure 19. Displacement response at point S1 computed by using partially discontinuous elements in terms of (a) the cases of unstable results and (b) the cases of stable results. Figure 19. Displacement response at point S1 computed by using partially discontinuous elements in terms of (a) the cases of unstable results and (b) the cases of stable results.
(a) (b) Figure 19. Displacement response at point S1 computed by using partially discontinuous elements in terms of (a) the cases of unstable results and (b) the cases of stable results.
(a) (b) Figure 20. Traction response at point S1 computed by using partially discontinuous elements in terms of (a) the cases of unstable results and (b) the cases of stable results.
Similarly, the displacement and traction responses at point S2, which is next to the corner of the square hole, are shown in Figures 21 and 22  From the figures of the results comparison in this section, it can be seen that the timedomain BEM can achieve high computational accuracy in the initial stage of computation. However, the numerical damping becomes larger and larger as the number of analysis steps increases, that is, the computational accuracy becomes lower and lower. There are two main reasons for this phenomenon: one is that the large time step adopted in the timedomain BEM computation leads to a lower accuracy of the result, while the result may not converge when the time step is small; the other is that the fundamental solutions of the time-domain BEM for elastodynamic problems are discontinuous piecewise functions, which leads to the element integral accuracy being low. The result error accumulates continuously as the number of analysis steps increases, resulting in lower and lower computational accuracy. The above two points are the difficulties to be overcome when using the From the figures of the results comparison in this section, it can be seen that the timedomain BEM can achieve high computational accuracy in the initial stage of computation. However, the numerical damping becomes larger and larger as the number of analysis steps increases, that is, the computational accuracy becomes lower and lower. There are two main reasons for this phenomenon: one is that the large time step adopted in the timedomain BEM computation leads to a lower accuracy of the result, while the result may not converge when the time step is small; the other is that the fundamental solutions of the time-domain BEM for elastodynamic problems are discontinuous piecewise functions, which leads to the element integral accuracy being low. The result error accumulates continuously as the number of analysis steps increases, resulting in lower and lower computational accuracy. The above two points are the difficulties to be overcome when using the From the figures of the results comparison in this section, it can be seen that the timedomain BEM can achieve high computational accuracy in the initial stage of computation. However, the numerical damping becomes larger and larger as the number of analysis steps increases, that is, the computational accuracy becomes lower and lower. There are two main reasons for this phenomenon: one is that the large time step adopted in the time-domain BEM computation leads to a lower accuracy of the result, while the result may not converge when the time step is small; the other is that the fundamental solutions of the time-domain BEM for elastodynamic problems are discontinuous piecewise functions, which leads to the element integral accuracy being low. The result error accumulates continuously as the number of analysis steps increases, resulting in lower and lower computational accuracy. The above two points are the difficulties to be overcome when using the time-domain BEM. They are also the reasons why the time-domain BEM is more suitable for calculating shortterm dynamic responses such as explosions and shocks rather than long-term dynamic responses. The algorithm proposed in this paper can improve the computational stability and accuracy of the time-domain BEM to a certain extent. However, further research is needed to thoroughly solve the above two problems.

Conclusions
To treat the corners of body surfaces in a simple manner while restricting the number of nodes, we adopted partially discontinuous elements, which are continuous on surfaces and discontinuous between adjacent surfaces, to discretize the boundary of the model for solving the elastodynamic problems. The effects of the numerical integration scheme of the partially discontinuous element and the positions of collocation points on the stability and accuracy of the computational results were studied. We found the following: (1) The computational accuracy of the singular integral is sensitive to the time step. For instance, the traditional singular integration scheme cannot obtain reliable results in the case where the value of the time step adopted is small. The proposed improved singular element subdivision method takes the time step into consideration, which significantly reduces the side effect of the time step on the accuracy of the singular integral. Meanwhile, it will not introduce too many Gauss points' numbers. The proposed singular element subdivision method can improve the stability and accuracy of the dynamic response results in a short time. However, the accumulated error becomes larger and larger as the number of analysis steps increases, which leads to the rapid damping of the amplitude of the result.
(2) The effects of the node position in the completely discontinuous element and the partially discontinuous element on the results were compared. In terms of the stability of the computation results, the partially discontinuous element has obvious advantages over the completely discontinuous element. In terms of computational accuracy, the accuracy of the results of the partially discontinuous elements is more stable than that of the completely discontinuous elements and insensitive to the node position.
(3) When using the discontinuous element, the inward shrink ratio of the element node should be carefully determined. Otherwise, it may cause unstable results. For partially discontinuous elements, the recommended value range of the inward shrink ratio is [0.15, 0.3] in the case of adopting linear quadrilateral elements. It should be noted that the more complex the model boundary, the smaller the available range of the inward shrink ratio.  Data Availability Statement: Not Applicable, the study does not report any data.

Conflicts of Interest:
The authors declare no conflict of interest.