Abstract
In the seismic analysis and design of the underground structure, the response displacement method, as a pseudostatic method, has been widely adopted for its solid theoretical background, clear physical concept, and ease of implementation. The subgrade modulus is an essential parameter to the response displacement method, and a few approaches are available to determine its value. However, the existing methods neglect the interaction between the radial and tangential subgrade modulus and the influence of actual ground deformation, resulting in an inaccurate estimation. This study presents a solution to overcome these defects for the response displacement method adopted in the transverse seismic analysis of the shield tunnel with a circular cross section. First, the analytical solutions of subgrade modulus for ground deformation modes described by the Fourier series are derived based on the theory of elasticity. The ratio of the radial displacement to tangential displacement is introduced to create a link between the radial and tangential subgrade modulus. Based on the solutions of subgrade modulus for different ground deformation modes, the displacement fitting method is proposed to derive the subgrade modulus corresponding to the actual ground deformation. With this method, the subgrade modulus would adjust according to the ground displacement. Finally, a case study is conducted to illustrate the validity of the displacement fitting method.
1. Introduction
During the last decades, progress has been made in studying seismic design and analysis of underground structures [1, 2]. Various methods related to seismic analysis of underground structure in transverse direction have been proposed, like the dynamic time-history analysis method, seismic coefficient method [3], BART method [4], free field deformation method [5], and response displacement method [6]. Among these methods, the response displacement method has been widely recommended and adopted for its solid theoretical background, clear physical concept, and ease of implementation [7–9]. However, the subgrade modulus, an essential parameter, has not been thoroughly studied, deteriorating the accuracy of the response displacement method.
Based on the observation of dynamic behaviors of an underground tank, undersea tunnels, and a rock tunnel during earthquakes, it has become clear that the deformation of underground structures is governed by strain in the surrounding ground. The response displacement method has been proposed based on this observational result [6, 10]. When adopting this method in the transverse seismic analysis of the shield tunnel, the ground springs are usually set up to represent the ground. One end of the ground spring is attached to the tunnel structure, and the ground displacement under the seismic load is applied at the other end of the ground spring. So, the ground spring stiffness reflects the deformability of the ground and the ability to transfer the ground load to the tunnel structure, which makes its value a critical parameter for the response displacement method.
The stiffness of the ground spring depends on the subgrade modulus, and a few approaches have been presented to estimate the subgrade modulus. Analytical solutions for the tunnel with a circular cross section are available. Kawashima [6] offered the solutions of radial and tangential subgrade modulus for different ground deformation modes. Using the elastic theory method, Huang and Cao [11] derived the radial subgrade modulus under uniform radial ground deformation and the tangential subgrade modulus under oval ground deformation. The static finite element method has also been put forward and recommended to estimate the subgrade modulus, suitable for the tunnel with a circular or a rectangular cross section [7]. Li et al. [12] optimized the strategy by which the unit force acts on the tunnel perimeter to improve the efficiency of the static finite element method. It can be found that the prescribed boundary condition is required for the analytical and numerical approach. However, the actual deformation or stress at the boundary may not be consistent with the predefined condition. Besides, the interaction between the radial and tangential subgrade modulus has not been studied. These defects may lead to an inaccurate estimation of the subgrade modulus.
In order to overcome the difficulty in estimating the subgrade modulus, some improved response displacement methods have been proposed. Liu et al. [13, 14] introduced an integral response deformation method in which the ground is modeled by 2D or 3D finite solid elements. Xu et al. [15] incorporated substructure analysis method into the response displacement method and regarded the tunnel and part of surrounding ground as the research object. In these methods, the ground is no longer simplified as the ground spring, which avoids the estimation of subgrade modulus and makes it suitable for the seismic analysis of tunnels with complex cross sections, like the horse shape, arch shape, or composite shape. However, these methods have complicated analysis procedures. Some essential issues have not been clarified, like the reasonable dimension of the substructure, which makes these methods hard to be adopted in engineering practice.
This study presents a general solution of the subgrade modulus for the response displacement method adopted in the seismic analysis and design of shield tunnels with a circular cross section. First, the analytical solutions of subgrade modulus for different ground deformation modes are obtained based on the theory of elasticity. Then, based on the solutions of subgrade modulus for ground deformation modes described by the Fourier series, the displacement fitting method is introduced to derive the subgrade modulus according to the actual deformation of the ground. Finally, the displacement fitting method is applied to a case study to demonstrate the validity.
2. Analytical Solutions of Subgrade Modulus for Different Ground Deformation Modes
Figure 1 presents the response displacement method. The seismic response of the ground, including the displacement, seismic shear stress, and acceleration at various depths, is first obtained through free-field analysis. Then, the relative displacement ∆uh at the tunnel perimeter with respect to the bottom of the tunnel is applied to the end of the ground springs; the seismic shear stress and inertia force of the structure are exerted on the tunnel structure. Ground springs are set up in both radial and tangential directions to represent the ground. The ground spring stiffness depends on the subgrade modulus (1). Kawashima [6] suggested the radial and tangential subgrade modulus solutions (2) for prescribed ground deformation modes which are described by the Fourier functions (3). The subgrade modulus depends on into which Fourier function the radial ground deformation would fall. However, the solutions from Kawashima [6] have the following drawbacks: (1) the subgrade modulus for uniform tangential movement at the tunnel perimeter is neglect; (2) the tangential subgrade modulus should not be supposed to be zero when the deformation parameter m is equal to 1 (m = 1); (3) the interaction between the radial and tangential subgrade modulus is not clarified when m > 1.where kspring is the ground spring stiffness, k is the subgrade modulus, and A is the area of the region the ground spring represents.

