Stability Analysis of an Explicit Integration Algorithm with 3D Viscoelastic Artificial Boundary Elements

Viscoelastic artificial boundary elements are one of the most commonly used artificial boundaries when solving dynamic soilstructure interactions or near-field wave propagation problems. However, due to the lack of clear and practical stability criteria for the explicit algorithm that considers the influence of viscoelastic artificial boundary elements, the determination of the stable time increment in such numerical analyses is still a challenge. In this study, we proposed a numerical stability analysis method for the explicit algorithm with a 3D viscoelastic artificial boundary element based on the idea of a subsystem. )rough this method, the artificial boundary subsystem that controls the stability of the overall numerical system is determined, and the analytical solution for the stability condition of the explicit integration algorithm with 3D viscoelastic artificial boundary elements is obtained. On this basis, the maximum time increment for solving dynamic problems with viscoelastic artificial boundary elements can be determined.


Introduction
During the dynamic analysis of a soil-structure interaction (SSI) system or the numerical simulation of a near-field wave propagation problem, a finite domain is usually intercepted from a semi-infinite medium to establish the near-field model, and the cutoff boundaries should be manipulated to simulate the wave radiation effect. Currently, artificial boundary (also known as an absorbing boundary or a nonreflecting boundary) techniques are the most commonly used methods for wave radiation problems. Typical artificial boundary techniques include the boundary element method (BEM) [1,2], the perfectly matched layers (PML) [3][4][5], the transmitting boundary [6,7], the viscous boundary [8], and the viscoelastic boundary [9][10][11] techniques. Among these methods, the viscoelastic boundary has received much attention because of the advantages of high accuracy and the spatial decoupling characteristics. However, the preprocessing operations of the viscoelastic boundary in the general finite element programs are still slightly complicated, since the equivalent spring-damper systems should be applied on each boundary node of the near-field model. To solve this problem, the viscoelastic artificial boundary element technique is developed [12][13][14], which is accomplished by setting a layer of elements with specified material properties along the normal direction of the cutoff boundaries to work identically to the traditional viscoelastic boundary. rough the viscoelastic artificial boundary element technique, the precision of the original viscoelastic boundary can be maintained, and the preprocessing process can be further simplified; thus, it has gained widespread use in recent years [15][16][17][18][19].
However, when the explicit integration method is adopted for a calculation, the stiffness, viscous damping, and geometrical dimensions of the viscoelastic artificial boundary elements will change the numerical stability of the original computational system. At present, due to the lack of clear and practical stability criteria of the explicit integration algorithm considering the influence of viscoelastic artificial boundary elements, the judgment and selection of the stable time increment in such a numerical analysis is still a challenge. erefore, when viscoelastic artificial boundary elements are adopted in dynamic analysis, unconditionally stable implicit integration algorithms are usually employed to avoid the possible numerical instability caused by viscoelastic artificial boundary elements [20,21]. Due to the iterative solving process in the implicit algorithm, the computational burden is heavy, and the efficiency is poor, especially for large-scale three-dimensional (3D) dynamic SSI problems. erefore, it is necessary to study the stability of the explicit integration algorithm in numerical systems containing viscoelastic artificial boundary elements.
A conventional method for the stability analysis of a numerical system including the artificial boundaries is to establish the transfer matrix of the overall system according to the time-domain integration algorithm, and then the stability condition can be obtained through the eigenvalue analysis of the transfer matrix [22]. However, it is usually difficult to derive the analytical eigenvalue solutions of the large-scale matrix, and obtaining the eigenvalues through numerical methods is a very time-consuming process.
erefore, it is difficult to determine the maximum stable time increment for the explicit calculation through this method. Guan and Liao [23] as well as Liao and Liu [24] proposed a stability analysis method based on a subsystem containing artificial boundaries to study the numerical stability of transmitting boundaries. e subsystems adopted in this method are composed of several rows and lines of nodes along the tangential direction and the normal direction of the transmitting boundary. However, the numbers of degrees of freedom in such subsystems are still relatively large; thus, the stability conditions of the artificial boundaries can only be analyzed through numerical methods. e practical analytical form of the stability criterion, however, has not been given.
In this study, we propose a numerical stability analysis method for explicit integration algorithms based on the idea of subsystems. en, three typical 3D subsystems containing the artificial boundary elements are established, and the spectral radiuses of the transfer matrixes of those subsystems are analyzed to determine the local stability condition. e artificial boundary subsystem that controls the stability of the overall system is determined through a comparison analysis, and the analytical solution for the stability condition of the explicit stepwise integration algorithm considering the 3D viscoelastic artificial boundary elements is proposed.

