Numerical Study of Influencing Factors of Safety and Stability of Tunnel Structure under Airport Runway

Abstract

A six-degree-of-freedom mathematical model and mechanical balance equation of a “five-point-contact” aircraft are established in this study. The model and equation are used to investigate the safety and stability of a tunnel structure under the runway of an airport, particularly when aircraft taxi or move on the runway. ABAQUS is used to construct a three-dimensional finite element model of the cooperative deformation of the airport runway–soil–tunnel structure. The analysis focuses on the response and evolution of structural safety mechanical indices from the perspective of three influencing factors: type of aircraft, road surface, and burial depth. The results show that the distribution position of the main landing gear wheel is more concentrated using the dynamic load equation of different aircraft. A rigid pavement is not easily deformed when subjected to aircraft loads, whereas a flexible pavement has an excellent attenuation effect on diffusing forces. The shear stresses on the upper and lower arches of the tunnel structure differ depending on the pavement material. The deformation of the arches under shear stress is more intense than that of other parts. With an increase in burial depth, the tunnel structure withstanding the aircraft load disturbance exhibits an attenuation trend. The disturbance caused by soil stress to the tunnel structure must not be ignored. When the burial depth of the tunnel exceeds 64 m, the tunnel structure ceases to be disturbed by aircraft loads. The research results can significantly guide airport construction and be used as a reference for investigating the safety and stability of substructures under airport runways.

1. Introduction

In recent years, many countries have not been sufficiently attentive to environmental protection due to their pursuit of economic development as it rapidly advances. Consequently, the emission of “greenhouse gases” has become an environmental problem that must be urgently resolved. Moreover, “carbon neutrality” has become a popular topic of discussion worldwide. Because of the severe impact of automobile exhaust on the environment, governments have urged people to travel using low-carbon transportation, thus promoting the rapid development of subways. Air travel has the advantages of speed and efficiency, and the number of travelers preferring this mode of transportation to save time and cost has continued to increase. Presently, most cities have initiated projects linking subways and airports to optimize urban transportation planning further. Accordingly, a form of travel that links transportation modes has been created, directly enabling the transfer of passengers to a corresponding flight without leaving the station. For example, the internal cargo lane of Beijing Capital International Airport crosses airport taxiway L, a shield tunnel under the runway of Taipei Songshan Airport. The underground service lane of Guangdong Zhuhai Airport crosses the apron; two railroad tunnels of Chengdu Metro Line 10 are under the apron and taxiway at Shuangliu Airport. In the UK, two railroad tunnels are under the taxiway and station apron of Heathrow Airport [1,2,3,4,5].

A subway tunnel crossing an airport runway is a large-scale project. To resolve the problems of high cost, time consumption, and resource consumption associated with this type of project, scholars generally employ numerical simulation in research. Aircraft load has a strong impact and can cause considerable damage. In view of the foregoing, the means for applying aircraft load to finite element numerical models has become a significant problem that requires solving. Research results are synthesized and summarized in the following three ways according to aircraft load application: (1) When the carload is simplified, pavement design is approximated as an elliptical ground area. In China, however, the current specification and engineering design mainly uses a circle for simplicity, and the aircraft load is typically simulated using this method in the preliminary study [6,7,8,9,10]. (2) The domestic “Specification for Airport Cement Concrete Pavement Design” (MHT5004-2010) prescribes the shape of a single wheel print of the main landing gear of the aircraft. According to this specification, the wheel print area is equated to a rectangle based on the principle of area equivalence and then converted to a nodal force or uniform load. This approach is employed by most scholars because it has been developed based on specifications [11,12,13]. (3) The aircraft load is simplified to a concentrated force that propagates to the bottom plate of the pavement at a 45° concrete punching tangent angle; thus, a load-spreading area is formed [14]. After solving the problem of aircraft load application through numerical simulations, the research of most scholars focuses on several areas, such as parking lots, airport runways, aircraft carriers, dynamization, and mechanization of system equipment subjected to pressure waves [15,16,17]. Such studies are implemented to analyze the airport runway material, dynamic response, and stability under aircraft load [18,19,20,21]. In addition, because airport construction starts late, although investigations on the dynamic response of the space under airport runways subjected to aircraft loads have been conducted, the number of these studies are limited compared with those performed on airport pavement structures. For example, Sun et al., (2018) [14] attributed aircraft loads to nodal and surface dynamic loads according to the finite element method. They also analyzed the dynamic response of the equivalent support section of a tunnel under four working conditions with the action of nodal dynamic loads, concluding that the structural dynamic effect of aircraft loads is evident. Chen et al. (2019) [7] used FLAC3D software to apply aircraft load as a concentrated force on the model. They compared the numerical simulation results with field monitoring data to analyze the influence law of aircraft load on the tunnel boring machine tube sheet at different locations and determine whether the tube sheet structure is safe. Wei et al. (2022) [22] established a six-degree-of-freedom mathematical model of the A380–800 aircraft and a calculation model of the airport runway–soil–tunnel structure collaborative deformation. They also analyzed the influencing factors of the mechanical response of tunnel lining considering different angles.