The analysis model to derive the analytical solutions of subgrade modulus for different deformation modes is shown in Figure 2. The shield tunnel with outer radius R1 is assumed to be constructed in a homogeneous isotropic infinite elastic ground with an elastic modulus of E1 and Poison’s ratio of ν1. The prescribed radial and tangential displacement described by Fourier functions is applied at the tunnel perimeter. The Airy stress function method is adopted to derive the stress change corresponding to the displacement at the tunnel perimeter. According to the definition, the subgrade modulus is the ratio of stress change to the displacement (4). The ground deformation modes are divided into three scenarios according to the deformation parameter m as listed in Table 1. As the ground displacement at the tunnel perimeter is predefined, there is no need to consider the original stress of the ground. The following assumptions are made in the derivation: (1) the ground is in the plane stress state, and the plane stress condition is in a direction perpendicular to the cross section of the tunnel; (2) there is no slip at the interface of the ground and tunnel lining; (3) the stress change and the displacement of the ground at infinity from the tunnel axis must be zero.

2.1. Subgrade Modulus for Deformation Mode I
As for the first case in the deformation mode I, the ground displacement only exits in the radial direction and distributes uniformly at the tunnel perimeter. According to Timoshenko and Goodier [16], the Airy stress function is
Based on the theory of elasticity [16], the stress change and the displacement of the ground are given bywhere ,
The following boundary conditions are applied to solve the constants in (6).(1)The boundary condition of stress at the infinity area for the ground is(2)The boundary condition of stress at r = R1 for the ground layer is
The constants and can be solved and expressed by .
According to the definition of subgrade modulus (4), the radial subgrade modulus is
As for the second case in the deformation mode I, the ground displacement only exits in the tangential direction and distributes uniformly around the tunnel perimeter. According to Timoshenko and Goodier [16], the Airy stress function is
Based on the theory of elasticity [16], the stress change and the displacement of the ground are given bywhere ;
The following boundary condition is applied to solve the solution of the constant in (12).(1)The boundary condition of stress at r = R1 for the ground is
The constants can be solved and expressed by .
According to the definition of subgrade modulus (4), the tangential subgrade modulus is
The negative sign enters the expression since the direction of the resulting displacements is opposite to the sense of the applied stress[17].
2.2. Subgrade Modulus for Deformation Mode II
Similar to the solution of subgrade modulus for deformation mode I, the Airy stress function for the deformation mode II is given by (16). The stress change and the ground displacement are given by (17):where ;
The following boundary conditions are applied to solve the constants in (17).(1)The boundary condition of stress at infinity area for the ground is(2)The boundary condition of stress at r = R1 for the ground is(3)The ground displacement is symmetrical to the vertical axis, so, for an arbitrary angle θ1, the following must be satisfied:(4)For an arbitrary angle θ2, the ground displacement meets the single-valued condition presented by the following: So, the coefficients of and in (17) are required to be zero; that is,(5)For this deformation model, the radial displacement of the ground ur is proportional to the logarithm of the radial coordinate, as shown by (17). Mathematically, this is a consequence of the infinitely long tunnel in a semi-infinite space, but it has no physical meaning [18, 19]. A characteristic dimension R3 should be incorporated into the solution to scale down the displacements with radial distance [18, 19]. In this solution, the radial displacement at r = R3 is arbitrarily set to zero; that is, The nonzero coefficients, , , and , are solved by the combination of (18)∼(20) and (22)∼(23) and are expressed by and . So, the boundary condition of ground displacement at r = R1 is where(6)The ratio of to is introduced to create a link between the radial and tangential subgrade modulus; that is,
With (25) and (27), the radial and tangential subgrade modulus can be expressed as follows:
2.3. Subgrade Modulus for Deformation Mode III
Similar to the previously mentioned derivation, the Airy stress function for the deformation mode III is given by (29), and the stress change and the displacement of the ground are given by (30):where ; .
The following boundary conditions are applied to solve the constants in (30).(1)The boundary condition of stress at infinity area for the ground is(2)The boundary condition of stress for the ground at r = R1 is The nonzero coefficients, and , are solved by the combination of (31) and (32) and are expressed by and . So, the boundary condition of ground displacement at r = R1 is where(3)The ratio of to is introduced to create a link between the radial and tangential subgrade modulus; that is,
With (34) and (36), the radial and tangential subgrade modulus can be expressed as (28):
The analytical solutions of radial and tangential subgrade modulus in the previously mentioned derivation are summarized in (38). It can be seen from this equation that the radial and tangential subgrade modulus are independent when m = 0, and the link between the radial and tangential subgrade modulus is built by introducing parameters and when m > 1. This analytical solution overcomes the defects in the solutions offered by Kawashima [6].
3. Solution of Subgrade Modulus Based on Displacement Fitting Method
The analytical solution of subgrade modulus deduced in the previous section is based on the fixed ground deformation modes. When this solution is adopted in the response displacement method, an error is likely to occur in the results as the ground deformation from the analysis result may not fall into the deformation shape upon which the subgrade modulus is derived. This conflict can be eliminated by employing the displacement fitting method. According to this method, the calculated ground displacement in the response displacement method is fitted by the Fourier functions (39). Hence, the ground displacement is decomposed to a series of simple deformation modes. The ground stress change introduced by the ground displacement is the sum of the products of the displacement and the corresponding subgrade modulus for each deformation mode (40). So, the subgrade modulus can be updated based on the definition (41), and the ground deformation can be recalculated with the newly estimated subgrade modulus. With an iterative process, the value of subgrade modulus will become more accurate, corresponding to actual ground deformation. The calculation flowchart of the displacement fitting method is shown in Figure 3.

