A method for calculating stability of rock mass considering fatigue and time effect

Abstract

The application provides a rock mass stability calculation method considering fatigue and time effect, comprising the following steps: dividing rock into a plurality of base elements to obtain mechanical parameters of the base elements; obtaining Hoek-Brown elastic-plastic damage model of the base elements; obtaining time-dependent damage variable of the base elements; obtaining fatigue damage variable of the base elements; obtaining coupling damage variable; obtaining t m+1 nominal stress tensor σ' of a time step m+1 ; obtaining coupling damage variable value of the t m+1 time step; and analyzing stability of the rock. The application establishes Hoek-Brown elastic-plastic damage model according to the nonlinearity and microscopic non-uniformity of rock mass strength, obtains coupling damage variable of the base elements based on the time effect characteristics of the rock mass and in combination with fatigue damage variable and load damage variable, finally obtains the value of the coupling damage variable changing with time, and realizes stability calculation of the non-homogeneous rock mass with time under the rheological effect and fatigue effect.

Description

Technical Field

[0001] This invention relates to the field of rock mass engineering calculation, and in particular to a method for calculating rock mass stability that takes into account fatigue and time effects. Background Technology

[0002] In recent years, rock mass engineering construction in many fields such as transportation, water conservancy, and mining in my country has flourished, including projects such as slopes, hydropower plants, tunnels, and foundations. Rock masses are geological structures composed of intact rock and structural planes, containing numerous faults, joints, and weak interlayers. Due to the complexity of the environment in which geological masses exist, compared to general media, rock masses exhibit heterogeneity and time effects, and are subject to environmental fatigue. The stability of rock mass engineering is becoming increasingly prominent. Therefore, researching rock mass stability calculation methods that consider rheological and fatigue effects is of great significance and provides important guidance for ensuring the safe construction and operation of rock mass engineering projects in complex environments.

[0003] Currently, stability calculations for rock mass engineering generally employ ideal elastoplastic macroscopic models. Traditional rock mass stability calculation methods, using constitutive models with linear strength criteria, fail to reflect the nonlinear characteristics of rock mass strength. They treat each stratum as homogeneous material and assign calculation parameters, failing to reflect the heterogeneity of the rock mass medium through macroscopic calculations. Furthermore, rheological models using component combinations cannot reflect the essence of rock mass deterioration under the combined effects of load, aging, and fatigue, nor the corresponding stability changes. Therefore, there are significant errors in considering the stability of rock masses that comprehensively consider both rheological and fatigue effects. Summary of the Invention This invention provides a method for calculating rock mass stability that takes into account fatigue and time effects, in order to overcome the above-mentioned technical problems.

[0004] To achieve the above objectives, the technical solution of the present invention is as follows: A method for calculating rock mass stability considering fatigue and time effects includes the following steps: S1: Divide the rock into several basic units. V In order to obtain the mechanical parameters of the basic element; The mechanical parameters of the basic elements include empirical parameters that reflect the characteristics of the rock mass. m b Parameters of rock mass fragmentation and integrity s Rock mass quality parameters a and the elastic modulus before damage E 0; S2: Based on the mechanical parameters of the basic element, obtain the Hoek-Brown elastoplastic damage model of the basic element; S3: Obtain the time-dependent damage variable of the primitive. ; S4: Obtain the fatigue damage variable of the primitive.D r ; S5: Based on load damage variables D m The aforementioned time-related damage variables and the fatigue damage variable D r Obtain the coupling damage variable Dc: S6: Based on the Hoek-Brown elastoplastic damage model of the basic element and the coupled damage variables, obtain... t m+1 nominal stress tensor at time step ; S7: According to the... t m+1 nominal stress tensor at time step , obtain the t m+1 The value of the coupling damage variable at the time step.

[0005] Beneficial Effects: This invention provides a method for calculating rock mass stability considering fatigue and time effects. Addressing the nonlinearity and microscopic heterogeneity of rock mass strength, the method divides the rock into several basic elements, obtaining empirical parameters reflecting rock mass characteristics, parameters of rock mass fragmentation and integrity, rock mass mass quality parameters, and the elastic modulus before damage. Based on these parameters, a Hoek-Brown elastoplastic damage model is established. Then, based on the time-effect characteristics of the rock mass, the time-dependent damage variables of the basic elements are obtained. Combined with fatigue damage variables and load damage variables, the coupled damage variables of the basic elements are obtained, ultimately yielding the values ​​of the coupled damage variables that change over time. This enables the calculation of the stability of heterogeneous rock masses under rheological and fatigue effects over time. Attached Figure Description

