Numerical Simulation Method for Stability Analysis of Jointed Tunnels Based on Discontinuous Deformation Analysis

Abstract

This invention discloses a numerical simulation method for jointed tunnel stability analysis based on discontinuous deformation analysis, comprising the following steps: Step 1: Based on the original DDA input program dl, the original DDA is modified in the C language environment to obtain different joint parameters from different probability density functions and generate joint lines and excavation chamber boundaries; Step 2: The joint lines generated by the dl program in Step 1 are imported into the dc program to realize the system program cutting the blocks; Step 3: The overall equilibrium equation of the block system considering damping is compiled into the df program to realize the block displacement solution; Step 4: The overall change and individual deformation of the block system after calculation are plotted, and the dg program provides a visualization interface; This method can replace the original DDA program, reflect more complex jointed block systems, explore the influence of joint physical and mechanical parameters on the deformation of tunnel surrounding rock, and has the characteristics of high simulation accuracy and strong practicality.

Description

Technical Field

[0001] This invention relates to the field of computer-aided design technology in geotechnical engineering, specifically to a numerical simulation method for stability analysis of jointed tunnels based on discontinuous deformation analysis. Background Technology

[0002] When studying geotechnical engineering problems, due to the special nature of geotechnical engineering, although field tests can obtain the necessary data, they require a lot of manpower and material resources, which is very uneconomical. At the same time, since the development of joints in rock mass is unknown, field tests can only obtain macroscopic data. Geotechnical engineering numerical simulation software is very effective for static and dynamic analysis of unknown joint development and other situations.

[0003] The stability of a tunnel is generally related to the rock strength and the integrity of the rock mass. In particular, under the long-term action of the original rock stress, geological interfaces of different scales, such as faults, joints and weak interlayers, develop inside the rock mass.

[0004] The deformation of jointed tunnel chambers during excavation often deviates significantly from the predictions of classical continuous medium theory. Because the continuous medium method assumes continuous materials and neglects the crucial influence of discontinuous interfaces within the rock mass, it is ill-suited for stability analysis of jointed rock tunnels. For discontinuous media such as engineering rock masses, more suitable numerical analysis methods include the discrete element method (DEM) and discontinuous deformation analysis. Discontinuous deformation analysis can reflect many behavioral characteristics of discrete media and effectively simulate the discontinuous deformation features of objects. In discontinuous deformation analysis, the block elements can both interact to form a unified system and move independently, thus simulating and solving the characteristics of discrete media.

[0005] The DDA program was developed by Chinese PhD Shi Genhua using the C language more than 30 years ago. Although its discontinuous deformation analysis method is very suitable for numerical simulation analysis of geotechnical engineering, the original DDA program has become inadequate as the research objects and problems have become increasingly complex, and it is urgent to improve it for actual engineering. Summary of the Invention

[0006] To address the aforementioned problems, this invention aims to provide a numerical simulation method for jointed tunnel stability analysis based on discontinuous deformation analysis. This method can replace the original DDA program, reflect more complex jointed block systems, and explore the influence of joint physical and mechanical parameters on tunnel surrounding rock deformation. It features high simulation accuracy and strong practicality.

[0007] To achieve the above objectives, the technical solution adopted by the present invention is as follows:

[0008] Numerical simulation method for stability analysis of jointed tunnels based on discontinuous deformation analysis, including

[0009] Step 1: Based on the original DDA input program dl, modify the original DDA in the C language environment to obtain different joint parameters from different probability density functions, and generate joint lines and excavation chamber boundaries;

[0010] Step 2: Import the joint lines generated by the dl program in Step 1 into the dc program to enable the system program to cut the block:

[0011] Step 3: Incorporate the overall equilibrium equations of the damped block system into the df program to solve for the block displacements;

[0012] Step three describes the process of incorporating the overall equilibrium equations of the damped block system into the df program.

[0013] S3.1 Based on the principle of stationary potential energy, the d'Alembert principle is used to consider dynamic equilibrium, and the effect of damping is added to the bulk system.

[0014] S3.2 The viscous damping force is represented by the time step and the displacement increment, and assembled into a damping matrix;

[0015] The calculation process of the damping matrix in step S3.2 includes:

[0016] S3.21 Let the potential energy change of the entire system be divided into:

[0017]

[0018] in: —The system's displacement matrix, velocity matrix, and acceleration matrix; —The system's inertial force array; —Mass surface density of the bulk system; —The system's damping force array; —The system's external load array; —The gravity array of the system; —The elastic strain energy of the system, where C is the viscous damping coefficient;

[0019] S3.22 Suppose the system has n independent displacement components, and...

[0020]

[0021] Therefore, the first variation of Π is zero.

[0022]

[0023] S3.23 Viscous resistance is represented as a matrix of time variation and displacement increment, as shown below:

[0024]

[0025] The potential energy of viscous damping force is expressed as

[0026]

[0027] in, , Indicates the time-step displacement increment; Indicates the viscous damping coefficient; , This represents the matrix consisting of the displacement increments at each time step; Δ is the time step size.

[0028] S3.24 Taking the second-order variational equation of the above equation, we obtain the damping matrix as follows:

[0029] ;

[0030] S3.3 Add the damping matrix from step S3.2 to the overall equilibrium equation and solve for it;

[0031] S3.4 The generated results are placed into the dgdt and dtat files;

[0032] Step 4: Draw and calculate the overall changes and individual deformations of the block system. The dg program provides a visual interface.

[0033] Preferably, the process of modifying the original DDA in step one includes:

[0034] In the modified program dl, the probability density function for the joint trace length is:

[0035]

[0036] in, denoted by , where w represents the mean of the negative exponential distribution;

[0037] By estimating the sample mean, a random variable following a negative exponential distribution is generated using a computer-generated random number sequence within the range (0,1).

[0038]

[0039] Where Z represents a uniform random sequence on (0, 1);

[0040] The dip and dip angle of rock joints in S1.2 follow a normal or log-normal distribution, with the following probability density function:

[0041]

[0042] in The variance represents the joint dip and dip angle; w represents the joint trace length;

[0043] Using two uniform random number sequences located on (0, 1) generated by computer , Then, the Box-Muller method is used to generate random variables that follow a normal distribution.

[0044]

[0045] Or for

[0046]

[0047] S1.3 Determine the joint spacing and length parameters and define the cavity based on the user input parameters;

[0048] S1.4 Put the results of the segmented joint lines into dcdt, dtat, and dlps.

[0049] Preferably, the process of importing the joint lines into the DC program in step two includes:

[0050] S2.1 Calculate all joint line intersections based on the joint lines obtained from the dl program in step one;

[0051] S2.2 Divide the complete joint line into independent joint segments using the intersection points of the joint lines to obtain the DC program for independent blocks;

[0052] S2.3 Input the divided joint line results into blck, dtat, and dcps for easy access and editing later.

[0053] The beneficial effects of this invention are: This invention discloses a numerical simulation method for jointed tunnel stability analysis based on discontinuous deformation analysis. Compared with the prior art, the improvement of this invention lies in:

[0054] This invention addresses the shortcomings of the original DDA program in numerical simulation of jointed rock masses. It designs a numerical simulation method for jointed tunnel stability analysis based on discontinuous deformation analysis (DDA), which mainly includes: (1) Improvement of the pre-processing program: taking into account the statistical characteristics of joint distribution in order to reflect the real situation; (2) Improvement of the tunnel excavation modeling method: modifying the "cutting" part in the original program so that it can cut out the blocks of the initial support of the tunnel according to the pre-set process, so as to consider the initial support effect; (3) Improvement of the calculation and solution program: introducing the concept of damping into the overall equilibrium equation of the block system to speed up the iterative solution; (4) Improvement of the post-processing program: using Matlab to further mine the program calculation data in order to show the situation of the entire system. It can be seen that the numerical simulation of jointed tunnel stability analysis using this method can replace the original DDA program, reflect the more complex jointed block system, explore the influence of the physical and mechanical parameters of joints on the deformation of the surrounding rock of the tunnel, and has the advantages of high simulation accuracy and strong practicality. Attached Figure Description

[0055] Figure 1 This is a schematic diagram of a partial improvement to the DDA program of the present invention.

[0056] Figure 2 The input parameter table for the dl program of this invention.

[0057] Figure 3 The following is an improved DDA program operation flowchart for this invention.

