## 3 Numerical Simulation of Tunnel Face Reinforcement Scheme  
### Model Setup  
As the launching section of the tunnel passes through muddy soil layer, reinforcement of the tunnel face is necessary prior to shield construction. To determine a suitable reinforcement plan, MIDAS/GTS NX 2022R1 software was used to conduct 3D numerical simulations of different reinforcement schemes and the optimal plan was selected for implementation. Given that the tunnel face is primarily located in muddy and clayey soils, two reinforcement schemes were modeled: Scheme 1: 4 m long \(\Phi\)850@600 triple-axis mixing piles +2 m thick backfill layer of cohesive soil as counter-pressure. Scheme 2: Replacement with cohesive soil down to the base of the muddy soil layer +4 m long \(\Phi\)600@450 high-pressure jet grouting piles +2 m thick cohesive soil backfill for counter-pressure. The reinforcement zone (with overburden fill, covering segments 1-8) is 9.6 m long; the transition zone between the reinforced and unreinforced areas (segment 9) is 1.2 m; and the unreinforced zone (without overburden fill, segments 10-30) is 25.2 m long. According to the Technical Code for Monitoring in Urban Rail Transit Projects (GB 50911-2013) [21], the early warning threshold for ground deformation in this project is 30 mm. MIDAS/GTS is a general finite element analysis software widely used in geotechnical engineering, particularly for shield tunnel simulations. The geometry model in GTS/NX supports a variety of element types, including scalar, 1D, 2D, solid, interface, and rigid link or interpolation elements. These can be applied in linear and nonlinear stress analysis, seepage analysis, and stress-seepage coupled analysis. In MIDAS/GTS, the virtual work principle is expressed using stress-deformation equations, commonly constrained by stress-strain relationships. The Hu-Washizu variational principle is applied. \[\delta G_{\text{ext}}=\int_{\Omega}\Bigl{[}(\
abla\delta\mathbf{u})^{T} \boldsymbol{\sigma}+\delta\boldsymbol{\varepsilon}^{T}(\boldsymbol{\sigma}( \boldsymbol{\varepsilon})-\boldsymbol{\sigma})+\delta\boldsymbol{\sigma}^{T}( \
abla\mathbf{u}-\boldsymbol{\varepsilon})\Bigr{]}\text{d}\Omega \tag{1}\]  
Assuming that the constitutive equation can always satisfy the strain-stress relationship and the coordination relationship between \(\boldsymbol{\varepsilon}\) and \(\
abla u\), it can be deduced as a general expression of the virtual work principle:  
\[\delta G_{ext}=\int_{\Omega}\left(\
abla\delta\mathbf{u}\right)^{T} \boldsymbol{\sigma}(\mathbf{u})\text{d}\Omega \tag{2}\]  
In the finite element method, the displacement \(\mathbf{u}\) represented by the shape function difference is:  
\[u^{h}=Nd^{e} \tag{3}\]  
Using \(\boldsymbol{\varepsilon}^{k}=\
abla u^{h}=B\boldsymbol{\varepsilon}^{k}\), the principle equation of virtual work in the unit can be obtained:  
\[\delta G_{\text{ext}}=\delta\mathbf{d}^{T}\mathbf{F}=\delta\mathbf{d}^{T} \left(\sum_{\boldsymbol{\varepsilon}}\int_{\Omega_{\boldsymbol{\varepsilon} }}\mathbf{B}^{T}\mathbf{D}\mathbf{B}\text{d}\Omega\right)\mathbf{d}=\delta \mathbf{d}^{T}\mathbf{K}\mathbf{d} \tag{4}\]  
The stiffness matrix \(K^{e}\) of the element is:  
\[k^{e}=\int_{\Omega_{\boldsymbol{\varepsilon}}}B^{T}DB\text{d}\Omega \tag{5}\]  
Considering that the reinforced tunnel section is 9.6 m long and the tunnel diameter is 6.2 m, the calculation model length is typically 3-5 times the tunnel diameter. Thus, a model 36 m long, 30 m wide, and 33.29 m high (sum of all strata thicknesses) was selected. Displacement constraints were applied using the software's automatic constraint function: the model ends were constrained in the x-direction, the other three side faces in both xand y directions, and the base in x, y, and z directions. The excavation step per cycle was set to 1.2 m, equivalent to one analysis step.