The foregoing studies show that aircraft load considerably affects the safety of tunnels under airport runways. However, no unified research system or reliable evaluation method for examining the safety and stability of an underground structure affected by aircraft load has ever been formulated locally or internationally. Moreover, the consideration of the aircraft load excitation equation is not comprehensive, and most of the studies ignore the influence of lift when establishing the aircraft dynamic equation. Most of the investigations only focus on one influencing factor, and other factors are insufficiently considered. Accordingly, a comprehensive numerical study is essential to analyze the factors affecting structural safety and stability of tunnels under airport runways. This study formulates a six-degree-of-freedom mathematical model of aircraft considering lift. Moreover, a dynamic equilibrium equation is derived, and a collaborative deformation finite element numerical model of an airport pavement–base layer–soil–tunnel lining is established. The model is employed to examine the factors influencing an existing shield tunnel. The foregoing can be used as a reference for the design of tunnel structures under pavements built during airport construction.

2. Aircraft Dynamic Model and Vibration Excitation Equation

An aircraft is in a complex state of force during takeoff, landing, and taxi; accordingly, its dynamics are modeled to derive vibration equilibrium equations to improve research. Professor Ling of Tongji University and Professor Cheng of the Civil Aviation University of China [23,24,25] developed the theoretical basis for a “three-point-contact” aircraft with six degrees of freedom by investigating the levelness of airport pavement. Based on the foregoing theory, Wei and Yang [22,26] derived a mathematical model appropriate for the A380–800 aircraft with six degrees of freedom. The foregoing studies only considered the force state of an aircraft at rest in the relevant dynamic model established, neglecting the force state of the aircraft in motion (takeoff, landing, and taxi). The dynamic model of the “three-point-contact” aircraft ignores this factor because passenger volume increases significantly. Moreover, the aircraft’s specifications have evolved, and the main landing gear configuration has become increasingly complex. Accordingly, based on previous research, the current study establishes a six-degree-of-freedom mathematical model and vibration equilibrium equations for a “five-point-contact” aircraft considering lift. The results of this study are anticipated to pave the way for subsequent numerical simulation works and be used as reference for airport construction projects.

2.1. Mathematical Model of “Five-Point-Contact” Aircraft with Six Degrees of Freedom

Most existing studies are based on mathematical models with one or two degrees of freedom that are relatively simple. However, owing to the complex motion of an aircraft, limiting the degree of freedom to the aforementioned extent as well as ignoring the aircraft’s vertical, pitch, and lateral rotations is not a comprehensive approach. Based on the theoretical basis of a classic 1/4 vehicle model, a dynamic model of a six-degree-of-freedom “five-point-contact” large aircraft considering the aircraft’s takeoff, landing, and taxiing motion states is proposed. The established model is shown in Figure 1.

With respect to vertical motion, the reed mass of the large aircraft is M (kg); N is the lift force acting on the large aircraft (kN); and Z is the vertical displacement of the entire aircraft (m). For the front wheels, left rear landing gear, right rear landing gear, left rear wing wheel landing gear, and right rear wing wheel landing gear, the unsprung masses are m1, m2, m3, m4, and m5 (kg), respectively; the stiffness coefficients are K1, K2, K3, K4, and K5, respectively; and the damping coefficients are C₁, C₂, C₃, C₄, and C₅, respectively, and incentives for the existence of uneven roads are q₁, q₂, q₃, q₄, and q₅, respectively. The vertical displacements (m) of the sprung masses and the unsprung masses of the front wheels, left rear landing gear, right rear landing gear, left rear wing wheel landing gear, and right rear wing wheel landing gear are Z_a, Z_b, Z_c, Z_d, and Z_e, respectively, and z₁, z₂, z₃, z₄, and z₅, respectively.

I_x is the rotational inertia of the aircraft model around the x axis (kg·m²); and a, b, and c are the vertical distances from the front landing gear, rear landing gear, and wing wheel landing gear, respectively, to the x axis (m).

I_y is the rotational inertia of the aircraft model around the y axis (kg·m²); and L₂, L₃, L₄, and L₅ are the vertical distances from the left rear landing gear, right rear landing gear, left wing wheel landing gear, and right wing wheel landing gear to the y-axis (m), respectively.

2.2. Balance Equation of “Five-Point-Contact” Aircraft Vibration

During takeoff and taxi, an aircraft is subjected to lift. Assume that the aircraft has uniform acceleration and a constant lift coefficient during takeoff. According to “Aerodynamics—Concepts, quantities and symbols—Part 4: Aerodynamic forces, moments, their coefficients and derivatives of aircraft (GB/T 16638.4-2008)” [27], the following can be used to analyze the force state of the aircraft:

$$\begin{array}{l}F=G-N\\ N=C_y\rho S{V}_0{}^2/2\end{array}$$

(1)

where F is the combined force on the aircraft (kN); G is the glide weight of the aircraft (kN); N is the lift force on the aircraft (kN); C_y is the lift coefficient; ρ is the air density (kg/m ³); S is the wing area (m²); and V₀ is the relative speed (m/s).

