# Multiscale Method for Seismic Response of Near-Source Sites

Authors: Shutao Li; Jingbo Liu; Zhou Yang; Xin Bao; Fei Wang; Xiaofeng Wang; Mohsen Saleh Asheghabadi
Journal: Advances in Civil Engineering
 (2020)

---

## Abstract

The traditional source-site-structure model for the calculation of seismic response of underground structures at near-source sites is restricted by the grid scale and the size of the structure. As a result, an excessive number of elements in the model make the numerical solving process difficult. To solve problems such as an inefficient computation and challenging nonlinear simulation, a multiscale analysis method for the calculation of the seismic response of underground structures at near-source sites is developed. The generalized free-field seismic response of the near-source region is obtained by establishing a large-scale calculation model of the source site and is used to simulate the fracture mechanism of faults and the process of seismic wave propagation. Then, using the method of seismic wave input based on artificial boundary substructures, the free-field motion of the wave is transformed into the equivalent seismic load, which is the seismic wave input data for the small-scale region of interest. Finally, with the help of local elements with special shapes to realize the grid transition of different scales, a small-scale model with reasonable soil-underground structure interaction is established, and the seismic response of the overall model can be effectively solved. The calculation and analysis of the seismic response of underground structures in irregular terrain are carried out. Compared with the results obtained directly from the source-site-structure model, the multiscale method has satisfactory accuracy and meets the needs of engineering design. Since the number of elements is fewer and the calculation time is much shorter than those required by the traditional model, the advantages in computational efficiency of the new method are highlighted. In addition, the reflected waves are too weak to have a considerable impact on structures because of the great energy loss at the reflection interface, which further proves the feasibility of the closed artificial boundary substructure method.


---

## Body

## 1. Introduction

Most scholars consider that near-source earthquakes correspond to ground motions that depend strongly on the fracture mechanism of faults, reflecting significant rupture directivity effects and permanent displacements [1–5]. The near-source area is generally considered to be less than 20 km from the fault fracture surface [2–5]. Due to the short distance from the source, the ground motion generated in this area by the source contains a relatively large number of near-field and midfield components [6]. The dynamic soil-structure interaction (SSI) of underground structures under the action of near-source earthquakes has always been a research focus in the field of seismic engineering [7–15]. An underground structure is in an approximate semi-infinite space, and the dynamic response of the structure is completely different from the dynamic response of a structure above ground [16]. At present, the research methods for the SSI problem mainly include theoretical analysis, model testing and numerical simulation [17]. With the rapid development of computational capabilities and the continuous improvement in finite element and finite difference methods, numerical simulation technology has become the main research means.According to the relative position of the earthquake source and the structure, the dynamic SSI problem can be divided into internal and external source problems. The SSI problem in near-source regions or far-field explosions is generally treated as an external source problem. Since the near-field and midfield terms of the near-source ground motion cannot be ignored, the seismic wave field in the near-source region is complicated, and it is difficult to find a suitable simplified model for describing the site response. The proposed response spectrum or the time-history data exhibit larger errors at some characteristic points. If the solution of the structural dynamic reaction is directly based on the seismic wave field, which is predicted by the near-source seismic analysis, it is necessary to solve the problem of seismic wave input. Many studies have been conducted on this subject; in an SSI model that does not include the earthquake source, the external load is input into the calculation area through the cutoff boundary. Considering the characteristics of the wave field, the wave method was first proposed by Liu and Lu and uses viscoelastic artificial boundaries [18]. The input motion can be transformed into an equivalent input load on the artificial boundary, aiming to keep the motion and stress of the nodes on the artificial boundary the same as those of the free-wave field. Several researchers have demonstrated that the wave method is effective and accurate [19–22]. However, the model and calculation process of the wave method is more complicated, and the direction and magnitude of the equivalent input seismic load of each node on the artificial boundary need to be determined individually. Based on the theory of the wave method, a new seismic wave input method based on an artificial boundary substructure is proposed [23]; boundary source technology is also the basis of this method. In the method proposed above, calculation of the equivalent input seismic load is realized by dynamic analysis of the artificial boundary substructure.With the research development of physics of earthquakes and ground motion simulation [2, 3, 24–29], it is possible to consider the dynamic SSI problem as an internal source problem. The problems such as the seismic wave input and the artificial boundary conditions can be solved by establishing a large-scale overall calculation model (source-site-structure model). Two main problems exist in the large-scale internal source method. First, for the analysis of the source region, the size of the elements ranges from several meters to several tens of meters. For the SSI problem, the size of the elements needs to be controlled to the decimeter level. If these two scales cannot be unified, a large number of elements will need to be used in the calculation model. For a wide range of sites, dynamic explicit methods are often used to improve the efficiency of computation. If local mesh encryption is performed in the area of interest, it does not improve the stability limit of the overall model, and the computation scale will still be too large or the model may be impossible to solve. Second, the intensity of the ground motion in the near-source region is relatively large, and the structure and soil easily reach the state of nonlinearity. The limitation of the calculation step will be strict, which increases the difficulty of calculation.Multiscale science studies the coupled phenomena at distinct length or time scales [30]. For finite element analysis in structural engineering, different parts of a structure can be modeled by using elements of different types and scales. For example, the beam element or the shell element is usually adopted for large-scale structural members, and the solid element is usually adopted for small-scale or key structural members. Researchers from various countries have conducted preliminary research on multiscale analysis. Dhia et al. used the Arlequin method to couple two different mechanical models [31, 32] in the coupled region by a coupling operator, and energy distribution can be realized between the two models. Based on the XFEM and global multiplex grid method, the multiscale three-dimensional crack propagation analysis of a simple structure was carried out by Rannou et al. [33, 34]. The multiscale method used in the analysis of bridges and underground structures was proposed by Deng et al. [35] as follows: first, establish the overall calculation model of large-scale coarse grids, perform dynamic analysis, determine the key components and dangerous parts of the structure, subdivide the local mesh, perform dynamic analysis of the small-scale mesh model, and obtain the stress and deformation of the key components and dangerous parts. The field reduction method (DRM) proposed by Bielak et al. is similar to the abovementioned approach [36–38]. The first step of the DRM is to simulate the source and propagation path through the large-scale free-field model, and the second step is to input the equivalent local force into the small-scale local model. The DRM overcomes the problem that the physical scale cannot be unified. In addition, Pan et al. also proposed mesh reduction procedures that can be applied in BEM approaches [39, 41].To solve the calculation problem of the large-scale source-site-structure model, this paper draws on the ideas of multiscale analysis. The seismic response of underground structures in near-source regions is proposed to be solved by the dynamic analysis of free fields with large-scale grids, while the analysis of soil-structure dynamic interactions is proposed to be solved by that with small-scale grids. The calculation efficiency can be greatly improved to an acceptable computational accuracy when the proposed modeling method is used. The computational stability of the viscoelastic artificial boundary element in the dynamic explicit calculation is preliminarily solved. This paper presents the detailed steps required to implement the proposed method, which is validated by modeling large-scale complex fields.


## 2. Multiscale Method Based on Artificial Boundary Substructure

### 2.1. Determination of Equivalent Seismic Load Input

