Simulation and analysis method for influence of small-spacing tunnel blasting on adjacent tunnel
By establishing a three-dimensional finite element model and using the fluid-structure interaction method, combined with fast Fourier transform and wavelet transform, the accuracy problem of tunnel blasting vibration simulation in existing technologies has been solved, enabling fine analysis and parameter optimization design of the dynamic response of tunnel structures.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- SOUTHWEST JIAOTONG UNIV
- Filing Date
- 2022-12-26
- Publication Date
- 2026-05-12
AI Technical Summary
Existing tunnel blasting design methods fail to accurately simulate the dynamic response of blasting vibrations to adjacent tunnels, neglect the propagation process of blasting shock waves in the fracture zone, and the boundary of the fracture zone is difficult to determine, resulting in significant discrepancies between the analysis results and the actual situation, and failing to reveal the dynamic response characteristics of the structure in depth.
A three-dimensional finite element model was established using the fluid-structure interaction method to construct rock mass and concrete material models, determine the explosive detonation time, conduct dynamic response analysis of the tunnel lining, and analyze the vibration spectrum under blasting load by fast Fourier transform and wavelet transform to simulate the velocity amplitude variation law of each monitoring point under blasting load.
It achieves precise simulation of blasting vibration, can guide the design of blasting parameters, quickly identify the most dangerous area of the structure, improve the accuracy and comprehensiveness of the analysis, is applicable to any type of blast hole, and the simulation results are highly reliable.
Smart Images

