## 4 Numerical Simulation  
### Numerical Model Establishment and Verification  
Based on the actual tunnel project, considering the influence of water pressure, the following numerical model was established by using the finite element software Midas-GTS. According to the Saint-Venant's Principle, the stress redistribution area generally appears within 3 to 5 times the width of the excavation. Therefore, in this model, the left and right sides of the model were four times the hole diameter, and the lower part was four times the hole diameter. The tunnel length was 100 m, the buried depth was 120 m, and the overall size was 120 m (X-direction) \(\times\) 100 m (Y-direction) \(\times\) 180 m (Z-direction), as shown in Figure 15.  
Model boundary conditions: The upper surface boundary conditions of the model are free, and the remaining surface boundary conditions are all normal constraints. Constitutive model: In previous studies, the Mohr-Coulomb (M-C) model, which does not fully consider the engineering characteristics of loess, has been widely used in the numerical analysis of loess tunnels. In this study, in order to improve the simulation accuracy, the Drucker-Prager (D-P) criterion was used. The D-P criterion considers the influence of the intermediate principal stress on the yield and failure; the yield surface is smooth, which can help to avoid the difficulty in the numerical calculation caused by the M-C criterion yield surface at the corners; that is, it avoids the so-called singularity [36; 37]. [PERSON] and [PERSON] put forward the Drucker-Prager criterion [38] in 1952. The expression is  
\[\alpha I_{1}+\sqrt{\beta_{2}}-k=0 \tag{2}\]  
Figure 15: Numerical simulation model. In the \(\pi\) plane, if the D-P criterion is the inscribed circle of the M-C, then \(\alpha\) and \(k\) can be expressed by the following formulas:  
\[\alpha=\frac{\sin\varphi}{\sqrt{3}\sqrt{3+\sin^{2}\varphi}} \tag{3}\]  
\[k=\frac{3c\cos\varphi}{\sqrt{3}\sqrt{3+\sin^{2}\varphi}} \tag{4}\]  
where \(I_{1}\) is the stress tensor first invariant, \(J_{2}\) is the second invariant of stress deviation, \(\varphi\) is the internal friction angle and \(c\) is cohesion. The loess was cored on-site, and the surrounding rock parameters were determined through laboratory tests (Table 1). The elastoplastic constitutive model was used for the initial support and secondary lining. Because the role of the advanced pre-grouting was to squeeze and strengthen the rock and soil through slurry veins, the simulation was carried out by improving the elastic modulus, cohesion and internal friction angle of the surrounding rock. The specific steps of the numerical simulation are as follows:  
1. Establish model grids, assign rock mass parameters, apply static boundary conditions, and then generate the initial ground stress field. Considering the influence of groundwater, apply hydrostatic pressure at the groundwater level. 2. Simulate the excavation process, using the three-step-seven-step method for excavation. The width of the upper step is 5 m, the height of the upper step from the vault is 2.2 m, the length of the upper step is 5 m, and the length of the lower step is 9 m. During the excavation process, each cycle of excavation footage is 1 steel arch with a spacing of 0.5 m. The upper arch, upper left pilot pit, upper right pilot pit, core soil and lower bench excavation are each staggered 3 m for parallel operation. The invert closed distance is 12 m from the tunnel face and the excavation step is 3 m.
3. The model grid is divided into excavation grid, primary support grid and secondary lining grid. When the tunnel is excavated, to which part of the excavation, the required part is passivated. After the excavation of each part of the steps is completed, the initial support grid is activated by the activation plate, the excavation module is passivated, the initial support is performed, before releasing 30% of the ground stress.