^{1}

^{2}

^{1}

^{1}

^{2}

^{2}

^{2}

^{2}

^{2}

^{1}

^{2}

Viscoelastic artificial boundary elements are one of the most commonly used artificial boundaries when solving dynamic soil-structure 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. Through 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.

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) [

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. Therefore, 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 [

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 [

In this study, we propose a numerical stability analysis method for explicit integration algorithms based on the idea of subsystems. Then, 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. The 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 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

3D viscoelastic artificial boundary element (hexahedron).

We assume that the hexahedron element has a length of 2

The element stiffness matrix formed by the 3D viscoelastic artificial boundary can be expressed as [

Then, the equivalent elastic modulus and Poisson ratio of the viscoelastic artificial boundary elements can be deduced:

The damping matrix is proportional to the stiffness matrix, which is given as follows [

Then, the damping coefficient can be derived:

Although the assumption that the thickness of the boundary element layer

3D viscoelastic artificial boundary elements.

The main purpose of the stability analysis in the explicit integration algorithm is to obtain the discretized time increment

Homogeneous model with a uniform discretization mesh.

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 [

In the following sections, two examples are adopted to illustrate the above deduction. First, we focus on the one-dimensional (1D) wave propagation problem. The discrete model of a 1D shear beam and the highest-order mode shape are shown in Figure

1D shear beam model and the highest-order vibration mode.

Then, two subsystems are intercepted from the shear beam model for analysis. Subsystem A is composed of two half elements between two adjacent nodes of the vibration mode with fixed constraints imposed on both ends, as shown in Figure

1D local subsystems: (a) subsystem A and (b) subsystem B.

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 _{S} is the velocity of the shear wave and

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 (

Equation (

Through similar analysis procedures, the shear stiffness

A comparison between (

Then, the previous deduction in the two-dimensional (2D) condition is verified. The 2D plane model and the corresponding highest-order mode are shown in Figure

The highest-order vibration mode of the 2D node system.

Two 2D subsystems are intercepted from the original computational model, as shown in Figure

2D local subsystems: (a) subsystem A and (b) subsystem B.

The motion equation for the 2D subsystem of subsystem A can be expressed as follows:

Then, the natural frequency and the stability condition in the central differential method of the 2D subsystem of subsystem A can be calculated as

Equation (

The other subsystem shown in Figure

Then, 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 (

The 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 highest-order mode with specified constraints applied on the boundary nodes of the subsystem to simulate the shape of the highest-order vibration mode. Then, 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 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. The 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

The following conditions should be met to ensure the stability of the numerical integration algorithm [

If transfer matrix

If transfer matrix

Therefore, the numerical stability analysis of the subsystem is based on the generation of the transfer matrix and the calculation of the spectral radius.

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

Different types of 3D artificial boundary subsystems. (a) Surface boundary subsystem. (b) Edge boundary subsystem. (c) Corner boundary subsystem.

The 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. The 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

For the viscoelastic artificial boundary elements, the elastic modulus

Then, each type of subsystem is focused on, and the corresponding stability analysis is conducted in the following sections.

Consider a regular hexahedral element from the surface boundary subsystem (Figure

Regular hexahedral element.

The concentrated mass matrix of the regular hexahedral element of the internal medium can be expressed as

The stiffness matrix of the regular hexahedral element of the internal medium can be expressed as

Node 1:

The internal medium element contains the mass matrix and the stiffness matrix, while the artificial boundary element contains the damping matrix in addition to the mass matrix and stiffness matrix. The damping matrix can be obtained by (

It can be seen from (

By combining (

The eigenvalues of transfer matrix

By substituting (

Equation (_{P} and _{S}, the mass density

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 (

The maximum eigenvalue of the transfer matrix

By substituting (

For the corner boundary subsystem, following the same procedure matrix assembly, the motion equation of the central node can be expressed as follows:

The motion equations of the corner boundary subsystem in different directions are also coupled with each other. By combining (

The maximum eigenvalue of the transfer matrix

By substituting (

A 3D homogeneous elastic half space model is established, as shown in Figure ^{3}, the shear wave velocity _{S} is 500 m/s, the compression wave velocity _{P} is 935 m/s, and the Poisson ratio is 0.3. The viscoelastic artificial boundary elements are applied on the cutoff boundaries to absorb the outgoing waves, and their material parameters are calculated through (

3D homogeneous model and viscoelastic artificial boundary elements. (a) 3D homogeneous model. (b) Viscoelastic artificial boundary elements.

The time history of the concentrated impulsive force.

According to (

It is important to note that the minimum distances from the geometric center of the model to the surfaces _{S}, to the edges _{E}, and to the corners _{C} gradually increase (_{S} < _{E} < _{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 (

Three kinds of instability states with time increments of (a) 0.0021 s, (b) 0.00176 s, and (c) 0.00071 s.

As we reduce the time increment to the maximum stable time increment of the surface boundary subsystem (

Then, the time increment is further reduced to the maximum stable time increment of the edge boundary subsystem (

Finally, the maximum stable time increment of the corner boundary subsystem is adopted for the explicit analysis (

Wave propagation calculated according to the stability conditions of the corner subsystem. (a) Time 0.05 s. (b) Time 0.10 s. (c) Time 0.15 s. (d) Time 0.20 s. (e) Time 0.25 s. (f) Time 0.30 s.

Vertical displacements of observation points A, B, C, and D.

The 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 ^{3} to 2500 kg/m^{3}, the velocity of the shear wave of the internal medium _{S} from 300 m/s to 1000 m/s, the element size

The maximum time increments under different (a) mass densities, (b) velocities of shear waves, (c) element sizes, (d) model sizes, and (e) Poisson ratios.

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

When the Poisson ratio

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. The stability condition of the corner subsystem given by (

In this study, a numerical stability analysis method for the explicit algorithm based on the idea of a subsystem is proposed. Through 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:

For the stability problem of large-scale numerical calculations, the stability analysis of the explicit time-domain 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. This 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. Thus, the sufficient condition can basically meet the needs of the stability analysis in practical applications.

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. The former has more stringent numerical stability conditions and requires smaller integration steps to meet overall system stability conditions compared to that of the latter.

The 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. Then, the stability conditions based on the other integral format can be obtained according to a similar process introduced in this paper.

The 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.

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

The authors declare that they have no conflicts of interest regarding the publication of this paper.

This study was supported by the National Natural Science Foundation of China (nos. 51878384 and U1839201) and the National Key Research and Development Program of China (no. 2018YFC1504305). Financial support from these organizations is gratefully acknowledged.