[0058] Figure 4 This is a geometric simulation diagram of the tunnel structure in Embodiment 2 of the present invention.

[0059] Figure 5 This is a diagram showing the deformation around the tunnel chamber in Embodiment 2 of the present invention.

[0060] Figure 6 This is a diagram showing the deformation of the surrounding rock in the horizontal and vertical directions in Embodiment 2 of the present invention.

[0061] Among them, Figure 4 In the figure, Figure (a) shows the tunnel structure and geometric parameters, and Figure (b) shows the tunnel structure and calculation model.

[0062] exist Figure 5In the figure, Figure (a) shows the deformation around the tunnel chamber when the calculation time step is 300, Figure (b) shows the deformation around the tunnel chamber when the calculation time step is 600, Figure (c) shows the deformation around the tunnel chamber when the calculation time step is 900, Figure (d) shows the deformation around the tunnel chamber when the calculation time step is 1200, and Figure (e) shows the deformation around the tunnel chamber when the calculation time step is 1500.

[0063] exist Figure 6 In the figure, Figure (a) shows the horizontal deformation of the surrounding rock, and Figure (b) shows the vertical deformation of the surrounding rock. Detailed Implementation

[0064] To enable those skilled in the art to better understand the technical solutions of the present invention, the technical solutions of the present invention will be further described below in conjunction with the accompanying drawings and embodiments.

[0065] Example 1: As Figure 1-6 As shown, when running the original DDA program for iterative solution, it is often difficult to converge due to the large number of blocks. If the program iteration process is observed, it can be found that some blocks collide and reflect in closed cavities after falling off, producing a situation similar to "oscillation", which leads to a long iterative calculation solution time.

[0066] This phenomenon is caused by the slow dissipation of system energy. In order to avoid encountering these long-term energy dissipation phenomena during iterative solutions;

[0067] Based on the above, this embodiment provides a numerical simulation method for jointed tunnel stability analysis based on discontinuous deformation analysis (DDA), the sequence of which is as follows: Figure 1 , Figure 3 As shown, the specific steps include:

[0068] Step 1: Based on the original DDA input program dl, modify it in the C language environment, such as... Figure 2 Users can smoothly input parameters to obtain different joint parameters from different probability density functions, and generate joint lines, excavation chamber boundaries, joint domains, etc.

[0069] The key to obtaining different joint parameters from different probability density functions lies in the fact that the geometric parameters of rock mass joints mostly conform to statistical distribution characteristics; the trace length of joints follows a negative exponential distribution; the dip and dip angle of joints follow a normal or log-normal distribution; the mean and variance of a certain distribution can be estimated through a computer-generated random sequence, thereby generating a more accurate joint line. Therefore, the specific process of step one above includes:

[0070] In the modified program dl, the probability density function of the joint trace length can be expressed as:

[0071]

[0072] In the formula, The mean of the negative exponential distribution can be considered as the mean of the joint trace lengths, where w represents the joint trace length.

[0073] By estimating the sample mean and using a computer-generated random number sequence within the range (0,1), a random variable following a negative exponential distribution can be generated.

[0074]

[0075] In the formula, Z represents a uniform random sequence on (0, 1);

[0076] The dip and dip angle of rock joints in S1.2 follow a normal or log-normal distribution, with the following probability density function:

[0077]

[0078] in The variance represents the joint dip and dip angle; w represents the joint trace length;

[0079] Using two uniform random number sequences located on (0, 1) generated by computer , Then, the Box-Muller method can be used to generate random variables that follow a normal distribution.

[0080]

[0081] Or for

[0082]

[0083] S1.3 Determine parameters such as joint spacing and length based on user input parameters and define the cavern;

[0084] S1.4 Place the results of the segmented joint lines into dcdt, dtat, and dlps for easy viewing;

[0085] Step 2: Import the joint lines generated by the dl program in Step 1 into the dc program to enable the system program to cut jointed surrounding rock and other blocks:

[0086] The key to the cutting jointed surrounding rock block lies in the definition of materials during the modeling process. In the original DDA program, user-defined material line parameters are used. When a block is passed through by a material line, it is assigned the parameters represented by that material line. In reality, a model cannot be accurately built with only a limited number of lines. This embodiment modifies the original DDA program so that multiple methods can be selected when defining material properties. Therefore, the specific process of step two above includes the following steps:

[0087] S2.1 Calculate all joint line intersections based on the joint lines obtained from the dl program in step one;

[0088] S2.2 Divide the complete joint line into independent joint segments using the intersection points of the joint lines to obtain the DC program for independent blocks;

[0089] S2.3 Input the divided joint line results into blck, dtat, and dcps for easy access and editing later;

[0090] Step 3: Incorporate the overall equilibrium equations of the damped block system into the df program to avoid the phenomenon of energy not dissipating easily for a long time during iteration, thereby accelerating the iteration process and realizing the solution of block displacement. This includes the following steps:

[0091] S3.1 Based on the principle of stationary potential energy, the d'Alembert principle is used to consider dynamic equilibrium, and the effect of damping is added to the bulk system.

[0092] S3.2 The viscous damping force is represented by the time step and the displacement increment, and assembled into a damping matrix;

[0093] S3.3 Add the damping matrix introduced in step S3.2 to the overall equilibrium equation and solve for it;

[0094] S3.4 The generated results are placed into the dgdt and dtat files for easy use later;

[0095] Step 4: Draw the overall changes and individual deformations of the block system after calculation. The dg program provides a drawing interface that can visually reflect the deformation and changes of the project, so as to view the deformation situation.

[0096] Preferably, the block division process described in step S2.2 can form irregular polygons, and calculations are performed using the coordinates of the polygon vertices. The core code is as follows:

[0097] void Area()

[0098] {

[0099] for (i=1; i<= n1; i++)

[0100] {

[0101] for (j=1; j<= 6; j++)

[0102] {

[0103] area[i][j] = 0;

[0104] }

[0105] / * x1 y1 to x2 y2 anti-clockwise around the block * /

[0106] for (j=k0[i][1]; j<=k0[i][2]; j++)

[0107] {

[0108] x1 = d[j][1];

[0109] y1 = d[j][2];

[0110] x2 = d[j+1][1];

[0111] y2 = d[j+1][2];

[0112] area0 = (x1*y2-x2*y1);

[0113] area[i][1] += area0 / 2;

[0114] area[i][2] += area0*(x1+x2) / 6;

[0115] area[i][3] += area0*(y1+y2) / 6;

[0116] area[i][4] += area0*(x1*x1+x2*x2+x1*x2) / 12;

[0117] area[i][5] += area0*(y1*y1+y2*y2+y1*y2) / 12;

[0118] area[i][6] += area0*(2*x1*y1+2*x2*y2+x1*y2+x2*y1) / 24;

[0119] }

[0120] }

[0121] }。

[0122] Preferably, the calculation process of the damping matrix in step S3.2 includes:

[0123] S3.21 Let the potential energy variation of the entire system be expressed as:

[0124]

[0125] in: —The system's displacement matrix, velocity matrix, and acceleration matrix; —The system's inertial force array; —Mass surface density of the bulk system; —The system's damping force array; —The system's external load array; —The gravity array of the system; —The elastic strain energy of the system, where C is the viscous damping coefficient, which is proportional to the mass, and the proportionality coefficient can be determined based on relevant empirical data;

[0126] S3.22 Suppose the system has n independent displacement components, that is,

[0127]

[0128] Therefore, the first variation of Π being zero can be expressed as:

[0129]

[0130] S3.23 Viscous resistance is represented as a matrix of time variation and displacement increment, as shown below:

[0131]

[0132] The potential energy of viscous damping force can be expressed as

[0133]

[0134] In the formula, , Indicates the time-step displacement increment; Indicates the viscous damping coefficient; , This represents the matrix consisting of the displacement increments at each time step; Δ is the time step size.

[0135] S3.24 To achieve system equilibrium and minimize the potential energy of the viscous damping force, a second-order variational calculus of the above equation yields the damping matrix as follows:

[0136] .

[0137] Preferably, the core code for adding the damping component in step three includes:

[0138] void damp()

