An Arc-Consistent Viscous-Spring Artificial Boundary for Numerical Analysis of Seismic Response of Underground Structures

. A new arc-consistent viscous-spring artiﬁcial boundary (ACVAB) was proposed by changing a traditional ﬂat artiﬁcial boundary based on the theory of viscous-spring artiﬁcial boundaries. Through examples, the concept underpinning the establishment and speciﬁc setting of the boundary in the ﬁnite element software were described. Through comparison with other commonly used artiﬁcial boundaries in an example for near-ﬁeld wave analysis using the two-dimensional (2D) half-space model, the reliability of the ACVAB was veriﬁed. Furthermore, the ACVAB was used in the numerical analysis of the eﬀects of an earthquake on underground structures. The results were compared with the shaking table test results on underground structures. On this basis, the applicability of the ACVAB to a numerical model of seismic response of underground structures was evaluated. The results show that the boundary is superior to common viscous-spring boundaries in terms of accuracy and stability, and therefore, it can be used to evaluate radiation damping eﬀects of seismic response of underground structures and is easier to use.


Introduction
Given increasing urbanisation, underground rail transit has developed rapidly.Numerical simulation, as one of the main methods to study the seismic performance of underground structures at present, can be used to reproduce dynamic characteristics of underground structures under seismic load.
e infinite domain (infinite foundation) can be simulated by introducing an artificial boundary [1][2][3][4][5][6][7][8], so the artificial boundary theory is one of the important bases for seismic analysis and numerical model study of underground structures.In the numerical simulation of near-field wave motion, artificial boundary conditions directly affect the accuracy of the calculations and a reasonable artificial boundary can simulate the propagation process of waves in an infinite domain medium.An ideal artificial boundary should be stable, accurate, efficient, and easy to realise, but there are often contradictions among the above characteristics.e generalised artificial boundaries are classified into global artificial boundaries and local artificial boundaries.Global artificial boundaries are mainly characterised by coupled motion of all boundary nodes in space and time and satisfaction of all field equations and mathematical and physical conditions in the infinite domain, which can accurately simulate the infinite domain [9][10][11].In the analysis of large-scale complex wave motion, especially nonlinear generalised structures, global artificial boundaries impose an onerous computational burden.
e artificial boundaries introduced above are mostly expressed in the form of differential equations of boundaries by an approximate treatment of wave equations.However, in the combination with the explicit finite element method, the discrete form is inconvenient to engineering design and there is no obvious advantage in its simulation accuracy.By referring to other studies [36][37][38], several boundaries, such as viscous boundary, viscous-spring boundary, consistent viscous-spring boundary elements, and remote artificial boundary, are discrete artificial boundary conditions established using approximate solutions to one-sided wave equations.
ey can be easily realised by the appropriate finite element analysis software.
Alterman et al. [36,37] performed numerical simulation on near-field wave motion and compared it with concise analytical solutions of integral transformation based on the remote artificial boundary.
e principle of the remote boundary is that for the problem of near-field wave motion, when the generalised structure is far enough from the artificial boundary, the artificial boundary problem can be ignored; that is, the artificial boundary can be set arbitrarily.
is processing method of the remote artificial boundary has been still used to provide accurate solutions hitherto.However, this method runs contrary to the goal of high efficiency and is invalid in current engineering design and calculation.e viscous boundary [33,34] was first studied by Lysmer et al., whose working principle [38,39] is derived based on the assumption that there is no energy reflection at the boundary.It is found from the research that the viscous boundary does not permit attenuation of waves in the medium and then does not consider the stiffness recovery of the medium, which is likely to raise problems of structural drift and low accuracy.Deek and Randolph (1994) [35] were the first to propose the expression of the two-dimensional (2D) plane strain problem of the viscous-spring boundary.To cope with the problems of structural drift and low accuracy in the use of viscous boundaries, Liu (2006) [38,40,41] proposed the concept of consistent viscousspring artificial boundary elements based on the cylindrical wave equation.A viscous-spring artificial boundary with good elastic recovery performance was established using discrete spring dampers at the boundary, which overcame the problem of low-frequency instability, indeed resulting in good stability [40][41][42].At present, the boundary has been widely utilised in finite element analysis software; however, when solving the dynamic problem of the complex 2D semiinfinite domain, the section of the viscous-spring boundary is not directly introduced into the most finite element analysis software.For this reason, researchers need to calculate the parameters of spring dampers and apply spring dampers one by one, requiring significant time in preprocessing.Moreover, in the derivation of the viscous-spring artificial boundary, the cylindrical wave assumption is adopted, but the boundary is simplified to a flat boundary, which is inconsistent with reality.
e concentrated viscous-spring dynamic artificial boundary requires onerous preprocessing and is of low accuracy and inconsistent with the cylindrical wave assumption.In view of these problems, a new artificial boundary, namely a 2D arc-consistent viscous-spring artificial boundary (ACVAB) element, was established based on the theory of consistent viscous-spring artificial boundary and viscous-spring artificial boundary element [38,40].Section 1 introduces the concept and realisation method of the ACVAB.Section 2 introduces the methods of setting equivalent stiffness and damping.In Section 3, the example in near-field wave analysis using the 2D half-space model and example in seismic response of underground structures are described.Section 4 discusses the results of examples and verifies the accuracy of the ACVAB and its applicability to the analysis of the seismic response of underground structures.e conclusions are summarised in Section 5.

