## 4 Millisecond Delay Determination  
To improve the contrast of the analysis, the five blast holes segments (3, 5, 7, 9 and 11, see Figure 5b) of the sixth part (other parts are labeled as 1, 2, 3, 4 and 5, see Figure 5a) of the double side drift method are selected. The numbers in Figure 5 indicate the sequence of excavation or detonation. The differential millisecond delay interval is set between the segments, and the charge arrangement is shown in Figure 5b. The transverse direction of the tunnel is the X-axis, and the Y-axis is the vertical direction. According to the actual engineering conditions, the size of the model is 18 \(\times\) 10 m, and the size of the tunnel is 13.56 \(\times\) 6.88 m. Gmsh software is used to build models and generate the computational mesh, as shown in Figure 6. In the FDEM, the numerical simulation results are affected by the mesh size. When the mesh size is small enough, the FDEM numerical simulation results tend to be stable [35; 43]. Therefore, using a small mesh size does not have a significant effect on the peak blasting vibration velocity. The normal direction around the model is fixed and set as absorbing boundaries. The values of the millisecond delay are 5, 10, 20, 30, 40, and 50 ms, respectively. By simulating these six working conditions, the peak blasting vibration velocity of the monitoring element on the right boundary of the model in Figure 6 is obtained. Then, the influence of a different millisecond delay on the peak blasting vibration velocity is analyzed to obtain a reasonable millisecond delay interval, which provides a reference for the selection and optimization of a reasonable millisecond delay. According to engineering experience and data, the meso-parameters input during calculation are shown in Table 1. \begin{table}
\begin{tabular}{c c} \hline \hline
**Parameters** & **Value** \\ \hline Triangular element & \\ \hline Density (kg/m\({}^{3}\)) & 2300 \\ Elastic Modulus \(E\) (GPa) & 20 \\ Friction angle (\({}^{\circ}\)) & 30 \\ Normal contact penalty parameter \(p_{n}\) (GPa) & 2000 \\ Tangential contact penalty parameter \(p_{t}\) (GPa/m) & 2000 \\ Poisson & 0.25 \\ \hline Joint element & \\ \hline Cohesive (MPa) & 1.5 \\ Tensile strength \(f_{t}\) (MPa) & 14 \\ Internal (\({}^{\circ}\)) & 30 \\ Tensile break energy release rate \(Gf_{I}\) (J/m\({}^{2}\)) & 30,000 \\ Shear fracture energy release rate \(Gf_{II}\) (J/m\({}^{2}\)) & 100,000 \\ Blasting load peak (surrounding holes/other holes) MPa & 104/316 \\ \hline \hline \end{tabular}
\end{table} 
Table 1: Meso input parameters of MultiFracS for tunnel calculation model. Figure 5: (a) Schematic diagram of blasting excavation construction process of the double side drift method; (b) schematic diagram of the blast hole arrangement and networked detonation in part VI. As shown in Figure 7, the exponential blasting load function \(P(t)\) is applied to the blast hole [44]:  
\[P(t)=P_{0}(e^{-at}-e^{-\beta t})/(e^{-at_{0}}-e^{-\beta t_{0}}) \tag{14}\]  
where \(t_{0}\) is the load rise time, which can be obtained by \(t_{0}=[1/(\beta-\alpha)]\log(\beta/\alpha)\), and \(\beta/\alpha\) is the control parameter of pressure decay. In this study, \(\beta/\alpha=100\) and \(t_{0}=10\) us are selected for simulation. \(P(t)\) only acts on the initial surface of the blast holes. Therefore, the gas flow into the crack is not considered. \(P_{0}\) is the peak load, which is obtained from the blasting technical parameters [45]. The peak load acting on the cutting holes, auxiliary holes, and bottom plate holes is all 316 MPa, and the peak load acting on the peripheral blast holes is 104 MPa. Figure 6: Schematic diagram of tunnel calculation model and monitoring unit for the double-side-drift method. Figure 7: Pressure–time curve of blasting load waveform. Figure 8 shows that the blasting vibration velocity of monitoring element A under different millisecond delays calculated by MultiFracS software.