To study the dynamic response of the SSI system during an earthquake, we first need to establish a model for the soil within a limited range from the layered foundation, and a reasonable artificial boundary needs to be set to absorb the scattering wave generated internally. According to the basic principle of the wave method [18], the seismic wave input can be converted into an equivalent seismic load applied on the cutoff boundary. The equation of motion of the SSI system under the equivalent seismic load can be written as follows:(1)Mu¨+Cu˙+Ku=F,where u, u, and u¨ are the displacement, velocity, and acceleration vectors of the finite element model nodes, respectively; F is the equivalent input seismic load vector applied to the cutoff boundary; M, C, and K are the mass, damping, and stiffness matrices of the SSI model, respectively; and the artificial boundary conditions can be integrated into M, C, and K. The incident seismic wave is determined by the source, propagation path, and local site condition, which is independent of the structure, so the determination of the equivalent input seismic load is only related to the free field. To obtain the equivalent input seismic load, we can establish a free-field finite element; the equivalent input seismic load is applied to the nodes on the cutoff boundary, and the motion equation of the model is as follows:(2)M0u¨+C0u˙+K0u=F,where M0, C0, and K0 are the mass, damping, and stiffness matrices of the free-field finite element model, respectively. Under a reasonable equivalent input seismic load, the motion of the free-wave field satisfies the following formula:(3)M0u¨0+C0u˙0+K0u0=F,where u0, u˙0, and u¨0 are the displacement, velocity, and acceleration vectors corresponding to the free-field model nodes, respectively. According to formula (3), if the free-wave field motion is applied to all the nodes of the free-field model, the constrained reaction force of the nodes on the cutoff boundary is the equivalent input seismic load F.


### 2.2. Seismic Wave Input Method Based on Artificial Boundary Substructure

The element nodes in the finite element model can be divided into three categories according to position: the artificial boundary nodes (), the internal nodes () adjacent to the artificial boundary nodes, and the remaining internal nodes. The artificial boundary substructure is defined as a structure consisting of elements containing artificial boundary nodes, as shown in Figure 1. According to the above node classification, the motion equations (formulas (1) and (2)) of the SSI model and the free-field finite element model can be written in the form of the following block matrix:(4)MIIMIC0MCIMCCMCB0MBCMBBu¨Iu¨Cu¨B+CIICIC0CCICCCCCB0CBCCBBu˙Iu˙Cu˙B+KIIKIC0KCIKCCKCB0KBCKBBuIuCuB=00FB,MII0MIC00MCI0MCC0MCB00MBC0MBB0u¨Iu¨Cu¨B+CII0CIC00CCI0CCC0CCB00CBC0CBB0u˙Iu˙Cu˙B+KII0KIC00KCI0KCC0KCB00KBC0KBB0uIuCuB=00FB.Figure 1
Artificial boundary substructure.The subscripts B, C, and I represent the artificial boundary nodes (), the internal nodes () adjacent to the artificial boundary, and the remaining internal nodes, respectively; the superscript 0 represents the free-field model; FB is the equivalent input seismic load acting on the artificial boundary nodes. If the element mesh on the artificial boundary of the free-field model is the same as that on the artificial boundary of the soil-structure model, then(5)MCB0=MCBMBC0=MBCMBB0=MBB,CCB0=CCBCBC0=CBCCBB0=CBB,KCB0=KCBKBC0=KBCKBB0=KBB.The equivalent input seismic load can be written as follows:(6)00FB=MII0MIC00MCI0MCC0MCB0MBCMBBu¨I0u¨C0u¨B0+CII0CIC00CCI0CCC0CCB0CBCCBBu˙I0u˙C0u˙B0+KII0KIC00KCI0KCC0KCB0KBCKBBuI0uC0uB0.The expression of the equivalent input seismic load is(7)FB=MBCu¨C0+MBBu¨B0+CBCu˙C0+CBBu˙B0+KBCuC0+KBBuB0.According to formula (7), the equivalent seismic load is independent of the internal nodes of the model. Therefore, by extracting the displacement history of all the nodes (() and ()) on the artificial boundary substructure in the free-wave field, the dynamic analysis of the artificial boundary substructure can be carried out to calculate the constrained reaction force of the artificial boundary nodes (). As a result, input of the equivalent seismic load can be realized, as shown in Figure 2.Figure 2
Seismic response analysis model of the SSI system.


### 2.3. Idea of Multiscale Method

The core idea of the multiscale method is to create a transition from the large-scale model with a coarse mesh to the small-scale model with a fine mesh based on the artificial boundary substructure and the seismic wave input method. The free-wave field motion calculated by the source-site model can be transformed into the equivalent seismic load input applied on the SSI model for the dynamic reaction analysis of the structure. The implementation steps of the multiscale method are as follows:Step 1: establish a large-scale calculation model of the source site with a coarse mesh, simulate the mechanism of the source fracture, carry out the dynamic calculation of the whole region, obtain the generalized free-site motion in the near-source area, and extract the displacement history of the corresponding nodes (() and ()) on the artificial boundary substructure in the free field.Step 2: load the displacement history of the free-field motion on the corresponding nodes (() and ()), carry out the dynamic calculation of the artificial boundary substructure, and obtain the constraint reaction force history of the nodes () on the artificial boundary.Step 3: establish the soil-structure interaction model, and set the history of the reaction force obtained in the previous step as the equivalent seismic load input to carry out the dynamic calculation of the SSI model and solve the dynamic response of the structure.A three-step flowchart of the multiscale method is shown in Figure3. It is important to note the following conditions: (1) The SSI model is surrounded by an area with a coarse mesh, but the mesh near the structure is a fine mesh; therefore, there exists a mesh transition area. If the lower calculation precision is accepted, the mesh transition area can be created directly via several triangular or irregular quadrilateral elements, as shown in Figure 4. If higher calculation precision is required to prevent the influence of the spectral difference between the meshes on different scales, the research results of other scholars can be used for the collaborative calculation of elements in the transition region via bridge domain coupling theory [42], the displacement coupling method [43], collaborative working interface technology [44], etc. (2) Linear elastic analysis is usually adopted for the large-scale source-site model, but nonlinear analysis can be used for the SSI model if necessary.Figure 3
Three-step flowchart of the multiscale method.Figure 4
Grid transition area diagram.


## 2.1. Determination of Equivalent Seismic Load Input

To study the dynamic response of the SSI system during an earthquake, we first need to establish a model for the soil within a limited range from the layered foundation, and a reasonable artificial boundary needs to be set to absorb the scattering wave generated internally. According to the basic principle of the wave method [18], the seismic wave input can be converted into an equivalent seismic load applied on the cutoff boundary. The equation of motion of the SSI system under the equivalent seismic load can be written as follows:(1)Mu¨+Cu˙+Ku=F,where u, u, and u¨ are the displacement, velocity, and acceleration vectors of the finite element model nodes, respectively; F is the equivalent input seismic load vector applied to the cutoff boundary; M, C, and K are the mass, damping, and stiffness matrices of the SSI model, respectively; and the artificial boundary conditions can be integrated into M, C, and K. The incident seismic wave is determined by the source, propagation path, and local site condition, which is independent of the structure, so the determination of the equivalent input seismic load is only related to the free field. To obtain the equivalent input seismic load, we can establish a free-field finite element; the equivalent input seismic load is applied to the nodes on the cutoff boundary, and the motion equation of the model is as follows:(2)M0u¨+C0u˙+K0u=F,where M0, C0, and K0 are the mass, damping, and stiffness matrices of the free-field finite element model, respectively. Under a reasonable equivalent input seismic load, the motion of the free-wave field satisfies the following formula:(3)M0u¨0+C0u˙0+K0u0=F,where u0, u˙0, and u¨0 are the displacement, velocity, and acceleration vectors corresponding to the free-field model nodes, respectively. According to formula (3), if the free-wave field motion is applied to all the nodes of the free-field model, the constrained reaction force of the nodes on the cutoff boundary is the equivalent input seismic load F.


## 2.2. Seismic Wave Input Method Based on Artificial Boundary Substructure