The Equivalent Physical Parameters of the 3D
Viscoelastic Artificial Boundary Elements e viscoelastic artificial boundary elements are made up of a layer of equivalent elements extending outward along the normal direction at the cutoff boundaries of the near-field model, as shown in Figure 1. e outermost nodes of the artificial boundary elements are fixed, and the material parameters are designated to simulate the wave absorption effect of the viscoelastic artificial boundaries.
We assume that the hexahedron element has a length of 2a, a width of 2b, and a height of h. e expression of the element matrix can be derived from the principle of virtual work [14]: where where E is the elastic modulus, G is the shear modulus, and μ is the Poisson ratio of the artificial boundary element. When h is much smaller than a and b, the stiffness matrix of the equivalent artificial boundary element is approximately e element stiffness matrix formed by the 3D viscoelastic artificial boundary can be expressed as [11,14] [K] � ab 9h Mathematical Problems in Engineering 3 where where K BT and K BN are the normal and tangential stiffness of the spring, respectively; R is the distance from the wave source to the artificial boundary node; G is the shear modulus of the internal medium; and α T and α N are the tangential and normal viscoelastic artificial boundary parameters, whose recommended values are 0.67 and 1.33 for 3D conditions, respectively. To make the equivalent element have the same rigidity as the viscoelastic artificial boundary, the following formula should be established: en, the equivalent elastic modulus and Poisson ratio of the viscoelastic artificial boundary elements can be deduced: e damping matrix is proportional to the stiffness matrix, which is given as follows [14]: en, the damping coefficient can be derived: where ρ, C S , and C P are the mass density, shear wave velocity, and compression wave velocity, respectively, of the internal medium. e mass density ρ of the viscoelastic artificial boundary element is zero.
Although the assumption that the thickness of the boundary element layer h is much smaller than the width L has been introduced during the derivation of the equivalent stiffness matrix of the viscoelastic artificial boundary element, further studies [13,14] have shown that the thickness of the boundary element h is quite robust and can be flexibly selected with little impact on the accuracy. Since the stability of the explicit algorithm is usually controlled by the minimum element size of the model, in practical applications, it is recommended to set the thickness of the boundary element h to be consistent with the internal element size, namely, h � L, to exclude the effect of the boundary element size on the overall stability, as shown in Figure 2.