The and in (40) are the and in (36), respectively. The and can also be derived based on (38). The radial and tangential displacement corresponding to the and are given by (42). By rotating this deformation shape degree, it turns out to be (43), which is consistent with the deformation mode listed in Table 1. At the same time, the parameter and are replaced by and in (44), respectively. The initial subgrade modulus is suggested to take subgrade modulus when m = 0 in (36) to avoid determining the value of and .
4. Application to a Design Case
4.1. Project Overview
A design case has been chosen to illustrate the validity of the displacement fitting method for the response displacement method. This case study is based on the Longquan tunnel belonging to the water diversion project in central Yunnan. The detailed information of the ground adopted in this study is listed in Table 2. The tunnel is buried in the gravelly soil layer and has a cover depth of 35.6 m. The estimated peak ground acceleration for the region of Longquan tunnel site is 0.2 g for a 10% probability of exceedance in 50 years as stated in the report of seismic ground motion parameter zonation [20]. According to the seismic ground motion parameters zonation map of China, the tunnel site has a seismic intensity of VIII degree [21].
This tunnel is built with the shield tunnel method and adopts the segmental lining as the tunnel structure. The segmental ring has an outer radius R1 of 3.1 m, a thickness of 0.3 m, and a width B of 1.2 m. Six segments, including three standard segments (block B1, B2, and B3), two adjacent segments (block L1 and L2), and one key segment (block F), compose the segmental ring. Two aligned steel bolts connect every two adjacent segments in one segmental ring at the longitudinal joint. The segmental rings are attached with sixteen bolts at the circumferential joint. The detailed layout of the segmental ring is shown in Figure 4. All the segments were precast with the reinforced high-performance concrete with a grade of C50.