[0139] {

[0140] for (ii=1; ii<= n1; ii++)

[0141] {

[0142] / * initialize matrix for damp matrix in the beginning * /

[0143] for (jj=1; jj<=6; jj++)

[0144] {

[0145] for (ll=1; ll<=6; ll++)

[0146] {

[0147] Ci[jj][ll]=0;

[0148] }

[0149] }

[0150] / *correct non-zero elements of the matrix * /

[0151] x0=area[ii][2] / area[ii][1];

[0152] y0=area[ii][3] / area[ii][1];

[0153] Ci[1][1]=area[ii][1];

[0154] Ci[2][2]=area[ii][1];

[0155] Ci[3][3]=area[ii][4]-x0*area[ii][2]+area[ii][5]-y0*area[ii][3];

[0156] Ci[3][4]=-area[ii][6]+x0*area[ii][3];

[0157] Ci[4][3]=Ci[3][4];

[0158] Ci[3][5]=area[ii][6]-x0*area[ii][3];

[0159] Ci[5][3]=area[ii][6]-x0*area[ii][3];

[0160] Ci[3][6]=0.5*(area[ii][4]-x0*area[ii][2]-area[ii][5]+y0*area[ii][3]);

[0161] Ci[6][3]=Ci[3][6];

[0162] Ci[4][4]=area[ii][4]-x0*area[ii][2];

[0163] Ci[4][6]=0.5*(area[ii][6]-x0*area[ii][3]);

[0164] Ci[6][4]=Ci[4][6];

[0165] Ci[5][5]=area[ii][5]-y0*area[ii][3];

[0166] Ci[5][6]=0.5*(area[ii][6]-x0*area[ii][3]);

[0167] Ci[6][5]=Ci[5][6];

[0168] Ci[6][6]=0.25*(area[ii][4]-x0*area[ii][2]+area[ii][5]-y0*area[ii][3]);

[0169] / * add damp matrix to a[][] * /

[0170] i1=k2[i];

[0171] if (pp==0) i2=n[i1][1]+n[i1][2]-1;

[0172] if (pp==1) i2=n[i1][1];

[0173] for (jj=1; jj<=6; jj++)

[0174] {

[0175] for (ll=1; ll<=6; ll++)

[0176] {

[0177] ji = 6*(jj-1) + ll;

[0178] a[i2][ji] += 2*0.02*o0 / tt*gg*Ci[jj][ll];

[0179] }

[0180] }

[0181] }

[0182] }

[0183] Example 2: Unlike Example 1 above, the following example is used to verify the feasibility and superiority of the numerical simulation method for jointed tunnel stability analysis based on discontinuous deformation analysis described in Example 1:

[0184] The southern route of the Central Asia-China Gas Pipeline adopts the Uzbekistan-Tajikistan-China (WTJZ) scheme: starting from Turkmenistan, passing through Uzbekistan, Tajikistan, and Kyrgyzstan, and reaching Xinjiang, China, with a total length of about 1,000 km. The key engineering projects include 41 tunnels that cross 71 km, with the longest tunnel being 3,150 m. The tunnels are distributed within a range of about 260 km in Tajikistan.

[0185] The No. 1 tunnel of the South Route of the Central Asia-China Gas Pipeline is the first mountain crossing project for the pipeline as it enters Tajikistan from west to east. The tunnel is located about 3 km southwest of Dushanbe. The tunnel entrance is relatively convenient to access, while the Tajik A384 national highway is about 3 km outside the exit. New roads need to be built to the tunnel entrance for both the entrance and exit. The tunnel has a horizontal length of 1860 m and a designed longitudinal slope of 10.48%.

[0186] Tunnel entrance boundary mileage coordinates: (X=4255568.027, Y=474073.396, H=714.59);

[0187] Design coordinates of the tunnel entrance: (X=4255564.4429, Y=474076.8837);

[0188] Design coordinates of the tunnel exit: (X=4254235.9448, Y=475378.6876);

[0189] The tunnel exit boundary mileage coordinates are: (X=4254310.042, Y=475379.742, H=910.17).

[0190] Bedrock is exposed on the surface south of the tunnel entrance, and the lithology is calcareous mudstone. According to the engineering geological survey results, joints are well-developed in the area through which the tunnel passes, with two main sets of joints:

[0191] ① Joint group 1, dip direction 215°~225°, dip angle 73°~83°;