Stability Analysis Method for the Explicit
Integration Algorithm Based on the Idea of a Subsystem e main purpose of the stability analysis in the explicit integration algorithm is to obtain the discretized time increment Δt, which satisfies the stability requirement in the calculation. For a clear theoretical derivation, in the stability analysis, a homogeneous medium and a uniform discretization mesh are usually employed, as shown in Figure 3, to obtain the general stability criterion. For an inhomogeneous model with nonuniform mesh generation, the stable time increment of the entire computational system can be comprehensively determined according to the smallest mesh size and the material properties.
From the previous studies on the numerical stability of integration algorithms, a common conclusion has been reached that the numerical stability is closely related to the cutoff frequency of the computational system [25] (the highest-order natural frequency of the model). Considering that the mode corresponding to the cutoff frequency is usually presented as the localized staggered vibration of adjacent nodes, it can be deduced that the stability of the overall system can be investigated through the numerical stability analysis of a local subsystem model.
In the following sections, two examples are adopted to illustrate the above deduction. First, we focus on the one-dimensional (1D) wave propagation problem. e discrete model of a 1D shear beam and the highest-order mode shape are shown in Figure 4, where L is the length of the beam element; the red circles represent the finite element nodes; the dashed lines represent the highest-order model; and the yellow squares represent the nodes of the vibration mode, which remain stationary during vibration. In addition, the shear modulus, the mass density, and the cross-sectional area of the beam model are designated as G, ρ, and A, respectively. en, two subsystems are intercepted from the shear beam model for analysis. Subsystem A is composed of two  Mathematical Problems in Engineering half elements between two adjacent nodes of the vibration mode with fixed constraints imposed on both ends, as shown in Figure 5(a). Since the nodes of the vibration mode are fixed, the vibration mode of subsystem A is identical to the localized shape of the highest-order mode of the overall beam model. erefore, the natural frequency of subsystem A is consistent with the highest-order natural frequency of the entire shear beam. Subsystem B is composed of two beam elements, with fixed constraints on the two endpoints. erefore, there is only one degree of freedom in subsystem B, as shown in Figure 5(b). Because of the different vibration modes between subsystem B and the overall beam model, the natural frequency of subsystem B is different from the cutoff frequency of the original system.
Since there is a classical solution for the stability condition in a 1D undamped linear elastic system, to compare the stability conditions of the subsystems with those of the overall model, the damping of the subsystem is not considered here. According to the concept of the shear beam, the shear stiffness k 1D A , mass m 1D A , and natural frequency ω 1D A of the 1D subsystem A can be obtained as follows: where C S is the velocity of the shear wave and C S � ��� G/ρ. In addition, according to the shear stiffness and mass of subsystem A, the single degree of freedom motion equation can be established. When the central difference method is adopted for the time-domain stepwise integration calculation, the following condition should be met to ensure the stability of the explicit calculation: By substituting (11) into (12), the numerical stability condition of subsystem A can be obtained as follows: Equation (13) is the exact stability condition of the central difference algorithm in the 1D wave propagation problem, indicating that the stability condition of the numerical integration algorithm in a 1D continuous medium can be obtained through the stability analysis of an A-type subsystem.
rough similar analysis procedures, the shear stiffness k 1D B , mass m 1D B , natural frequency ω 1D B , and stability condition in the center differential method of subsystem B can be expressed as follows: A comparison between (13) and (14) shows that the stability condition of subsystem B is loose compared to that of subsystem A, which means that the stability condition of subsystem B provides an upper limit estimation of the maximum stable time increment of the entire calculation system. en, the previous deduction in the two-dimensional (2D) condition is verified. e 2D plane model and the corresponding highest-order mode are shown in Figure 6, in which the 4-node isoparametric element with a side length of L is adopted for model discretization. e mass density and the elastic modulus are designated as ρ and E. ere is a classical solution of the stability condition in the 2D undamped linear elastic system, namely, Δt ≤ (L/C P ) [25], where C P is the velocity of the compression wave. e premise of the above classical solution is that the Poisson ratio is 0, leading to C P � ��� E/ρ. To compare the numerical stability conditions of the subsystems with the classical solution, the Poisson ratio of the system is set to zero. Similarly, the red circles represent the finite element nodes, and the yellow squares represent the nodes of the vibration mode.