The vertical vibration balance equation for the spring-loaded mass, M, under vibration excitation owing to the unevenness of the road surface during aircraft taxi is established based on D’Alembert’s principle (the center of mass of the aircraft model is considered as the origin) [22,23,24,25,26,28]:

$$\begin{array}{ll}M\ddot{Z}+& C_1\left(\dot{Z_a}-\dot{z_1}\right)+K_1\left(Z_a-z_1\right)+C_2\left(\dot{Z_b}-\dot{z_2}\right)+K_2\left(Z_b-z_2\right)+\\ & C_3\left(\dot{Z_c}-\dot{z_3}\right)+K_3\left(Z_c-z_3\right)+C_4\left(\dot{Z_d}-\dot{z_4}\right)+K_4\left(Z_d-z_4\right)+\\ & C_5\left(\dot{Z_e}-\dot{z_5}\right)+K_5\left(Z_e-z_5\right)+Mg-N=0\end{array}$$

(2)

The pitch rotation equation of the “five-point-contact” model of an aircraft with the spring-loaded mass, M, is

$$\begin{array}{ll}{I}_x\ddot{\theta }+& C_1\left(\dot{Z_a}-\dot{z_1}\right)·a+K_1\left(Z_a-z_1\right)·a-C_2\left(\dot{Z_b}-\dot{z_2}\right)·b-K_2\left(Z_b-z_2\right)·b-\\ & C_3\left(\dot{Z_c}-\dot{z_3}\right)·b-K_3\left(Z_c-z_3\right)·b-C_4\left(\dot{Z_d}-\dot{z_4}\right)·c-K_4\left(Z_d-z_4\right)·c-\\ & C_5\left(\dot{Z_e}-\dot{z_5}\right)·c-K_5\left(Z_e-z_5\right)·c=0\end{array}$$

(3)

The equation for the lateral rotation equilibrium of the “five-point-contact” aircraft model with a spring-loaded mass, M, is

$$\begin{array}{ll}{I}_y\ddot{\psi }+& C_3\left(\dot{Z_c}-\dot{z_3}\right)\cdot L_3+K_3\left(Z_c-z_3\right)\cdot L_3+C_5\left(\dot{Z_e}-\dot{z_5}\right)\cdot L_5+K_5\left(Z_e-z_5\right)\cdot L_5-\\ & C_2\left(\dot{Z_b}-\dot{z_2}\right)\cdot L_2-K_2\left(Z_b-z_2\right)\cdot L_2-C_4\left(\dot{Z_d}-\dot{z_4}\right)\cdot L_4-K_4\left(Z_d-z_4\right)\cdot L_4=0\end{array}$$

(4)

The equilibrium equations of the unsprung masses of the front and rear wing wheel landing gears in the vertical direction for the “five-point-contact” aircraft model are as follows:

$$\begin{array}{l}{m}_1\ddot{z_1}-C_1\left(\dot{Z_a}-\dot{z_1}\right)-K_1\left(Z_a-z_1\right)+c_1\dot{z_1}+k_1{z}_1-c_1\dot{q_1}-k_1{q}_1=0\\ m_2\ddot{z_2}-C_2\left(\dot{Z_b}-\dot{z_2}\right)-K_2\left(Z_b-z_2\right)+c_2\dot{z_2}+k_2{z}_2-c_2\dot{q_2}-k_2{q}_2=0\\ m_3\ddot{z_3}-C_3\left(\dot{Z_c}-\dot{z_1}\right)-K_3\left(Z_c-z_3\right)+c_3\dot{z_3}+k_3{z}_3-c_3\dot{q_3}-k_3{q}_3=0\\ m_4\ddot{z_4}-C_4\left(\dot{Z_d}-\dot{z_4}\right)-K_4\left(Z_d-z_4\right)+c_4\dot{z_4}+k_4{z}_4-c_4\dot{q_4}-k_4{q}_4=0\\ m_5\ddot{z_5}-C_5\left(\dot{Z_e}-z_5\right)-K_5\left(Z_e-z_5\right)+c_5\dot{z_5}+k_5{z}_5-c_5\dot{q_5}-k_5{q}_5=0\end{array}$$

(5)

A corresponding equivalence exists between the vertical displacement of the spring-loaded mass and vertical displacement of the unsprung mass of the “five-point-contact” aircraft model. The distribution of the vertical displacement at each point of the spring-loaded mass is as follows:

$$\begin{array}{c}{Z}_a=Z+a\theta \\ Z_b=Z-b\theta -L_2\psi \\ Z_c=Z-b\theta +L_3\psi \\ Z_d=Z-c\theta -L_4\psi \\ Z_e=Z-c\theta +L_5\psi \end{array}$$

(6)

2.3. Excitation Equation of Aircraft Load Vibration

Because the tunnel is under an airport runway, the impact of the aircraft load is considerable. Hence, the means for establishing the aircraft dynamic load equation becomes the primary consideration. Based on the provisions of Article 5.0.2.1, “Specification for Airport Cement Concrete Pavement Design (MHT 5004-2010)” [29], an equivalent area method is used to simplify the shape of the grounded wheel print of the aircraft wheel. The equivalent area method of this simplification technique is shown in Figure 2.

