Time-Domain Analysis of Underground Station-Layered Soil Interaction Based on High-Order Doubly Asymptotic Transmitting Boundary

Based on the modified scale boundary finite element method and continued fraction solution, a high-order doubly asymptotic transmitting boundary (DATB) is derived and extended to the simulation of vector wave propagation in complex layered soils. The high-order DATB converges rapidly to the exact solution throughout the entire frequency range and its formulation is local in the time domain, possessing high accuracy and good efficiency. Combining with finite element method, a coupled model is constructed for time-domain analysis of underground station-layered soil interaction. The coupled model is divided into the near and far field by the truncated boundary, of which the near field is modelled by FEM while the far field is modelled by the high-order DATB. The coupled model is implemented in an open source finite element software, OpenSees, in which the DATB is employed as a super element. Numerical examples demonstrate that results of the coupled model are stable, accurate and efficient compared with those of the extended mesh model and the viscous-spring boundary model. Besides, it has also shown the fitness for long-time seismic response analysis of underground station-layered soil interaction. Therefore, it is believed that the coupled model could provide a new approach for seismic analysis of underground station-layered soil interaction and could be further developed for engineering.


Introduction
The fast development of urbanization has catalyzed for underground space exploitation, especially urban rail transit. The dynamic response of embedded structures in layered soil such as metro station and tunnels is greatly affected by soil-structure interaction (SSI) [Wolf (1989)], of which the radiation damping of the semi-infinite foundation is the main difficulty and key point of simulation. Extensive numerical approaches have been developed over the past 40 years to model the unbounded domain. Most of them could be summarized as approximate and rigorous methods. Approximate methods like the artificial boundary conditions [Dan (1999[Dan ( , 2004] are spatially and temporally local, which are computationally efficient but must be set with a quite distance from the interested region to obtain required accuracy [Lysmer (1969)]. Besides, most of the artificial boundaries are only singly asymptotic at high-frequency, resulting in a very high order of approximation for reasonable accuracy [Hagstrom and Hariharan (1998)] when confronting the evanescent wave decayed promptly in semi-infinite layered medium. While rigorous methods like the boundary element method (BEM) [Estorff and Hagen (2005); Luco and Barros (2010)], the thin layer method (TLM) [Kausel and Peek (1982); Kausel (2010)] and the scaled boundary finite element method (SBFEM) [Song and Wolf (1997); Wolf (2004)] are spatially and temporally global, thus computationally expensive. The SBFEM is a semi-analytical method, of which only the boundary is discretized and no fundamental solution is required to satisfy the damping condition of unbounded domain. Therefore, it is an appealing method for researching the unbounded domain. The application of the SBFEM could be found in Song (2010a, 2010b); Bransch and Lehmann (2011); Yaseri, Bazyar and Hataf (2014)], but it is still lack in efficiency because of globalization in the time domain, restricting the broad application of SBFEM. To balance well between accuracy and computational efficiency, there are many researches on the application of SBFEM in layered soils. Birk et al. [Birk and Behnke (2012)] modify the SBFEM in the case of layered soils by replacing the scale center with line and successfully apply it to the three-dimension problems. Chen et al. [Chen, Birk and Song (2015)] apply SBFEM to elastic wave propagation simulation in layered soil in the time domain. However, it is time-consuming because of convolutional integration solution. Based on the SBFEM, Bazyar et al. [Bazyar and Song (2008)] construct a highorder transmitting boundary using continued fraction solution, realizing time localization. Extra factor matrices are introduced by Birk et al. [Birk and Song (2009)] to improve the numerical stability, making it suitable for large-scale problems. However, it could not model evanescent waves because of singly asymptotic at high-frequency and only the propagating waves is absorbed. Hence, a high-order doubly asymptotic transmitting boundary (DATB) is proposed by Prempramote et al. [Prempramote, Song, Tin-Loi et al. (2010)] for scalar wave propagation in semi-infinite layered medium, in which continued fraction solution is expanded at both high and low frequency, obtaining dynamic stiffness solution closed to the exact solution. Wang et al. [Wang, Jin, Prempramote et al. (2011)] and Gao et al. [Gao, Jin, Wang et al. (2011)] expand this boundary to unbounded reservoir and apply it with FEM to the seismic response analysis of dam-reservoir systems. Furthermore, the high-order DATB is applied for problem of vector wave propagation in single layered soil [Prempramote (2011)]. However, it does not go further application on complicated layers situation. In all, lots of researches have concentrated on the solution in the frequency domain and it is still relied on time-consuming convolutional integration solution in the time domain, resulting in low efficiency. Besides, it is rarely reported the application of high-order DATB on time-domain analysis of underground station-layered soil interaction, especially in multiple layered soils. In this paper, the excellent high-order DATB is extended to simulation of the vector wave propagation in the complicated layered soils. Combining the high-order DATB with the finite element method (FEM), a coupled model is constructed for the seismic response analysis of the realistic complicated underground station-layered soil systems. The coupled model is numerically implemented in OpenSees to study the seismic response of the underground station-layered soil systems. OpenSees, the open system for earthquake engineering simulation (http://opensees.berkeley.edu), is an open source general purpose finite element software for modeling structural and geotechnical systems. The accuracy and efficiency of the coupled model are validated by numerical examples.

Summary of SBFEM of semi-infinite layers with constant depth
Underground station is usually buried in the complicated layered soils. Thus, the key equations of the SBFEM for vector wave problem in layered soils [Prempramote, Song, Tin-Loi et al. (2010)] are briefly presented in this section. As shown in Fig. 1, elastic semi-infinite layers with constant depth are assumed to be isotropic. As a special case, the scaling center of semi-infinite layers in the SBFEM is located at infinity [Li, Cheng, Deeks et al. (2005)]. The model is divided into the near field and far field by the truncated boundary. One-dimensional line element is brought in to discretize the truncated boundary B v shown in Fig. 1(a) and a typical line element is shown in Fig. 1(b). Besides, free boundary condition is imposed on the upper boundary Bu and fixed boundary condition is imposed on the bottom boundary Bb.
(a) geometry and discretization (b) a typical line element where the coefficient matrices are independent of ξ and defined in the Prempramote [Prempramote (2011)].

The dynamic stiffness matrix of an unbounded domain [ ( )]
S ω ∞ (superscript ∞ denotes for unbounded domain) is defined by the excitation force-displacement relationship as where { } R is nodal forces on the truncated boundary.
(2) can be transformed into the formulation in dynamic stiffness as where ω is the excitation frequency. Based on Eq. (4), the high-order DATB could be constructed with the help of continued fraction solution and would be described in Section 3.

High-order doubly asymptotic transmitting boundary
In this section, continued fraction solution is applied in Eq. (4) to approximate the dynamic stiffness matrix of unbounded domain at both high-frequency and lowfrequency range. The transmitting boundary condition is constructed using the continued fraction solution of dynamic stiffness matrix and the auxiliary variable technique. The solution process of the continued fraction expansion at both high frequency and low frequency range can be found in Prempramote [Prempramote (2011)], only the key equations are presented.
To start with, the continued fraction solution is firstly determined at the high-frequency limit ( ω → ∞ ) and is expressed as To determine a valid solution over the whole frequency range, the continued fraction solution for the residual of the high frequency range is constructed at the low-frequency limit ( 0 Y ω at the low-frequency limit is expressed as are coefficient matrices to be calculated recursively. The factor matrix ( ) [ ] i L X is introduced to improve numerical stability of the continued fraction solution. L M is the order of the continued fraction solution at low frequency range. Substituting continued fraction solution of dynamic stiffness matrix (Eqs. (5)-(6)) into Eq.
(3) as well as introducing auxiliary variables, the high-order DATB in the frequency domain is constructed. Apply the inverse Fourier transform, this transmitting boundary is formulated in the time domain as: where [ ]    .
where subscript "A" represents the auxiliary variables. Eq. (7) is a first order ordinary differential equation and can be solved using the standard time integration scheme, such as Newmark method. Thus, the time-consuming convolutional integration is avoided in the transient analysis and the computational efficiency is improved.

Coupled model and implementation in OpenSees
Employing the same grids and shape functions at the common boundary (i.e., Bv), the high-order doubly asymptotic transmitting boundary can be seamlessly coupled with FEM.  Combining Eqs. (7) and (10), the coupled model could be built. As the coefficient matrices of Eq. (8) result from the repeated inversion in the continued fraction solution, possible spurious modes may exist, resulting in the instability of Eq. (10). To guarantee the numerical stability, a spectrum shifting technique is proposed by Birk et al. [Birk, Prempramote and Song (2012) Eq. (11) has the same characteristics as the finite element equation, i.e., the coefficient matrices are symmetric and sparse. This coupled model is implemented in the open source finite element software OpenSees. The high-order DATB is employed as a super element in OpenSees. As the coeffect matrices of the super element are constant, they can be solved only once by MATLAB and read into OpenSees during the analysis. The main solution procedure of the coefficient matrices in MATLAB are summarized as Fig. 2. Moreover, the TCL (Tool Command Language) script for the input of the DATB super element is also generated by MATLAB. As the coefficient matrices of the DATB super element is solved and stored in binary file before the transient analysis, the DATB super element in OpenSees simply reads in the coefficient matrices according to the analysis procedure and thus computationally efficient.  [Liu, Du, Du et al. (2006)] (VSB) is also employed to make a comparison with DATB. All the methods are conducted through using the open source finite element software OpenSees. In Section 5.1, surface traction is applied at double semi-infinite layered soil. And in Section 5.2, a metro station buried in double layered soil subjected to El-Centro seismic wave is analyzed. Moreover, in Section 5.3, the calculation efficiency is discussed among numerical examples.

Semi-infinite layered soil under surface traction of triangular impulse
In this example, the double semi-infinite layered soil subjected to surface traction in the inclined direction of 45° is illustrated in Fig. 3. The top of the soil is free boundary condition and the bottom of the soil is fixed. Both have the height of 0.5. The material properties of soils are presented in Fig. 3. Standard four-node quadrilateral element in OpenSees is selected to model near and far field for plane-strain condition. The near field needs to be discretized with 100 quadrilateral elements (size of element is 0.1×0.1). As for the boundary between near and far field, 10 two-node line elements of equal length are used for discretization in the DATB model. While for the EMM, total discretization area is 80×1 so that responses are not affected by the waves reflected at the truncated boundary under correspondingly analyzing time. 8000 quadrilateral elements (size of element is 0.1×0.1) are used for discretization. The time history of triangle impulse ( ) p t is plotted in Fig. 4 and the maximum surface traction is denoted as P . The total analyzing time is 30 and the dimensionless time step t ∆ is 0.1. . The displacement responses at Point A and B are plotted in Figs. 5 and 6, respectively. Excellent agreement is shown between the response of the DATB and the EMM through the whole calculating time, demonstrating the high accuracy of the DATB. As for the VSB, unsatisfactory results can be seen in both Point A and B. Besides, the results of Point B deviate much more from the EMM than those of Point A, manifesting a significant error at the boundary of the VSB. It is indicated the deficiency of the VSB in the long-time dynamic analysis.

Metro station embedded in semi-infinite layered soil
In Section 5.2, a metro station embedded in semi-infinite layered soil is studied under El-Centro wave. Two layered soils are selected and the metro station is chosen as doublelayered island platform station. As shown in Fig. 7, the physical dimension of the near field model is 66 m×55 m and the metro station is horizontally located in the middle of near field. The top of the soil (the ground surface) is free boundary condition and the bottom of the soil is fixed. Three observant Point A, B and C are selected for displacement. As depicted in Fig. 8, the metro station is mainly consisted of the roof slab (RS), the middle slab (MS), the base slab (BS), the lateral wall (LW) and the columns. The thickness of RS, MS, BS and LW is 0.8 m, 0.5 m, 0.9 m and 0.8 m, respectively. The columns are spatial in practical and the principle of equivalent stiffness is introduced considering the plane-strain condition. Therefore, the width of the columns is equivalent to 0.4 m. The near field soils are modeled by 3300 quadrilateral elements (size of elements is 1 m×1 m). Another 55 two-node line elements are used for discretizing the vertical boundary to model far field in the DATB. As for the EMM, the entire discretization area is determined as 3000 m×55 m. Totally 84171 elements are used to model near field soil as well as underground station in the EMM. 186 quadrilateral elements are applied to model the metro station. The metro station is assumed to deform together with the soils during the seismic as one kind of the soil-structure interaction (SSI). Therefore, the common girds are shared by station and soils in modeling. The mesh model of the near field is illustrated in Fig. 9. . It can be observed from three figures that excellent agreement exists between the DATB and the EMM in both horizontal and vertical displacement among three observant points. However, the results of the VSB are less accurate compared with the DATB. Similar as the observations from Section 5.1, the results of the boundary (Point C) as plotted in Fig. 13 have a much more deviation than those of the inside points (A and B).
Through the examples, it can be concluded that the DATB is much more accurate in modeling semi-infinite layered medium as well as solving time-domain problem than the VSB. The excellent agreement in Section 5.2 between the DATB and the EMM indicates that the DATB might be applicable to practical engineering project.   (1) and (2) represent the example in Section 5.1 and metro station example in Section 5.2. For numerical example Ⅰ as well as Ⅱ, the calculation time of the DATB is close to that of the VSB, showing a good efficiency while ensuring a high accuracy. Therefore, it is shown that high-order doubly transmitting boundary balances well between accuracy with efficiency as it can greatly improve the accuracy while does not spend too much time. So far, the high-order doubly asymptotic boundary has been proved accurate and efficient. It is demonstrated that this coupled model is well suitable for long-time analysis of underground station-layered soil interaction in the time domain.

Conclusions
The high-order DATB is extended to the simulation of elastic wave propagation in semiinfinite multiply layered soils. This transmitting boundary converges to the exact solution rapidly throughout the entire frequency range. It is also demonstrated that the coupled model is suitable for long-time analysis of underground station-layered soil interaction. Besides, the coupled model is achieved through OpenSees and can be further developed for engineering application.