[0006] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0007] Figure 1 This is a flowchart of the rock mass stability calculation method considering fatigue and time effects of the present invention. Figure 2a This is a probability density distribution diagram of the RVE mechanical properties when x0 is 10 in an embodiment of the present invention; Figure 2b This is a probability density distribution diagram of the RVE mechanical properties when ζ is 0.5 in an embodiment of the present invention; Figure 3This is a schematic diagram of the Poytin-Thomson model in an embodiment of the present invention; Figure 4 This is a graph showing the variation of rock mass damage variables with fatigue cycles in an embodiment of the present invention; Figure 5 This is a calculation model diagram of a rock sample under mesoscopic distribution in an embodiment of the present invention; Figure 6a In an embodiment of the present invention, m is when ζ is 0.5 b Value distribution cloud map; Figure 6b In an embodiment of the present invention, m is when ζ is 2. b Value distribution cloud map; Figure 6c In an embodiment of the present invention, m is when ζ is 4. b Value distribution cloud map; Figure 6d In an embodiment of the present invention, m is when ζ is 6. b Value distribution cloud map; Figure 6e In an embodiment of the present invention, m is when ζ is 8. b Value distribution cloud map; Figure 6f In an embodiment of the present invention, m is when ζ is 10. b Value distribution cloud map; Figure 7a shows the m when ζ is 0.5 in an embodiment of the present invention. b Schematic diagram of RVE statistics for value regions; Figure 7b shows m when ζ is 2 in an embodiment of the present invention. b Schematic diagram of RVE statistics for the value region; Figure 7c shows m when ζ is 4 in an embodiment of the present invention. b Schematic diagram of RVE statistics for value regions; Figure 7d shows the m when ζ is 6 in an embodiment of the present invention. b Schematic diagram of RVE statistics for value regions; Figure 7e In an embodiment of the present invention, m is when ζ is 8. b Schematic diagram of RVE statistics for value regions; Figure 7f In an embodiment of the present invention, m is when ζ is 10. b Schematic diagram of RVE statistics for the value region; Figure 8a In an embodiment of the present invention, m is when ζ is 0.5 b And a diagram showing the evolution of damage values; Figure 8b In an embodiment of the present invention, m is when ζ is 4. b And a diagram showing the evolution of damage values; Figure 9 The creep curves of rocks under different stress levels are shown in the embodiments of the present invention. Figure 10a This is a curve showing the change of the rock sample damage zone over time when ζ is 6 in an embodiment of the present invention; Figure 10b This is a curve showing the change of the rock sample damage zone over time when ζ is 6 in an embodiment of the present invention; Figure 11 This is a finite element model diagram for slope calculation in an embodiment of the present invention; Figure 12a This is a diagram showing the variation of the slope's elastic modulus and the damaged area when ζ is 2 in an embodiment of the present invention. Figure 12b This is a diagram showing the variation of the slope's elastic modulus and the damaged area when ζ is 6 in an embodiment of the present invention. Figure 12c This is a diagram showing the variation of the slope's elastic modulus and the damaged area when ζ is 10 in an embodiment of the present invention. Figure 13a This is a schematic diagram of the slope damage zone when the month is 0 (month) in an embodiment of the present invention; Figure 13b This is a schematic diagram of the slope damage zone when the month is 4 (month) in an embodiment of the present invention; Figure 13c This is a schematic diagram of the slope damage zone in an embodiment of the present invention when the month is August. Figure 13d This is a schematic diagram of the slope damage zone when the month is December in an embodiment of the present invention; Figure 13e This is a schematic diagram of the slope damage zone when the month is 16 (month) in an embodiment of the present invention; Figure 13f This is a schematic diagram of the slope damage zone in an embodiment of the present invention when the month is 20 (month). Figure 14 This is a graph showing the variation of the slope safety factor over time in an embodiment of the present invention; Figure 15 This is a graph showing the variation of maximum coupled damage of the slope over time (t / month) under different freeze-thaw cycles in an embodiment of the present invention. Figure 16a This is a schematic diagram of slope displacement after 80 freeze-thaw cycles with a time of 0 (months) in an embodiment of the present invention; Figure 16b This is a schematic diagram of slope displacement after 80 freeze-thaw cycles over a period of 4 months, as shown in an embodiment of the present invention. Figure 16c This is a schematic diagram of slope displacement after 80 freeze-thaw cycles over a period of 8 months, as shown in an embodiment of the present invention. Figure 16dThis is a schematic diagram of slope displacement after 80 freeze-thaw cycles over a period of 12 months, as shown in an embodiment of the present invention. Figure 16e This is a schematic diagram of slope displacement after 80 freeze-thaw cycles over a period of 16 months, as shown in an embodiment of the present invention. Figure 16f This is a schematic diagram of slope displacement after 80 freeze-thaw cycles over a period of 20 months, as shown in an embodiment of the present invention. Figure 17 This is a graph showing the change in displacement of slope monitoring points over time in an embodiment of the present invention; Figure 18a This is a schematic diagram of the slope damage zone when the anchor cable prestress is 200KN in an embodiment of the present invention; Figure 18b This is a schematic diagram of the slope damage zone when the anchor cable prestress is 300KN in an embodiment of the present invention; Figure 18c This is a schematic diagram of the slope damage zone when the anchor cable prestress is 400KN in an embodiment of the present invention; Figure 18d This is a schematic diagram of the slope damage zone when the anchor cable prestress is 500KN in an embodiment of the present invention. Detailed Implementation

[0008] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0009] This invention provides a method for calculating rock mass stability considering fatigue and time effects, such as... Figure 1 As shown; including the following steps: S1: To study microscopic discontinuities caused by the effects of voids and crystals from the perspective of continuous medium damage mechanics, it is necessary to homogenize the microscopic mechanical properties of the rock material. First, the rock is divided into several basic units. V Secondly, statistical methods are used to analyze the primitives. V The physical and mechanical properties are described to obtain the mechanical parameters of the basic element; The mechanical parameters of the basic elements include, but are not limited to, empirical parameters reflecting the characteristics of the rock mass. m b Parameters of rock mass fragmentation and integrity s Rock mass quality parameters a and the elastic modulus before damage E 0, etc.

