Improved triangular mesh numerical manifold method for simulating progressive failure of jointed rock slopes

Abstract

The present invention discloses an improved triangular grid numerical manifold method for simulating the progressive failure of a jointed rock slope, which includes the following steps: establishing an initial geometric model of the jointed rock slope; establishing a triangular grid discrete calculation model of the initial geometric model; constructing a global displacement field of a singular manifold unit according to the triangular grid discrete calculation model, and calculating the singular function integral through an enhanced generalized Duffy transform formula; calculating the contact force according to a cone complementary contact force calculation model that simultaneously considers friction and adhesion; calculating the stress field near the tip of the crack according to the singular function integral and the contact force, and determining the cracking condition of the crack contained in the jointed rock slope in combination with the critical length cracking criterion. The improved method proposed by the present invention can obtain a more accurate crack propagation path and progressive failure process, which is mainly reflected in: being able to reasonably describe the end displacement discontinuity and stress singularity near the crack tip; the numerical integration accuracy of the singular function is higher; improving the accuracy and stability of the contact force calculation; and being able to simultaneously consider the tensile shear and compressive shear fractures of the rock mass.

Description

Technical Field

[0001] The invention belongs to the field of numerical calculation methods and simulations of geotechnical engineering, engineering geology and mountain disasters, and specifically relates to an improved triangular grid numerical manifold method for simulating progressive failure of jointed rock slopes. Background Art

[0002] Natural rock slopes often contain discontinuities such as joints, faults, and cracks. The progressive destruction process from continuous small deformation to discrete large movement is very common in mountain disasters such as collapse and landslide. From the theoretical analysis of uniaxial tensile cracks, it can be seen that in addition to the discontinuity of the displacement field near the crack tip, the stress field also has singularity. Four basis functions can be introduced on the basis of the classical numerical manifold method: Describe this singularity. This involves two singular terms and 1 / r. However, traditional numerical integration methods such as the Gauss integration rule cannot achieve satisfactory integration accuracy. Secondly, the progressive failure process of jointed rock slopes involves a large number of contact interactions between crack surfaces and between fractured bodies. The contact force calculation of the classical numerical manifold method is problematic and empirically dependent, and its stability is insufficient. In addition, the cracking criterion used in the classical numerical manifold method cannot accurately reflect the tensile shear and compressive shear failure of rock materials at the same time. Therefore, the classical triangular mesh numerical manifold method often lacks accuracy, has poor stability, and its computational efficiency needs to be improved when simulating the progressive failure process of jointed rock slopes.

[0003] Accurate prediction and simulation of the progressive failure process of jointed rock slopes are critical to understanding the failure mechanism, evolution process, and designing preventive measures. In view of the shortcomings of the classical triangular mesh numerical manifold method, the following four improvements are needed: first, it should be able to reflect the discontinuity of the displacement field and the singularity of the stress field near the crack tip; second, the numerical integration accuracy of the singular function should meet the requirements; third, the calculation of the contact force between crack surfaces or discrete bodies should be accurate and stable; fourth, under complex stress conditions, it should be able to describe the tensile shear and compressive shear failure of geotechnical materials. Summary of the invention

[0004] In view of the above-mentioned deficiencies in the prior art, the improved triangular grid numerical manifold method for simulating the progressive failure of jointed rock slopes provided by the present invention solves the problems that the classical triangular grid numerical manifold method often has insufficient accuracy, poor stability and needs to improve calculation efficiency when simulating the progressive failure process of jointed rock slopes.

[0005] In order to achieve the above-mentioned invention object, the technical solution adopted by the present invention is: an improved triangular grid numerical manifold method for simulating the progressive failure of jointed rock slopes, comprising the following steps:

[0006] S1. Establish the initial geometric model of jointed rock slope;

[0007] S2, establishing a triangular mesh discrete calculation model of the initial geometric model;

[0008] S3, constructing the global displacement field of the singular manifold unit according to the triangular mesh discrete calculation model, and calculating the singular function integral by the enhanced generalized Duffy transform formula;

[0009] S4, calculating the contact force according to the cone complementary contact force calculation model that considers both friction and adhesion;

[0010] S5. The stress field near the tip of the crack is calculated based on the singular function integral and the contact force, and the cracking condition of the cracks in the jointed rock slope is determined in combination with the critical length cracking criterion.

[0011] Further: S2 is specifically:

[0012] The initial geometric model of the jointed rock slope is covered with a three-node triangular mesh, and geometric operations of cutting and intersecting the three-node triangular mesh and the initial geometric model are performed to generate manifold units and obtain a triangular mesh discrete calculation model.

[0013] Further: In S3, the method for calculating the singular function integral is specifically:

[0014] The singular manifold unit is triangulated with the crack tip as the vertex to generate several sub-triangles. The integral of each sub-triangle is calculated using the enhanced generalized Duffy transform formula, and the integrals of all sub-triangles are added together to obtain the singular function integral.

[0015] Further: In S3, the global displacement field u of the singular manifold unit S The expression for (x,y) is:

[0016]

[0017] Where the subscript S represents any singular manifold unit, is the undetermined degree of freedom in the x and y directions of the i-th node of the triangle of the singular manifold unit, and its specific expression is:

[0018]

[0019] In the formula, u i and v i are the displacements of node i in the x and y directions respectively; All are undetermined coefficients;

[0020] T S i is the node interpolation matrix of the singular manifold unit, and its specific expression is:

[0021]

[0022] In the formula, w i is the shape function of the i-th node of the triangle, {ψ1, ψ2, ψ3, ψ4} are the four trigonometric functions of the first two terms of the Williams displacement series introduced near the crack tip, and their specific expressions are:

[0023]

[0024] Where r and θ are the polar diameter and polar angle of the local polar coordinates with the crack tip as the pole and the crack direction as the polar axis.

[0025] Further: In S3, the expression of the enhanced generalized Duffy transform formula is specifically:

[0026]

[0027] Where I is the numerical integral of the singular function, and β is the minimum integer that ensures that 2β-1-αβ is a positive integer, which is used to simplify the integral term u 2β-1-αβ , α is a real number, f(u,w) is an arbitrary function of the integral variables u and w, e w and e -w are all exponential functions of the integral variable w.

[0028] Further: In S4, the expression of the cone complementary contact force calculation model considering both friction and adhesion is specifically:

[0029] p j =(p n ,p s ,p t ) T

[0030]

[0031] In the formula, p n is the normal contact force, p s and p t is the contact force in two mutually perpendicular tangential directions on the tangent plane, g n is the normal contact distance, g s and g t is the contact distance between two mutually perpendicular tangential directions on the tangent plane, μ * is the equivalent friction coefficient, is the contact spacing vector, p j is the contact force vector, Υ j is the friction cone, It is the dual cone corresponding to the friction cone.

[0032] Further: S5 includes the following sub-steps:

[0033] S51, calculating the principal stress at the crack tip, including the first to third principal stresses;

[0034] S52, judging whether the first principal stress is greater than the material tensile stress, if so, proceeding to S53; if not, proceeding to S54;

[0035] S53, finding the point closest to the tensile stress of the material in the direction perpendicular to the first principal stress, and calculating the first distance from the crack tip, judging whether the first distance is not less than the crack extension length in the direction of the crack surface during tensile failure, if so, the current crack tip is extended by the length of the first distance in the direction perpendicular to the maximum principal stress, and entering S55, if not, the current crack does not extend, and entering S55;

[0036] S54, calculating the Mohr-Coulomb strength force according to the first to third principal stresses, determining whether the Mohr-Coulomb strength force is greater than 0, and if so, searching for the point whose Mohr-Coulomb strength force is closest to 0 along the first set direction, and calculating the second distance from the crack tip, searching for the point whose Mohr-Coulomb strength force is closest to 0 along the second set direction, and calculating the third distance from the crack tip, determining the cracking condition of the current crack according to the relationship between the second distance, the third distance and the crack extension length along the direction of the crack surface during tensile failure, and entering S55;

[0037] S55. Based on the obtained cracking condition of the current crack, determine the cracking condition of the cracks contained in the jointed rock slope.

[0038] Further: In S52, the first principal stress is specifically the maximum principal stress perpendicular to the crack surface near the crack tip, and its expression is specifically:

[0039]

[0040] In the formula, K I is the I-mode stress intensity factor, L T is the crack extension length along the direction of the crack surface during tensile failure, and its specific expression is:

[0041]

[0042] In the formula, K IC is the critical I-mode stress intensity factor for material instability failure, σ T The tensile stress of the material.

[0043] Further: In S54, the expression for calculating the Mohr-Coulomb strength force is specifically:

[0044] F(σ1,σ3)=σ1(1-sinφ)-σ3(1+sinφ)-2ccosφ=0

[0045] In the formula, c is the cohesion of the material, φ is the internal friction angle of the material, and σ3 is the third principal stress;

[0046] First set direction The specific expression is:

[0047]

[0048] In the formula, β3 is the direction of the minimum principal stress, β f is the set stress direction, and β f =π / 4+φ / 2;

[0049] Second setting direction The specific expression is:

[0050]

[0051] The method for determining the cracking condition of the current crack is specifically as follows:

[0052] If max(L A ,L B )<L T , then the current crack does not extend, where L A is the second distance, L B is the second distance;

[0053] If L A ≥L T >L B , the current crack extends along the first set direction to a second distance;