In this study, the safety and stability of the underground structure of airport pavement subjected to the loads of different aircraft types are compared, and the weak and damage-prone areas of the tunnel structure under the most unfavorable conditions are observed. For this purpose, A380–800, B-747–400, and B-737–800 aircraft with different landing gear configurations based on the data of MHT 5004–2010 are selected and used in the study. These aircraft are the most massive in terms of load and specifications. The main landing gear configurations are shown in Figure 3; the load application points (center-of-mass points) of the main landing gear of each type are marked.

The “Dload” subroutine implements aircraft loading in the numerical simulation. Hence, this section is based on Equation (3.0.2) of MHT 5004–2010. The aircraft load excitation equations for the following study are

$$\mathrm{A}\mathrm{380-800} \mathrm{aircraft} \mathrm{single} \mathrm{wheel} \mathrm{load}: F_{s_1}=\frac{\rho G}{n_{c_2}·n_{w_2}}=\frac{0.97·5620}{2·10}=272.6 \mathrm{kN},$$

$$\mathrm{B-747-400} \mathrm{aircraft} \mathrm{single} \mathrm{wheel} \mathrm{load}: F_{s_2}=\frac{\rho G}{n_{c_2}·n_{w_2}}=\frac{0.952·3978}{4·4}=236.7 \mathrm{kN},$$

$$\mathrm{B-737-800} \mathrm{aircraft} \mathrm{single} \mathrm{wheel} \mathrm{load}: F_{s_3}=\frac{\rho G}{n_{c_2}·n_{w_2}}=\frac{0.95·792.6}{2·2}=188.2 \mathrm{kN},$$

where is the aircraft’s main landing gear load distribution factor; G is the maximum weight of the aircraft (kN); $N_{c_2}$ is the number of main landing gears; and $N_{w_2}$ is the number of wheels in the main landing gear.

The aircraft dynamic load equation can be transformed into a cycle-varying simple harmonic vibration load. The aircraft load equation can be scaled up by considering the effect of aircraft vibration [9,11]. Accordingly, the aircraft’s main landing gear dynamic load equation is

$$\begin{array}{c}{F}_M=F_s+0.1F_{s}sin\left(\omega t\right)\\ f=\frac{v}{\pi D}\\ \omega =2\pi f\end{array}$$

(7)

where $\overline{\omega }$ is the rotational circular frequency of the wheel (known as angular frequency); t is the time variable (s); and $f$ is the frequency (Hz). The wheel diameter is 1.5 m, and the aircraft’s takeoff speed is considered as 80 m/s. According to the foregoing formula, the circular frequency of wheel rotation is $f=\frac{80}{1.5\pi }=16.97 \mathrm{Hz}, \omega =2\pi f=106.66 \mathrm{rad}/\mathrm{s}$. Then, the main landing gear dynamic load equations for the three types of aircraft are as follows:

$$\mathrm{A}\mathrm{380-800}, F_{M_1}=272.6+27.26sin\left(106.66t\right) \left(\mathrm{kN}\right);$$

$$\mathrm{B-747-400}, F_{M_2}=236.7+23.67sin\left(106.66t\right) \left(\mathrm{kN}\right);$$

$$\mathrm{and} \mathrm{B-737-800}, F_{M_3}=188.2+18.82sin\left(106.66t\right) \left(\mathrm{kN}\right).$$

Because the Dload subroutine is loaded in the form of model-applied pressure, the main landing gear dynamic loading equation is transformed into a pressure–time equation using $P_M=\frac{F_M}{S}$, where S is the simplified wheel impression area of the machine wheel; the specific dimensions of the wheel are shown in Figure 3.

The main landing gear pressure–time equations for A380–800, B-747–400, and B-737–800 (unit: kPa) are

$$P_{M_1}=1456.19+145.62sin\left(106.66t\right) \left(\mathrm{kPa}\right);$$

$$P_{M_2}=1420.77+142.08sin\left(106.66t\right) \left(\mathrm{kPa}\right);$$

$$P_{M_3}=1458.91+145.89sin\left(106.66t\right) \left(\mathrm{kPa}\right).$$

The main landing gear pressure–timescale curves for the three aircraft types using MATLAB are shown in Figure 4.

3. Numerical Simulation

ABAQUS is widely regarded as the most powerful finite element software, which can analyze complex structural mechanical systems of solid mechanics, especially manage very large and complex problems and simulate highly nonlinear problems. In this study, ABAQUS software was employed to establish a numerical calculation model of the airport runway–base layer–soil–tunnel structure. The basic model size was defined as 80 m × 80 m × 80.8 m. The airport road surface was designed according to the code, and the mechanical response of the tunnel concrete lining and structurally weak parts under aircraft load was observed. A schematic of the numerical model is shown in Figure 5.

3.1. Model Meshing