[0010] Specifically, a power law with a threshold value is used to describe the extreme value distribution law of rock mass material strength. Assume the mechanical parameters in the primitive... P z It can be described using the Weibull statistical distribution function. The method for obtaining the mechanical parameters of the basic element is as follows: (1) In the formula, P 0 represents the macroscopic properties of the rock sample; x z Use Weibull random numbers; P z These are the mechanical parameters of the basic unit; Weibull random numbers x z The method to obtain it is as follows: First, based on the cumulative distribution function and probability density function of the Weibull distribution, as follows: (2) (3) In the formula: W ( x ) is the cumulative distribution function (CDF) of the Weibull distribution; w ( x ) is the probability density function (PDF) of the Weibull distribution. >0 is the scaling parameter; ζ For shape parameters; x It is a random variable; In equation (3), as ζ The increase in the rock (4) In the formula: Use Weibull random numbers. These are uniformly distributed random numbers; S2: Based on the mechanical parameters of the basic element, obtain the Hoek-Brown elastoplastic damage model of the basic element to reflect the nonlinear strength and plastic softening characteristics of the rock mass. Specifically, load-damage variables are typically introduced based on the existing rock mechanics yield function. D m This is used to characterize the weakening effect of damage on strength. This embodiment assumes that damage affects the parameters during loading.m b and s To address the weakening effect, a Hoek-Brown elastoplastic damage model considering the basic elements with plastic damage is proposed, and its expression is as follows: (5) (6) In the formula, f Let be the yield function; J 2 represents the second deviatoric stress invariant; θ σ For Rhodes angle; D m For load damage variables; σ ci The uniaxial compressive strength of intact rock; p It is the hydrostatic pressure; m bg for m b The corresponding plastic potential function parameters; s g for s The corresponding plastic potential function parameters; a g for a The corresponding plastic potential function parameters; g Let be the plastic potential function; when m bg = m b , s g = s , a g = a When the flow is related, it is a flow rule; otherwise, it is a flow rule unrelated.

[0011] Among them, load damage variables D m The following is obtained: Due to the stress redistribution effect caused by engineering excavation, plastic damage will occur inside the rock mass. The damage increases with the equivalent plastic strain. Accumulated and exacerbated, load damage variables D m It can be represented as The power form: (7) In the formula, α To determine the initial slope of the softening curve of the rock material after damage, the value range is (0, +∞); β The parameter used to determine the maximum damage value of the rock has a value range of [0,1). Equivalent plastic strain; (8) In the formula, ε p1 Let x be the principal plastic strain along the x-axis of the spatial coordinate system (i.e., principal stress space); ε p2 Let y be the principal plastic strain along the y-axis in the spatial coordinate system (i.e., principal stress space); ε p3 Let Z be the principal plastic strain along the z-axis of the spatial coordinate system (i.e., principal stress space); S3: Obtain Time-Related Damage Variables The method is as follows; When the axial compression is less than the long-term strength value, the mechanical model of the element can be used. Figure 3 The Poytin-Thomson model is shown, and the creep equation of the primitives is established as follows: (9) In the formula, In response to the situation; For axial stress; K 1. K 2 are the mechanical parameters of the elastic element; η The viscosity coefficient of the sticky pot; For the rheological time; In this embodiment, the Poytin-Thomson model is a prior art technique that includes two elastic elements, therefore K 1. K 2 represents the mechanical parameters of the two elastic elements in the model.