[0054] If L B ≥L T >L A , the current crack extends along the second set direction by a third distance;

[0055] If L A ≥L T And L B ≥L T , the crack extends to the length of the second distance and the third distance respectively along the first set direction and the second set direction at the same time.

[0056] The beneficial effects of the present invention are:

[0057] (1) The present invention first adds the four trigonometric functions contained in the first two terms of the Williams displacement series near the crack tip to the global displacement approximation of the manifold unit where the crack tip is located to reflect the properties of the field variables at the crack tip; then proposes an enhanced generalized Duffy transform method to improve the numerical integration accuracy of singular functions; and then adopts a cone complementary contact force calculation model that can simultaneously consider friction and adhesion to accurately calculate the contact force during the progressive failure process of the rock slope; finally, the critical length cracking criterion is introduced to determine the direction and distance of crack propagation. Compared with the classic triangular mesh numerical manifold method, the improved method proposed in the present invention can obtain a more accurate crack propagation path and progressive failure process, which is mainly reflected in:

[0058] 1) It can reasonably describe the displacement discontinuity and stress singularity near the crack tip;

[0059] 2) The numerical integration accuracy of singular functions is higher;

[0060] 3) Improved the accuracy and stability of contact force calculation;

[0061] 4) It can simultaneously consider the tensile shear and compressive shear failure of rock mass.

[0062] (2) The present invention implements an improved triangular grid numerical manifold method, which can obtain a reasonable prediction of the progressive failure evolution process of jointed rock slopes, improve the accuracy and stability of numerical calculations, and can be adopted by relevant scientific research, commercial computer codes and programs. BRIEF DESCRIPTION OF THE DRAWINGS

[0063] Figure 1 The present invention is a flow chart of the improved triangular grid numerical manifold method for simulating the progressive failure of jointed rock slopes.

[0064] Figure 2 The manifold unit containing the crack tip and its triangulation.

[0065] Figure 3 Verify the geometric model for the enhanced generalized Duffy transform method.

[0066] Figure 4 The calculation results are verified for the enhanced generalized Duffy transform method.

[0067] Figure 5 Validate the geometry for contact force calculations.

[0068] Figure 6 Schematic diagram of the tensile shear failure criterion.

[0069] Figure 7 Schematic diagram of compression shear failure criterion.

[0070] Figure 8 This is the geometric model for the progressive failure analysis of jointed rock slope.

[0071] Fig. 9 These are the calculation results of the progressive failure example of jointed rock slope. DETAILED DESCRIPTION

[0072] The specific implementation modes of the present invention are described below so that those skilled in the art can understand the present invention. However, it should be clear that the present invention is not limited to the scope of the specific implementation modes. For those of ordinary skill in the art, as long as various changes are within the spirit and scope of the present invention as defined and determined by the attached claims, these changes are obvious, and all inventions and creations utilizing the concept of the present invention are protected.

[0073] Embodiment 1:

[0074] like Figure 1 As shown, in one embodiment of the present invention, the improved triangular grid numerical manifold method for simulating the progressive failure of jointed rock slope includes the following sub-steps:

[0075] S1. Establish the initial geometric model of jointed rock slope;

[0076] S2, establishing a triangular mesh discrete calculation model of the initial geometric model;

[0077] S3, constructing the global displacement field of the singular manifold unit according to the triangular mesh discrete calculation model, and calculating the singular function integral by the enhanced generalized Duffy transform formula;

[0078] S4, calculating the contact force according to the cone complementary contact force calculation model that considers both friction and adhesion;

[0079] S5. The stress field near the tip of the crack is calculated based on the singular function integral and the contact force, and the cracking condition of the cracks in the jointed rock slope is determined in combination with the critical length cracking criterion.

[0080] In S1, the initial geometric model is established by commercial or open source pre-processing software, such as AutoCAD and FreeCAD.

[0081] The S2 is specifically:

[0082] The initial geometric model of the jointed rock slope is covered with a three-node triangular mesh, and geometric operations of cutting and intersecting the three-node triangular mesh and the initial geometric model are performed to generate manifold units and obtain a triangular mesh discrete calculation model.

[0083] In this embodiment, the manifold unit includes a non-singular manifold unit and a singular manifold unit. The non-singular manifold unit is a manifold unit without a crack tip, and the singular manifold unit is a manifold unit with a crack tip.

[0084] In S3, the global displacement field u of the singular manifold unit S The expression for (x,y) is:

[0085]

[0086] Where the subscript S represents any singular manifold unit, is the undetermined degree of freedom in the x and y directions of the i-th node of the triangle of the singular manifold unit, and its specific expression is:

[0087]