Description of the Equivalent Viscous-Spring Artificial
Boundary.In the 2D arc-consistent viscous-spring artificial boundary element, a global stiffness matrix and a damping matrix of the model are derived, as shown in formulas ( 1) and (2) [35].In the finite element analysis software, the consistent viscous-spring artificial boundary is realised by extending a layer of the same type of elements along the normal to the outer surface of the established model area, fixing nodes on the outer element surface, and setting springs and dampers.Based on the equivalence principle of the matrix, the consistent stiffness and damping matrixes of the element are transformed into an equivalent stiffness matrix and an equivalent damping matrix of the boundary element, thus meeting the requirements of resilience and damping on the boundary (Figure 1).
where K BT , K BN , C BT , and C BN represent the tangential stiffness, normal spring stiffness, tangential damping coefficient, and normal damping coefficient, respectively.

Description of the ACVAB.
When processing the artificial boundary in a semi-infinite domain, to simplify the 2 Shock and Vibration model, previous researchers have mostly adopted the flat artificial boundary (Figure 1).In the derivation of the viscous-spring artificial boundary, the cylindrical wave assumption is used, so it is inconsistent with actual situations when simplifying the boundary to be flat.erefore, a new viscous-spring boundary processing method, that is, the ACVAB, was proposed in this study.
In the specific setting, firstly, a finite element model of an arc boundary of radius R is established taking the point of action of the scattered wave source as the centre, so that the distance from the wave source to the arc boundary remains unchanged.Secondly, a layer of ordinary elements can be extended along the normal direction on the boundary of the established finite element analysis model and the outer boundary is fixed.An appropriate material coefficient is given through calculation to replace the original discrete spring dampers.e calculation and setting of the coefficient are described in the following section, and the model is displayed in Figure 2.
is method does not need to calculate and arrange spring dampers one by one, which greatly simplifies the preprocessing work.Meanwhile, because the radius of the model is fixed, the error caused by the approximation of the coefficient of material property is avoided, and the accuracy of the dynamic finite element analysis is thus improved.