[0012] Simplifying equation (9) yields (10) In the formula, P, Q, and R are all intermediate parameters, where, , , ; From formula (10), we can obtain that the elastic modulus changes with time. t The expression that changes is: (11) E ( t () represents the elastic modulus as a function of time. t A changing function; According to the theory of continuous damage mechanics, the degree of damage to a material modulus can be determined by the elastic modulus method, that is: (12) In the formula, These are time-dependent damage variables, i.e., damage variables caused by time effects; E 0 represents the elastic modulus before damage; S4: Obtain fatigue damage variables D r ,as follows: Assume the first The damage value of the cross section under secondary fatigue disturbance is D t It can be measured that: (13) In the formula: D t For the first Damage values ​​of rock cross-sections under secondary fatigue disturbance; D r For fatigue damage variables D r That is, the first N fatigue damage variables prior to the first fatigue disturbance; This is the number of fatigue disturbances; N For the first Total number of fatigue disturbances before the current fatigue disturbance; W and U All are material parameters; After the first Damage value of rock cross-section under secondary fatigue disturbance D t It can be obtained from the rock porosity and longitudinal wave velocity, and the calculation formula is shown in Equation (14): (14) In the formula, n 0 represents the porosity of the rock before fatigue. n Porosity of the rock after fatigue. The longitudinal wave velocity after rock fatigue is given; the relationship between the number of fatigue cycles and fatigue damage is fitted using equation (13), and the results are as follows: Figure 4 As shown.

[0013] Step 5: Based on load damage variables D m The aforementioned time-related damage variables and the fatigue damage variable D r Obtain the coupling damage variable Dc: For engineering projects in cold northern regions during service, rock masses undergo aging damage and fatigue damage under long-term loading, as well as freeze-thaw fatigue damage. This paper couples these various damage types to analyze the rock mass stability under the influence of multiple factors. Fatigue damage variables are comprehensively considered. Dr Load damage variables D m and time-related damage variables D (t), then the coupling damage variable is obtained as follows: (15) In the formula, Dc is the coupling damage variable; D m For load damage variables; S6: Based on the Hoek-Brown elastoplastic damage model of the basic element and the coupled damage variables, obtain... t m+1 nominal stress tensor at time step ; S61: The stress-strain constitutive relation considering time-fatigue-load coupled damage is as follows: (16) In the formula, Stress to account for time-fatigue-load coupled damage; For elastic strain; Stiffness matrix considering coupled damage; , It is a fourth-order symmetric tensor. For Kronecker product; G ( D c () represents the damage shear modulus; K ( D c ) represents the bulk modulus; where, , i , j , k , l All are subscripts of fourth-order symmetric tensors; , , , All are components of a fourth-order symmetric tensor; among them, G ( D c ), K ( D c The initial shear modulus can be used. G 0 and initial bulk modulus K 0 represents: (17) In the formula, G 0 represents the initial shear modulus; K 0 represents the initial bulk modulus; Poisson's ratio, For the damage variable is D C The elastic modulus at time; Equation (17) reflects the weakening effect of time-coupled damage on the elastic modulus of the rock mass.

[0014] Damage not only affects the elastic modulus but also the mechanical parameters of the primitives in the Hoek-Brown elastoplastic damage model. The variation of Hoek-Brown elastoplasticity over time is obtained from the following formula: (18a) (18b) In the formula: for t Moment value; For the initial value; for t Moment s value; For the initial s value; The Hoek-Brown elastoplastic damage model of formula (5) uses a fully implicit stress reflection algorithm. The entire solution process is divided into three steps: elastic prediction, plastic correction, and damage correction. Normally, when rock begins to deform, it is in the elastic stage. When the deformation reaches a certain level, it enters the plastic deformation stage. Therefore, we first assume that the rock is in the elastic deformation stage and solve it according to the following steps: S62: Obtain the... t m+1 Predicted stress at time steps to determine whether the rock has entered the plastic correction stage; Assumption Then, given the strain increment The expression for the predicted stress is obtained as follows: (19) In the formula, For the first t m+1 Predicted stress at time step; Stiffness matrix considering coupled damage; For the first t m Response under time steps; For the first t m The given strain increment at the time step; For the first t m+1 Predicting coupled damage variables at time steps; For the first t m Coupling damage variables at time step; The symbol for double dot product is used in tensor operations to represent the multiplication of components in a specific order. S63: The first t m+1 Substituting the predicted stress at the time step into the yield function of formula (5) f In the process of judgment, if: (20) Then the first t m+1 Stress at time step equal to the t m+1 Predicted stress at time step, i.e. If the condition is not met, proceed to step S62; otherwise, proceed to step S64 and enter the plasticity correction stage.

[0015] S64: Based on the plastic potential function, the Hoek-Brown elastoplastic damage model is plastically modified to obtain the modified equivalent plastic strain. To avoid t n+1 If the calculated stress at the end of the time step deviates from the yield surface, this embodiment employs a fully implicit stress regression algorithm to correct the elastically predicted stress back to the yield surface. During the plastic correction process, the first... t m Strain increment at time step With damage variables remaining constant, in the plastic state, the strain increment includes the elastic strain increment. and plastic strain increment The expression for the plastic strain increment is: (twenty one) In the formula: This represents the increment of plastic strain. To correspond to the plasticity factor of the bth plastic surface, since formulas (5) and (6) are hexagonal pyramids with 6 surfaces in the principal stress space, in formula (21), q The number of stresses exceeding the surface; The plasticity function corresponding to the bth plastic surface is formula (6); when the return stress is located at the yield surface, edge, or cusp, one, two, or three plasticity factors need to be solved respectively. Because the Hoek-Brown elastoplastic damage model has numerical singularities, it cannot be directly differentiated. Furthermore, using a smoothing method would reduce the solution accuracy. Therefore, this embodiment is based on... The location of the stress in space is discussed in different cases to determine whether the final updated stress is located at a cusp, yield surface, or edge, thus avoiding the difficulties caused by global differentiation. This step is a well-known finite element algorithm technique in the field and will not be discussed in detail here.

[0016] S65: Obtain from the corrected equivalent plastic strain t m+1 nominal stress tensor at time step The equivalent plastic strain obtained from the plastic correction step is calculated according to equations (7) and (15). t m+1 Coupled damage variables at time steps And use equation (22) to correct the stress obtained in the previous step again.

[0017] therefore, t m+1 nominal stress tensor at time step The following information was obtained: (twenty two) In the formula: and They are the first t m and t m+1 Damage elasticity matrix at time step; For the first t m+1 The plastic strain increment at the time step; It is the first t m+1 Coupling damage variables at time step; For the first t m Coupling damage variables at time step; S7: According to the... t m+1 nominal stress tensor at time step , obtain the t m+1 The value of the coupling damage variable at the time step.

[0018] Specifically, S6 is executed iteratively to solve the problem, and the result is the output of the first iteration. t m+1 The values ​​of coupled damage variables at each time step. To ensure the quadratic convergence rate during the overall iterative solution process, the consistent tangent modulus of the model should also be provided after the stress solution is completed. This embodiment uses Fortran to write the above numerical calculation process, realizing the finite element solution of the Hoek-Brown elastoplastic damage model of rock. The damage values ​​of each element are then output. This is prior art and will not be elaborated here.

[0019] S8: According to the... t m+1 The stability of the rock is analyzed by examining the values ​​of the coupled damage variable at each time step. Specifically, based on the values ​​of the coupled damage variable at each time step, the trend of the coupled damage variable over time can be obtained. When the values ​​of the coupled damage variable at two adjacent time steps are less than an empirical threshold, it is considered that the trend of the coupled damage variable over time tends to be stable, and the rock stability is relatively good. This is a well-known technique in the field and will not be elaborated on here.

[0020] This section uses a two-dimensional uniaxial compression numerical test of rock as an example to introduce the numerical solution process of the Hoek-Brown elastoplastic damage model, which considers the time effect of the rock HB microstructure elastoplastic damage model. First, a two-dimensional planar calculation model is established according to the standard specimen dimensions, as follows: Figure 5 As shown. The model is 0.01m high and 0.005m wide, with fixed boundary conditions applied to the bottom and displacement boundary conditions applied to the top. The model was meshed using quadrilateral Solid elements, resulting in a total of 1250 elements.

[0021] by m b Taking the average value as an example, we will explain the parameter assignment process for each primitive. m b0 With the shape parameters remaining unchanged, ζ The change m b The distribution pattern within the studied space is completely different, exhibiting microscopic disorder. This disorder reflects the unique discreteness of the rock. The Monte Carlo method is used to analyze the intensity parameters that exhibit disorder in physical space. m b Perform the assignment. First, generate values ​​within the range (0,1). c indivual( c (The total number of primitives) are random numbers that follow a uniform distribution. u L Then, the inverse method is used to map this set of random numbers to random numbers that follow the Weibull distribution, and the mapping function is shown in equation (23): (twenty three) Random attribute values ​​are assigned to the materials (primaries) of each partitioning unit. First, the method of this application is used to generate the properties by writing a program. c A set of random numbers following a Weibull distribution are generated, and then a loop statement in Python is used to generate... c A material array, finally... c The primitives are iterated and assigned values.

[0022] The calculation object is gneiss. In addition to the given calculation parameters, the time-effect calculation parameters are obtained using the parameter identification method proposed by Xu Hongfa. P =0.7524e-3, Q =-0.5416e-3, R =0.028. Table 1 shows the calculation parameters of the HB microscopic elastoplastic damage model considering the time effect.

[0023] Table 1. Detailed calculation parameters of the HB-D model

[0024] Model microparameter influence analysis: Let m b0 If the value is fixed at 6.36, then the difference is... ζ Value m b The distribution pattern is as follows Figures 6a-6f As shown. It can be seen that, ζ The smaller the value, m b The greater the dispersion, the greater the deviation of its maximum value from the average value. When ζ When the value is 0.5, in the RVE element (primary element) of the specimen m b The maximum value is 1.97, which is 9.85 times the average; when ζ When =10, in the RVE element of the specimen m b The maximum value is 0.22, which is 1.1 times the average. It can also be seen that... ζ The higher the value, the more RVE units within the rock... m b The values ​​fluctuate only slightly around the average value, indicating that the rock samples are relatively homogeneous. ζ The smaller the value, m b The values ​​fluctuate significantly around the average value, indicating that the rock sample has strong heterogeneity.

[0025] Figure 7a- Figure 7f Statistics were collected on different ζ Under the value, each m b The number of cells between value regions. (From Figure 7a-) Figure 7f It can be seen that when ζ When =0.5, in the overall model m b Only 133 RVEs had values ​​close to the mean of 0.2, accounting for approximately 10% of the total units. Meanwhile, m bThe large range of fluctuation in the values ​​indicates that the overall model is significantly heterogeneous. ζ The larger the value, the more RVE units there are. m b The closer the value is to the mean, the better. ζ When =4, m b The number of RVE units with values ​​less than 0.01 and greater than 0.3 is already zero. When ζ When =10, on each RVE unit m b The values ​​are only distributed between 0.1 and 0.3, and the overall mechanical properties of the rock specimens exhibit obvious homogeneity. Figure 7a- Figure 7f This was explained from the perspective of the number of RVE units in different regions. ζ The influence of the value on the homogeneity of rocks, i.e. ζ The higher the value, the more homogeneous the rock is; conversely, the lower the value, the stronger the heterogeneity of the rock.

[0026] Apply displacement boundary conditions at the top of the model and record the load under different load steps. m b The evolution of values ​​and damage values, the results are as follows: Figure 8a and Figure 8b As shown in the figure, damage gradually accumulates during loading, and the RVE elements with weaker mechanical parameters enter the plastic zone earlier, while their strength parameters... m b Under the influence of damage, it continuously weakens, eventually forming a fracture damage zone. Due to shape parameters... ζ The difference, m b The distribution density is different. When ζ When smaller ( ζ= 0.5), the simulated sample exhibits strong heterogeneity, and the resulting damage bands also differ significantly from those under homogeneous conditions. With ζ As the value increases, the RVE units in each sample... m b Value at average m b0 As the fluctuation range decreases, the heterogeneity of the simulated sample gradually weakens. At this point, the morphology of the damage zone generated by the sample under load is closer to that under completely homogeneous conditions. The calculation results show that the time-dependent rock Hoek-Brown elastoplastic damage model proposed in this application can effectively describe the heterogeneous characteristics of rocks and reproduce the fracture process of loaded rocks.

[0027] This embodiment analyzes the influence of time and shape parameters on localization: triaxial creep testing is used to evaluate the shape parameters. ζThe aging model established under condition 10 (highly homogeneous specimens) was validated. During the triaxial creep test, a constant confining pressure of 5 MPa was used, with an initial axial pressure of 10 MPa. When the creep deformation rate was less than 0.001 mm / d, the next loading stage was initiated, with each stage increasing by 20 MPa. Numerical calculations used the same loading method. Based on the parameter identification method, using gneiss creep test data as an example, the aging calculation parameters were obtained. P =0.8245e -3 , Q =-0.4731e -3 , R =0.025. The method for obtaining the parameters of the HB model is the same as described above, as detailed in Table 1.

[0028] The final comparison results between the calculated and experimental values ​​are as follows: Figure 9 As shown. From Figure 9 It is evident that when creep enters the steady-state creep stage, no significant damage appears. At lower stress levels, the accumulation of microscopic damage within the internal structure is insufficient to cause macroscopic failure of the specimen; therefore, it gradually approaches a stable state. At higher stress levels, the macroscopic mechanical properties of the specimen gradually deteriorate. When the microscopic damage accumulates to a certain extent, the accelerated creep stage begins, and the macroscopic mechanical properties of the specimen manifest as accelerated creep. The calculated results are close to the experimental results, demonstrating the correctness of the established model.

[0029] Figure 10a and Figure 10b ,exist ζ The variation of rock sample damage values ​​over time under conditions =6 and =8. As shown in the figure, under the final stage of load, load damage has already occurred inside the rock sample. ζ When the value is 6, the maximum internal damage value of the rock sample is 8.457e. -2 As time changes, the mechanical parameters of the rock sample weaken, leading to a continuous accumulation and increase in damage value under constant load (gradually increasing to 1.115e). -1 and 1.504 e -1 ), and eventually, destruction will occur. When ζ When the value is 8, the variation pattern is similar, but the shape of the damage zone changes, and the more homogeneous the sample ( ζ The larger the value, the more obvious the symmetry of the damage zone.

[0030] Long-term stability evaluation of slope engineering in cold regions: This application presents a numerical simulation study on the long-term stability of a slope in northern China affected by freeze-thaw cycles. The selected section K400+184 is located along the Dandong-Altay Expressway at the border of Ji'an and Baishan cities. The slope is adjacent to the Linfeng Reservoir. The main rock composition of the study section is strongly weathered quartzite, with well-developed joints and fissures, relatively fragmented rock mass, good structural bonding, massive structure, and a basic rock mass quality grade of Class IV. The geographical location of the slope is as follows: Figure 11 As shown.

[0031] A finite element calculation model of the slope is established based on the actual site conditions, such as... Figure 11 As shown in the figure. The model is 36m long and 22.6m high, divided into 4786 nodes and 4610 elements. Horizontal constraints are applied to both sides of the model, and fixed constraints are applied to the bottom. The mean values ​​of the calculation parameters, preliminarily determined based on the geological survey report and laboratory tests, are shown in Table 2. Assume the rock's elastic modulus, uniaxial compressive strength, etc. σ ci , m b and s The values ​​follow a Weibull distribution. The effects of shape parameter variations and time variations on the slope safety factor and damage mode are studied. Specifically, when... ζ When the value is 10, the corresponding parameter can be considered equivalent to the mean.

[0032] Table 2. Calculation parameters of the microscopic HB-D model

[0033] Analysis of the influence of shape parameters on slope stability: Figures 12a-12c Different shape parameters are given ζ The distribution map of the elastic modulus (for simplicity, only the distribution law of the elastic modulus is given) and the distribution cloud map of the slope damage zone are shown. As can be seen from the figure, within the selected parameter range ( ζ =2 to ζ =10), shape parameters ζ The impact on the maximum damage value of the slope is relatively small, with the magnitude of damage variation being on the order of 1×10. -4 Left and right. But shape parameters ζ It has a significant impact on the morphology of damage distribution. ζ The smaller the value, the more "branched" structures appear in the slope damage zone compared to homogeneous conditions. ζ As the value increases, the slope damage zone gradually becomes more concentrated and smoother, more closely resembling the morphology under homogeneous conditions. Calculation results show that, within a certain range, the Weibull parameter value has a relatively small impact on the overall safety of the slope, but a significant influence on its failure mode.

[0034] Analysis of the impact of aging on slope stability: Without considering freeze-thaw action, we only analyze the impact of aging on slope stability. Figures 13a-13f The variation of the maximum damage value of the slope after excavation over time is presented. The calculation results show that as time progresses, the rock mass parameters continuously weaken, approaching their long-term strength values, thus causing a significant increase in the damage value. Initially, the maximum damage value within the slope is 0.27, and this maximum damage value increases continuously over time. t In October, the maximum damage value inside the slope reached 0.2855, an increase of 0.0155 from the initial value. The rate of change of the damage value was more obvious in the early stage and gradually approached zero in the later stage.

[0035] Figure 14 This shows the variation of the slope safety factor over time. Similar to the variation of displacement and damage values, the safety factor decreases significantly in the initial stage, dropping from 1.59 to 1.27 in the first 12 months. After 12 months, the safety factor gradually stabilizes at 1.26 and no longer shows significant changes.

[0036] Analysis of the Impact of Aging and Freeze-Thaw Fatigue on Slope Engineering Stability: Based on field investigation, the selected section has been subjected to long-term disturbances caused by freeze-thaw fatigue cycles. Therefore, analyzing its long-term safety reserve under freeze-thaw fatigue cycle action is essential. Since rock undergoes coupled damage under the combined effects of freeze-thaw cycles and loads, and the freeze-thaw damage value is related to the number of freeze-thaw cycles, the relationship can be expressed as: (twenty four) In the formula, D r This is a fatigue damage variable, which in this embodiment is a freeze-thaw fatigue damage variable; N In this embodiment, the number of freeze-thaw fatigue cycles is used.

[0037] During the calculation process, the parameters in the model are considered to be E , m b and s The value changes over time (months), and the time-effect related parameters are the same as those mentioned above. Figure 15 This study investigates the variation of maximum slope damage over time after 10, 40, and 80 freeze-thaw cycles. The calculations assume the entire slope model is affected by freeze-thaw fatigue damage. The results show that when the number of freeze-thaw cycles... N When =10, the initial stage of excavation is... t At month 0, the maximum coupled damage value within the slope was 0.283. Over time, the rock mass parameters continuously weakened, approaching their long-term strength values, thus causing a significant increase in damage values. When time... tIn October, the maximum coupled damage value within the slope reached 0.3035, an increase of 0.0205 from the initial value. This was achieved after 20 freeze-thaw cycles. N =40, at t =At month 0, the maximum coupled damage value inside the slope is 0.3031. t =In 20 months, the maximum coupled damage value inside the slope increased by 0.3136. When the number of freeze-thaw cycles... N =80, at t =0 months, the maximum coupled damage value inside the slope is 0.3214, when t =In 20 months, the maximum coupled damage value inside the slope increased to 0.328. The rate of change of the coupled damage value was more obvious in the early stage and gradually approached zero in the later stage.

[0038] Figures 16a-16f This shows the change in slope displacement over time after 80 freeze-thaw cycles. The calculation results indicate that the slope displacement changes over time after excavation. Using the maximum displacement value as the standard, when... t When the slope reaches 0 degrees Celsius, the maximum displacement is 0.016 m. t In 20 months, the maximum displacement of the slope was 0.023m, an increase of approximately 41.8%. Figure 17 The figure shows the displacement variation over time at various monitoring points on the slope. As can be seen from the figure, the slope displacement changes significantly in the initial stage of excavation, but gradually stabilizes and no longer changes over time. Taking monitoring point No. 2, which has the largest displacement value, as an example, when... t When the months are 0, 4, 8, 12, 16 and 20, the displacement of this point is 0.015m, 0.018m, 0.02m, 0.021m, 0.0214m and 0.216m, respectively. The magnitude of the change gradually decreases and eventually approaches zero.

[0039] Prestressed anchor cables were used on-site to reinforce the slope and ensure its long-term stability. The anchor cable lengths in this project ranged from 15 to 18 meters, with an anchorage length of 8 meters. Each anchor cable consisted of 6 bundles of steel strands, and the anchor cable inclination angle was 20°. The prestressing of the anchor cables was designed to be 600 kN. To reasonably calculate the prestressing of the anchor cables, the contact problem between the soil and the anchor cables was calculated using the Embedded region method based on the above model. The prestressing of the anchor cables was set using the temperature field method, and the prestressing calculation formula is as follows: (25) In the formula, ΔT The preset temperature reduction value is in °C. F N This represents the prestress value of the anchor cable, in Newtons (N). E This refers to the elastic modulus of the material, expressed in Pa. AThe cross-sectional area is expressed in meters (m²). 2 ; is the coefficient of thermal expansion of the anchor cable, expressed in 1 / °C.

[0040] With anchor cable prestresses set at 200kN, 300kN, 400kN, 500kN, and 600kN respectively, the final damage zone distribution of the slope was calculated as follows: Figures 18a-18d As shown.

[0041] The calculation results show that after applying prestressed anchor cables, the slope damage zone shifts inward, reducing the potential risk of slope collapse. As the prestress increases, the slope damage value gradually decreases. Compared to when no prestressed anchor cables are applied, the slope damage value decreases significantly when 200kN anchor cables are applied, with the maximum damage value decreasing from 0.328 to 0.202. When the prestress exceeds 400kN, the slope damage value no longer changes significantly. Based on the calculation results and considering construction costs, it is recommended that the anchor cable prestress be set at 400kN for this slope.

[0042] This invention provides a rock mass stability calculation method considering rheological and fatigue effects, which can evaluate the long-term stability of geotechnical engineering under freeze-thaw fatigue disturbance, and provide a methodological basis for the design of construction parameters for engineering projects in special areas. Based on the failure zones of rock samples obtained from numerical compression tests under different shape parameters, and the distribution of damage zones in slope engineering stability calculations, this method demonstrates its rationality.

[0043] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for calculating rock mass stability considering fatigue and time effects, characterized in that, Includes the following steps: S1: Divide the rock into several basic units. V In order to obtain the mechanical parameters of the basic element; The mechanical parameters of the basic elements include empirical parameters that reflect the characteristics of the rock mass. m b Parameters of rock mass fragmentation and integrity s Rock mass quality parameters a and elastic modulus before damage E 0; S2: Based on the mechanical parameters of the basic element, obtain the Hoek-Brown elastoplastic damage model of the basic element; S3: Obtain the time-dependent damage variable of the primitive. ; S4: Obtain the fatigue damage variable of the primitive. D r ; S5: Based on load damage variables D m The aforementioned time-related damage variables and the fatigue damage variable D r Obtain the coupling damage variable Dc: in, (15) In the formula, Dc is the coupling damage variable; D m For load damage variables; For time-related damage variables; D r For fatigue damage variables; D t For the first Damage values ​​of rock cross-sections under secondary fatigue disturbance; This is the number of fatigue disturbances; S6: Based on the Hoek-Brown elastoplastic damage model of the basic element and the coupled damage variables, obtain... t m+1 nominal stress tensor at time step ; S7: According to the... t m+1 nominal stress tensor at time step , obtain the t m+1 Values ​​of coupling damage variables at each time step; S8: According to the... t m+1 The stability of the rock is analyzed by taking the values ​​of the coupled damage variables at each time step.

2. The method for calculating rock mass stability considering fatigue and time effects according to claim 1, characterized in that, In S1, the mechanical parameters of the basic element are obtained as follows: (1) In the formula, P 0 represents the macroscopic properties of the rock sample; x z Use Weibull random numbers; P z These are the mechanical parameters of the basic unit; Among them, Weibull random numbers x z The method to obtain it is as follows: (2) (3) In the formula: W ( x ) is the cumulative distribution function of the Weibull distribution; w ( x ) is the probability density function of the Weibull distribution; >0 is the scaling parameter; ζ For shape parameters; x It is a random variable; (4) In the formula: These are uniformly distributed random numbers.

3. The method for calculating rock mass stability considering fatigue and time effects according to claim 1, characterized in that, In S2, the Hoek-Brown elastoplastic damage model of the basic element is obtained as follows: (5) (6) In the formula, f Let be the yield function; J 2 represents the second deviatoric stress invariant; θ σ For Rhodes angle; D m For load damage variables; σ ci The uniaxial compressive strength of intact rock; p It is the hydrostatic pressure; m bg for m b The corresponding plastic potential function parameters; s g for s The corresponding plastic potential function parameters; a g for a The corresponding plastic potential function parameters; g Let be the plastic potential function; (7) In the formula, α To determine the initial slope of the softening curve of the rock material after damage, the value range is (0, +∞); β The parameter used to determine the maximum damage value of the rock has a value range of [0,1). Equivalent plastic strain; (8) In the formula, ε p1 The principal plastic strain along the x-axis in the spatial coordinate system; ε p2 The principal plastic strain along the y-axis in the spatial coordinate system; ε p3 The value represents the principal plastic strain along the z-axis in the spatial coordinate system.

4. The method for calculating rock mass stability considering fatigue and time effects according to claim 1, characterized in that, In step S3, the time-related damage variable of the primitive is obtained. The method is as follows; (9) In the formula, In response to the situation; For axial stress; K 1. K 2 are the mechanical parameters of the elastic element; η The viscosity coefficient of the sticky pot; For the time of rheometry; Simplifying equation (9) yields (10) In the formula, P, Q, and R are all intermediate parameters, where, , , ; From formula (10), we can obtain: (11) E ( t () represents the elastic modulus as a function of time. t A changing function; (12) In the formula, These are time-dependent damage variables, i.e., damage variables caused by time effects; E 0 represents the elastic modulus before damage.

5. The method for calculating rock mass stability considering fatigue and time effects according to claim 1, characterized in that, In S4, the fatigue damage variable of the basic element D r The following information was obtained: (13) In the formula: D t For the first Damage values ​​of rock cross-sections under secondary fatigue disturbance; D r For fatigue damage variables D r That is, the first N fatigue damage variables prior to the first fatigue disturbance; This is the number of fatigue disturbances; N For the first Total number of fatigue disturbances before the current fatigue disturbance; W and U All are material parameters; (14) In the formula, n 0 represents the porosity of the rock before fatigue. n Porosity of the rock after fatigue. The longitudinal wave velocity is the velocity of the rock after fatigue.

6. The method for calculating rock mass stability considering fatigue and time effects according to claim 1, characterized in that, In S6, obtain t m+1 nominal stress tensor at time step The method is as follows: S61: The stress-strain constitutive relation considering time-fatigue-load coupled damage is as follows: (16) In the formula, Stress to account for time-fatigue-load coupled damage; For elastic strain; Stiffness matrix considering coupled damage; , It is a fourth-order symmetric tensor. For Kronecker product; G ( D c () represents the damage shear modulus; K ( D c () represents the bulk modulus; in, (17) In the formula, G 0 represents the initial shear modulus; K 0 represents the initial bulk modulus; Poisson's ratio, For the damage variable is D C The elastic modulus at that time; (18a) (18b) In the formula: for t Moment value; For the initial value; for t Moment s value; For the initial s value; S62: Obtain the... t m+1 Predicted stress at time steps to determine whether the rock has entered the plastic correction stage; Assumption Then, given the strain increment The expression for the predicted stress is obtained as follows: (19) In the formula, For the first t m+1 Predicted stress at time step; Stiffness matrix considering coupled damage; For the first t m Response under time steps; For the first t m The given strain increment at the time step; For the first t m+1 Predicting coupled damage variables at time steps; For the first t m Coupling damage variables at time step; This represents a double dot product; S63: Based on the yield function f ,like: (20) Then the first t m+1 Stress at time step equal to the t m+1 Predicted stress at time step, i.e. If necessary, execute S62; otherwise, execute S64. S64: Based on the plastic potential function, the Hoek-Brown elastoplastic damage model is plastically modified to obtain the modified equivalent plastic strain. S65: Obtain from the corrected equivalent plastic strain t m+1 nominal stress tensor at time step as follows: (22) In the formula: and They are the first t m and t m+1 Damage elasticity matrix at time step; For the first t m+1 The plastic strain increment at the time step; It is the first t m+1 Coupling damage variables at time step; For the first t m Coupling damage variables at time steps.

Images