[0088] In the formula, u i and v i are the displacements of node i in the x and y directions respectively; All are undetermined coefficients;

[0089] T S i is the node interpolation matrix of the singular manifold unit, and its specific expression is:

[0090]

[0091] In the formula, w i is the shape function of the i-th node of the triangle, {ψ1, ψ2, ψ3, ψ4} are the four trigonometric functions of the first two terms of the Williams displacement series introduced near the crack tip, and their specific expressions are:

[0092]

[0093] Where r and θ are the polar diameter and polar angle of the local polar coordinates with the crack tip as the pole and the crack direction as the polar axis.

[0094] In S3, the method for calculating the singular function integral is specifically:

[0095] The singular manifold unit is triangulated with the crack tip as the vertex to generate several sub-triangles. The integral of each sub-triangle is calculated using the enhanced generalized Duffy transform formula, and the integrals of all sub-triangles are added together to obtain the singular function integral.

[0096] In S3, the expression of the enhanced generalized Duffy transform formula is specifically:

[0097]

[0098] Where I is the numerical integral of the singular function, and β is the minimum integer that ensures that 2β-1-αβ is a positive integer, which is used to simplify the integral term u 2β-1-αβ , α is a real number, f(u,w) is an arbitrary function of the integral variables u and w, ew and e -w are all exponential functions of the integral variable w.

[0099] In S4, the expression of the cone complementary contact force calculation model considering both friction and adhesion is specifically:

[0100] p j =(p n ,p s ,p t ) T

[0101]

[0102] In the formula, p n is the normal contact force, p s and p t is the contact force in two mutually perpendicular tangential directions on the tangent plane, g n is the normal contact distance, g s and g t is the contact distance between two mutually perpendicular tangential directions on the tangent plane, μ * is the equivalent friction coefficient, is the contact spacing vector, p j is the contact force vector, Υ j is the friction cone, It is the dual cone corresponding to the friction cone.

[0103] The method for calculating the contact force according to the cone complementary contact force calculation model that considers both friction and adhesion is as follows:

[0104]

[0105] In the formula, the superscript k is the calculation step, and the subscript r is the number of iterations. is the contact force vector, ω is the algorithm parameter, and 0<ω<1. In this embodiment, ω=0.2, N, b and They are all temporary matrices and are calculated as follows:

[0106]

[0107] Among them, H is the transformation matrix that converts the contact force in the local coordinate system at the contact point into the global contact force vector, which is determined by the contact normal and tangent directions, M is the equivalent mass matrix, F is the equivalent load vector, and I is the unit matrix. is the lower triangular matrix of H, △t is the time step, Trace() means to find the sum of the diagonal elements of the matrix, Proj{} is the projection operator, which is calculated as follows:

[0108]

[0109] Where μ is the friction coefficient and c is the bonding strength;

[0110] The S5 comprises the following sub-steps:

[0111] S51, calculating the principal stress at the crack tip, including the first to third principal stresses;

[0112] S52, judging whether the first principal stress is greater than the material tensile stress, if so, proceeding to S53; if not, proceeding to S54;

[0113] S53, finding the point closest to the tensile stress of the material in the direction perpendicular to the first principal stress, and calculating the first distance from the crack tip, judging whether the first distance is not less than the crack extension length in the direction of the crack surface during tensile failure, if so, the current crack tip is extended by the length of the first distance in the direction perpendicular to the maximum principal stress, and entering S55, if not, the current crack does not extend, and entering S55;

[0114] S54, calculating the Mohr-Coulomb strength force according to the first to third principal stresses, determining whether the Mohr-Coulomb strength force is greater than 0, and if so, searching for the point whose Mohr-Coulomb strength force is closest to 0 along the first set direction, and calculating the second distance from the crack tip, searching for the point whose Mohr-Coulomb strength force is closest to 0 along the second set direction, and calculating the third distance from the crack tip, determining the cracking condition of the current crack according to the relationship between the second distance, the third distance and the crack extension length along the direction of the crack surface during tensile failure, and entering S55;

[0115] S55. Based on the obtained cracking condition of the current crack, determine the cracking condition of the cracks contained in the jointed rock slope.

[0116] In S52, the first principal stress is specifically the maximum principal stress perpendicular to the crack surface near the crack tip, and its expression is specifically:

[0117]

[0118] In the formula, K I is the I-mode stress intensity factor, L T is the crack extension length along the direction of the crack surface during tensile failure, and its specific expression is:

[0119]

[0120] In the formula, K IC is the critical I-mode stress intensity factor for material instability failure, σ T The tensile stress of the material.

[0121] In S54, the expression for calculating the Mohr-Coulomb strength force is specifically:

[0122] F(σ1,σ3)=σ1(1-sinφ)-σ3(1+sinφ)-2ccosφ=0