Setting of Equivalent Stiffness and
Damping of the ACVAB e artificial boundary introduced in the analysis of dynamic interaction between structures and foundation is generally deduced based on the assumption that there is no energy reflected at the boundary, such as the viscous boundary [33].Based on the assumption that 2D scattered waves are cylindrical, Deeks [35] deduced the artificial boundary condition as follows: where r b represents the coordinates of the artificial boundary in a polar coordinate system; C b and K b are the viscous damper and linear spring applied on the artificial boundary r � r b , respectively; c s � ��� G/ρ  denotes the shear wave velocity; and G and ρ separately denote the shear modulus and mass density of the medium.
Based on the cylindrical wave equation, Liu et al. [38,[40][41][42] established a 2D viscous-spring artificial boundary, in which the spring coefficient K B and damping coefficient C B of the equivalent physical system separately are expressed as follows: Tangential boundary: Normal boundary: where K BN and K BT refer to the normal and tangential stiffness of the spring; R represents the distance from the wave source to the artificial boundary point; c s and c P denote the S-wave and P-wave velocities; and α T and α N denote the parameters of the tangential and normal viscous-spring artificial boundaries, respectively.As shown in Figure 2, the artificial boundary element B consists of continuously and uniformly distributed springs and dampers.As illustrated in Figure 3, the displacement shape function of the two-node artificial boundary element can be expressed by the linear interpolation function.e springs distributed in the normal direction of element B and the corresponding displacement shape function are shown in Figure 4.
Assuming u and v represent the tangential and normal displacements of the artificial boundary, the displacements of nodes i and j of the artificial boundary element are expressed as follows:

Shock and Vibration
As shown in Figure 4, the normal displacement shape function of the artificial boundary element is expressed as follows: e tangential displacement shape function is the same as the normal displacement shape function.
e geometry matrix of the artificial boundary element corresponding to formula (3) is shown as follows: e elastic stiffness matrix of springs distributed on the viscous-spring artificial boundary is expressed as follows: e stiffness matrix of the viscous-spring artificial boundary element can be derived based on the principle of virtual work as follows: By substituting formulas ( 8) and ( 9) into formula (10), the stiffness matrix of the consistent viscous-spring artificial boundary element is obtained: rough the same steps to form the stiffness matrix of the consistent viscous-spring boundary element, the damping matrix of the 2D arc-consistent viscous-spring boundary element can be attained: Formulas ( 11) and ( 12) demonstrate the stiffness matrix and damping matrix of the 2D arc-consistent viscous-spring artificial boundary element.When the arc viscous-spring boundary is segmented into infinite elements with the minimum thickness (extremely dense meshing), each element approximates a rectangular element.A four-node rectangular element is used to simulate the 2D ACVAB.As shown in Figure 5, the displacement shape function of the element is described as follows: In formula (12), l and h represent the side length and the extended thickness of the element, respectively; i, j, k, and l represent the nodes of the element, of which i and j indicate the nodes of the aforementioned consistent viscous-spring artificial boundary: because nodes l and k of the element are fixed, the explicit expression of the stiffness matrix of plane 4 Shock and Vibration strain can be found, thus giving the equivalent stiffness matrix of the element corresponding to two nodes i and j： where  ρ represents the mass density of the equivalent element;  c s and  c p indicate the S-wave and P-wave velocities of the medium of the equivalent element.
When the thickness of the artificial boundary zone is h ≈ 0, the stiffness matrix is approximated as follows: By comparing formulas (11) and (15), according to the following definitions of S-wave and P-wave velocities, the equivalent shear modulus and elastic modulus of the boundary element are obtained as follows: where  G,  E, and  v denote the equivalent shear modulus, equivalent elastic modulus, and equivalent Poisson's ratio of the equivalent viscous-spring boundary element, respectively; h, R, and G refer to the thickness of the equivalent element, distance from source of the wave to the artificial boundary point, and shear modulus of the medium, respectively.Shock and Vibration e damping matrix of the equivalent boundary element is formed using damping proportional to the stiffness, and then, let en, where  η BT and  η BN separately represent the scale factors related to tangential and normal stiffness.By substituting formulas (11) and (12) into formula (19), the following formula is obtained: To realise the boundary element in a convenient manner, the damping scale coefficient of the equivalent element is set as the average of the parameters in two directions.

Example for Near-Field Wave Analysis Using the 2D Half-Space Model.
To verify the accuracy of the ACVAB proposed in this study, an example of the near-field wave analysis using the 2D half-space model was compared with other commonly used artificial boundaries.When modelling the ACVAB, a semicircular near-field observation area with a radius of 50 m was taken, a layer of elements with the unit thickness was extended around the model, and the nodes on the outer layer were fixed.e medium of elastic semiinfinite domain was selected for research, and the density of the medium, Young's modulus, and Poisson's ratio are 2000 kg/m 3 , 20 MPa, and 0.3, respectively.Considering the Rayleigh damping of the medium, a � 0.616 and β � 0.000312.e coefficient of material property of the arcconsistent viscous-spring element is calculated by referring to formulas (19) and (21).Half of a sinusoidal pulse at the frequency of 1 Hz was applied at the centre, and surface point A at the axis of symmetry and surface point B at a burial depth of 10 m were taken as points for observation.
e 2D half-space model for near-field wave analysis is shown in Figure 6.To verify the effects of the ACVAB, multiple models including the concentrated viscous boundary, ACVAB, concentrated viscous-spring boundary, fixed boundary, and remote boundary were also established for comparisons.

Example for Seismic Response of Underground Structures.
To assess the applicability of the ACVAB proposed in this research to the seismic analysis of underground structures, a numerical model for seismic response of underground structures was built based on the ACVAB.e model results were compared with data from shaking table tests [43], thus allowing the analyst to determine whether the boundary is applicable to the analysis of radiation damping effects in the seismic response of underground structures.

Description of Project Background.
e geomorphic unit of Xi'an Feitian Road Station belongs to the second-and third-grade loess tableland.e lithologic characteristics of the site strata are summarised in Tables 1 and 2. e table is arranged from top to bottom according to the order of soil layers from shallow to deep.
e underground subway station is a reinforced concrete structure with a total height of 14.01 m and a total width of 19.2 m. e longitudinal spacing of the centre pillar is 9 m. e cross section of the centre pillar measures 0.8 m × 1.2 m. e depth of soil above the roof is 3.459 m. e density of concrete ρ is 2.5 g/cm 3 , its modulus of elasticity E is 35 GPa, and Poisson's ratio ] is 0.15.A typical cross section is shown in Figure 7.
e station is located in Chang'an District of Xi'an City.e site type in this area is class II, the basic seismic intensity is 8 degrees, and the characteristic period of seismic response spectrum is 0.4 s.

Description of Calculation Model.
Due to the length of Feitian Station, the three-dimensional dynamic interaction system of loess subway underground structure was considered as a two-dimensional plane strain problem in the seismic analysis by the time history analysis method.e typical section of loess site and subway underground structure was numerically simulated, and the model was established as shown in Figure 8. e ACVAB proposed in this study is adopted as the artificial boundary of calculation model.TAFT waves are chosen as the incident seismic waves, and the seismic fluctuation input problem is transformed into a wave source problem to realise the seismic wave input [41].e plastic damage model of concrete in ABAQUS finite element analysis software was adopted to simulate the mechanical behavior of the prototype subway station concrete.e CPE4R element was used for on-site soil, and the CPE4 element was used to model the subway station structure.e initial stress on the soil was calculated using the geostatic module in ABAQUS.To simulate the in situ stress on the soil and how it affects the adjacent underground structure, the model states of excavation, support, construction of underground structure, and backfilling were established to simulate the construction process of this subway station, and the stress state in the soil after construction was taken as the initial stress field for the later dynamic analysis.

Comparison and Discussion of Results of the Example for Near-Field Wave Analysis Using the 2D Half-Space Model.
To validate the accuracy of the 2D ACVAB, the displacements of observation points on each boundary model were compared and analysed.e remote boundary with complex modelling, low calculation efficiency, and optimal simulation effects was taken as the standard solution to seek for the artificial boundary with accuracy approaching the remote boundary and high modelling efficiency and calculation efficiency.As shown in Figure 9, due to reflection-induced interference, the fixed boundary model differs greatly from the standard solution obtained by the remote boundary and cannot reflect the actual propagation of waves in the semiinfinite domain medium.In the model of the concentrated viscous-spring boundary, the reflection of waves is reduced, but there is a significant downward drift in the middle and rear sections, so it is considered that its accuracy cannot meet the requirements of engineering practice.For the models of viscous boundary, ACVAB, and remote boundary, after inputting pulse waves, energy rapidly decreases, showing only small fluctuations and being almost completely absorbed by the boundary element, which meets the accuracy requirements imposed by engineering practice.rough comparison, the advantages and disadvantages of the commonly used artificial boundaries and the ACVAB proposed in this study are summarised in Table 3.
To compare the accuracy of artificial boundaries, Table 4 is obtained by dividing the displacement of observation points of each boundary model by the displacement of observation points on the remote boundary (standard solutions).After processing the data in Table 4, the factors of accuracy improvement using the proposed ACVAB compared with other commonly used artificial boundaries can be determined, as listed in Table 5.By combining this with the data in Tables 4 and 5, the accuracy of the model of the ACVAB increases with the increase in the burial depth.For surface observation point A (used for monitoring displacement), the accuracy of the ACVAB is improved by 0.6 to 3.4% compared with other commonly used artificial boundaries.In terms of observation point B at the deep part, its accuracy increases by 3.1 to 16.7% in comparison with other commonly used artificial boundaries.From the perspective of the mean value, its accuracy increases by 3.5 to 15.4% on average.In conclusion, it is considered that the calculation accuracy and stability of the ACVAB are superior to those of the viscous boundary and concentrated viscousspring boundary.e accuracy is closer to that of the remote boundary, and the setting method is simpler.It is easier to realise this in the finite element analysis software.

Comparison and Discussion of Results of Example for Seismic Response of Underground Structures.
e numerical simulation results of loess subway stations based on the ACVAB were compared with shaking table test results [43].Under the action of seismic waves, the acceleration response time history of the foundation in the model and corresponding Fourier spectra in the numerical simulation and shaking table test were compared, as shown in Figure 10.e comparison of acceleration time history of the structure and corresponding Fourier spectra is shown in Figure 11.
As illustrated in Figure 10, the waveforms and amplitudes of acceleration response time history of the foundation in the model and composition characteristics of Fourier spectra recorded in the numerical simulation and shaking table test are coincided.At a low acceleration, the numerical simulation and test results are consistent in terms of waveforms and amplitudes of acceleration response time history of the foundation in the model and composition characteristics of Fourier spectra; however, as the acceleration increases, the simulated accelerations of the foundation show certain deviation from the test result, with the error within 15%. is result indicates that the proposed ACVAB can be used to simulate the propagation of seismic waves in the soil medium of the semi-infinite domain, which well solves the transmission of ground motion energy from near-field to far-field areas.
As demonstrated in Figure 11, the waveforms and amplitudes of acceleration response time history of the structure in the model and composition characteristics of Fourier spectra recorded in the numerical simulation and shaking table test are similar.When the acceleration is low, consistent waveforms and amplitudes of acceleration response time history of the structure in the model and composition characteristics of Fourier spectra are found in the numerical simulation and test data.As the acceleration rises, the numerical simulation shows a certain error (within 25%) from the test results in terms of acceleration, which is related to the model of structural materials and accuracy of processing method for contact between the structure and foundation.Shock and Vibration

Conclusions
e concept of a new ACVAB was proposed based on the theory of the consistent viscous-spring artificial boundary and viscous-spring artificial boundary element.e establishment ideas of the artificial boundary and the implementation method in the finite element analysis software were described, and its applicability was studied by analysing numerical examples.
(1) Considering the error caused by the approximation in the setting of the concentrated viscous-spring boundary and the cumbersome preprocessing work, the new artificial boundary, namely the ACVAB, was proposed in this research.e operation of this element model is simpler than that of the concentrated viscous-spring boundary, and its high accuracy and good stability were proven by example.(2) rough the example of near-field wave analysis using the 2D half-space model, the ACVAB was compared with other commonly used artificial boundaries.e results demonstrate that the boundary with high accuracy and simple setting method is easy to realise in the finite element analysis software and improves both modelling efficiency and calculation speed.

Figure 3 :Figure 4 :
Figure 3: Details of the element (the length l of the element is minimised).

Figure 8 :
Figure 8: Schematic representation of the time history analysis calculation model.

Figure 9 :
Figure 9: Displacement time history of observation points: (a) Point A and (b) Point B.

Table 2 :
Soil properties from shallow to deep.

Table 1 :
Physical and mechanical properties of the soil layer.

Table 3 :
Advantages and disadvantages of artificial boundary.

Table 4 :
Comparison of fitting accuracy of time history curve.

Table 5 :
Factors of accuracy improvement of ACVAB.