The element nodes in the finite element model can be divided into three categories according to position: the artificial boundary nodes (), the internal nodes () adjacent to the artificial boundary nodes, and the remaining internal nodes. The artificial boundary substructure is defined as a structure consisting of elements containing artificial boundary nodes, as shown in Figure 1. According to the above node classification, the motion equations (formulas (1) and (2)) of the SSI model and the free-field finite element model can be written in the form of the following block matrix:(4)MIIMIC0MCIMCCMCB0MBCMBBu¨Iu¨Cu¨B+CIICIC0CCICCCCCB0CBCCBBu˙Iu˙Cu˙B+KIIKIC0KCIKCCKCB0KBCKBBuIuCuB=00FB,MII0MIC00MCI0MCC0MCB00MBC0MBB0u¨Iu¨Cu¨B+CII0CIC00CCI0CCC0CCB00CBC0CBB0u˙Iu˙Cu˙B+KII0KIC00KCI0KCC0KCB00KBC0KBB0uIuCuB=00FB.Figure 1
Artificial boundary substructure.The subscripts B, C, and I represent the artificial boundary nodes (), the internal nodes () adjacent to the artificial boundary, and the remaining internal nodes, respectively; the superscript 0 represents the free-field model; FB is the equivalent input seismic load acting on the artificial boundary nodes. If the element mesh on the artificial boundary of the free-field model is the same as that on the artificial boundary of the soil-structure model, then(5)MCB0=MCBMBC0=MBCMBB0=MBB,CCB0=CCBCBC0=CBCCBB0=CBB,KCB0=KCBKBC0=KBCKBB0=KBB.The equivalent input seismic load can be written as follows:(6)00FB=MII0MIC00MCI0MCC0MCB0MBCMBBu¨I0u¨C0u¨B0+CII0CIC00CCI0CCC0CCB0CBCCBBu˙I0u˙C0u˙B0+KII0KIC00KCI0KCC0KCB0KBCKBBuI0uC0uB0.The expression of the equivalent input seismic load is(7)FB=MBCu¨C0+MBBu¨B0+CBCu˙C0+CBBu˙B0+KBCuC0+KBBuB0.According to formula (7), the equivalent seismic load is independent of the internal nodes of the model. Therefore, by extracting the displacement history of all the nodes (() and ()) on the artificial boundary substructure in the free-wave field, the dynamic analysis of the artificial boundary substructure can be carried out to calculate the constrained reaction force of the artificial boundary nodes (). As a result, input of the equivalent seismic load can be realized, as shown in Figure 2.Figure 2
Seismic response analysis model of the SSI system.


## 2.3. Idea of Multiscale Method

The core idea of the multiscale method is to create a transition from the large-scale model with a coarse mesh to the small-scale model with a fine mesh based on the artificial boundary substructure and the seismic wave input method. The free-wave field motion calculated by the source-site model can be transformed into the equivalent seismic load input applied on the SSI model for the dynamic reaction analysis of the structure. The implementation steps of the multiscale method are as follows:Step 1: establish a large-scale calculation model of the source site with a coarse mesh, simulate the mechanism of the source fracture, carry out the dynamic calculation of the whole region, obtain the generalized free-site motion in the near-source area, and extract the displacement history of the corresponding nodes (() and ()) on the artificial boundary substructure in the free field.Step 2: load the displacement history of the free-field motion on the corresponding nodes (() and ()), carry out the dynamic calculation of the artificial boundary substructure, and obtain the constraint reaction force history of the nodes () on the artificial boundary.Step 3: establish the soil-structure interaction model, and set the history of the reaction force obtained in the previous step as the equivalent seismic load input to carry out the dynamic calculation of the SSI model and solve the dynamic response of the structure.A three-step flowchart of the multiscale method is shown in Figure3. It is important to note the following conditions: (1) The SSI model is surrounded by an area with a coarse mesh, but the mesh near the structure is a fine mesh; therefore, there exists a mesh transition area. If the lower calculation precision is accepted, the mesh transition area can be created directly via several triangular or irregular quadrilateral elements, as shown in Figure 4. If higher calculation precision is required to prevent the influence of the spectral difference between the meshes on different scales, the research results of other scholars can be used for the collaborative calculation of elements in the transition region via bridge domain coupling theory [42], the displacement coupling method [43], collaborative working interface technology [44], etc. (2) Linear elastic analysis is usually adopted for the large-scale source-site model, but nonlinear analysis can be used for the SSI model if necessary.Figure 3
Three-step flowchart of the multiscale method.Figure 4
Grid transition area diagram.


## 3. Application of Viscoelastic Artificial Boundary Elements in the Dynamic Explicit Calculation

In this paper, the size of the source-site model is large, and the selection of 4-node planar strain elements will cause the number of elements to be millions to tens of millions. The dynamic implicit calculation method uses more resources and results in a slower calculation speed than those of the dynamic explicit calculation method. As a result, the dynamic explicit calculation method should be adopted. For the treatment of the cutoff boundary in the soil, the use of equivalent viscoelastic artificial boundary elements can significantly improve the preprocessing efficiency and result in acceptable calculation accuracy [45]. However, when the artificial boundary elements participate in the calculation of the whole model, the physical size and the material parameters of the elements become important factors that would restrict the step size of the incremental dynamic explicit calculation and affect the calculation stability of the whole model. When the wave reaches the elements on the artificial boundary, calculation instability easily occurs. This issue is discussed in this section.


### 3.1. Viscoelastic Artificial Boundary Elements in Two Dimensions

The elasticity coefficient and damping coefficient of a two-dimensional equivalent physical system of a viscoelastic artificial boundary are described below [46]:Tangent direction:(8)KBT=αTGR,CBT=ρCS.Normal direction:(9)KBN=αNGR,CBN=ρCP,where KBN and KBT are the stiffnesses in the normal and tangent directions, respectively, R is the distance from the wave source to the nodes on the artificial boundary, CS and CP are the velocities of the S wave and P wave, respectively, G is the shear modulus of the medium, ρ is the mass density of the medium, and aT and aN are the parameters of the viscoelastic artificial boundary in the tangent and normal directions, respectively.The method of the equivalent viscoelastic artificial boundary elements is used to extend a layer of the same type of elements along the normal direction on the boundary of the finite element model. The outer nodes of the model are fixed, and the material properties of the equivalent elements are defined to be equivalent to the viscoelastic artificial boundary, as shown in Figure5. The equivalent stiffness, Poisson’s ratio, and equivalent damping of the artificial boundary elements are [45](10)E˜=αNhGR⋅1+μ˜1−2μ˜1−μ˜,μ˜=αN/αT−22αN/αT−1,η˜=ρR2GCSαT+CPαN,where h is the thickness of the equivalent elements and the mass density of the equivalent elements is zero. In the two-dimensional model, it is recommended that aT=0.5 and aN=1 [45] so that μ˜=0. Thus, all the equivalent physical parameters of the artificial boundary elements are determined.Figure 5
Artificial boundary element size.


### 3.2. Stability of Dynamic Explicit Calculation