[0123] In the formula, c is the cohesion of the material, φ is the internal friction angle of the material, and σ3 is the third principal stress;

[0124] First set direction The specific expression is:

[0125]

[0126] In the formula, β3 is the direction of the minimum principal stress, β f is the set stress direction, and β f =π / 4+φ / 2;

[0127] Second setting direction The specific expression is:

[0128]

[0129] The method for determining the cracking condition of the current crack is specifically as follows:

[0130] If max(L A ,L B )<L T , then the current crack does not extend, where L A is the second distance, L B is the second distance;

[0131] If L A ≥L T >L B , the current crack extends along the first set direction to a second distance;

[0132] If L B ≥L T >L A , the current crack extends along the second set direction by a third distance;

[0133] If L A ≥L T And L B ≥L T , the crack extends to the length of the second distance and the third distance respectively along the first set direction and the second set direction at the same time.

[0134] Embodiment 2:

[0135] This example provides experimental data for Example 1.

[0136] Step 1: Establish the initial geometric model and calculation parameters of the jointed rock slope;

[0137] The initial geometric model of the jointed rock slope is established using commercial or open source pre-processing software, and the calculation parameters are given at the same time.

[0138] Step 2, establishing a triangular mesh discrete calculation model of the initial geometric model;

[0139] The initial geometric model of the entire jointed rock slope is covered with a three-node triangular mesh. The mesh and the initial geometric model are subjected to geometric operations such as cutting and intersection to generate manifold units and form a mesh discrete calculation model.

[0140] Step 3, construct the global displacement field of the singular manifold unit;

[0141] like Figure 2 As shown, the triangle contains the crack tip T(0.193118,1.0), and its three node coordinates are A(0.13535,1.0548), B(0.180931,0.98), and C(0.226512,1.0548), so w 1 =-4.566083+14.101604x-8.535521y,w 2 =-10.969422+0x+10.969422y,w 3 =6.684484-13.368983x+6.684499y. Assume that there is a point D with coordinates D(-0.137146,1.052645). Figure 2 The local polar coordinate system shown in the figure shows that the local polar coordinates corresponding to point D are (r,θ)=(0.07684,2.386809), and we can get

[0142]

[0143] ={0.102147, 0.257693, -0.250960, -0.117724}, and further we can get T S i (x,y).

[0144] is the degree of freedom to be solved.

[0145] Step 4: Triangulate the manifold unit containing the crack tip with the crack tip as the vertex, see Figure 2The manifold unit △123 is divided into four sub-triangles △T12, △T2F, △TF3 and △T31. The following enhanced generalized Duffy transform formula is applied to each sub-triangle to calculate the integral, and then the singular function integral is calculated by adding them together.

[0146] The expression of the enhanced generalized Duffy transform formula is as follows:

[0147]

[0148] Where β is the smallest integer that ensures that 2β-1-αβ is a positive integer to simplify the integral term u 2β-1-αβ .

[0149] like Figure 3 In the enhanced generalized Duffy transform method verification geometric model shown in the figure, it is assumed that the crack tip is at the vertex A(0,0) of the triangle △ABC (integral domain), the coordinates of point B are B(0,1), and point C is on the ray that is 30° from the positive x-axis. By changing the distance l∈[0.1,10] between points A and C, different points C(lcos30°,lsin30°) can be obtained. Accordingly, the triangle △ABC may be an acute triangle, a right triangle, or an obtuse triangle, so as to test the influence of the triangle shape on the calculation accuracy of the enhanced generalized Duffy transform method.

[0150] For the singular integral term 1 / r = 1 / r α (α=1,β=1), the integral ∫1 / rd△ is calculated by using the generalized Duffy transform (GDT) and enhanced generalized Duffy transform (MGDT) based on 8×8 Gaussian integration points, and the reference solution is obtained by using the generalized Duffy transform method based on 50×50 Gaussian integration points. The relative error is shown in Figure 4 As shown in the curve on the left side of : when the distance l between points A and C changes from 0.1 to 10, the relative errors of the GDT and MGDT methods are both less than 0.01%. At the same time, it can be seen that from l>1.3, the calculation accuracy of the MGDT method is improved by about 1 order of magnitude compared with the GDT method.

[0151] For the singular integral Similarly, the integral is calculated using the generalized Duffy transform (GDT) and enhanced generalized Duffy transform (MGDT) based on 8×8 Gaussian integration points. △, the reference solution is obtained by using the generalized Duffy transform method and 50×50 Gaussian integration points. The relative error is as follows Figure 4As shown in the curve on the left side of the graph, when l changes from 0.1 to 10, the relative errors of the GDT and MGDT methods are both less than 0.01%. At the same time, it can be seen that from l>1.0, the calculation accuracy of the MGDT method is about 10 times higher than that of the GDT method. It can be seen that the enhanced generalized Duffy transform method proposed in the present invention effectively improves the singular terms 1 / r and The numerical integration accuracy of .