In this study, Visual Studio programming software was chosen to write the Dload subroutine, which was combined with the subroutine interface packaged with ABAQUS; the subroutine enables the imposition of aircraft moving load. This approach has high calculation accuracy and can effectively save computational resources. Because the overall model of the structure is larger than the aircraft wheel print area, the accuracy of the software post-processing module results is affected if the mesh division is large. Therefore, the key study area is divided into meshes by local encryption. The considered global grid size of the structure was 4 m × 4 m × 4 m, the grid size of the runway was 1 m × 1 m × 1 m, and the tunnel lining structure was cut into 8 parts circumferentially. Each part was divided into 12 grids by limiting the number of grids; that is, the number of circumferential grids is 96. The specific size of the grid division is shown in Figure 6. The number of pavement–soil–tunnel structure grids was 158,928, and that of the tunnel lining structure grids was 38,016. The mesh check option in ABAQUS software checked the quality of the meshing and made the simulation results more precise.

3.2. Boundary Conditions and Element Types

This study focuses on the numerical modeling of tunnels under airport pavements. Because of the high excitation impact effect of aircraft loads, the problem of stress wave reflection in the numerical simulation is encountered; therefore, the boundary conditions applied to the soil must be considered. Many types of soil boundaries, such as viscous, transmissive, and viscoelastic boundaries (each with its own characteristics and application ranges), are used in finite element calculations [30]. Based on previous studies, the problem of stress wave reflection inside the soil can be solved by extending the model size. In addition, because of the large scale of the numerical model built in this study, the stress wave reflection of the soil can be handled well; hence, an artificial static boundary can be set in the soil model. Hence, in the numerical simulation, the normal displacements, U, in the X and Y directions are all set to 0 for the front, rear, left, and right of the model; the bottom is completely fixed (U1 = U2 = U3 = UR1 = UR2 = UR3 = 0).

You et al. [31] conducted an in-depth study on the establishment of a finite element model of an asphalt pavement structure. They compared and analyzed the effects of various contact conditions and mesh cell types on the calculation and accuracy of the model. Based on their conclusions, the C3D8R (8-node linear unit) cell was chosen to simulate each component of the numerical model to represent the system characteristics as accurately as possible and improve computational accuracy and efficiency.

3.3. Material Properties and Contact Settings

The three-dimensional (3D) finite element model simulates the theory of elastic laminar system and satisfies the following basic assumptions [31]: (1) All layers, except for the soil, are composed of isotropic homogeneous linear elastic materials. (2) The soil layer is infinite in both horizontal and depth directions, whereas the structural layer of the pavement structure is limited in depth and remains infinite in the horizontal direction. (3) Each layer connection is continuous. Because this study is mainly concerned with the stability of the tunnel structure, the airport pavement structure is modeled using the linear elastic material properties considering the foregoing assumptions. The elastic–plastic Mohr–Coulomb failure criterion is used as the soil constitutive model to simulate the nonlinear behavior of soil under aircraft loads with Rayleigh damping and small strain damping based on the Rayleigh damping equation:

$$\left[C\right]=\alpha \left[M\right]+\beta \left[K\right]$$

(8)

where $\left[M\right]$ is the mass matrix, and $\left[K\right]$ is the stiffness matrix. The Rayleigh damping factors, a and β, can be calculated using the following equations:

$$\alpha =\zeta \frac{2{\overline{\omega }}_1{\overline{\omega }}_2}{{\overline{\omega }}_1+{\overline{\omega }}_2}$$

(9)

$$\beta =\zeta \frac{2}{{\overline{\omega }}_1+{\overline{\omega }}_2}$$

(10)

where $\zeta$ is the damping ratio, and ${\overline{\omega }}_1$ and ${\overline{\omega }}_2$ are the first and second intrinsic frequencies of the system, respectively. The main material properties of the numerical model selected according to the MH/T 5004–2010 specification for airport cement concrete pavement design and MH/T 5010–2017 specification for asphalt pavement design of civil airports [29,32] are listed in

Table 1. The main material parameters.

Item	Material	Thickness (m)	Unit Weight (kg·m−3)	Elastic Modulus (MPa)	Poisson’s Ratio	Damping Constant $\mathit{\alpha }\left({\mathit{s}}^{-1}\right)$	Damping Constant $\mathit{\beta }\left(\mathit{s}\right)$	Mohr–Coulomb Model
								C (MPa)	$\mathit{\phi }(°)$
Rigid pavement	Cement concrete pavement	0.4	2440	36,000	0.2	0.1	0.001	/	/
	Cement-stabilized gravel base	0.4	2100	1500	0.25	0.1	0.002	/	/
Flexible pavement	Asphalt concrete pavement	0.2	2380	1800	0.25	0.1	0.001	/	/
	Cement-stabilized gravel base	0.3	2100	1500	0.25	0.1	0.002	/	/
	Broken stone base course	0.3	2000	200	0.2	0.1	0.002	/	/
Foundation course	Soil	80	1800	30	0.4	1	0.01	0.001	25
Tunnel lining	Concrete	1.5	2440	36,000	0.2	0.1	0.001	/	/

.

Since this study focuses on the dynamic response of the tunnel structure, the contact between the soil layers is simplified to be fully contacted, and the excavation layer model is divided by direct cutting. The interaction between the soil and tunnel structure has been considered in modeling. Hard contact is set as the normal contact, and the “penalty function” is set as the tangential contact between the two; the friction coefficient used is 0.577 [33].

3.4. Setup Analysis Step