4.2. Seismic Response of the Site
The seismic response of the site is obtained based on an equivalent linear analysis performed using DEEPSOIL software [22]. Variation of normalized shear modulus reduction (G/Gmax) and material damping (%) with cyclic shear strain for each ground layer are given by the field survey and are shown in Figure 5. The seismic bedrock interface is set to 70.0 m below the ground surface and loaded with an artificial seismic wave with the peak acceleration 0.25 g for a 5% probability of 50-year exceedance (Figure 6). The site response along the height of the tunnel is shown in Figure 7 when the maximum relative displacement between the tunnel vault and bottom is achieved.



4.3. Beam-Spring Model
The beam-spring model adopted to simulate the segmental lining in this study is shown in Figure 8. To consider the coupling effect between segmental rings when the staggered assemble strategy is adopted, this model contains three segmental rings, where the middle one with the entire width of segment B is the target ring. In comparison, the other two rings have a half width of the segment. Beam elements with stiffness equal to the segments are adopted to simulate the segments, and the bending stiffness is EcI for the middle ring and EcI/2 for the other two. The mechanical behaviors of the segment joint include compression along the tangential direction, shear along the radial direction, and rotation around the joint, which are simulated by axial, radial, and rotational springs with the stiffness of kA, kR, and kθ, respectively. The stiffness of the axial and radial springs has a negligible effect on segmental lining behaviors [23]. In this study, the axial and radial spring stiffness kA and kR are empirically set to ten times the compression stiffness and shear stiffness of segments, respectively. The rotation stiffness of the segment joint kθ is nonlinear and depends on the axial force and rotation angle. A simplified bilinear relationship proposed by Janssen [24] is adopted (45). Shear springs are set up along radial and tangential directions to transmit forces between adjacent rings under the coupling effect. All the joint parameters are summarized in Table 3. Both radial and tangential ground springs are set up. The ground spring stiffness is estimated using (1) and the displacement fitting method.where Ec is the elastic modulus of concrete and is 34.5 × 103 MPa, N is the axial force under static load from ground and is set to 2150 kN/m, M is bending moment at the joint, and h′ is the effective joint height and is 0.3 m in this study.

4.4. Results and Discussion
The initial radial and tangential subgrade modulus is set according to (38) when m = 0. The characteristic dimension R3 is three times R1. This case study also conducted the analysis in which the subgrade modulus is set according to Kawashima [6] to investigate the error induced when determining subgrade modulus based on fixed ground deformation modes.
4.4.1. Iteration Process
In total, ten iterative steps are conducted for this case, in which the parameter n in (39) is set to 3. Figure 9 presents the fitting results in the first and tenth steps. It can be found that the Fourier functions with n = 3 are enough to describe the deformation shape of the ground. Besides, the coefficient of determination (R2) for the fitting functions is 1.00 in each step.

(a)

(b)
The convergence of the error between two successive steps is shown in Figure 10. The relative error for bending moment M, radial displacement ur, and radial subgrade modulus kr at 30° of the target segmental ring are plotted against the step number. The relative error decreases fast during the iteration procedure. In the third step, the relative error has fallen below 1.05 × 10−2, acceptable in engineering practice. If the iteration termination condition is that the relative error of M is less than 0.1%, only five iterative steps are needed. The relative error indicates that the displacement fitting method efficiently determines the subgrade modulus according to the ground displacement for the response displacement method.