[0152] Step 5: Accurate calculation of contact force is very important in the simulation of progressive failure of jointed rock slopes. The present invention proposes a cone complementary contact force calculation model that considers both friction and adhesion using the Lagrange multiplier method. Figure 5 The verification example shown in the figure uses three rigid disks (radius r = 10m, density ρ = 2500kg / m 3 , gravitational acceleration g = -10m / s 2 ) are stacked on a fixed foundation. Without considering friction, the analytical solutions of the three contact forces are F1 = 7853981.63397N, F2 = 15707963.26795N, and F3 = 23561944.90192N, as shown in Table 1.

[0153] Table 1 Calculation results of contact forces F1, F2 and F3

[0154]

[0155] It can be seen from Table 1 that, for this verification example, the contact force calculation accuracy of the improved numerical manifold method proposed in the present invention is improved by more than 2 orders of magnitude compared with the classical numerical manifold method, which verifies the effectiveness of the cone complementary contact force calculation model and calculation process established in the present invention.

[0156] Step 6, using the critical length cracking criterion to determine the tensile shear and compressive shear fracture of the rock material;

[0157] For the tensile shear fracture of rock materials, Figure 6 As shown in , a local polar coordinate system is established along the crack line. Under tensile shear stress conditions, the crack will extend in a direction that is θ-angled with the crack line, where θ is the angle between the vertical direction of the maximum principal stress and the crack line (counterclockwise is positive). For the compression shear fracture of rock materials, as shown in Figure 7 As shown, a local polar coordinate system is established along the crack line. Under compressive shear stress conditions, there are two crack extension directions, located on the upper and lower sides of the crack surface. The following steps are used to determine whether a crack has occurred:

[0158] (1) Calculate the three principal stresses σ1, σ2 and σ3 at the crack tip;

[0159] (2) When σ1>σ TWhen σ1≈σ T point and calculate its distance L from the crack tip θ , if L θ ≥L T Then the crack tip expands along the direction θ θ length; otherwise, the crack does not extend.

[0160] (3) When F(σ1,σ2,σ3)>0, along the direction Find the point where F≈0 and calculate its distance L from the crack tip A ; Along direction Find the point where F≈0 and calculate its distance L from the crack tip B .further:

[0161] ① If L M =max(L A ,L B )<L T , then the current crack does not expand;

[0162] ②If L A ≥L T >L B , then the current crack is along the direction Extended L A length;

[0163] ③If L B ≥L T >L A , then the current crack is along the direction Extended L B length;

[0164] ④ If L A ≥L T And L B ≥L T , then the crack along and Expand L A and L B length.

[0165] Specific application example 1: Progressive failure of jointed rock slope:

[0166] 1.1. Establish the initial geometric model of a jointed rock slope Figure 8 As shown: The entire geometric model is 7.2m wide and 7.6m high, with very obvious vertical cracks and nearly horizontal faults. The slope will experience typical failure modes such as collapse, rockfall, dumping, and sliding. The main discontinuity surface is represented by a red solid line. The calculation parameters are set as follows: The rock mass density is 2500kg / m 3, Young's modulus is E = 5 GPa, Poisson's ratio is v = 0.25, and the tensile strength is σ T =1MPa, internal friction angle φ = 25°, cohesion c = 1MPa, gravitational acceleration g = -20m / s 2 The time step is 0.0005s, the total time is 500 calculation steps, a total of 0.25s. The contact surface friction coefficient μ = 0.5773 (equivalent to the friction angle ψ = 30°).

