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

Authors: Shutao Li; Jingbo Liu; Xin Bao; Yifan Jia; Lan Xiao; Xigui Tao; Yicun Chen; Weifeng Zhang
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

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.


---

## Body

## 1. 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–5], the transmitting boundary [6, 7], the viscous boundary [8], and the viscoelastic boundary [9–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–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. Through 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–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. 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 [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. Therefore, 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. Therefore, 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. The 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. The 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. 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.


## 2. The Equivalent Physical Parameters of the 3D Viscoelastic Artificial Boundary Elements

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 Figure1. The 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.Figure 1
3D viscoelastic artificial boundary element (hexahedron).We assume that the hexahedron element has a length of 2a, a width of 2b, and a height of h. The expression of the element matrix can be derived from the principle of virtual work [14]:(1)K˜=∫BTDBdΩ,where K˜ is the stiffness matrix, D is the elastic modulus matrix, and B is the geometric matrix of the equivalent element. By fixing nodes 5–8 of the element, the stiffness matrix of the equivalent element can be obtained:(2)K˜1=29ahG˜2b+2abG˜h+bhp2ahq12+hG˜12bq12−bG˜12hq12+hG˜1229ahG˜2b+2abG˜h+bhpaaq12+aG˜12bq12−bG˜12aq12+aG˜1229ahG˜2b+2abph+bhG˜2a,where(3)p=E˜1−μ˜1+μ˜1−2μ˜,q=E˜μ˜1+μ˜1−2μ˜,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(4)K˜=ab9h4G˜002G˜00G˜002G˜004G˜002G˜00G˜002G˜04m002m00m002m4G˜002G˜00G˜004G˜002G˜00G˜0sym.4m002m00m4G˜002G˜00sym.4G˜002G˜04m002msym.4G˜004G˜04m.The element stiffness matrix formed by the 3D viscoelastic artificial boundary can be expressed as [11, 14](5)K=ab9h4KBT002KBT00KBT002KBT004KBT002KBT00KBT002KBT04KBN002KBN00KBN002KBN4KBT002KBT00KBT004KBT002KBT00KBT0sym.4KBN002KBN00KBN4KBT002KBT00sym.4KBT002KBT04KBN002KBNsym.4KBT004KBT04KBN,where(6)KBT=aTGR,KBN=aNGR,where KBT and KBN 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:(7)abhG˜=abKBT,abhm=abKBN.Then, the equivalent elastic modulus and Poisson ratio of the viscoelastic artificial boundary elements can be deduced:(8)E˜=αNhGR⋅1+μ˜1−2μ˜1−μ˜,μ˜=α−22α−1,α=αNαT.The damping matrix is proportional to the stiffness matrix, which is given as follows [14]:(9)C=η˜K.Then, the damping coefficient can be derived:(10)η˜=ρR3G2CSαT+CPαN,where ρ, CS, and CP are the mass density, shear wave velocity, and compression wave velocity, respectively, of the internal medium. The mass density ρ˜ of the viscoelastic artificial boundary element is zero.Although the assumption that the thickness of the boundary element layerh 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.Figure 2
3D viscoelastic artificial boundary elements.


## 3. Stability Analysis Method for the Explicit Integration Algorithm Based on the Idea of a Subsystem