Figure CN115935753B_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of tunnel engineering technology, specifically to a method for simulating and analyzing the impact of blasting in small-clearance tunnels on adjacent tunnels. Background Technology
[0002] With the experience accumulated in tunnel and underground engineering design and construction in recent years, highway tunnel construction in my country is trending towards larger cross-sections. Large-section tunnels involve larger excavation cross-sections and larger amounts of explosives, which can significantly impact the structures of adjacent tunnels. The impact of tunnel excavation blasting vibrations on the structures of adjacent tunnels is a key research issue in tunnel engineering construction. Subsequent tunnel blasting can have a certain impact on the safety and stability of the lining structure and surrounding rock of preceding tunnels.
[0003] In existing tunnel blasting design methods, the selection of blasting parameters is guided by the goal of breaking the rock mass to be excavated and reducing disturbance to the surrounding rock. Blasting parameters are determined based on semi-empirical formulas, and then the blasting parameters are continuously adjusted according to the over-excavation and under-excavation of the surrounding rock and the degree of fragmentation after blasting. The dynamic impact of blasting vibration on adjacent tunnels is rarely considered, which may have an adverse effect on the safety and stability of adjacent tunnels.
[0004] Currently, there is a lack of practical theoretical formulas for calculating the dynamic response of nearby tunnels under blasting loads. Existing numerical simulation methods for tunnel blasting generally simplify the blasting load as a triangular load in time history curves, using empirical values for the loading and unloading periods, and determining the peak load using semi-empirical formulas. Ignoring the propagation process of the blast shock wave in the fractured zone, the equivalent elastic boundary is simplified to a circular envelope of the fractured zone of each blast hole. The blasting load is applied to the equivalent elastic boundary, and then, based on the principle of force equivalence, the blasting load is evenly distributed on the excavation profile surface according to a certain proportion before finite element calculation. While this method is simple to operate, the loading and unloading times and peak load in the blasting load time history curves are all empirical or semi-empirical values, which are somewhat arbitrary. The simplification of the equivalent elastic boundary deviates from the actual situation, and the boundary of the fractured zone is difficult to determine accurately. The accuracy of the results is affected by empirical parameters, leading to significant discrepancies with reality. Existing analysis methods generally only focus on the time history curves of structural velocities and peak stresses, failing to deeply reveal the dynamic response characteristics of the structure, thus hindering better optimization of design parameters and construction guidance. Summary of the Invention
[0005] To address the aforementioned shortcomings in existing technologies, this invention provides a simulation and analysis method for the impact of small-clearance tunnel blasting on adjacent tunnels. This method solves the problems of existing technologies being unable to simulate structural dynamic response, neglecting the propagation of blast shock waves in the fractured zone, and having difficulty accurately determining the boundary of the fractured zone.
[0006] To achieve the aforementioned objectives, the present invention employs the following technical solution: a method for simulating and analyzing the impact of small-clearance tunnel blasting on adjacent tunnels, comprising the following steps:
[0007] S1. Establish a three-dimensional finite element model for blasting construction of tunnels with small clearance;
[0008] S2. Construct the rock mass material model and concrete material model of the three-dimensional finite element model to obtain the complete three-dimensional finite element model;
[0009] S3. Free boundaries are used for the top of the complete three-dimensional finite element model, the surface of the preceding tunnel lining, the soil of the subsequent tunnel, and the surface of the initial support; non-reflective boundaries are used for the remaining boundaries.
[0010] S4. Determine the detonation time of the explosives in different sections and submit it to the solver for analysis.
[0011] S5. Perform dynamic response analysis of the preliminary tunnel lining to obtain the peak distribution map of the dynamic response of the monitoring section;
[0012] S6. Based on the peak distribution map of the dynamic response of the monitoring section, a preliminary tunnel vibration spectrum analysis was conducted to obtain the variation law of the velocity amplitude of each monitoring point in the time domain and frequency domain under the blasting load.
[0013] Furthermore, the specific implementation of step S1 is as follows:
[0014] S1-1. Establish finite element model components based on the actual tunnel cross-sectional geometry, left and right line clearance, tunnel depth, and construction sequence; wherein, the distance between the left and right boundaries and the lower boundary of the model and the tunnel is 1.5D, where D is the tunnel diameter;
[0015] S1-2. Establish blast holes at the corresponding positions on the tunnel face of the subsequent tunnel according to the actual working conditions;
[0016] S1-3, According to the formula:
[0017]
[0018] P1=C0+C1μ+C2μ 2 +C3μ 3 +(C4+C5μ+C6μ 2 E1
[0019] The pressure P of the explosive gas and the pressure P1 of the air are obtained; where V is the relative volume; E0 is the initial internal energy density; A, B, R1, R2, and ω are parameters of the equation of state; μ = ρ / ρ0 - 1, where ρ is the current air density and ρ0 is the air density at the initial moment; C0, C1, C3, C4, C5, and C6 are coefficients of the polynomial equation; E1 is the internal energy density; and V0 is the initial relative volume.
[0020] S1-4. Based on the pressure P of the explosive gas and the pressure P1 of the air, complete the setting of air and explosive, and complete the construction of a three-dimensional finite element model for tunnel blasting construction with small clearance.
[0021] Furthermore, the specific implementation of step S2 is as follows:
[0022] S2-1, According to the formula:
[0023]
[0024] A rock mass material model simulated using a bilinear elastoplastic constitutive model was obtained; where σ0 is the initial yield stress of the material under external forces; σ y ε is the final yield stress; ε is the strain rate under a certain external load; C and p are the strain rate parameters of the material itself; E p E represents the hardening modulus of plastics. p =E tan E / (EE tan E is Young's modulus. tan β is the tangent modulus; β is the hardening parameter. When it is 1, the material exhibits isotropic hardening properties, and when it is 0, it exhibits kinematic hardening properties. For effective plastic strain;
[0025] S2-2. For the initial support and secondary lining of the concrete material model, a bilinear elastoplastic material constitutive model is used.
[0026] Furthermore, the specific implementation of step S5 is as follows:
[0027] S5-1. The face of the tunnel face in the blasting zone of the subsequent tunnel is taken as the cross section of the preceding tunnel as the ZDM-0-0 section.
[0028] S5-2. Along the longitudinal direction of the tunnel, select a section every 1m within ±2 of the ZDM-0-0 section. The direction of tunnel excavation for the subsequent tunnel is negative, and the direction of tunnel excavation for the preceding tunnel is positive. This will result in monitoring sections ZDM-(-2)-(-2), ZDM-(-1)-(-1), ZDM-1-1, and ZDM-2-2. Monitoring points will be selected at key locations on each monitoring section.
[0029] S5-3. Obtain the comprehensive vibration velocity cloud map of the secondary lining structure of the tunnel at different times;
[0030] S5-4. Extract the dynamic response-time relationship curves of each monitoring point from the comprehensive vibration velocity cloud map to obtain the peak distribution map of the dynamic response of the monitoring section.
[0031] Furthermore, the specific implementation of step S5-4 is as follows:
[0032] S5-4-1. Input the distance coordinates of each monitoring point on the same cross section, the coordinates of the points at equal intervals on the outline of the preceding tunnel, and the coordinates of the monitoring points into MATLAB.
[0033] S5-4-2. Use the max(A(2,:)) function to obtain the peak value of the dynamic response; where A is a matrix composed of dynamic response-time model data, the first row is time data, and the second row is dynamic response data;
[0034] S5-4-3. Based on the time data and dynamic response data, construct the dynamic response-time relationship curve using the cross-sectional contour coordinate points, and obtain the derivative of the tunnel cross-sectional contour curve to obtain the position derivative of each monitoring point; obtain the orthogonal direction between the monitoring point and the cross-sectional contour.
[0035] S5-4-4. Multiply the peak value of the dynamic response at the monitoring point by a coefficient to obtain the distance from the data point to the monitoring point;
[0036] S5-4-5. Based on the distance from the data point to the monitoring point and the orthogonal direction between the monitoring point and the cross-sectional profile, obtain the coordinates of the data point location;
[0037] S5-4-6. Repeat steps S5-4-2 to S5-4-5 until the coordinates of all monitoring points on the same cross section are obtained; connect the coordinates of the data points corresponding to different monitoring points to obtain the peak distribution map of the dynamic response of the monitoring cross section.
[0038] Furthermore, the specific implementation of step S6 is as follows:
[0039] S6-1. Based on the peak distribution map of the dynamic response of the monitoring section, perform a fast Fourier transform on the dynamic response time history curve of each monitoring point to obtain the dynamic response spectrum at each monitoring point.
[0040] S6-2, According to the formula:
[0041]
[0042]
[0043]
[0044]
[0045]
[0046] Obtain the wavelet basis function ψ a,b (t); where a is the scale factor; b is the translation amount; and t is time; Let w be the Fourier transform of ψ(t), where w is the frequency; ψ(t) is a finite-energy signal, and L is the frequency.2 L(R) is a square-integrable space; L(R) is an integrable space.
[0047] S6-3, According to the formula:
[0048]
[0049] Perform continuous wavelet transform; where f(t) is the signal to be analyzed. yes The conjugate function of WT f (a,b) are wavelet transform coefficients; R is the real number field;
[0050] S6-4. The optimal wavelet basis function is determined by using the modulus average method based on the continuous wavelet transform coefficients and the wavelet basis function.
[0051] S6-5. Obtain the three-dimensional time-frequency diagram of each monitoring point by performing continuous wavelet transform based on the optimal wavelet basis function;
[0052] S6-6. Analyze the wavelet three-dimensional time-frequency diagrams of each monitoring point to obtain the variation law of the velocity amplitude of each monitoring point in the time domain and frequency domain under the blasting load.
[0053] The beneficial effects of this invention are as follows:
[0054] 1. Using fluid-structure interaction to simulate blasting construction avoids the tedious calculation of equivalent loads and the errors caused by using some ideal assumptions. It is applicable to any type of blast hole and can effectively simulate the micro-delay effect and group hole effect of blasting. The simulation results can guide the design of blasting parameters such as micro-delay blasting time interval and charge coefficient.
[0055] 2. The peak curve of the dynamic response on the cross section can be obtained intuitively, and the most dangerous area of the structure can be quickly identified. At the same time, the fast Fourier transform and wavelet transform are used to analyze the vibration spectrum. Compared with conventional analysis methods, the analysis has great advantages in terms of precision, accuracy and comprehensiveness. Attached Figure Description
[0056] Figure 1 This is a flowchart of the present invention;
[0057] Figure 2 This is a diagram illustrating the steps of drilling holes in an embodiment of the present invention;
[0058] Figure 3 This is a finite element three-dimensional model diagram of an embodiment of the present invention;
[0059] Figure 4 This is a top view of the secondary lining monitoring section layout according to an embodiment of the present invention;
[0060] Figure 5 This is a diagram showing the layout of monitoring points for the secondary lining in an embodiment of the present invention.
[0061] Figure 6 This is the composite vibration velocity cloud map at 17.3ms in an embodiment of the present invention;
[0062] Figure 7 This is the time history curve of the comprehensive vibration velocity at the arch monitoring point of the ZDM-0-0 monitoring section in an embodiment of the present invention;
[0063] Figure 8 This is a peak distribution diagram of dynamic vibration velocity at the monitoring section of ZDM-0-0 according to an embodiment of the present invention;
[0064] Figure 9 The first embodiment of the present invention is the vibration velocity time history curve in the X direction of the monitoring point on the right arch waist of the tunnel.
[0065] Figure 10 The Fourier spectrum of vibration velocity in the X direction at the monitoring point of the right arch waist of the tunnel is shown in this embodiment of the invention.
[0066] Figure 11 This is the process for finding the optimal scale of the wavelet basis function in an embodiment of the present invention;
[0067] Figure 12 This is a sym5 wavelet scale-modulus average value diagram according to an embodiment of the present invention;
[0068] Figure 13 This is a contour map of wavelet coefficients in sym5, as described in an embodiment of the present invention.
[0069] Figure 14 The X-axis vibration velocity wavelet three-dimensional time-frequency diagram of the monitoring point on the right arch waist of the secondary lining of the tunnel in this embodiment of the invention is shown. Detailed Implementation
[0070] The specific embodiments of the present invention are described below to enable those skilled in the art to understand the present invention. However, it should be understood that the present invention is not limited to the scope of the specific embodiments. For those skilled in the art, various changes are obvious as long as they are within the spirit and scope of the present invention as defined and determined by the appended claims. All inventions utilizing the concept of the present invention are protected.
[0071] like Figure 1 As shown, a method for simulating and analyzing the impact of small-clearance tunnel blasting on adjacent tunnels includes the following steps:
[0072] S1. Establish a three-dimensional finite element model for blasting construction of tunnels with small clearance;
[0073] S2. Construct the rock mass material model and concrete material model of the three-dimensional finite element model to obtain the complete three-dimensional finite element model;
[0074] S3. Free boundaries are used for the top of the complete three-dimensional finite element model, the surface of the preceding tunnel lining, the soil of the subsequent tunnel, and the surface of the initial support; non-reflective boundaries are used for the remaining boundaries.
[0075] S4. Determine the detonation time of the explosives in different sections and submit it to the solver for analysis.
[0076] S5. Perform dynamic response analysis of the preliminary tunnel lining to obtain the peak distribution map of the dynamic response of the monitoring section;
[0077] S6. Based on the peak distribution map of the dynamic response of the monitoring section, a preliminary tunnel vibration spectrum analysis was conducted to obtain the variation law of the velocity amplitude of each monitoring point in the time domain and frequency domain under the blasting load.
[0078] The specific implementation method of step S1 is as follows:
[0079] S1-1. Establish finite element model components based on the actual tunnel cross-sectional geometry, left and right line clearance, tunnel depth, and construction sequence; wherein, the distance between the left and right boundaries and the lower boundary of the model and the tunnel is 1.5D, where D is the tunnel diameter;
[0080] S1-2. Establish blast holes at the corresponding positions on the tunnel face of the subsequent tunnel according to the actual working conditions;
[0081] S1-3, According to the formula:
[0082]
[0083] P1=C0+C1μ+C2μ 2 +C3μ 3 +(C4+C5μ+C6μ 2 E1
[0084] The pressure P of the explosive gas and the pressure P1 of the air are obtained; where V is the relative volume; E0 is the initial internal energy density; A, B, R1, R2, and ω are parameters of the equation of state; μ = ρ / ρ0 - 1, where ρ is the current air density and ρ0 is the air density at the initial moment; C0, C1, C3, C4, C5, and C6 are coefficients of the polynomial equation; E1 is the internal energy density; and V0 is the initial relative volume.
[0085] S1-4. Based on the pressure P of the explosive gas and the pressure P1 of the air, complete the setting of air and explosive, and complete the construction of a three-dimensional finite element model for tunnel blasting construction with small clearance.
[0086] The specific implementation method of step S2 is as follows:
[0087] S2-1, According to the formula:
[0088]
[0089] A rock mass material model simulated using a bilinear elastoplastic constitutive model was obtained; where σ0 is the initial yield stress of the material under external forces; σ y ε is the final yield stress; ε is the strain rate under a certain external load; C and p are the strain rate parameters of the material itself; E p E represents the hardening modulus of plastics. p =E tan E / (EE tan E is Young's modulus. tan β is the tangent modulus; β is the hardening parameter. When it is 1, the material exhibits isotropic hardening properties, and when it is 0, it exhibits kinematic hardening properties. For effective plastic strain;
[0090] S2-2. For the initial support and secondary lining of the concrete material model, a bilinear elastoplastic material constitutive model is used.
[0091] The specific implementation method of step S5 is as follows:
[0092] S5-1. The face of the tunnel face in the blasting zone of the subsequent tunnel is taken as the cross section of the preceding tunnel as the ZDM-0-0 section.
[0093] S5-2. Along the longitudinal direction of the tunnel, select a section every 1m within ±2 of the ZDM-0-0 section. The direction of tunnel excavation for the subsequent tunnel is negative, and the direction of tunnel excavation for the preceding tunnel is positive. This will result in monitoring sections ZDM-(-2)-(-2), ZDM-(-1)-(-1), ZDM-1-1, and ZDM-2-2. Monitoring points will be selected at key locations on each monitoring section.
[0094] S5-3. Obtain the comprehensive vibration velocity cloud map of the secondary lining structure of the tunnel at different times;
[0095] S5-4. Extract the dynamic response-time relationship curves of each monitoring point from the comprehensive vibration velocity cloud map to obtain the peak distribution map of the dynamic response of the monitoring section.
[0096] The specific implementation method of step S5-4 is as follows:
[0097] S5-4-1. Input the distance coordinates of each monitoring point on the same cross section, the coordinates of the points at equal intervals on the outline of the preceding tunnel, and the coordinates of the monitoring points into MATLAB.
[0098] S5-4-2. Use the max(A(2,:)) function to obtain the peak value of the dynamic response; where A is a matrix composed of dynamic response-time model data, the first row is time data, and the second row is dynamic response data;
[0099] S5-4-3. Based on the time data and dynamic response data, construct the dynamic response-time relationship curve using the cross-sectional contour coordinate points, and obtain the derivative of the tunnel cross-sectional contour curve to obtain the position derivative of each monitoring point; obtain the orthogonal direction between the monitoring point and the cross-sectional contour.
[0100] S5-4-4. Multiply the peak value of the dynamic response at the monitoring point by a coefficient to obtain the distance from the data point to the monitoring point;
[0101] S5-4-5. Based on the distance from the data point to the monitoring point and the orthogonal direction between the monitoring point and the cross-sectional profile, obtain the coordinates of the data point location;
[0102] S5-4-6. Repeat steps S5-4-2 to S5-4-5 until the coordinates of all monitoring points on the same cross section are obtained; connect the coordinates of the data points corresponding to different monitoring points to obtain the peak distribution map of the dynamic response of the monitoring cross section.
[0103] The specific implementation method of step S6 is as follows:
[0104] S6-1. Based on the peak distribution map of the dynamic response of the monitoring section, perform a fast Fourier transform on the dynamic response time history curve of each monitoring point to obtain the dynamic response spectrum at each monitoring point.
[0105] S6-2, According to the formula:
[0106]
[0107]
[0108]
[0109]
[0110]
[0111] Obtain the wavelet basis function ψ a,b (t); where a is the scale factor; b is the translation amount; t is time; ψ(w) is the Fourier transform of ψ(t), w is the frequency; ψ(t) is a finite-energy signal, L 2 L(R) is a square-integrable space; L(R) is an integrable space.
[0112] S6-3, According to the formula:
[0113]
[0114] Perform continuous wavelet transform; where f(t) is the signal to be analyzed. yes The conjugate function of WT f(a,b) are wavelet transform coefficients; R is the real number field;
[0115] S6-4. The optimal wavelet basis function is determined by using the modulus average method based on the continuous wavelet transform coefficients and the wavelet basis function.
[0116] S6-5. Obtain the three-dimensional time-frequency diagram of each monitoring point by performing continuous wavelet transform based on the optimal wavelet basis function;
[0117] S6-6. Analyze the wavelet three-dimensional time-frequency diagrams of each monitoring point to obtain the variation law of the velocity amplitude of each monitoring point in the time domain and frequency domain under the blasting load.
[0118] In one embodiment of the present invention, the Chengdu Longquanshan No. 1 Tunnel has a maximum excavation face area of approximately 230㎡, which is a large-section highway tunnel. The net distance between the left and right lines of the tunnel is 22-44m. The tunnel body mainly passes through mudstone strata with a basic rock mass quality grade of V. It also passes through silty mudstone with a basic rock mass quality grade of IV. The overall surrounding rock conditions of the tunnel are poor.
[0119] The tunnel was excavated using the CD method, and the blast holes were arranged as follows: Figure 2 As shown in Table 1, the blasting parameters are as follows.
[0120] Table 1
[0121]
[0122] To investigate the dynamic response of the preceding tunnel during blasting excavation of the subsequent tunnel under different clearance conditions, and to reveal the dynamic response characteristics of the preceding tunnel under blasting load, a three-dimensional numerical model of the Longquanshan Tunnel under different clearance conditions (22m, 27m, 32m) was established to analyze the dynamic response of the preceding tunnel under the blasting action of the subsequent tunnel.
[0123] A three-dimensional finite element model of a tunnel with small clearance for blasting construction was established. The maximum excavation width of the tunnel is 21.25m, and the model dimensions are 120×71.45×60m. All elements in the model are Solid164 solid elements, and the model mainly includes rock mass, initial support, and secondary lining. The clearances between the tunnel sections are 22m, 27m, and 32m. The left and right boundaries of the model are 30m from the tunnel, and the bottom boundary is 30m from the tunnel. The tunnel depth is set to 30m. The three-dimensional model is shown below. Figure 3 As shown.
[0124] According to the actual working conditions, blast holes are established at the corresponding positions on the tunnel face, and air and explosives are established separately. Considering that air and explosives cannot overlap, the position of the explosives must be taken into account for the air. The pressure P of the explosive gas and the pressure P1 of the air are obtained; the values of state parameters A, B, R1, R2, and ω are shown in Table 2, and the values of coefficients C0, C1, C3, C4, C5, and C6 of the polynomial equation are shown in Table 3.
[0125] Table 2
[0126]
[0127] Table 3
[0128]
[0129] The rock mass material model and concrete material model of the three-dimensional finite element model were constructed to obtain the complete three-dimensional finite element model; the surrounding rock parameters are shown in Table 4, and the concrete strength parameters are shown in Table 5.
[0130] Table 4
[0131]
[0132] Table 5
[0133]
[0134] The delay time for the blasting load in each adjacent segment is set to 3ms, and the solution is submitted for analysis.
[0135] Dynamic response analysis of the lining of the preceding tunnel was conducted, and the peak distribution map of the dynamic response at the monitoring sections was obtained. To study the variation law of vibration velocity of each mass point in the adjacent lining structure during the blasting of the subsequent tunnel, the ZDM-0-0 section was selected as the cross section of the preceding tunnel corresponding to the working face of the subsequent tunnel blasting zone. A section was selected every 1m within ±2 of this section along the tunnel longitudinal direction, with the direction of subsequent tunnel excavation being negative and the opposite being positive, resulting in monitoring sections ZDM-(-2)-(-2), ZDM-(-1)-(-1), ZDM-1-1, and ZDM-2-2. Monitoring points were selected at key locations on each monitoring section. More key monitoring points were selected on the blast-facing side of each monitoring section of the preceding tunnel, and relatively fewer on the blast-back side. The monitoring sections and monitoring points are as follows: Figure 4 and Figure 5 As shown.
[0136] Obtain comprehensive vibration velocity contour maps of the secondary lining structure of the tunnel at different times. Comprehensive vibration velocity can significantly reflect the vibration acceleration response of the lining structure, which is beneficial for investigating the propagation state of blasting vibration elastic waves in the lining structure. The comprehensive vibration velocity contour map at 17.3 ms is shown below. Figure 6 As shown.
[0137] Taking the ZDM-0-0 monitoring section as an example, the comprehensive vibration velocity time history curve of the arch monitoring point of the ZDM-0-0 monitoring section is as follows: Figure 7As shown, the curve exhibits a clear "five-peak" pattern, with the peak times occurring approximately between 0.5–2 ms, 4–5 ms, 8–9.5 ms, 12–13.5 ms, and 16–17.5 ms, respectively. This "five-peak" pattern effectively reflects the "micro-delay effect" of the blasting load. Based on this result, the micro-delay blasting time interval can be rationally determined to avoid energy concentration and reduce the impact of blasting on the tunnel.
[0138] The peak-time vibration velocity signals from each monitoring point on the same cross-section were imported into a self-written MATLAB program in matrix form. The distance coordinates of equally spaced points on the tunnel outline and the coordinates of the monitoring points were also imported into the MATLAB program to automatically plot the distribution map of the peak dynamic vibration velocity of the monitoring cross-section. Figure 8 As shown. Figure 8 The vibration can clearly reflect the dynamic response state of the structure, which is of great significance for studying the safety of blast vibration response. The maximum value of the comprehensive vibration velocity of each monitoring section appears on the right side wall. Therefore, key monitoring should be carried out on the right arch foot, right side wall and right arch waist on the blast-facing side where the vibration velocity peak is larger. If necessary, reinforcement measures should be taken in time to prevent damage and ensure structural safety.
[0139] Vibration spectrum analysis of the pilot tunnel: Taking the "right arch waist monitoring point" in the ZDM-0-0 monitoring section with the largest dynamic response peak as an example, a fast Fourier transform of its dynamic response time history curve is derived to obtain the dynamic response spectrum at each monitoring point. The X-direction vibration velocity time history curves of the pilot tunnel's right arch waist monitoring point under different clearances are shown below. Figure 9 As shown, the Fourier spectrum of the vibration velocity in the X direction at the monitoring point on the right arch waist of the tunnel after Fast Fourier Transform is as follows: Figure 10 As shown, the Fourier spectrum amplitudes in all directions at different clearances are relatively large in the range of 0–100 Hz, with the dominant Fourier frequency appearing in the range of 10–50 Hz. The monitoring points on the right wall at clearances of 27 m and 32 m show attenuations of 27.27% and 48.37% respectively compared to the monitoring point at clearance of 22 m. This indicates that as the clearance increases, the Fourier spectrum amplitude of the structure near the right wall continuously decreases, revealing that the blast stress wave undergoes dissipation during propagation in the structure.
[0140] Select the optimal wavelet basis function. The optimal scale is determined using the modulus-average method combined with wavelet coefficient contour maps as evidence; the process is as follows: Figure 11 Taking the sym5 wavelet basis as an example, wavelet coefficients at different scales can be obtained using continuous wavelet transform. Scale-mode mean plots and contour maps can then be plotted. (See...) Figure 12 and Figure 13 Since the scale-mode average reaches its maximum at a scale of 30, 30 can be initially selected as the optimal scale. From the contour map, it can be found that the changes in brightness are most obvious at a scale of 30. Therefore, the scale of 30 can be selected as the optimal scale for analyzing the vibration velocity curve using the sym5 wavelet.
[0141] Using the selected optimal wavelet basis function, a continuous wavelet transform is performed to obtain the wavelet three-dimensional time-frequency plots for each monitoring point, as shown below. Figure 14 As shown, under the blasting load, the maximum vibration velocity in the X direction of the right arch waist of the tunnel is distributed between 0 and 0.03 s, spanning five time intervals. The maximum value occurs around 0.0145 s, within the fourth time interval, with an amplitude of 10.69. This reflects, to some extent, that the "micro-delay" and "time-delay" blasting considered when applying the blasting load effectively dispersed the blasting vibration energy. In terms of the energy distribution with frequency, the energy is mainly distributed in the frequency range of 5–300 Hz, with the dominant vibration frequency around 17.61 Hz.
[0142] This invention uses fluid-structure interaction to simulate blasting operations, avoiding cumbersome equivalent load calculations and errors caused by ideal assumptions. It is applicable to any type of borehole and can effectively simulate blasting differential effects and multi-hole effects. The simulation results can guide the design of blasting parameters such as differential blasting time intervals and charge coefficients. The analysis method provided by this invention can intuitively derive the peak dynamic response curve on the cross-section, quickly identifying the most dangerous area of the structure. Furthermore, it uses Fast Fourier Transform and Wavelet Transform to analyze the vibration spectrum, offering significant advantages in terms of detail, accuracy, and comprehensiveness compared to conventional analysis methods.
Claims
1. A method for simulating and analyzing the impact of small-clearance tunnel blasting on adjacent tunnels, characterized in that, Includes the following steps: S1. Based on the pressure of the explosive gas and the air pressure, complete the setting of air and explosive, and complete the construction of a three-dimensional finite element model for tunnel blasting construction with small clearance. S2. Construct the rock mass material model and concrete material model of the three-dimensional finite element model to obtain the complete three-dimensional finite element model; S3. Free boundaries are used for the top of the complete three-dimensional finite element model, the surface of the preceding tunnel lining, the soil of the subsequent tunnel, and the surface of the initial support; non-reflective boundaries are used for the remaining boundaries. S4. Determine the detonation time of the explosives in different sections and submit it to the solver for analysis. S5. Perform dynamic response analysis of the preliminary tunnel lining to obtain the peak distribution map of the dynamic response of the monitoring section; S6. Based on the peak distribution map of dynamic response of the monitoring section, conduct preliminary tunnel vibration spectrum analysis to obtain the variation law of velocity amplitude of each monitoring point in the time domain and frequency domain under blasting load; The specific implementation method of step S5 is as follows: S5-1. The face of the tunnel face in the blasting zone of the subsequent tunnel is taken as the cross section of the preceding tunnel as the ZDM-0-0 section. S5-2. Along the longitudinal direction of the tunnel, select a section every 1m within ±2 of the ZDM-0-0 section. The direction of tunnel excavation for the subsequent tunnel is negative, and the direction of tunnel excavation for the preceding tunnel is positive. This will result in monitoring sections ZDM-(-2)-(-2), ZDM-(-1)-(-1), ZDM-1-1, and ZDM-2-2. Monitoring points will be selected at key locations on each monitoring section. S5-3. Obtain the comprehensive vibration velocity cloud map of the secondary lining structure of the tunnel at different times; S5-4. Extract the dynamic response-time relationship curves of each monitoring point from the comprehensive vibration velocity cloud map to obtain the peak distribution map of the dynamic response of the monitoring section; The specific implementation method of step S6 is as follows: S6-1. Based on the peak distribution map of the dynamic response of the monitoring section, perform a fast Fourier transform on the dynamic response time history curve of each monitoring point to obtain the dynamic response spectrum at each monitoring point. S6-2, According to the formula: Obtain wavelet basis functions ;in, a Scale factor; b This is the translation amount; t For time; for Fourier transform, For frequency; For energy-limited signals, It is a square-integrable space; It is an integrable space; S6-3, According to the formula: Perform continuous wavelet transform; where, The signal to be analyzed, yes The conjugate function, These are the wavelet transform coefficients; R For the real number field; S6-4. The optimal wavelet basis function is determined by using the modulus average method based on the continuous wavelet transform coefficients and the wavelet basis function. S6-5. Obtain the three-dimensional time-frequency diagram of each monitoring point by performing continuous wavelet transform based on the optimal wavelet basis function; S6-6. Analyze the wavelet three-dimensional time-frequency diagrams of each monitoring point to obtain the variation law of the velocity amplitude of each monitoring point in the time domain and frequency domain under the blasting load.
2. The method for simulating and analyzing the impact of small-clearance tunnel blasting on adjacent tunnels according to claim 1, characterized in that, The specific implementation method of step S1 is as follows: S1-1. Establish finite element model components based on the actual tunnel cross-sectional geometry, left and right line clearance, tunnel depth, and construction sequence; wherein, the distance between the left and right boundaries and the lower boundary of the model and the tunnel is 1.5D, where D is the tunnel diameter; S1-2. Establish blast holes at the corresponding positions on the tunnel face of the subsequent tunnel according to the actual working conditions; S1-3, According to the formula: The pressure of the explosive gas is obtained P and air pressure ;in, Relative volume; The initial internal energy density; These are the parameters of the state equation; , The current density of air. The air density at the initial moment; The coefficients of the polynomial equation; Internal energy density; This represents the initial relative volume; e It is a natural constant; S1-4. Based on the pressure of the explosive gas P and air pressure Complete the setup of air and explosives, and construct a three-dimensional finite element model for tunnel blasting with small clearance.
3. The method for simulating and analyzing the impact of small-clearance tunnel blasting on adjacent tunnels according to claim 2, characterized in that, The specific implementation method of step S2 is as follows: S2-1, According to the formula: A rock mass material model simulated using a bilinear elastoplastic material constitutive model was obtained; among which, The initial yield stress of a material when subjected to external forces; This is the final yield stress; The strain rate under a given external load; C and The strain rate parameter is the material's own strain rate. For the hardening modulus of plastic, , For Young's modulus, Tangent modulus; This is the hardening parameter. When it is set to 1, the material exhibits isotropic hardening properties, and when it is set to 0, it exhibits kinematic hardening properties. For effective plastic strain; S2-2. For the initial support and secondary lining of the concrete material model, a bilinear elastoplastic material constitutive model is used.
4. The method for simulating and analyzing the impact of small-clearance tunnel blasting on adjacent tunnels according to claim 3, characterized in that, The specific implementation method of step S5-4 is as follows: S5-4-1. Input the distance coordinates of each monitoring point on the same cross section, the coordinates of the points at equal intervals on the outline of the preceding tunnel, and the coordinates of the monitoring points into MATLAB. S5-4-2. Use the max(A(2,:)) function to obtain the peak value of the dynamic response; where A is a matrix composed of dynamic response-time model data, the first row is time data, and the second row is dynamic response data; S5-4-3. Based on the time data and dynamic response data, construct the dynamic response-time relationship curve using the cross-sectional contour coordinate points, and obtain the derivative of the tunnel cross-sectional contour curve to obtain the position derivative of each monitoring point; obtain the orthogonal direction between the monitoring point and the cross-sectional contour. S5-4-4. Multiply the peak value of the dynamic response at the monitoring point by a coefficient to obtain the distance from the data point to the monitoring point; S5-4-5. Based on the distance from the data point to the monitoring point and the orthogonal direction between the monitoring point and the cross-sectional profile, obtain the coordinates of the data point location; S5-4-6. Repeat steps S5-4-2 to S5-4-5 until the coordinates of all monitoring points on the same cross section are obtained; connect the coordinates of the data points corresponding to different monitoring points to obtain the peak distribution map of the dynamic response of the monitoring cross section.