In the process of numerical simulation, two analysis steps are set to simulate the stress state of the tunnel structure in the soil more realistically. The static general analysis step and the implicit dynamic analysis step were carried out in which the gravity was applied in the first step and the aircraft load was applied in the second step, respectively. The numerical simulations in this study were pre-equilibrated for the entire model. Moreover, initial stresses were applied to eliminate large model deformation and displacement due to the self-weight of the structure and the soil stress.

After setting the initial conditions, the Dload subroutine was executed on the UMAT interface included in ABAQUS to load the aircraft and dynamically analyze this load. Pressure was applied to the runway area surface of the model in the Load module in ABAQUS. Selected User was defined to call the Umat interface at the input pressure and select the Dload subroutine file in the Job module to complete the application of aircraft load. Owing to the large size of the model, the analysis is performed at a 0.005-s time step to conserve resources; however, the calculation results are output to 100 frames.

3.5. Verification of Numerical Model Accuracy

The numerical simulation technique is popular among scholars because of its advantages, such as conserving resources and cost savings; the accuracy of its results is also a considerable benefit. To verify whether the Dload subroutine can accurately impose the aircraft moving load, a dynamic analysis is conducted on the model based on the calculation example reported in the literature [34].

A uniform infinitely long thin plate is 5 m wide and 0.1 m thick. It has a mass density of 2500 kg/m³, a modulus of elasticity of 10 MPa and a Poisson’s ratio of 0. The Rayleigh damping factors are $=0.1s^{-1},\beta =0.01s$, respectively. The bottom of the slab is vertically restrained, and a sinusoidal excitation load ($P=\sin \left(20\pi t\right)\frac{\mathrm{N}}{{\mathrm{m}}^2}, t\in \left[0, 0.5\right]$) is applied to within 1 m of the center of the slab. The C3D8R cell is used to divide the mesh of the model. The mesh size is 0.1 m × 0.1 m × 0.1 m, and the analysis step time is 0.5 s. The schematics of the specific calculation example and finite element model are shown in Figure 7.

The vertical displacement time curves for observation points A, B, and C are shown in Figure 8. The general trend of the graphs indicates that the application of the excitation load by the Dload subroutine developed using Visual Studio is generally in agreement with the results reported in the literature [34]. The aircraft excitation load acting on the model was effectively simulated by the Dload subroutine.

4. Analysis of Numerical Simulation Results

Developments in airport construction are rapidly increasing. This section considers the influence of different factors on the structural mechanical response of tunnel structures under airport pavements from different perspectives. Because the use of shield tunnels is increasing, the summarized laws provide a research basis for airport and shield tunnel construction.

4.1. Influence of Different Aircraft on Tunnel Structure

Aircraft types vary; hence, the main landing gear structures and load distribution coefficients differ. Consequently, the dynamic load excitation effect generated by aircraft taxi also varies. Because different types of aircraft takeoff, land, and taxi in airports, the runway is subjected to daily load excitations. Hence, this section discusses the analysis of the influence of different aircraft types on the dynamic response of the tunnel structure under the runway. It also presents the dynamic instability and weak positions of the structure.

The radial stress and peak displacement of the tunnel lining structure of the three models are shown in Figure 9; different aircraft dynamic load equations are applied to the models. Moreover, the mechanical response of the tunnel structure varies. As indicated by the pressure–time equation for each aircraft (Section 2.3), the B-737–800 model has its highest load excitation, followed by the A380–800 model and then the B-747–400 model. As shown in Figure 9, the B-747–400 and B-737–800 aircraft loads have the greatest and least effects on the tunnel structure, respectively; the effect of the A380–800 aircraft load is between those of the two. This phenomenon occurs due to the different configurations of the main landing gear. The main landing of the B-737–800 aircraft is in the form of “double wheels” having a small contact area with the pavement. Moreover, because its fuselage weight is considerably less than those of the other two types, its disturbance effect on the structure is not considerable. The A380–800 and B-747–400 main landing gears are “three-axle two-wheel” and “two-axle two-wheel” configurations, respectively. Although the taxiing weight of B-747–400 is considerably less than that of A380–800, in calculating the dynamic load for each model, the main landing gear distribution coefficients of the two aircraft are not the same. Moreover, the width and footprint of the main landing gear of B-747–400 are smaller than those of A380–800.

The numerical segment shown in Figure 9 indicates that the peak radial stress at the top part of the tunnel occurs in the positive direction, indicating that the top plate is subjected to tensile stress. In contrast, the stress on the rest of the tunnel is in the negative direction; hence, the stress is compressive. The radial displacement generated at the top plate position shown in Figure 9b is positive, whereas the peak displacement at the bottom plate position is negative, conforming to the trend of the tunnel lining structure stresses. Concrete has an outstanding compressive capacity and relatively weak tensile capacity; hence, the tunnel structure can sustain aircraft loads. However, the top slab area tends to form weak areas; hence, it must be the focus of structural design and reinforcement.