The stability limit has a great influence on the reliability and accuracy of the dynamic explicit calculation. To improve the computation efficiency, the increment must be determined by analysis. The adopted increment should be as small as possible and not greater than the stability limit. The undamped and damped stability limits in the dynamic explicit calculations are as follows [47]:Undamped stability limit:(11)Δtstable=2ωmax.Damped stability limit:(12)Δtstable=2ωmax1+ξ2−ξ,where ωmax is the highest frequency of the material in the system and ξ is the critical damping portion of the highest frequency mode of the material. In this study, the highest frequency is based on a complex set of interaction factors, and it is difficult to calculate the exact value. An alternative method is to apply an effective and conservative estimation, which is defined by the length Le of the elements and the compression wave velocity CP of the materials [47]:(13)Δtstable=LeCP.Two important assumptions are adopted in deriving the stiffness matrix of the equivalent artificial boundary elements: one is that the thickness of the boundary element is much smaller than the width, as shown in Figure5; the other is that the mass density of any artificial boundary element is zero, such that the compression wave velocity in these elements is infinitely large. If the dynamic explicit calculation method is used, the artificial boundary elements will become the most unstable part in the model. To ensure the smooth progress of the calculation, it is necessary to study the size and the compression wave velocity of the artificial boundary elements. According to formula (13), two factors (Le and CP) are important for determining the stability of the boundary elements. Lu et al. [44] studied the effect of the thickness-to-width ratio of the two-dimensional viscoelastic artificial boundary elements on the calculation accuracy based on the dynamic implicit algorithm. Their conclusion was that the thickness of the artificial boundary elements has little influence on the calculation accuracy. Even if the thickness increases by 500%, the error in the results is still not obvious. Referring to the above conclusions, to exclude the influence of the element size on the calculation stability, it is recommended to set the thickness and width of the artificial boundary elements to be equal in the dynamic explicit calculation, as shown in Figure 6.Figure 6
Artificial boundary element size for dynamic explicit calculation.To study the influence of the compression wave velocity on the stability limit, a two-dimensional homogeneous elastic half-space model is established as an example. Using dimensionless units, a plane with a width of 300 and a height of 150 is taken as the area for calculation, and the size of every element is 10 × 10. A layer of two-dimensional elements with a thickness of 10 is extended outward along the normal of the cutoff boundary to create the equivalent artificial boundary elements, and the outermost nodes of the boundary elements are fixed. The shear wave velocity in the model is 3500, the compression wave velocity is 6548, the density is 2500, and Poisson’s ratio is 0.3. For the artificial boundary elements, both the thickness and width are 10, the equivalent elastic modulus is 2.04×109, the mass density is 1 (the material density of the elements in the general computing software cannot be set to 0, so the value 1 is adopted as an approximation), aT is 0.5, and aN is 1.0. The example model is shown in Figure 7, in which the positions of the observation points A, B, C, and D are also shown. The load history for the model is shown in Figure 8.Figure 7
Two-dimensional homogeneous elastic half-space model.Figure 8
Load history.According to formula (13), the stability limit of the calculated region is 0.00152, and that of the artificial boundary elements is 0.000221. To ensure that the computation efficiency is not restricted by the artificial boundary elements, the mass density of the artificial boundary elements is increased to reduce the compression wave velocity of the boundary elements. As a result, the stability limit of the artificial boundary elements is not less than that of the calculation region. The increment step is set to be fixed and is equal to 0.0001 in the dynamic explicit calculation. The calculation results show that the increase in the mass density of the artificial boundary elements has no obvious effect on the calculation accuracy, which is reflected in the calculation results of the four observation points in Figure 9.Figure 9
Effects of different mass densities of the boundary elements on theY-direction displacement at the observation points: (a) node (A), (b) node (B), (c) node (C), and (d) node D.
(a)
(b)
(c)
(d)According to the principle that the stability limit of the boundary elements should not be less than the stability limit of the material elements, the mass density value of the boundary elements can be estimated by formula (13). In this example, if the mass density of the boundary elements increases, the stability limit will become longer. To reduce the total time of calculation, the constant increment step can be increased appropriately.Figure10 compares the calculation results by increasing the constant increment step by 10 times the initial step size considered. The results show that this method can significantly improve the efficiency of the dynamic explicit calculation (the former calculation time is 126, and the latter calculation time is 15). The proposed method can satisfy the needs of large-scale numerical calculations with acceptable calculation accuracy.Figure 10
Comparison of theY-direction displacement at the observation points after adjusting the density of the artificial boundary elements and the calculating step size: (a) points A and B, (b) points C and D.
(a)
(b)


## 3.1. Viscoelastic Artificial Boundary Elements in Two Dimensions

The elasticity coefficient and damping coefficient of a two-dimensional equivalent physical system of a viscoelastic artificial boundary are described below [46]:Tangent direction:(8)KBT=αTGR,CBT=ρCS.Normal direction:(9)KBN=αNGR,CBN=ρCP,where KBN and KBT are the stiffnesses in the normal and tangent directions, respectively, R is the distance from the wave source to the nodes on the artificial boundary, CS and CP are the velocities of the S wave and P wave, respectively, G is the shear modulus of the medium, ρ is the mass density of the medium, and aT and aN are the parameters of the viscoelastic artificial boundary in the tangent and normal directions, respectively.The method of the equivalent viscoelastic artificial boundary elements is used to extend a layer of the same type of elements along the normal direction on the boundary of the finite element model. The outer nodes of the model are fixed, and the material properties of the equivalent elements are defined to be equivalent to the viscoelastic artificial boundary, as shown in Figure5. The equivalent stiffness, Poisson’s ratio, and equivalent damping of the artificial boundary elements are [45](10)E˜=αNhGR⋅1+μ˜1−2μ˜1−μ˜,μ˜=αN/αT−22αN/αT−1,η˜=ρR2GCSαT+CPαN,where h is the thickness of the equivalent elements and the mass density of the equivalent elements is zero. In the two-dimensional model, it is recommended that aT=0.5 and aN=1 [45] so that μ˜=0. Thus, all the equivalent physical parameters of the artificial boundary elements are determined.Figure 5
Artificial boundary element size.


## 3.2. Stability of Dynamic Explicit Calculation

The stability limit has a great influence on the reliability and accuracy of the dynamic explicit calculation. To improve the computation efficiency, the increment must be determined by analysis. The adopted increment should be as small as possible and not greater than the stability limit. The undamped and damped stability limits in the dynamic explicit calculations are as follows [47]:Undamped stability limit:(11)Δtstable=2ωmax.Damped stability limit:(12)Δtstable=2ωmax1+ξ2−ξ,where ωmax is the highest frequency of the material in the system and ξ is the critical damping portion of the highest frequency mode of the material. In this study, the highest frequency is based on a complex set of interaction factors, and it is difficult to calculate the exact value. An alternative method is to apply an effective and conservative estimation, which is defined by the length Le of the elements and the compression wave velocity CP of the materials [47]:(13)Δtstable=LeCP.Two important assumptions are adopted in deriving the stiffness matrix of the equivalent artificial boundary elements: one is that the thickness of the boundary element is much smaller than the width, as shown in Figure5; the other is that the mass density of any artificial boundary element is zero, such that the compression wave velocity in these elements is infinitely large. If the dynamic explicit calculation method is used, the artificial boundary elements will become the most unstable part in the model. To ensure the smooth progress of the calculation, it is necessary to study the size and the compression wave velocity of the artificial boundary elements. According to formula (13), two factors (Le and CP) are important for determining the stability of the boundary elements. Lu et al. [44] studied the effect of the thickness-to-width ratio of the two-dimensional viscoelastic artificial boundary elements on the calculation accuracy based on the dynamic implicit algorithm. Their conclusion was that the thickness of the artificial boundary elements has little influence on the calculation accuracy. Even if the thickness increases by 500%, the error in the results is still not obvious. Referring to the above conclusions, to exclude the influence of the element size on the calculation stability, it is recommended to set the thickness and width of the artificial boundary elements to be equal in the dynamic explicit calculation, as shown in Figure 6.Figure 6
Artificial boundary element size for dynamic explicit calculation.To study the influence of the compression wave velocity on the stability limit, a two-dimensional homogeneous elastic half-space model is established as an example. Using dimensionless units, a plane with a width of 300 and a height of 150 is taken as the area for calculation, and the size of every element is 10 × 10. A layer of two-dimensional elements with a thickness of 10 is extended outward along the normal of the cutoff boundary to create the equivalent artificial boundary elements, and the outermost nodes of the boundary elements are fixed. The shear wave velocity in the model is 3500, the compression wave velocity is 6548, the density is 2500, and Poisson’s ratio is 0.3. For the artificial boundary elements, both the thickness and width are 10, the equivalent elastic modulus is 2.04×109, the mass density is 1 (the material density of the elements in the general computing software cannot be set to 0, so the value 1 is adopted as an approximation), aT is 0.5, and aN is 1.0. The example model is shown in Figure 7, in which the positions of the observation points A, B, C, and D are also shown. The load history for the model is shown in Figure 8.Figure 7
Two-dimensional homogeneous elastic half-space model.Figure 8
Load history.According to formula (13), the stability limit of the calculated region is 0.00152, and that of the artificial boundary elements is 0.000221. To ensure that the computation efficiency is not restricted by the artificial boundary elements, the mass density of the artificial boundary elements is increased to reduce the compression wave velocity of the boundary elements. As a result, the stability limit of the artificial boundary elements is not less than that of the calculation region. The increment step is set to be fixed and is equal to 0.0001 in the dynamic explicit calculation. The calculation results show that the increase in the mass density of the artificial boundary elements has no obvious effect on the calculation accuracy, which is reflected in the calculation results of the four observation points in Figure 9.Figure 9
Effects of different mass densities of the boundary elements on theY-direction displacement at the observation points: (a) node (A), (b) node (B), (c) node (C), and (d) node D.
(a)
(b)
(c)
(d)According to the principle that the stability limit of the boundary elements should not be less than the stability limit of the material elements, the mass density value of the boundary elements can be estimated by formula (13). In this example, if the mass density of the boundary elements increases, the stability limit will become longer. To reduce the total time of calculation, the constant increment step can be increased appropriately.Figure10 compares the calculation results by increasing the constant increment step by 10 times the initial step size considered. The results show that this method can significantly improve the efficiency of the dynamic explicit calculation (the former calculation time is 126, and the latter calculation time is 15). The proposed method can satisfy the needs of large-scale numerical calculations with acceptable calculation accuracy.Figure 10
Comparison of theY-direction displacement at the observation points after adjusting the density of the artificial boundary elements and the calculating step size: (a) points A and B, (b) points C and D.
(a)
(b)