Internal medium elements
Viscoelastic artificial boundary elements  e normal constraints are applied on the four side nodes (nodes 2-5), which can accurately simulate the highestorder vibration mode shape. e above boundary conditions of subsystem A can be expressed as follows: where u is the nodal displacement, the subscripts x and y denote the x-and y-directions, respectively, and the subscripts 1-6 denote the node numbers in Figure 7. e motion equation for the 2D subsystem of subsystem A can be expressed as follows: en, the natural frequency and the stability condition in the central differential method of the 2D subsystem of subsystem A can be calculated as Equation (17) is the exact stability condition of the central difference algorithm for 2D wave problems, indicating that the stability condition of the numerical integration algorithm in 2D continuous medium can be obtained through the stability analysis of the 2D subsystem of subsystem A. e other subsystem shown in Figure 7(b) (subsystem B) consists of four adjacent elements, with fixed constraints imposed on all the outer nodes (nodes 2-9). en, the motion equation of subsystem B can be expressed as en, the natural frequency and the stability condition in the central differential method of subsystem B can be calculated as It can be seen from the comparison between (17) and (19) that, for 2D conditions, the maximum stable incremental time step of subsystem B also provides an upper limit estimation of the numerical stability condition of the entire calculation system. e results of the 1D and 2D subsystem analyses show that (1) if the cutoff frequency of the entire system and the corresponding vibration mode can be accurately determined, the partial subsystem can be intercepted from the entire system according to the characteristics of the highestorder mode with specified constraints applied on the boundary nodes of the subsystem to simulate the shape of the highest-order vibration mode. en, the maximum stable time increment of the entire system can be obtained through the stability analysis of the subsystem. (2) If it is difficult to determine the highest-order natural frequency and the corresponding mode of the entire system, a  Mathematical Problems in Engineering subsystem composed of all the elements connected to a single node can be intercepted to perform the stability analysis, whose maximum time increment provides an upper approximation of the stability condition of the overall system. When the viscoelastic artificial boundary elements are applied in the computational system, the numerically stable condition of the internal domain can be obtained. e cutoff frequency in the artificial boundary region and the corresponding vibration mode, however, is difficult to determine. In this situation, a B-type subsystem can be employed for the numerical stability analysis.
Furthermore, for the stability analysis of an inhomogeneous subsystem, the motion equation can be written in the following form according to the explicit stepwise integration algorithm:r where the vector U consists of the nodal acceleration, velocity, and displacement of the subsystem; Q is the external force vector; A is the transfer matrix, which is related to the discrete steps in the space domain and time domain, as well as the material parameters of the subsystem; and the superscript p denotes the natural number, t � p · Δt. e following conditions should be met to ensure the stability of the numerical integration algorithm [26]: (1) If transfer matrix A has no repeated eigenvalues, then ρ(A) ≤ 1 (ρ(A) � max|λ i |), where ρ(A) is the spectral radius and λ i denotes the ith eigenvalue of transfer matrix A (2) If transfer matrix A has repeated eigenvalues, the moduli of the repeated eigenvalues should be less than 1 erefore, the numerical stability analysis of the subsystem is based on the generation of the transfer matrix and the calculation of the spectral radius.

Stability Analysis of 3D Artificial Boundary Subsystems
When the 3D viscoelastic artificial boundary elements are employed in the numerical models, the B-type subsystem can be classified into three types according to the interception region, as shown in Figure 8: (a) the surface boundary subsystem, which is intercepted from the lateral or bottom surface of the model and consists of four artificial boundary elements and four internal elements; (b) the edge boundary subsystem, which is intercepted from the lateral or bottom edge regions and consists of six artificial boundary elements and two internal elements; and (c) the corner boundary subsystem is intercepted from the corner region and consists of seven artificial boundary elements and one internal element. e meshes in the model are uniformly generated, and the side length is L. e boundary condition of these three kinds of subsystems is obtained by fixing the 26 surface nodes of the system and letting the remaining central node be free. e subsystems with all fixed surface nodes may have a higher natural vibration frequency, which can ensure that the numerical stability conditions of the subsystem are as close as possible to the overall system. In addition, there is only one free node in such a subsystem, which is helpful for deriving the analytical solution of the stability condition.
For the internal medium, the elastic modulus, shear modulus, mass density, and Poisson ratio are designated E, G, ρ, and μ. Damping is not considered in the internal region. en, the stability condition of the internal domain can be calculated by [25] Δt int ernal ≤ L C P .
For the viscoelastic artificial boundary elements, the elastic modulus E, the damping coefficient η, the Poisson ratio μ, and the mass density ρ can be determined according to (8) and (10). Mathematical Problems in Engineering en, each type of subsystem is focused on, and the corresponding stability analysis is conducted in the following sections.