The vertical acceleration—time curve corresponding to each monitoring point—is shown in Figure 10. The figure indicates a relatively gentle range of variation in vertical acceleration when the structure is subjected to the B-737–800 aircraft load. When the structure is subjected to the A380–800 and B-747–400 aircraft loads, the corresponding mechanical responses were comparable. This is because the aircraft taxiing speed, number and configuration of main landing gears, and other factors are the same. However, the single wheel load and equivalent simplified area of the wheels depend on the model; consequently, the aircraft dynamic load equation also varies. The pressure–time equations of the three aircraft have the highest and least effects on A380–800 and B-737–800 aircraft, respectively; this observation corresponds to structural acceleration.

The vibration frequency of the tunnel structure varies with the taxiing position of the aircraft, as clearly shown in Figure 10. After the aircraft load is directly applied to the tunnel, the vertical acceleration abruptly changes, causing a greater impact on the structure. This indicates that the effect of aircraft taxiing on the runway substructure cannot be neglected. The top plate of the tunnel structure has the highest vertical acceleration variation of up to 259%, indicating that this section is more susceptible to impact and damage than the other parts. Because the taxiing weight of the A380–800 aircraft is maximum, and considering the most unfavorable conditions, the load of this aircraft is selected for succeeding studies.

4.2. Influence of Different Pavement Types on Tunnel Structure

Airport pavements are currently classified into rigid and flexible pavements. Before the technology of modified asphalt was perfected, asphalt materials with high-temperature stability are strictly required for flexible pavements. Consequently, rigid cement concrete pavements are mainly used in airports. This section presents two pavement material models that have been established in this study to analyze the mechanical response of the lower tunnel structure under the two pavements subjected to the A380–800 aircraft load. The material parameters and properties are listed in

Table 2. Aircraft parameter list.

Aircraft	Maximum Sliding Weight (kN)	Maximum Takeoff Weight (kN)	Maximum Landing Weight (kN)	Main Landing Gear Distribution Factor	Main Landing Gear Tire Pressure (MPa)
A380–800	5620.00	5600.00	3860.00	0.97	1.47
B-747–400	3978.00	3968.93	2857.63	0.952	1.38
B-737–800	792.60	790.04	663.80	0.95	1.47

.

To study the influence of different pavement materials on the structural force, the peak vertical displacement deformation of the pavement structure under the action of aircraft load is shown in Figure 11. Because the downward Z direction is specified as the positive direction in the modeling process, the deformation in this direction is positive. As shown in Figure 11, the degree of deformation of the flexible pavement significantly exceeds that of the rigid pavement. The peak vertical displacements of the rigid and flexible pavements are 14.3 and 19.6 mm, respectively; the difference between the two can reach 37.1%. The reason for this phenomenon is that the elastic modulus of the cement concrete rigid pavement is higher than that of the asphalt concrete flexible pavement. Hence, the rigid pavement does not easily deform under load and has a more prominent compressive capacity than the flexible pavement. However, because concrete is a brittle material, its tensile capacity is considerably less than its own compressive capacity. If the bending tensile stress due to the aircraft load action exceeds the ultimate tensile strength, then the concrete is prone to fracture. Therefore, both the compressive and tensile strengths of the runway must be considered in the design of airport rigid pavements. Because this study focuses on the structure under the runway, the stress state of the runway is not discussed here.

The relationship between the vertical displacement and time of each monitoring point on the central section of the tunnel is shown in Figure 12. This figure indicates that the vertical displacement of each monitoring point in the tunnel structure in the rigid pavement exceeds that in the flexible pavement. The application of the load to the flexible pavement results in satisfactory attenuation effects. Although the rigid pavement has satisfactory deformation resistance, it is not effective in diffusing force. The degree of deformation of each monitoring point in the tunnel varies with the position of the aircraft load. The vertical displacement of each monitoring point reaches the maximum when the aircraft taxis at the top of each point, further emphasizing that the aircraft load considerably impacts the safety and stability of the tunnel structure. The interference effect cannot be ignored. The curve in the figure shows that monitoring points S1 and S5 (top and bottom, respectively) undergo the largest deformation, verifying that the top and bottom of the tunnel are prone to damage after sustaining stress and forming weak areas. The peak displacement values of each monitoring point in different pavement structures are summarized and compared in

Table 3. Peak displacement values of the monitoring points.

Monitoring Points	Rigid Pavement	Flexible Pavement	Differential Value	Difference (%)
S1	1.156	0.972	0.184	15.9
S2	1.085	0.867	0.218	20.1
S3	1.077	0.837	0.240	22.3
S4	1.106	0.901	0.205	18.5
S5	1.187	1.031	0.156	13.1
S6	1.347	1.203	0.144	10.7
S7	1.456	1.316	0.140	9.6
S8	1.349	1.191	0.158	11.7

.

The table indicates that the largest difference in the overall displacement change of the two structures under the action of aircraft loads occurs at the left and right arch parts; the difference is less than that at the top and bottom plates. As shown in Figure 12b, the shear stress on the tunnel structure is considerably high, and extreme values in the positive and negative directions appear in the upper and lower arches, respectively. Because the pavement structures differ, the propagation attenuation effect on the force varies; consequently, the shear stress variation is also considerable. When subjected to shear stress, the tunnel structure undergoes relative dislocation deformation in the direction of the external force. The shear stresses on the upper and lower arches of the central section are relatively large; therefore, the changes in the vertical displacement at these locations are more drastic than those at other areas. Thus, the displacement of the tunnel structure arch, considering the two pavement types, exceeds those of the other parts.