## 4. Examples

To verify the calculation accuracy and efficiency of the multiscale method, the homogeneous flat field (4.1), complex stratified field (4.2), and large-scale complex field (4.3) were selected as examples to calculate. The numerical simulation results of the corresponding source-site-structure overall model are used as the measurement standard to verify the accuracy and efficiency of the multiscale method.


### 4.1. Examples of a Homogeneous Flat Field

#### 4.1.1. Example 1

The two-dimensional site is 8 km long and 4 km wide, the mass density is 2000 kg/m3, the shear wave velocity is 2000 m/s, and Poisson’s ratio is 0.3. The burial depth of the underground structure is 20 m. The mass density of C30 concrete is 2490 kg/m3, its elastic modulus is 30 GPa, and its Poisson’s ratio is 0.15. The earthquake source is 3.9 km from the surface. The mechanism of the source fracture is simulated by a vertical concentrated force. The model and load history are shown in Figure 11. The large-scale source-site model is composed of square plane strain elements of 10 m in width. The elements of the SSI model are refined into square plane strain elements 5 m in width. According to the beam element theory in finite elements, if the focus is on the overall force of the structure rather than the local force, beam elements can be used to simulate slender parts. In our two-dimensional models, the dimensions in out-of-plane direction can be ignored, so, the structure can be simulated by beam elements; it can also obtain the results such as axial force and bending moment conveniently which benefit the process in data analysis.Figure 11
Model of Example 4.1.1 (burial depth is 20 m): (a) the large-scale source-site model, (b) load history, (c) artificial boundary substructure model, (d) SSI model, and (e) underground structure model.
(a)
(b)
(c)
(d)
(e)Figures12 and 13 show that the node displacement and the bending moment of sections (A, B, C, and D) calculated by the multiscale method are in good agreement with those calculated by the source-site-structure overall model.Figure 12
Displacement history of the nodes of the artificial boundary elements: (a) node 1, (b) node 2, (c) node 3, and (d) node 4.
(a)
(b)
(c)
(d)Figure 13
Comparison of bending moments of the structural section (burial depth is 20 m): (a) section (A), (b) section (B), (c) section (C), and (d) section D.
(a)
(b)
(c)
(d)


#### 4.1.2. Example 2

The structure is buried at a depth of 200 m, and the remaining parameters are the same as those in Section 4.1. The source-site model, the artificial boundary substructure, and the SSI model are shown in Figure14. As shown in Figure 15, the bending moments of sections A and C calculated by the multiscale method are in good agreement with those calculated by the source-site-structure overall model.Figure 14
Model of Example 4.1.2 (burial depth is 200 m): (a) the source-site model, (b) artificial boundary substructure, and (c) SSI model.
(a)
(b)
(c)Figure 15
Comparison of bending moment values of structural section (burial depth is 200 m): (a) section (A), (b) section C.
(a)
(b)With the increase in the burial depth of the structure, the size of the artificial boundary substructure in the vertical direction and the number of nodes and elements participating in the calculation increase, which is inconvenient to the modeling and calculation.


#### 4.1.3. Example 3

An optimization is carried out on the artificial boundary substructure of the deeply buried structure. A closed substructure is adopted instead of the original open substructure, and the other calculation parameters remain unchanged. The closed substructure model is shown in Figure16. As shown in Figure 17, in the case of the closed artificial boundary substructure, the bending moments of sections A and C calculated by the multiscale method are in good agreement with those calculated by the source-site-structure overall model. At the same time, the number of elements and nodes participating in the calculation is significantly reduced, and the calculation efficiency significantly increases. When the structure is deeply buried, after the complex input wave field is reflected on the soil-structure interface and then the ground surface, the reflection wave transmits to the structure again but has a small effect on the structure.Figure 16
Model of Example 4.1.3 (burial depth is 200 m): (a) closed artificial boundary substructure, (b) SSI model.
(a)
(b)Figure 17
Comparison of the bending moment values of the structural section (burial depth is 200 m): (a) section (A), (b) section C.
(a)
(b)


### 4.2. Examples of a Complex Layered Field

According to the complex site conditions, the layered field with mountains and rivers is selected as a calculation example. The calculation field is 10 km long and 5 km wide. The material parameters of medium 1 in the layered field are the same as those used in Section 4.1.1. For the material parameters of medium 2, the shear wave velocity is 3500 m/s, and the mass density is 2500 kg/m3. The burial depth of the structure is 300 m, and the remaining parameters are the same as those in Section 4.1.1. The calculation model and load history are shown in Figure 18. It should be noted that the artificial boundary substructures of Examples 4.2.1 and 4.2.3 do not contain medium 2. The artificial boundary substructures of Examples 4.2.2 and 4.2.4 contain medium 2. The four example models are shown in Figure 19.Figure 18
The complex layered source-site model and load history: (a) the layered source-site model, (b) load history.
(a)
(b)Figure 19
Four example models of complex layered fields: (a) Example 4.2.1, (b) Example 4.2.2, (c) Example 4.2.3, and (d) Example 4.2.4.
(a)
(b)
(c)
(d)Figure20 shows the comparison of the calculation results of the bending moments of the four examples with the dynamic calculation results of the source-site-structure overall model; under the condition of the complex layered field, the bending moments of sections A, B, C, and D calculated by the multiscale method are in good agreement with those calculated by the source-site-structure overall model whether the substructure is open or closed.Figure 20
Section bending moment history of the four examples of a complex layered field: (a) section (A), (b) section (B), (c) section (C), and (d) section D.
(a)
(b)
(c)
(d)Table1 shows the maximum bending moment values and differences calculated by the source-site-structure overall model and the four examples. The calculation accuracy of the model containing medium 2 (Examples 4.2.1 and 4.2.3) is slightly higher than that of the model that does not contain medium 2 (Examples 4.2.2 and 4.2.4). The reason for the error may be that Examples 4.4.2 and 4.2.4 consider the effect of the scattered wave energy on the structure after reflection by the layered interface, while Examples 4.2.1 and 4.2.3 are equivalent to extending medium 1 into infinite space. Since the scattered wave energy is much weaker than the incident wave energy, the error caused by ignoring the layered interface reflection is very small.Table 1
Maximum bending moments values and errors of the observed sections in 4.2.1–4.2.4 (N·m).
A
B
C
D
Source-site-structure calculation model
−47640
−56223
−55742
−47787
Example 4.2.1
Moment
−46029
−55968
−56387
−46313
Error
−3.38%
0.40%
0.29%
−3.08%
Example 4.2.2
Moment
−46758
−55475
−55770
−46985
Error
−1.85%
−0.48%
−0.81%
−1.68%
Example 4.2.3
Moment
−46710
−56538
−56975
−46897
Error
−1.95%
1.43%
1.34%
−1.86%
Example 4.2.4
Moment
−46904
−56706
−57000
−47135
Error
−1.55%
1.73%
1.38%
−1.36%