The 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.Figure 3
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 [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. The discrete model of a 1D shear beam and the highest-order mode shape are shown in Figure4, 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.Figure 4
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 Figure5(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. Therefore, 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. Therefore, 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.Figure 5
1D local subsystems: (a) subsystem A and (b) subsystem B.
(a)(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 stiffnesskA1D, mass mA1D, and natural frequency ωA1D of the 1D subsystem A can be obtained as follows:(11)kA1D=4GAL,mA1D=ρAL,ωA1D=kA1DmA1D=2CSL,where CS is the velocity of the shear wave and CS=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:(12)ΔtA1D≤2ωA1D.By substituting (11) into (12), the numerical stability condition of subsystem A can be obtained as follows:(13)ΔtA1D≤LCS.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.Through similar analysis procedures, the shear stiffnesskB1D, mass mB1D, natural frequency ωB1D, and stability condition in the center differential method of subsystem B can be expressed as follows:(14)kB1D=2GAL,mB1D=ρAL,ωB1D=2CSL,ΔtB1D≤2⋅LCS.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.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 Figure6, in which the 4-node isoparametric element with a side length of L is adopted for model discretization. The mass density and the elastic modulus are designated as ρ and E. There is a classical solution of the stability condition in the 2D undamped linear elastic system, namely, Δt≤L/CP [25], where CP is the velocity of the compression wave. The premise of the above classical solution is that the Poisson ratio is 0, leading to CP=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.Figure 6
The highest-order vibration mode of the 2D node system.Two 2D subsystems are intercepted from the original computational model, as shown in Figure7. Subsystem A is composed of four quarters of an element. The four corner nodes of subsystem A (nodes 6–9) are the nodes of the vibration mode and remain stationary during the analysis. The normal constraints are applied on the four side nodes (nodes 2–5), which can accurately simulate the highest-order vibration mode shape. The above boundary conditions of subsystem A can be expressed as follows:(15)ux1=ux4=ux5,uy1=uy2=uy3,ux2=ux3=ux6=ux7=ux8=ux9=0,uy4=uy5=uy6=uy7=uy8=uy9=0,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.Figure 7
2D local subsystems: (a) subsystem A and (b) subsystem B.
(a)(b)The motion equation for the 2D subsystem of subsystem A can be expressed as follows:(16)ρL21001⋅u¨x1u¨y1+E4004⋅ux1uy1=0.Then, the natural frequency and the stability condition in the central differential method of the 2D subsystem of subsystem A can be calculated as(17)ωA2D=2CPL,ΔtA2D≤2ωA2D=LCP.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.The other subsystem shown in Figure7(b) (subsystem B) consists of four adjacent elements, with fixed constraints imposed on all the outer nodes (nodes 2–9). Then, the motion equation of subsystem B can be expressed as(18)ρL21001⋅u¨x1u¨y1+E2002⋅ux1uy1=0.Then, the natural frequency and the stability condition in the central differential method of subsystem B can be calculated as(19)ωB2D=2CPL,Δt≤2ωB2D=2⋅LCP.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.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(20)Up+1=A⋅Up+Qp,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.The following conditions should be met to ensure the stability of the numerical integration algorithm [26]:(1)
If transfer matrixA 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 matrixA has repeated eigenvalues, the moduli of the repeated eigenvalues should be less than 1Therefore, the numerical stability analysis of the subsystem is based on the generation of the transfer matrix and the calculation of the spectral radius.


## 4. 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 Figure8: (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. The meshes in the model are uniformly generated, and the side length is L.Figure 8
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 designatedE, G, ρ, and μ. Damping is not considered in the internal region. Then, the stability condition of the internal domain can be calculated by [25](21)Δtint ernal≤LCP.For the viscoelastic artificial boundary elements, the elastic modulusE˜, the damping coefficient η˜, the Poisson ratio μ˜, and the mass density ρ˜ can be determined according to (8) and (10).Then, each type of subsystem is focused on, and the corresponding stability analysis is conducted in the following sections.


### 4.1. Stability Analysis of the Surface Boundary Subsystem

Consider a regular hexahedral element from the surface boundary subsystem (Figure8(a)), as shown in Figure 9.Figure 9
Regular hexahedral element.The concentrated mass matrix of the regular hexahedral element of the internal medium can be expressed as(22)Me=ρL318⋯00180⋱⋮sym.⋱sym.181824×24.The stiffness matrix of the regular hexahedral element of the internal medium can be expressed as(23)Ke=k11k12⋯k18k21⋱⋮⋮⋱⋮k81⋯⋯k8824×24,kij=EL3161+μ1−2μa11a12a13a21a22a23a31a32a33,where(24)a11=1−μL2⋅ξiξj1+ηiηj31+ζiζj3+1−2μ2L2ηiηj1+ζiζj31+ξiξj3+ζiζj1+ξiξj31+ηiηj3,a22=1−μL2⋅ηiηj1+ζiζj31+ξiξj3+1−2μ2L2ξiξj1+ηiηj31+ζiζj3+ζiζj1+ξiξj31+ηiηj3,a33=1−μL2⋅ζiζj1+ξiξj31+ηiηj3+1−2μ2L2ξiξj1+ηiηj31+ζiζj3+ηiηj1+ξiξj31+ζiζj3,a12=1L21+ζiζj3μξiηj+1−2μ2ηiξj,a13=1L21+ηiηj3μξiζj+1−2μ2ζiξj,a21=1L21+ζiζj3μηiξj+1−2μ2ξiηj,a23=1L21+ξiξj3μηiζj+1−2μ2ζiηj,a31=1L21+ηiηj3μζiξj+1−2μ2ξiζj,a32=1L21+ξiξj3μζiηj+1−2μ2ηiζj.Node 1:ξ,η,ζ1=1,1,−1; Node 2: ξ,η,ζ2=−1,1,−1; Node 3: ξ,η,ζ3=−1,−1,−1; Node 4: ξ,η,ζ4=1,−1,−1; Node 5: ξ,η,ζ5=1,1,1; Node 6: ξ,η,ζ6=−1,1,1; Node 7: ξ,η,ζ7=−1,−1,1; and Node 8: ξ,η,ζ8=1,−1,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 (9) and (10). According to this method, the mass, stiffness, and damping matrix of a single internal medium element and artificial boundary element can be obtained. Then, the eight elements in the subsystem can be assembled according to the node number (the node number is shown in Figure 9). The fixed boundary conditions are applied to the surface nodes of the subsystem. Then, the motion equation of the central node can be expressed as follows:(25)ρ+ρ˜L32100010001u¨xu¨yu¨z+8ρL2aN272CSaT+CPaN100010001u˙xu˙yu˙z+ρCS2L32−48μ181−2μ+aN8ρCS2L29R100010001uxuyuz=0,where the subscripts x, y, and z denote the x-, y-, and z-directions, respectively.It can be seen from (25) that the motion equations in the x-, y-, and z-directions are decoupled and identical to each other. Therefore, the stability analysis can be performed in only one direction. Taking the x-direction as an example, the explicit time-domain stepwise integration format based on the central difference method can be expressed as follows [27]:(26)u˙xt+Δt2=u˙xt−Δt2+Δtu¨xt,uxt+Δt=uxt+Δtu˙xt+Δt2.By combining (25) and (26), the following expression can be obtained:(27)uxt+ΔtΔtu˙xt+Δt2=AuxtΔtu˙xt−Δt2,where(28)A=1−kΔt2m1−cΔtm−kΔt2m1−cΔtm,m=ρ+ρ˜L32,c=8ρL2aN272CSaT+CPaN,k=ρCS2L32−48μ181−2μ+aN8ρCS2L29R.The eigenvalues of transfer matrixA are(29)λ1,2=1−kΔt22m−cΔt2m±Δt2m−4km+k2Δt2+2ckΔt+c2.By substituting (29) into the stability condition ρA=maxλi≤1, the following expression can be obtained:(30)Δtsurface≤1k4km+c2−c.Equation (30) provides the analytical solution for the stability condition of the surface boundary subsystem. The stability condition is related to the wave velocities CP and CS, the mass density ρ, the Poisson ratio μ of the internal media, and the element size L.


### 4.2. 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:(31)ρ+3ρ˜L34100010001u¨xu¨yu¨z+ρL2aN362CSaT+CPaN16−10−11600016u˙xu˙yu˙z+13ρCS2L32−48μ121−2μ+aN4LR121−2μ−aNL4R0121−2μ−aNL4R32−48μ121−2μ+aN4LR00032−48μ121−2μ+aN4LRuxuyuz=0.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:(32)uxt+Δtuyt+Δtuzt+ΔtΔtu˙xt+Δt2Δtu˙yt+Δt2Δtu˙zt+Δt2=A⋅uxtuytuztΔtu˙xt−Δt2Δtu˙yt−Δt2Δtu˙zt−Δt2,where(33)A=1−k1Δt2m−k2Δt2m01−16cΔtmcΔtm0−k2Δt2m1−k1Δt2m0cΔtm1−16cΔtm0001−k1Δt2m001−16cΔtm−k1Δt2m−k2Δt2m01−16cΔtmcΔtm0−k2Δt2m−k1Δt2m0cΔtm1−16cΔtm000−k1Δt2m001−16cΔtm,m=ρ+3ρ˜L34,c=ρL2aN362CSaT+CPaN,k1=ρCS2L32−48μ361−2μ+4aNρCS2L23R,k2=ρCS2L61−2μ−aNρCS2L212R.The maximum eigenvalue of the transfer matrixA can be solved as follows:(34)maxλi=122−k1Δt2m+k2Δt2m−17cΔtm−−41−17cΔtm+−2+k1Δt2m−k2Δt2m+17cΔtm2.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:(35)Δtedge≤4mk1−k2+289c2−17ck1−k2.


### 4.3. 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:(36)ρ+7ρ˜L38100010001u¨xu¨yu¨z+ρL2aN722CSaT+CPaN1123−11−111231111123u˙xu˙yu˙z+13ρCS2L32−48μ241−2μ+aN14L3R141−2μ−aNL8R−141−2μ+aNL8R141−2μ−aNL8R32−48μ241−2μ+aN14L3R−141−2μ+aNL8R−141−2μ+aNL8R−141−2μ+aNL8R32−48μ241−2μ+aN14L3Ruxuyuz=0.The 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:(37)uxt+Δtuyt+Δtuzt+ΔtΔtu˙xt+Δt2Δtu˙yt+Δt2Δtu˙zt+Δt2=A⋅uxtuytuztΔtu˙xt−Δt2Δtu˙yt−Δt2Δtu˙zt−Δt2,where(38)A=1−k1Δt2m−k2Δt2m−k3Δt2m1−112cΔt3mcΔtm−cΔtm−k2Δt2m1−k1Δt2m−k3Δt2mcΔtm1−112cΔt3m−cΔtm−k3Δt2m−k3Δt2m1−k1Δt2m−cΔtm−cΔtm1−112cΔt3m−k1Δt2m−k2Δt2m−k3Δt2m1−112cΔt3mcΔtm−cΔtm−k2Δt2m−k1Δt2m−k3Δt2mcΔtm1−112cΔt3m−cΔtm−k3Δt2m−k3Δt2m−k1Δt2m−cΔtm−cΔtm1−112cΔt3m,m=ρ+7ρ˜L38,c=ρL2aN722CSaT+CPaN,k1=ρCS2L32−48μ721−2μ+14aNρCS2L29R,k2=ρCS2L121−2μ−aNρCS2L224R.The maximum eigenvalue of the transfer matrixA can be solved as follows:(39)maxλi=166−3k1Δt2m+3k2Δt2m−115cΔtm−−49−345cΔtm+−6+3k1Δt2m−3k2Δt2m+115cΔtm2.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:(40)Δtcorner≤36mk1−k2+13225c2−115c3k1−k2.


## 4.1. Stability Analysis of the Surface Boundary Subsystem

Consider a regular hexahedral element from the surface boundary subsystem (Figure8(a)), as shown in Figure 9.Figure 9
Regular hexahedral element.The concentrated mass matrix of the regular hexahedral element of the internal medium can be expressed as(22)Me=ρL318⋯00180⋱⋮sym.⋱sym.181824×24.The stiffness matrix of the regular hexahedral element of the internal medium can be expressed as(23)Ke=k11k12⋯k18k21⋱⋮⋮⋱⋮k81⋯⋯k8824×24,kij=EL3161+μ1−2μa11a12a13a21a22a23a31a32a33,where(24)a11=1−μL2⋅ξiξj1+ηiηj31+ζiζj3+1−2μ2L2ηiηj1+ζiζj31+ξiξj3+ζiζj1+ξiξj31+ηiηj3,a22=1−μL2⋅ηiηj1+ζiζj31+ξiξj3+1−2μ2L2ξiξj1+ηiηj31+ζiζj3+ζiζj1+ξiξj31+ηiηj3,a33=1−μL2⋅ζiζj1+ξiξj31+ηiηj3+1−2μ2L2ξiξj1+ηiηj31+ζiζj3+ηiηj1+ξiξj31+ζiζj3,a12=1L21+ζiζj3μξiηj+1−2μ2ηiξj,a13=1L21+ηiηj3μξiζj+1−2μ2ζiξj,a21=1L21+ζiζj3μηiξj+1−2μ2ξiηj,a23=1L21+ξiξj3μηiζj+1−2μ2ζiηj,a31=1L21+ηiηj3μζiξj+1−2μ2ξiζj,a32=1L21+ξiξj3μζiηj+1−2μ2ηiζj.Node 1:ξ,η,ζ1=1,1,−1; Node 2: ξ,η,ζ2=−1,1,−1; Node 3: ξ,η,ζ3=−1,−1,−1; Node 4: ξ,η,ζ4=1,−1,−1; Node 5: ξ,η,ζ5=1,1,1; Node 6: ξ,η,ζ6=−1,1,1; Node 7: ξ,η,ζ7=−1,−1,1; and Node 8: ξ,η,ζ8=1,−1,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 (9) and (10). According to this method, the mass, stiffness, and damping matrix of a single internal medium element and artificial boundary element can be obtained. Then, the eight elements in the subsystem can be assembled according to the node number (the node number is shown in Figure 9). The fixed boundary conditions are applied to the surface nodes of the subsystem. Then, the motion equation of the central node can be expressed as follows:(25)ρ+ρ˜L32100010001u¨xu¨yu¨z+8ρL2aN272CSaT+CPaN100010001u˙xu˙yu˙z+ρCS2L32−48μ181−2μ+aN8ρCS2L29R100010001uxuyuz=0,where the subscripts x, y, and z denote the x-, y-, and z-directions, respectively.It can be seen from (25) that the motion equations in the x-, y-, and z-directions are decoupled and identical to each other. Therefore, the stability analysis can be performed in only one direction. Taking the x-direction as an example, the explicit time-domain stepwise integration format based on the central difference method can be expressed as follows [27]:(26)u˙xt+Δt2=u˙xt−Δt2+Δtu¨xt,uxt+Δt=uxt+Δtu˙xt+Δt2.By combining (25) and (26), the following expression can be obtained:(27)uxt+ΔtΔtu˙xt+Δt2=AuxtΔtu˙xt−Δt2,where(28)A=1−kΔt2m1−cΔtm−kΔt2m1−cΔtm,m=ρ+ρ˜L32,c=8ρL2aN272CSaT+CPaN,k=ρCS2L32−48μ181−2μ+aN8ρCS2L29R.The eigenvalues of transfer matrixA are(29)λ1,2=1−kΔt22m−cΔt2m±Δt2m−4km+k2Δt2+2ckΔt+c2.By substituting (29) into the stability condition ρA=maxλi≤1, the following expression can be obtained:(30)Δtsurface≤1k4km+c2−c.Equation (30) provides the analytical solution for the stability condition of the surface boundary subsystem. The stability condition is related to the wave velocities CP and CS, the mass density ρ, the Poisson ratio μ of the internal media, and the element size L.


## 4.2. 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:(31)ρ+3ρ˜L34100010001u¨xu¨yu¨z+ρL2aN362CSaT+CPaN16−10−11600016u˙xu˙yu˙z+13ρCS2L32−48μ121−2μ+aN4LR121−2μ−aNL4R0121−2μ−aNL4R32−48μ121−2μ+aN4LR00032−48μ121−2μ+aN4LRuxuyuz=0.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:(32)uxt+Δtuyt+Δtuzt+ΔtΔtu˙xt+Δt2Δtu˙yt+Δt2Δtu˙zt+Δt2=A⋅uxtuytuztΔtu˙xt−Δt2Δtu˙yt−Δt2Δtu˙zt−Δt2,where(33)A=1−k1Δt2m−k2Δt2m01−16cΔtmcΔtm0−k2Δt2m1−k1Δt2m0cΔtm1−16cΔtm0001−k1Δt2m001−16cΔtm−k1Δt2m−k2Δt2m01−16cΔtmcΔtm0−k2Δt2m−k1Δt2m0cΔtm1−16cΔtm000−k1Δt2m001−16cΔtm,m=ρ+3ρ˜L34,c=ρL2aN362CSaT+CPaN,k1=ρCS2L32−48μ361−2μ+4aNρCS2L23R,k2=ρCS2L61−2μ−aNρCS2L212R.The maximum eigenvalue of the transfer matrixA can be solved as follows:(34)maxλi=122−k1Δt2m+k2Δt2m−17cΔtm−−41−17cΔtm+−2+k1Δt2m−k2Δt2m+17cΔtm2.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:(35)Δtedge≤4mk1−k2+289c2−17ck1−k2.


## 4.3. 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:(36)ρ+7ρ˜L38100010001u¨xu¨yu¨z+ρL2aN722CSaT+CPaN1123−11−111231111123u˙xu˙yu˙z+13ρCS2L32−48μ241−2μ+aN14L3R141−2μ−aNL8R−141−2μ+aNL8R141−2μ−aNL8R32−48μ241−2μ+aN14L3R−141−2μ+aNL8R−141−2μ+aNL8R−141−2μ+aNL8R32−48μ241−2μ+aN14L3Ruxuyuz=0.The 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:(37)uxt+Δtuyt+Δtuzt+ΔtΔtu˙xt+Δt2Δtu˙yt+Δt2Δtu˙zt+Δt2=A⋅uxtuytuztΔtu˙xt−Δt2Δtu˙yt−Δt2Δtu˙zt−Δt2,where(38)A=1−k1Δt2m−k2Δt2m−k3Δt2m1−112cΔt3mcΔtm−cΔtm−k2Δt2m1−k1Δt2m−k3Δt2mcΔtm1−112cΔt3m−cΔtm−k3Δt2m−k3Δt2m1−k1Δt2m−cΔtm−cΔtm1−112cΔt3m−k1Δt2m−k2Δt2m−k3Δt2m1−112cΔt3mcΔtm−cΔtm−k2Δt2m−k1Δt2m−k3Δt2mcΔtm1−112cΔt3m−cΔtm−k3Δt2m−k3Δt2m−k1Δt2m−cΔtm−cΔtm1−112cΔt3m,m=ρ+7ρ˜L38,c=ρL2aN722CSaT+CPaN,k1=ρCS2L32−48μ721−2μ+14aNρCS2L29R,k2=ρCS2L121−2μ−aNρCS2L224R.The maximum eigenvalue of the transfer matrixA can be solved as follows:(39)maxλi=166−3k1Δt2m+3k2Δt2m−115cΔtm−−49−345cΔtm+−6+3k1Δt2m−3k2Δt2m+115cΔtm2.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:(40)Δtcorner≤36mk1−k2+13225c2−115c3k1−k2.


## 5. Validation

A 3D homogeneous elastic half space model is established, as shown in Figure10. The dimensions of the model are 100 m × 100 m × 50 m, and the mesh size is 2 m × 2 m × 2 m. The mass density of the internal medium ρ is 2500 kg/m3, the shear wave velocity CS is 500 m/s, the compression wave velocity CP 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 (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.Figure 10
3D homogeneous model and viscoelastic artificial boundary elements. (a) 3D homogeneous model. (b) Viscoelastic artificial boundary elements.Figure 11
The time history of the concentrated impulsive force.According to (21), (30), (35), and (40), the maximum stable time increments for the internal region Δtinternal, the surface boundary subsystem Δtsurface, the edge boundary subsystem Δtedge, and the corner boundary subsystem Δtcorner in this model can be calculated, whose values are 0.0021 s, 0.00176 s, 0.00071 s, and 0.00032 s, respectively. The 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 surfacesLS, to the edges LE, and to the corners LC gradually increase (LS < LE < LC); 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=Δtinternal = 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>Δtsurface); thus, instability appears near the center of the bottom surface, as shown in Figure 12(a).Figure 12
Three kinds of instability states with time increments of (a) 0.0021 s, (b) 0.00176 s, and (c) 0.00071 s.
(a)(b)(c)As we reduce the time increment to the maximum stable time increment of the surface boundary subsystem (Δt=Δtsurface = 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>Δtedge).Then, the time increment is further reduced to the maximum stable time increment of the edge boundary subsystem (Δt=Δtedge = 0.00071 s). This time the simulation of the wave propagation can be accurately conducted until the waves arrive at the bottom corners of the model. The 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>Δtcorner).Finally, the maximum stable time increment of the corner boundary subsystem is adopted for the explicit analysis (Δt=Δtcorner = 0.00032 s). The 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.Figure 13
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.
(a)(b)(c)(d)(e)(f)Figure 14
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ρ from 1000 kg/m3 to 2500 kg/m3, the velocity of the shear wave of the internal medium CS 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. The results are shown in Figure 15.Figure 15
The maximum time increments under different (a) mass densities, (b) velocities of shear waves, (c) element sizes, (d) model sizes, and (e) Poisson ratios.
(a)(b)(c)(d)(e)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 CS is constant, the compression wave velocity CP can be expressed as follows:(41)CP=CS21−μ1−2μ.When the Poisson ratioμ increases, the compression wave velocity of the internal medium CP 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. The stability condition of the corner subsystem given by (40) is a sufficient condition for the numerical stability of the entire model.


## 6. Conclusion

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:(1)
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.(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. 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.(3)
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.(4)
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.



---

*Source: 1848169-2020-08-29.xml*