[0167] 1.2. The initial geometric model of the entire jointed rock slope is covered with a three-node triangular mesh. Figure 8 As shown in Figure 2. The mesh and the initial geometry model are cut and intersected, and finally 2589 manifold elements are generated. We will focus on seven crack tips ( Figure 8 The progressive expansion process (marked with ①, ②, ③, …, ⑦, etc.) demonstrates the simulation capability of the improved numerical manifold method proposed in the present invention.

[0168] 1.3. Construct the global displacement field for the manifold unit without crack tip and the manifold unit with crack tip respectively.

[0169] 1.4. Divide the manifold unit containing the crack tip into triangles with the crack tip as the vertex, and apply the enhanced generalized Duffy transform formula to calculate the integral for each sub-triangle, and then add them together.

[0170] 1.5. All contact forces are calculated using the cone complementary contact force calculation model proposed by the present invention that can simultaneously consider friction and adhesion.

[0171] 1.6. Calculate the seven crack tips ( Figure 8 The three principal stresses σ1, σ2 and σ3 (marked with ①, ②, ③, …, ⑦) are then used to determine whether each crack tip will extend and the extension length and direction according to the critical length cracking criterion.

[0172] 1.7. As the calculation steps progress and steps 1.3 to 1.6 above are repeated, the seven crack tips ( Figure 8 The gradual expansion process (marked with ①, ②, ③, …, ⑦) is shown in Fig. 9The upper two rows. It can be seen that from step 13 to step 100, the crack gradually expands under the action of gravity: at step 13, the tip of crack ① expands first; at step 16, the tip of crack ③ begins to expand; followed by cracks ⑤ (step 17) and ④ (step 18); at step 46, the tips of cracks ③ and ⑤ expand for the second time at the same time; from step 48 to step 30, the tip of crack ③ gradually expands to the lower left surface of the slope; at step 49, the expansion paths of the tips of cracks ② and ③ intersect; at step 63, the tip of crack ⑥ expands to the crack line where the tip of crack ⑦ is located; at step 100, the tip of crack ⑦ expands to the upper surface of the slope. At this point, all crack tips have completed their expansion, and the initially relatively complete slope has been cut into 10 discrete blocks by the expanded cracks. Fig. 9 The bottom row clearly shows the movement and collapse process of 10 discrete blocks from 200 steps to 500 steps. It can be seen that this application example fully demonstrates that the improved numerical manifold method proposed in this invention can obtain the entire progressive failure evolution process of jointed rock slope from continuous small deformation to discrete large movement.

[0173] The beneficial effects of the present invention are as follows: the present invention first adds the four trigonometric functions contained in the first two terms of the Williams displacement series near the crack tip to the global displacement approximation of the manifold unit where the crack tip is located to reflect the properties of the field variables at the crack tip; then proposes an enhanced generalized Duffy transform method to improve the numerical integration accuracy of singular functions; and then adopts a cone complementary contact force calculation model that can simultaneously consider friction and adhesion to accurately calculate the contact force during the progressive failure of the rock slope; finally, introduces the critical length cracking criterion to determine the direction and distance of crack propagation. Compared with the classic triangular mesh numerical manifold method, the improved method proposed by the present invention can obtain a more accurate crack propagation path and progressive failure process, which is mainly reflected in:

[0174] 1) It can reasonably describe the displacement discontinuity and stress singularity near the crack tip;

[0175] 2) The numerical integration accuracy of singular functions is higher;

[0176] 3) Improved the accuracy and stability of contact force calculation;

[0177] 4) It can simultaneously consider the tensile shear and compressive shear failure of rock mass.

[0178] The present invention realizes an improved triangular grid numerical manifold method, which can obtain a reasonable prediction of the progressive failure evolution process of jointed rock slopes, increases the accuracy and stability of numerical calculations, and can be adopted by relevant scientific research, commercial computer codes and programs.

[0179] In the description of the present invention, it is necessary to understand that the orientation or positional relationship indicated by the terms "center", "thickness", "upper", "lower", "horizontal", "top", "bottom", "inner", "outer", "radial", etc. is based on the orientation or positional relationship shown in the drawings, and is only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the device or element referred to must have a specific orientation, be constructed and operated in a specific orientation, and therefore cannot be understood as a limitation on the present invention. In addition, the terms "first", "second", and "third" are used only for descriptive purposes, and cannot be understood as indicating or implying the relative importance or the number of implicitly specified technical features. Therefore, the features defined by "first", "second", and "third" may explicitly or implicitly include one or more of the features.

Claims

1. An improved triangular grid numerical manifold method for simulating the progressive failure of jointed rock slopes is characterized by: The following steps are included: S1. Establish the initial geometric model of jointed rock slope; S2, establishing a triangular mesh discrete calculation model of the initial geometric model; S3, constructing the global displacement field of the singular manifold unit according to the triangular mesh discrete calculation model, and calculating the singular function integral by the enhanced generalized Duffy transform formula; S4, calculating the contact force according to the cone complementary contact force calculation model that considers both friction and adhesion; S5. Calculate the stress field near the tip of the crack according to the singular function integral and the contact force and determine the cracking condition of the cracks in the jointed rock slope in combination with the critical length cracking criterion; In S3, the expression of the enhanced generalized Duffy transform formula is specifically: Where I is the numerical integral of the singular function, and β is the minimum integer that ensures that 2β-1-αβ is a positive integer, which is used to simplify the integral term u 2β-1-αβ , α is a real number, is any function of the integral variables u and w, e w and e -w They are all exponential functions of the integral variable w; In S4, the expression of the cone complementary contact force calculation model considering both friction and adhesion is specifically: In the formula, p n is the normal contact force, p s and p t is the contact force in two mutually perpendicular tangential directions on the tangent plane, g n is the normal contact distance, g s and g t is the contact distance between two mutually perpendicular tangential directions on the tangent plane, is the equivalent friction coefficient, is the contact spacing vector, is the contact force vector, is the friction cone, It is the dual cone corresponding to the friction cone.