### 4.3. Example of a Large-Range and Complex Field with a Length of 100 km

According to the actual engineering problems, a large-range and complex site with a length of 100 km and a width of 20 km is selected for calculation. The shear wave velocity of the soil material is 3500 m/s, and the mass density is 2500 kg/m3. The earthquake source is located 10 km below the ground surface and 20 km from the structure. The mode of the field is shown in Figure 21, and the mesh size ratio between the field and the structure is 100 : 1. The underground structure model and load history are shown in Figures 11 and 18. In this section, the multiscale method is used to study the influence of the burial depth of the structure on the dynamic reaction of the structure. The accuracy and efficiency of the multiscale method are also further verified by the calculation in this section. The investigated burial depths of the structure are 50 m, 100 m, 200 m, 300 m, 400 m, 500 m, and 600 m, respectively. A closed artificial boundary substructure is used when the burial depth is at least 100 m.Figure 21
Large-range complex field with a length of 100 km, and the underground structure.Figure22 shows the comparison of the section bending moment calculated by the multiscale method and those calculated by the source-site-structure overall method for different burial depths. The results show that the values calculated by the multiscale method are in good agreement with those calculated by the source-site-structure overall model of a large range of complex fields.Figure 22
Moment history of the section for different depths: (a) section A (depth is 50 m), (b) section B (depth is 100 m), (c) section C (depth is 200 m), (d) section D (depth is 300 m), (e) section A (depth is 400 m), (f) section B (depth is 500 m), and (f) section A (depth is 600 m).
(a)
(b)
(c)
(d)
(e)
(f)
(g)Figure23 shows the relationship between the error and burial depth. The error refers to the difference between the maximum average displacement of the four corner points in the structure calculated by the multiscale model and those calculated by the source-site-structure overall model. For the closed substructure, the displacement error decreases gradually with increasing burial depth. Therefore, the wave energy is reduced after the complex wave field is reflected at the soil-structure interface and the ground surface. The deeper the underground structure is, the smaller the influence of the reflected wave on the structure is. The error is not obvious in the case of a burial depth greater than 100 m, so the closed substructure can be used safely in this situation.Figure 23
Relationship between the error and burial depth.Table2 shows a comparison of the calculation efficiency between the multiscale method and the source-site-structure overall method. The number of elements and the calculation time of the multiscale method are far lower than those of the source-site-structure overall method, and the advantages of using the multiscale method are obvious.Table 2
Comparison of the calculation efficiencies between two methods.
Example of large-scale site
Method of this paper
Source-site-structure overall calculation model
Source-site model
Artificial boundary substructure model
Soil-structure interaction model
Number of elements
2.001E + 5
160
1274
2.01E + 9
Calculation time
38 min
1 min
3 min
1860 minTable3 shows the specific data of the maximum node displacement and section bending moment at the four corners of the structure calculated by the multiscale method and the source-site-structure overall method under different burial depths. The error of the node displacement is controlled within 2%, and the error of the maximum section bending moment is not more than 3%, but in a few sections (A and C), the moment error is approximately 10% in the cases of 500 m and 600 m burial depths. Since the precision of the bending moment is second-order data, the error of the moment is larger than that of the displacement. Table 3 also shows that the displacement and the bending moment calculated by the multiscale method are generally larger than those calculated by the source-site-structure method. The reason for this difference may be that the artificial boundary elements of the small-scale SSI model cannot completely simulate the infinite domain space. A small reflection wave may exist in the scattering wave field, and the superposition of the reflection wave and the incident wave results in a larger structural reaction.Table 3
Maximum node displacement and bending moments of 4 positions of structure at different depths (N·m).
Depth(m)
Substructure type
Position
Node displacement (m)
Error (displacement) (%)
Section moment (N·m)
Error (moment) (%)
Source-site-structure calculation model
Method of this paper
Source-site-structure calculation model
Method of this paper
50
Open
A
9.786
9.801
0.16
−44844
−45302
1.02
B
9.805
9.818
0.13
−42397
−41847
−1.30
C
9.782
9.804
0.23
−46152
−45728
−0.92
D
9.801
9.821
0.21
−41177
−41597
1.02
100
Closed
A
9.827
9.934
1.08
−39796
−41856
1.40
B
9.822
9.927
1.07
−44640
−46171
0.73
C
9.824
9.932
1.10
−46400
−47758
0.26
D
9.818
9.925
1.09
−37440
-38777
2.71
200
Closed
A
9.702
9.833
1.35
−47532
−49242
3.60
B
9.654
9.785
1.36
−65236
−66380
1.76
C
9.699
9.838
1.43
−65344
−69613
6.53
D
9.651
9.790
1.44
−46495
−47766
2.73
300
Closed
A
9.338
9.475
1.47
60056
60706
1.08
B
9.255
9.390
1.45
−80187
−82536
2.93
C
9.337
9.480
1.52
−79990
−82029
2.55
D
9.256
9.395
1.50
60629
61573
1.56
400
Closed
A
8.783
8.900
1.34
75814
76362
0.72
B
8.673
8.787
1.31
−89635
−92216
2.88
C
8.786
8.906
1.37
−89308
−91961
2.97
D
8.678
8.794
1.34
76405
77297
1.17
500
Closed
A
8.087
8.169
1.02
87633
95427
8.90
B
7.959
8.037
0.98
−96458
−96644
0.19
C
8.090
8.174
1.04
−96151
−106515
10.78
D
7.962
8.041
0.99
88115
84800
3.76
600
Closed
A
7.303
7.348
0.62
97320
106536
9.47
B
7.184
7.211
0.38
−102589
−103490
0.88
C
7.304
7.349
0.62
−102144
−112777
10.41
D
7.181
7.212
0.42
97801
94402
3.48


## 4.1. Examples of a Homogeneous Flat Field

### 4.1.1. Example 1

The two-dimensional site is 8 km long and 4 km wide, the mass density is 2000 kg/m3, the shear wave velocity is 2000 m/s, and Poisson’s ratio is 0.3. The burial depth of the underground structure is 20 m. The mass density of C30 concrete is 2490 kg/m3, its elastic modulus is 30 GPa, and its Poisson’s ratio is 0.15. The earthquake source is 3.9 km from the surface. The mechanism of the source fracture is simulated by a vertical concentrated force. The model and load history are shown in Figure 11. The large-scale source-site model is composed of square plane strain elements of 10 m in width. The elements of the SSI model are refined into square plane strain elements 5 m in width. According to the beam element theory in finite elements, if the focus is on the overall force of the structure rather than the local force, beam elements can be used to simulate slender parts. In our two-dimensional models, the dimensions in out-of-plane direction can be ignored, so, the structure can be simulated by beam elements; it can also obtain the results such as axial force and bending moment conveniently which benefit the process in data analysis.Figure 11
Model of Example 4.1.1 (burial depth is 20 m): (a) the large-scale source-site model, (b) load history, (c) artificial boundary substructure model, (d) SSI model, and (e) underground structure model.
(a)
(b)
(c)
(d)
(e)Figures12 and 13 show that the node displacement and the bending moment of sections (A, B, C, and D) calculated by the multiscale method are in good agreement with those calculated by the source-site-structure overall model.Figure 12
Displacement history of the nodes of the artificial boundary elements: (a) node 1, (b) node 2, (c) node 3, and (d) node 4.
(a)
(b)
(c)
(d)Figure 13
Comparison of bending moments of the structural section (burial depth is 20 m): (a) section (A), (b) section (B), (c) section (C), and (d) section D.
(a)
(b)
(c)
(d)