Stability Analysis of the Surface Boundary Subsystem.
Consider a regular hexahedral element from the surface boundary subsystem (Figure 8(a)), as shown in Figure 9.
(28) e eigenvalues of transfer matrix A are By substituting (29) into the stability condition ρ(A) � max|λ i | ≤ 1, the following expression can be obtained: Equation (30) provides the analytical solution for the stability condition of the surface boundary subsystem. e stability condition is related to the wave velocities C P and C S , the mass density ρ, the Poisson ratio μ of the internal media, and the element size L.

Stability Analysis of the Edge Boundary Subsystem.
For the edge boundary subsystem, following the same procedure of matrix assembly, the motion equation of the central node can be expressed as follows: It can be seen from (31) that the motion equations in different directions are coupled with each other, and the stability analysis of the overall subsystem is required. By combining (26) and (31), the following expression can be obtained: where , Mathematical Problems in Engineering e maximum eigenvalue of the transfer matrix A can be solved as follows: (34) By substituting (34) into the stability condition ρ(A) � max|λ i | ≤ 1, the analytical solution of the stability condition of the edge boundary subsystem can be obtained:

Stability Analysis of the Corner Boundary Subsystem.
For the corner boundary subsystem, following the same procedure matrix assembly, the motion equation of the central node can be expressed as follows: (ρ + 7ρ)L 3 8 e motion equations of the corner boundary subsystem in different directions are also coupled with each other. By combining (26) and (36), the following expression can be obtained: where e maximum eigenvalue of the transfer matrix A can be solved as follows: By substituting (39) into the stability condition ρ(A) � max|λ i | ≤ 1, the analytical solution of the stability condition of the corner boundary subsystem can be obtained:

Validation
A 3D homogeneous elastic half space model is established, as shown in Figure 10. e dimensions of the model are 100 m × 100 m × 50 m, and the mesh size is 2 m × 2 m × 2 m. e mass density of the internal medium ρ is 2500 kg/m 3 , the shear wave velocity C S is 500 m/s, the compression wave velocity C P is 935 m/s, and the Poisson ratio is 0.3. e viscoelastic artificial boundary elements are applied on the cutoff boundaries to absorb the outgoing waves, and their material parameters are calculated through (8) and (10). A concentrated impulsive force is imposed on the geometric center of the model along the z-direction, and its time history is shown in Figure 11.
According to (21), (30), (35), and (40), the maximum stable time increments for the internal region Δt internal , the surface boundary subsystem Δt surface , the edge boundary subsystem Δt edge , and the corner boundary subsystem Δt corner in this model can be calculated, whose values are 0.0021 s, 0.00176 s, 0.00071 s, and 0.00032 s, respectively. e central difference method is adopted, and the accuracy of the above stability analysis can be verified by specifying different fixed time increments and conducting explicit dynamic analyses.
It is important to note that the minimum distances from the geometric center of the model to the surfaces L S , to the edges L E , and to the corners L C gradually increase (L S < L E < L C ); thus, the time for the wave crest arriving at the surfaces, the edges, and the corners also increases one by one.
First, the maximum stable time increment of the internal domain is adopted for the explicit analysis (Δt � Δt internal � 0.0021 s). When the pulse wave propagates within the internal region, the calculation can be executed accurately. However, once the waves reach the bottom surface of the model, the stability condition of the surface boundary subsystem is no longer satisfied (Δt > Δt surface ); thus, instability appears near the center of the bottom surface, as shown in Figure 12(a).
As we reduce the time increment to the maximum stable time increment of the surface boundary subsystem (Δt � Δt surface � 0.00176 s), the process of wave propagation in the internal region and around the bottom surface can be precisely simulated. However, computational instability occurs when the waves reach the bottom edges, as shown in Figure 12(b), since the stability condition of the edge boundary subsystem is not met (Δt > Δt edge ). en, the time increment is further reduced to the maximum stable time increment of the edge boundary subsystem (Δt � Δt edge � 0.00071 s). is time the simulation of the wave propagation can be accurately conducted until the waves arrive at the bottom corners of the model. e instability appears in the corner point, as shown in Figure 12(c), on account of the unsatisfied stability condition of the corner boundary subsystem (Δt > Δt corner ).
Finally, the maximum stable time increment of the corner boundary subsystem is adopted for the explicit analysis (Δt � Δt corner � 0.00032 s). e wave propagation process in the calculation model is shown in Figure 13, and the vertical displacements of observation points A, B, C, and D are shown in Figure 14. It can be seen that the dynamic calculation of the entire model can be successfully executed, and the outgoing waves are effectively absorbed by the artificial boundaries. e above dynamic analyses verify the precision of the stability conditions proposed in this study. From the comparison between the stability conditions of the internal domain and different substructures, it can be deduced that the numerical stability of the computational system is controlled by the corner region of the model. To validate this deduction, we successively change the important model parameters that would affect the numerical stability, namely, continuously changing the mass density of the internal medium ρ from 1000 kg/m 3 to 2500 kg/m 3 , the velocity of the shear wave of the internal medium C S from 300 m/s to 1000 m/s, the element size L from 1 m to 10 m, the model size     R from 1 m to 1000 m, and the Poisson ratio of the internal medium μ from 0.3 to 0.5, and the maximum time increments are investigated under different model parameters. e results are shown in Figure 15. Under different model parameters, the corner boundary subsystem provides the strictest stability condition, followed by the edge boundary subsystem. In addition, the surface boundary subsystem provides stricter stability conditions than the internal domain under different mass densities, velocities of shear waves, element sizes, and model sizes. However, when the Poisson ratio μ is relatively large, the stability condition of the internal medium is stricter than that of the surface subsystem. When the shear wave velocity C S is constant, the compression wave velocity C P can be expressed as follows: When the Poisson ratio μ increases, the compression wave velocity of the internal medium C P increases rapidly. According to (21), a larger compression wave velocity will lead to stricter stability conditions in the internal medium.
From the above analysis, it can be demonstrated that the numerical stability of the system containing the viscoelastic artificial boundary elements is controlled by the corner region of the model. e stability condition of the corner subsystem given by (40) is a sufficient condition for the numerical stability of the entire model.