[0192] ② Joint group 2, dips 305°~315°, dip angle 53°~63°; at the tunnel exit end 1K1+690~1K1+700, the surrounding rock is calcareous mudstone (containing silt) and there is a small amount of water seepage. The initial support concrete has local cracking and spalling, and the local longitudinal connecting steel bars have deformed.

[0193] Longitudinal cracks (approximately 20 meters in length and 2-3 mm in width) appeared in the arch wall near 1K1+663 to 1K1+683. Monitoring and measurement data showed that the settlement of the arch crown and the horizontal convergence rate of the arch waist were relatively large.

[0194] The numerical model was established based on the tunnel's structural characteristics and geological conditions. The tunnel cross-section is from K1+625 to K1+688. The tunnel's structural form and geometric parameters are shown in [reference needed]. Figure 4 (a); Computational model such as Figure 4 As shown in (b), the calculation area is 50m×50m. The model boundary of this area has double-sided constraint displacement in the x direction and upper and lower constraint displacement in the z direction. In the tunnel section with a large burial depth (about 100m), the depth exceeding the calculation area is converted into the initial stress.

[0195] The mechanical parameters of the rock mass used in the calculation are shown in Table 1. The initial support structure has a thickness of 0.2m and its mechanical parameters are the same as those of the rock mass.

[0196] Table 1: Rock mass mechanical parameters

[0197]

[0198] The computational model uses two sets of joints: the rock strata dip at 37° and the dip angle at 39.5°, and the tunnel cross-section dip at 136° and the dip angle at 90°. The physical parameters of the joints are shown in Table 2. To reflect the uncertainty of joint distribution in the rock mass, the joints in the program are generated statistically using joint parameters with a certain degree of randomness. The dip and dip angle of the joints are both normally or log-normally distributed, while the spacing and length of the joints are distributed negatively exponentially. The simulation is performed using computer-generated random numbers.

[0199] Table 2: Rock Mass Mechanical Parameters

[0200]

[0201] The tunnel was simulated using the Discontinuous Deformation Analysis (DDA) method, and the deformation around the tunnel chamber was as follows: Figure 5 As shown;

[0202] When the calculation time step is around 300 steps, the displacement of the block system has not yet developed, and the overall deformation of the tunnel is relatively small. When the calculation time step increases to 600 to 900 steps, the blocks around the tunnel chamber show rapid deformation development, with the horizontal displacement of the right sidewall developing rapidly, and some blocks of the floor slab and arch beginning to show vertical displacement. As the calculation time step continues to increase, the deformation of the tunnel chamber undergoes a rapid development process, and the displacement development tends to be stable. During this period, the horizontal displacement of the right sidewall is always the largest, and the vertical displacement of the arch is also relatively large. During the displacement development of the tunnel surrounding rock, the entire right sidewall experienced large-scale block instability, while the displacement of the left sidewall was very small, and the surrounding rock blocks did not move significantly. This is because the large displacement of the right sidewall and arch releases the stress around the chamber, and the left sidewall can maintain stability under the support.

[0203] Figure 6 The deformation of the surrounding rock in the horizontal and vertical directions is given, from Figure 6 It can be seen that the maximum displacement occurred at the junction of the side wall and the arch, where the block instability was caused by toppling failure; at the same time, the blocks at the bottom plate and the top of the arch showed large vertical displacement, and the overall displacement was large.

[0204] Figure 6 (a) shows the horizontal deformation of the surrounding rock, from Figure 6 As shown in (a), under the influence of joint surfaces, the tunnel sidewalls exhibit asymmetrical convergent deformation. The deformation of the surrounding rock at the left sidewall is very small, while the horizontal deformation of the surrounding rock at the right sidewall is larger. The reason for this phenomenon is that the large displacement of the right sidewall and the arch releases the stress around the tunnel, while the left sidewall remains stable under the support. The deformation locations and instability modes are consistent with the measured results from K1+663 to the tunnel face.