### 4.1.2. Example 2

The structure is buried at a depth of 200 m, and the remaining parameters are the same as those in Section 4.1. The source-site model, the artificial boundary substructure, and the SSI model are shown in Figure14. As shown in Figure 15, the bending moments of sections A and C calculated by the multiscale method are in good agreement with those calculated by the source-site-structure overall model.Figure 14
Model of Example 4.1.2 (burial depth is 200 m): (a) the source-site model, (b) artificial boundary substructure, and (c) SSI model.
(a)
(b)
(c)Figure 15
Comparison of bending moment values of structural section (burial depth is 200 m): (a) section (A), (b) section C.
(a)
(b)With the increase in the burial depth of the structure, the size of the artificial boundary substructure in the vertical direction and the number of nodes and elements participating in the calculation increase, which is inconvenient to the modeling and calculation.


### 4.1.3. Example 3

An optimization is carried out on the artificial boundary substructure of the deeply buried structure. A closed substructure is adopted instead of the original open substructure, and the other calculation parameters remain unchanged. The closed substructure model is shown in Figure16. As shown in Figure 17, in the case of the closed artificial boundary substructure, the bending moments of sections A and C calculated by the multiscale method are in good agreement with those calculated by the source-site-structure overall model. At the same time, the number of elements and nodes participating in the calculation is significantly reduced, and the calculation efficiency significantly increases. When the structure is deeply buried, after the complex input wave field is reflected on the soil-structure interface and then the ground surface, the reflection wave transmits to the structure again but has a small effect on the structure.Figure 16
Model of Example 4.1.3 (burial depth is 200 m): (a) closed artificial boundary substructure, (b) SSI model.
(a)
(b)Figure 17
Comparison of the bending moment values of the structural section (burial depth is 200 m): (a) section (A), (b) section C.
(a)
(b)


## 4.1.1. Example 1

The two-dimensional site is 8 km long and 4 km wide, the mass density is 2000 kg/m3, the shear wave velocity is 2000 m/s, and Poisson’s ratio is 0.3. The burial depth of the underground structure is 20 m. The mass density of C30 concrete is 2490 kg/m3, its elastic modulus is 30 GPa, and its Poisson’s ratio is 0.15. The earthquake source is 3.9 km from the surface. The mechanism of the source fracture is simulated by a vertical concentrated force. The model and load history are shown in Figure 11. The large-scale source-site model is composed of square plane strain elements of 10 m in width. The elements of the SSI model are refined into square plane strain elements 5 m in width. According to the beam element theory in finite elements, if the focus is on the overall force of the structure rather than the local force, beam elements can be used to simulate slender parts. In our two-dimensional models, the dimensions in out-of-plane direction can be ignored, so, the structure can be simulated by beam elements; it can also obtain the results such as axial force and bending moment conveniently which benefit the process in data analysis.Figure 11
Model of Example 4.1.1 (burial depth is 20 m): (a) the large-scale source-site model, (b) load history, (c) artificial boundary substructure model, (d) SSI model, and (e) underground structure model.
(a)
(b)
(c)
(d)
(e)Figures12 and 13 show that the node displacement and the bending moment of sections (A, B, C, and D) calculated by the multiscale method are in good agreement with those calculated by the source-site-structure overall model.Figure 12
Displacement history of the nodes of the artificial boundary elements: (a) node 1, (b) node 2, (c) node 3, and (d) node 4.
(a)
(b)
(c)
(d)Figure 13
Comparison of bending moments of the structural section (burial depth is 20 m): (a) section (A), (b) section (B), (c) section (C), and (d) section D.
(a)
(b)
(c)
(d)


## 4.1.2. Example 2

The structure is buried at a depth of 200 m, and the remaining parameters are the same as those in Section 4.1. The source-site model, the artificial boundary substructure, and the SSI model are shown in Figure14. As shown in Figure 15, the bending moments of sections A and C calculated by the multiscale method are in good agreement with those calculated by the source-site-structure overall model.Figure 14
Model of Example 4.1.2 (burial depth is 200 m): (a) the source-site model, (b) artificial boundary substructure, and (c) SSI model.
(a)
(b)
(c)Figure 15
Comparison of bending moment values of structural section (burial depth is 200 m): (a) section (A), (b) section C.
(a)
(b)With the increase in the burial depth of the structure, the size of the artificial boundary substructure in the vertical direction and the number of nodes and elements participating in the calculation increase, which is inconvenient to the modeling and calculation.


## 4.1.3. Example 3

An optimization is carried out on the artificial boundary substructure of the deeply buried structure. A closed substructure is adopted instead of the original open substructure, and the other calculation parameters remain unchanged. The closed substructure model is shown in Figure16. As shown in Figure 17, in the case of the closed artificial boundary substructure, the bending moments of sections A and C calculated by the multiscale method are in good agreement with those calculated by the source-site-structure overall model. At the same time, the number of elements and nodes participating in the calculation is significantly reduced, and the calculation efficiency significantly increases. When the structure is deeply buried, after the complex input wave field is reflected on the soil-structure interface and then the ground surface, the reflection wave transmits to the structure again but has a small effect on the structure.Figure 16
Model of Example 4.1.3 (burial depth is 200 m): (a) closed artificial boundary substructure, (b) SSI model.
(a)
(b)Figure 17
Comparison of the bending moment values of the structural section (burial depth is 200 m): (a) section (A), (b) section C.
(a)
(b)


## 4.2. Examples of a Complex Layered Field

According to the complex site conditions, the layered field with mountains and rivers is selected as a calculation example. The calculation field is 10 km long and 5 km wide. The material parameters of medium 1 in the layered field are the same as those used in Section 4.1.1. For the material parameters of medium 2, the shear wave velocity is 3500 m/s, and the mass density is 2500 kg/m3. The burial depth of the structure is 300 m, and the remaining parameters are the same as those in Section 4.1.1. The calculation model and load history are shown in Figure 18. It should be noted that the artificial boundary substructures of Examples 4.2.1 and 4.2.3 do not contain medium 2. The artificial boundary substructures of Examples 4.2.2 and 4.2.4 contain medium 2. The four example models are shown in Figure 19.Figure 18
The complex layered source-site model and load history: (a) the layered source-site model, (b) load history.
(a)
(b)Figure 19
Four example models of complex layered fields: (a) Example 4.2.1, (b) Example 4.2.2, (c) Example 4.2.3, and (d) Example 4.2.4.
(a)
(b)
(c)
(d)Figure20 shows the comparison of the calculation results of the bending moments of the four examples with the dynamic calculation results of the source-site-structure overall model; under the condition of the complex layered field, the bending moments of sections A, B, C, and D calculated by the multiscale method are in good agreement with those calculated by the source-site-structure overall model whether the substructure is open or closed.Figure 20
Section bending moment history of the four examples of a complex layered field: (a) section (A), (b) section (B), (c) section (C), and (d) section D.
(a)
(b)
(c)
(d)Table1 shows the maximum bending moment values and differences calculated by the source-site-structure overall model and the four examples. The calculation accuracy of the model containing medium 2 (Examples 4.2.1 and 4.2.3) is slightly higher than that of the model that does not contain medium 2 (Examples 4.2.2 and 4.2.4). The reason for the error may be that Examples 4.4.2 and 4.2.4 consider the effect of the scattered wave energy on the structure after reflection by the layered interface, while Examples 4.2.1 and 4.2.3 are equivalent to extending medium 1 into infinite space. Since the scattered wave energy is much weaker than the incident wave energy, the error caused by ignoring the layered interface reflection is very small.Table 1
Maximum bending moments values and errors of the observed sections in 4.2.1–4.2.4 (N·m).
A
B
C
D
Source-site-structure calculation model
−47640
−56223
−55742
−47787
Example 4.2.1
Moment
−46029
−55968
−56387
−46313
Error
−3.38%
0.40%
0.29%
−3.08%
Example 4.2.2
Moment
−46758
−55475
−55770
−46985
Error
−1.85%
−0.48%
−0.81%
−1.68%
Example 4.2.3
Moment
−46710
−56538
−56975
−46897
Error
−1.95%
1.43%
1.34%
−1.86%
Example 4.2.4
Moment
−46904
−56706
−57000
−47135
Error
−1.55%
1.73%
1.38%
−1.36%