2. The improved triangular grid numerical manifold method for simulating progressive failure of jointed rock slopes according to claim 1 is characterized in that: The S2 is specifically: The initial geometric model of the jointed rock slope is covered with a three-node triangular mesh, and geometric operations of cutting and intersecting the three-node triangular mesh and the initial geometric model are performed to generate manifold units and obtain a triangular mesh discrete calculation model.

3. The improved triangular grid numerical manifold method for simulating progressive failure of jointed rock slope according to claim 1, characterized in that: In S3, the method for calculating the singular function integral is specifically: The singular manifold unit is triangulated with the crack tip as the vertex to generate several sub-triangles. The integral of each sub-triangle is calculated using the enhanced generalized Duffy transform formula, and the integrals of all sub-triangles are added together to obtain the singular function integral.

4. The improved triangular grid numerical manifold method for simulating progressive failure of jointed rock slopes according to claim 3 is characterized in that: In S3, the global displacement field of the singular manifold unit The expression is: Where the subscript S represents any singular manifold unit, is the undetermined degree of freedom in the x and y directions of the i-th node of the triangle of the singular manifold unit, and its specific expression is: In the formula, and are the displacements of node i in the x and y directions respectively; All are undetermined coefficients; is the node interpolation matrix of the singular manifold unit, and its specific expression is: In the formula, w i is the shape function of the i-th node of the triangle, In order to introduce the four trigonometric functions of the first two terms of the Williams displacement series near the crack tip, the specific expressions are: Where r and θ are the polar diameter and polar angle of the local polar coordinates with the crack tip as the pole and the crack direction as the polar axis.

5. The improved triangular grid numerical manifold method for simulating progressive failure of jointed rock slopes according to claim 1 is characterized in that: The S5 comprises the following sub-steps: S51. Calculate the principal stress at the crack tip, including the first to third principal stresses; S52, judging whether the first principal stress is greater than the material tensile stress, if so, proceeding to S53; if not, proceeding to S54; S53, finding the point closest to the tensile stress of the material in the direction perpendicular to the first principal stress, and calculating the first distance from the crack tip, judging whether the first distance is not less than the crack extension length in the direction of the crack surface during tensile failure, if so, the current crack tip is extended by the length of the first distance in the direction perpendicular to the maximum principal stress, and entering S55, if not, the current crack does not extend, and entering S55; S54, calculating the Mohr-Coulomb strength force according to the first to third principal stresses, determining whether the Mohr-Coulomb strength force is greater than 0, and if so, searching for the point whose Mohr-Coulomb strength force is closest to 0 along the first set direction, and calculating the second distance from the crack tip, searching for the point whose Mohr-Coulomb strength force is closest to 0 along the second set direction, and calculating the third distance from the crack tip, determining the cracking condition of the current crack according to the relationship between the second distance, the third distance and the crack extension length along the direction of the crack surface during tensile failure, and entering S55; S55. Based on the obtained cracking condition of the current crack, determine the cracking condition of the cracks contained in the jointed rock slope.

6. The improved triangular grid numerical manifold method for simulating progressive failure of jointed rock slopes according to claim 5, characterized in that: In S52, the first principal stress is specifically the maximum principal stress perpendicular to the crack surface near the crack tip, and its expression is specifically: In the formula, K I is the I-mode stress intensity factor, L T is the crack extension length along the direction of the crack surface during tensile failure, and its specific expression is: In the formula, K IC is the critical I-mode stress intensity factor for material instability failure, σ T The tensile stress of the material.

7. The improved triangular grid numerical manifold method for simulating progressive failure of jointed rock slopes according to claim 6, characterized in that: In S54, the expression for calculating the Mohr-Coulomb strength force is specifically: In the formula, c is the cohesion of the material, φ is the internal friction angle of the material, is the third principal stress; First set direction The specific expression is: In the formula, is the direction of minimum principal stress, is to set the stress direction, and =π / 4+ϕ / 2; Second setting direction The specific expression is: The method for determining the cracking condition of the current crack is specifically as follows: If max(L A , L B )<L T , then the current crack does not extend, where L A is the second distance, L B is the third distance; If L A ≥ L T > L B , the current crack extends along the first set direction to a second distance; If L B ≥ L T > L A , the current crack extends along the second set direction by a third distance; If L A ≥ L T And L B ≥ L T , the crack extends to the length of the second distance and the third distance respectively along the first set direction and the second set direction at the same time.

Images