[0205] Figure 6 (b) gives the vertical deformation of the surrounding rock, from Figure 6 As shown in (b), the deformation of the right sidewall and the tunnel crown is relatively large under the influence of joint surfaces. The vertical deformation of the surrounding rock is the largest at the right sidewall, while the vertical displacement of the surrounding rock at the crown is relatively small. The failure modes of the right sidewall and the crown under the influence of joints are different. The "critical block" at the right sidewall moves first, and the failure mode is block toppling and instability with a large displacement. After the block at the sidewall moves, some blocks at the crown are suspended and then move, and the failure mode is block falling with a slightly smaller displacement. These deformed locations and instability modes are consistent with the measured results from K1+663 to the tunnel face.

[0206] The foregoing has shown and described the basic principles, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The embodiments and descriptions in the specification are merely illustrative of the principles of the invention. Various changes and modifications can be made to the invention without departing from its spirit and scope, and all such changes and modifications fall within the scope of the present invention as claimed. The scope of protection of this invention is defined by the appended claims and their equivalents.

Claims

1. A numerical simulation method for stability analysis of jointed tunnels based on discontinuous deformation analysis, characterized in that: include Step 1: Based on the original DDA input program dl, modify the original DDA in the C language environment to obtain different joint parameters from different probability density functions, and generate joint lines and excavation chamber boundaries; Step 2: Import the joint lines generated by the dl program in Step 1 into the dc program to enable the system program to cut the block: Step 3: Incorporate the overall equilibrium equations of the damped block system into the df program to solve for the block displacements; Step three describes the process of incorporating the overall equilibrium equations of the damped block system into the df program. S3.1 Based on the principle of stationary potential energy, the d'Alembert principle is used to consider dynamic equilibrium, and the effect of damping is added to the bulk system. S3.2 The viscous damping force is represented by the time step and the displacement increment, and assembled into a damping matrix; The calculation process of the damping matrix in step S3.2 includes: S3.21 Let the potential energy change of the entire system be divided into: ； in: —The system's displacement matrix, velocity matrix, and acceleration matrix; —The system's inertial force array; —Mass surface density of the bulk system; —The system's damping force array; —The system's external load array; —The gravity array of the system; —The elastic strain energy of the system, where C is the viscous damping coefficient; S3.22 Suppose the system has n independent displacement components, and... ； Therefore, the first variation of Π is zero. ； S3.23 Viscous resistance is represented as a matrix of time variation and displacement increment, as shown below: ； The potential energy of viscous damping force is expressed as ； in, , Indicates the time-step displacement increment; Indicates the viscous damping coefficient; , This represents the matrix consisting of the displacement increments at each time step; Δ is the time step size. S3.24 Taking the second-order variational equation of the above equation, we obtain the damping matrix as follows: ； S3.3 Add the damping matrix from step S3.2 to the overall equilibrium equation and solve for it; S3.4 The generated results are placed into the dgdt and dtat files; Step 4: Draw and calculate the overall changes and individual deformations of the block system. The dg program provides a visual interface.

2. The numerical simulation method for jointed tunnel stability analysis based on discontinuous deformation analysis according to claim 1, characterized in that: The process of modifying the original DDA as described in step one includes: In the modified program dl, the probability density function for the joint trace length is: ； in, denoted by , where w represents the mean of the negative exponential distribution; By estimating the sample mean, a random variable following a negative exponential distribution is generated using a computer-generated random number sequence within the range (0,1). ； Where Z represents a uniform random sequence on (0, 1); The dip and dip angle of rock joints in S1.2 follow a normal or log-normal distribution, with the following probability density function: ； in The variance represents the joint dip and dip angle; w represents the joint trace length; Using two uniform random number sequences located on (0, 1) generated by computer , Then, the Box-Muller method is used to generate random variables that follow a normal distribution. ； Or for ； S1.3 Determine the joint spacing and length parameters and define the cavity based on the user input parameters; S1.4 Put the results of the segmented joint lines into dcdt, dtat, and dlps.

3. The numerical simulation method for jointed tunnel stability analysis based on discontinuous deformation analysis according to claim 1, characterized in that: Step two describes the process of importing joint lines into the DC program, which includes... S2.1 Calculate all joint line intersections based on the joint lines obtained from the dl program in step one; S2.2 Divide the complete joint line into independent joint segments using the intersection points of the joint lines to obtain the DC program for independent blocks; S2.3 Input the divided joint line results into blck, dtat, and dcps for easy access and editing later.

Images