## 4.3. Example of a Large-Range and Complex Field with a Length of 100 km

According to the actual engineering problems, a large-range and complex site with a length of 100 km and a width of 20 km is selected for calculation. The shear wave velocity of the soil material is 3500 m/s, and the mass density is 2500 kg/m3. The earthquake source is located 10 km below the ground surface and 20 km from the structure. The mode of the field is shown in Figure 21, and the mesh size ratio between the field and the structure is 100 : 1. The underground structure model and load history are shown in Figures 11 and 18. In this section, the multiscale method is used to study the influence of the burial depth of the structure on the dynamic reaction of the structure. The accuracy and efficiency of the multiscale method are also further verified by the calculation in this section. The investigated burial depths of the structure are 50 m, 100 m, 200 m, 300 m, 400 m, 500 m, and 600 m, respectively. A closed artificial boundary substructure is used when the burial depth is at least 100 m.Figure 21
Large-range complex field with a length of 100 km, and the underground structure.Figure22 shows the comparison of the section bending moment calculated by the multiscale method and those calculated by the source-site-structure overall method for different burial depths. The results show that the values calculated by the multiscale method are in good agreement with those calculated by the source-site-structure overall model of a large range of complex fields.Figure 22
Moment history of the section for different depths: (a) section A (depth is 50 m), (b) section B (depth is 100 m), (c) section C (depth is 200 m), (d) section D (depth is 300 m), (e) section A (depth is 400 m), (f) section B (depth is 500 m), and (f) section A (depth is 600 m).
(a)
(b)
(c)
(d)
(e)
(f)
(g)Figure23 shows the relationship between the error and burial depth. The error refers to the difference between the maximum average displacement of the four corner points in the structure calculated by the multiscale model and those calculated by the source-site-structure overall model. For the closed substructure, the displacement error decreases gradually with increasing burial depth. Therefore, the wave energy is reduced after the complex wave field is reflected at the soil-structure interface and the ground surface. The deeper the underground structure is, the smaller the influence of the reflected wave on the structure is. The error is not obvious in the case of a burial depth greater than 100 m, so the closed substructure can be used safely in this situation.Figure 23
Relationship between the error and burial depth.Table2 shows a comparison of the calculation efficiency between the multiscale method and the source-site-structure overall method. The number of elements and the calculation time of the multiscale method are far lower than those of the source-site-structure overall method, and the advantages of using the multiscale method are obvious.Table 2
Comparison of the calculation efficiencies between two methods.
Example of large-scale site
Method of this paper
Source-site-structure overall calculation model
Source-site model
Artificial boundary substructure model
Soil-structure interaction model
Number of elements
2.001E + 5
160
1274
2.01E + 9
Calculation time
38 min
1 min
3 min
1860 minTable3 shows the specific data of the maximum node displacement and section bending moment at the four corners of the structure calculated by the multiscale method and the source-site-structure overall method under different burial depths. The error of the node displacement is controlled within 2%, and the error of the maximum section bending moment is not more than 3%, but in a few sections (A and C), the moment error is approximately 10% in the cases of 500 m and 600 m burial depths. Since the precision of the bending moment is second-order data, the error of the moment is larger than that of the displacement. Table 3 also shows that the displacement and the bending moment calculated by the multiscale method are generally larger than those calculated by the source-site-structure method. The reason for this difference may be that the artificial boundary elements of the small-scale SSI model cannot completely simulate the infinite domain space. A small reflection wave may exist in the scattering wave field, and the superposition of the reflection wave and the incident wave results in a larger structural reaction.Table 3
Maximum node displacement and bending moments of 4 positions of structure at different depths (N·m).
Depth(m)
Substructure type
Position
Node displacement (m)
Error (displacement) (%)
Section moment (N·m)
Error (moment) (%)
Source-site-structure calculation model
Method of this paper
Source-site-structure calculation model
Method of this paper
50
Open
A
9.786
9.801
0.16
−44844
−45302
1.02
B
9.805
9.818
0.13
−42397
−41847
−1.30
C
9.782
9.804
0.23
−46152
−45728
−0.92
D
9.801
9.821
0.21
−41177
−41597
1.02
100
Closed
A
9.827
9.934
1.08
−39796
−41856
1.40
B
9.822
9.927
1.07
−44640
−46171
0.73
C
9.824
9.932
1.10
−46400
−47758
0.26
D
9.818
9.925
1.09
−37440
-38777
2.71
200
Closed
A
9.702
9.833
1.35
−47532
−49242
3.60
B
9.654
9.785
1.36
−65236
−66380
1.76
C
9.699
9.838
1.43
−65344
−69613
6.53
D
9.651
9.790
1.44
−46495
−47766
2.73
300
Closed
A
9.338
9.475
1.47
60056
60706
1.08
B
9.255
9.390
1.45
−80187
−82536
2.93
C
9.337
9.480
1.52
−79990
−82029
2.55
D
9.256
9.395
1.50
60629
61573
1.56
400
Closed
A
8.783
8.900
1.34
75814
76362
0.72
B
8.673
8.787
1.31
−89635
−92216
2.88
C
8.786
8.906
1.37
−89308
−91961
2.97
D
8.678
8.794
1.34
76405
77297
1.17
500
Closed
A
8.087
8.169
1.02
87633
95427
8.90
B
7.959
8.037
0.98
−96458
−96644
0.19
C
8.090
8.174
1.04
−96151
−106515
10.78
D
7.962
8.041
0.99
88115
84800
3.76
600
Closed
A
7.303
7.348
0.62
97320
106536
9.47
B
7.184
7.211
0.38
−102589
−103490
0.88
C
7.304
7.349
0.62
−102144
−112777
10.41
D
7.181
7.212
0.42
97801
94402
3.48


## 5. Conclusion

In this study, based on the seismic wave input method for an artificial boundary, the efficient calculation of a large-scale complex field model is carried out by multiscale analysis with transitions created by using different meshes. The theoretical principle and implementation steps of the multiscale analysis method are presented. The efficiency and accuracy of the proposed method are verified by examples of homogeneous fields, complex layered fields, and large-scale actual engineering fields. The main advantages of this method are as follows:(1)
The multiscale method greatly reduces the number of elements of the source-site-structure model and the running time, significantly improving the work efficiency with good calculation accuracy. In this paper, the operation flow of the proposed method is simple. The displacement extraction from the nodes of a large-scale free field and the application of an equivalent seismic load can be realized automatically by using an automated program, and the process of modeling and calculation is completed quickly.(2)
By simulating the process of source fracturing, the complex wave field in the near-source region can be transformed into the ground motion field around the region where the structure is located. After recording the ground motion field, the SSI problem under the same topographic and geomorphic conditions can be solved iteratively. If the form or parameters of the structure are adjusted during the process of designation, it is not necessary to calculate the equivalent seismic load repeatedly. Only the SSI model needs to be analyzed, which can improve the work efficiency.(3)
In this paper, the processing method of the internal source can be used to consider the vertical and oblique incidences of a complex wave field, as well as the influence of special topography and multilayered media on the transmission of the complex wave field.



---

*Source: 8183272-2020-02-20.xml*