4.4.2. Subgrade Modulus Corresponding to the Ground Deformation
Fitting results for the displacement parameters in (37) (Table 4) reveal that the deformation mode with m = 2 is the dominant mode. The deformation mode with m = 1 or m = 3 also influences the ground deformation shape. So, a single deformation mode does not describe the actual ground deformation shape accurately.
The radial and tangential subgrade modulus corresponding to actual ground deformation is shown in Figure 9. Four singular points occur for the radial and tangential subgrade modulus, where the ground displacement is close to zero. The subgrade modulus changes violently near the singular points, which has a negligible effect on tunnel lining due to the small magnitude of the displacement. The radial and tangential subgrade modulus no longer distribute uniformly around the tunnel despite the variation near the singular points. Both the radial and tangential subgrade modulus gradually increase from tunnel vault to bottom. Comparing the subgrade modulus estimated by the displacement fitting method with the solutions from Kawashima [6] (Figure 11), there lies an apparent difference between them. The relationship for the subgrade modulus is > > > and > > = .

(a)

(b)
4.4.3. Internal Forces and Displacement of Tunnel Lining
The internal forces and displacement of the tunnel lining based on subgrade modulus for fixed deformation modes and actual deformation shape are presented in Figures 12 and 13. Table 5 lists the maximum and minimum values of internal forces and displacement. The percentage in Table 5 is the increase by comparing the results based on subgrade modulus for fixed deformation modes with those based on the displacement fitting method. It can be figured out that the seismic response of the segmental lining diverges when determining the subgrade modulus according to different approaches. The internal forces and displacement resulting from subgrade modulus for m = 2 are closest to those adopting the displacement fitting method and share the same distribution characteristic. However, the discrepancy in the results is notable. For example, the minimum axial force Nmin obtained based on the displacement fitting method is 11.4% smaller than that adopting subgrade modulus for m = 2. What is more, when the subgrade modulus is determined according to the fixed deformation mode with m = 0 or 1, the tangential subgrade modulus is ignored, leading to a different axial force distribution (Figure 12(b)) and integral structure rotation (Figure 13(b)). It can be concluded that when the subgrade modulus is determined based on a fixed deformation mode, the error will inevitably arise. The displacement fitting method offers the subgrade modulus corresponding to the actual ground deformation, which helps avoid evaluating the rationality of the results with subgrade modulus derived upon fixed deformation modes.

(a)

(b)

(c)

(a)

(b)
5. Conclusions
This study presents a solution to determine the subgrade modulus in the response displacement method adopted in the transverse seismic analysis of the shield tunnel with a circular cross section. Analytical solutions of subgrade modulus under different deformation modes are firstly deduced based on the theory of elasticity. The displacement fitting method is introduced to derive the subgrade modulus corresponding to ground displacement, and this method is applied to a design case. Some conclusions can be drawn as follows.(1)The subgrade modulus for the fixed deformation modes described by the Fourier series is independent of θ. The radial and tangential subgrade modulus are independent when the deformation parameter m is equal to zero (m = 0). For these deformation modes with m > 1, the interaction between the radial and tangential subgrade modulus should be considered. By introducing parameter Ω1 or Ωm, the link is built between the radial and tangential subgrade modulus.(2)The subgrade modulus is no longer constant around the tunnel when the ground deformation contains a series of deformation modes described by the Fourier series.(3)The displacement fitting method is proved to be efficient in deriving the subgrade modulus corresponding to the ground deformation and helps avoid the error introduced by determining subgrade modulus based on a fixed deformation mode.
Data Availability
The datasets generated and analyzed in the current study are available from the corresponding author upon reasonable request.
Conflicts of Interest
The authors declare no conflicts of interest.
Acknowledgments
The authors appreciate the support from the National Natural Science Foundation of China (nos. 51991394, 51878569, and 52078430) and Yunnan Provincial Science and Technology Department of China (no. 202002AF080003).