4.3. Influence of Different Buried Depth on Tunnel Structure

Shallow buried tunnel projects are large-scale projects because a shallow depth is not conducive to the safety and stability of a structure. In airport construction, because the aircraft load has considerable vibration frequency and impact, the safety assessment of a shallow buried tunnel is necessary. This section presents eight burial depths (0.5D, 1D, 1.5D, 2D, 2.5D, 3D, 3.5D, and 4D), as well as the development of eight numerical calculation models considering shallow tunnel setting conditions. A380–800 aircraft loads were applied, and various indicators of the mechanical performance of the structure were analyzed to provide reference values for airport construction.

The vertical acceleration–time history curve of each monitoring point in the tunnel structure is shown in Figure 13. The figure indicates that the change in the frequency of each monitoring point is consistent with the taxiing law and varies with the movement of the aircraft load. The vibration amplitude of S1 is the largest, verifying that the top plate position is primarily disturbed. The curve further indicates that the disturbance effect of aircraft loads on the structure is inversely proportional to its burial depth, i.e., the vibration caused by the aircraft load decreases with the increase in depth. The top plate (S1) and bottom plate (S5) are selected to derive the first principal stress curve for further investigating the effect of burial depth on the tunnel structure.

As shown in Figure 14, the variation range of each curve differs, and the first principal stress curve of the two monitoring points tends to be gentle with increasing burial depth. This indicates that the effect of aircraft loads on the structure gradually weakens, and the deeper the embedment depth, the lesser the structure disturbance due to the load; consequently, safety is considerably improved. To improve the numerical simulation of the long-term state of the tunnel considering the influence of the self-weight stress field of the soil on the internal force state of the tunnel, geo-stress balance was implemented using the finite element software. Variations among the initial values of the first principal stress at the measuring points were observed. In terms of the numerical results, the soil stress increases with the burial depth, and its influence on the tunnel structure also increases. Hence, the influence of soil stress on the structure cannot be ignored in the design.

The peak points of vertical displacement and their fitting curves, considering the burial depth of the S1 (roof) position, are plotted in Figure 15. The disturbance of the tunnel structure caused by the aircraft load attenuates with the increase in burial depth. Because the displacements of the monitoring points are all positive in the Z direction, the peak vertical displacement of the roof is obtained to predict the undisturbed burial position of the tunnel structure. The fitting curve relationship is as follows:

$$\begin{array}{l}y=1.843-0.487exp\left(\frac{x}{47.618}\right)\\ R^2=0.996\end{array}$$

(11)

where y is the peak vertical displacement (mm); x is the burial depth of the tunnel (m); and R² is the coefficient of determination of the fit. The fitting curve obtained using MATLAB is shown in Figure 15b. The curve indicates that the peak vertical displacement of the tunnel structure buried at a depth of 63.37 m is 0. Considering the calculation formula of the aircraft dynamic load equation and the relevant clauses in the “Code for Design of Cement Concrete Pavement for Civil Airports (MH/T 5004–2010)” [29], the aircraft dynamic load equation of the A380–800 aircraft yields the largest value. Moreover, when the burial depth exceeds 64 m, the tunnel structure under the airport runway is not affected by the aircraft load.

5. Conclusions

In this study, a six-degree-of-freedom mathematical model of a “five-point-contact” aircraft was established, and the vibration balance equation of the mathematical model was deduced. Using ABAQUS software as a platform, a 3D finite element calculation model of the airport pavement–soil layer–tunnel was formulated to study the variation law of tunnel structures under airport runways due to different influencing factors. The main conclusions are as follows.

(1)

Different types of main landing gears have various forms and dynamic load equations because of the variations in tire pressure, number of wheels, and single wheel load. The distribution position of the main landing gear wheels also affects the structural force; a more concentrated wheel distribution causes greater interference to the structure. The roof position of the tunnel structure is stretched, and the vertical acceleration change range can reach 259% under the action of aircraft loads. This position is a damage-prone section that must be considered in the design.

(2)

Rigid pavements have large elastic moduli and are not easily deformed when subjected to aircraft loads. For the pavement substructure, flexible pavements have a satisfactory attenuation effect on force transmission. In the airport substructure design, the flexible pavement can be selected as the airport runway to ensure safety. The shear stress at the upper and lower arches of the tunnel structure varies depending on the road surface material. Moreover, relative dislocation and deformation occur at the arch. Owing to the differences in pavement materials, the shear stress on the structure varies, and the deformation at the arch is more severe than at other positions.

(3)

The interference effect of aircraft loads on the structure gradually weakens with increasing embedment depth. However, the soil stress on the tunnel structure increases with the burial depth, and the degree of interference of the soil stress on the structure is greater than that of the aircraft load. Therefore, during construction, the selection of corresponding support methods to cope with various soil depths and environments is necessary. The fitting curve indicates that when the burial depth exceeds 64 m, the movement of aircraft on the runway above does not affect the buried tunnel.