Conclusion
In this study, a numerical stability analysis method for the explicit algorithm based on the idea of a subsystem is proposed.
rough this method, the artificial boundary subsystem that controls the stability of the overall system is determined, and the analytical solution for the stability condition of the explicit integration algorithm considering the 3D viscoelastic artificial boundary elements is derived. Based on the results, the following conclusions can be obtained: (1) For the stability problem of large-scale numerical calculations, the stability analysis of the explicit timedomain stepwise integration algorithm can be performed on the local subsystem. When the internal medium satisfies the stability condition of the algorithm, the numerical stability of the overall system is controlled by the corner region. is stability condition is a sufficient condition for the stability of the overall system, not a necessary and sufficient condition. However, the results of the example verification show that the difference between the sufficient condition and the necessary and sufficient condition is very small. us, the sufficient condition can basically meet the needs of the stability analysis in practical applications.
(2) When the viscoelastic artificial boundary elements are adopted, the stability of the overall system is different from the internal medium stability due to the influence of the quality, stiffness, and damping of the artificial boundary elements. e former has more stringent numerical stability conditions and requires smaller integration steps to meet overall system stability conditions compared to that of the latter.
(3) e stability conditions of the local subsystem are derived based on the format of the central difference, which is commonly used in dynamic numerical analysis. It is important to note that the stability condition is closely related to the integral format. If another time domain step integration format is applied, the method proposed in this paper can also be used by intercepting the representative subsystems and deriving the corresponding transfer matrix and spectral radius. en, the stability conditions based on the other integral format can be obtained according to a similar process introduced in this paper. (4) e stability conditions given in this paper are derived from regular hexahedron elements and are equally applicable to rectangular hexahedron elements. Since the stability condition is affected by the minimum size of the elements, the minimum side length of the rectangular elements can be used as a parameter to derive the stability condition.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest regarding the publication of this paper.