Method and system for constructing stratified rock constitutive model

Abstract

The invention provides a stratified rock mass constitutive model construction method and system. The stratified rock mass constitutive model construction method comprises the steps that a basic hypothesis is set according to basic mechanical characteristics of a stratified rock mass; based on the basic assumption, a stratified rock mass constitutive model is constructed in the equivalent continuum frame according to the mechanical property of gradual deterioration of the stratified rock mass; the method further comprises the step of verifying the stratified rock mass constitutive model through theoretical solution comparison, numerical test effect analysis, indoor test result comparison and round tunnel excavation effect simulation. According to the embodiment of the invention, the contribution of the layer in the deformation and failure process of the rock mass is highlighted, the influence of the stress state change on the mechanical parameters of the rock mass is reflected, and the gradual degradation of the mechanical parameters along with the equivalent plastic strain in the deformation and failure process of the layered rock mass is realized; and the reasonability and the accuracy of the numerical calculation simulation result in the stratified rock mass engineering construction process are improved.

Description

Technical Field

[0001] The invention relates to the field of rock mechanics and engineering technology, and in particular to a method and system for constructing a layered rock mass constitutive model. Background Art

[0002] Layered rock mass is a widely distributed engineering rock mass, commonly found in water conservancy and hydropower, tunnel transportation, rock slope and other projects. Its mechanical properties and engineering impacts are very complex. Constructing a reasonable and effective mechanical model to describe the mechanical behavior of layered rock mass is the key to ensuring the safe construction of rock engineering. The mechanical model of layered rock mass should not only reflect the essential properties of layered rock mass, but also reflect the deformation characteristics, stress characteristics and characteristics of the distribution of the failure area of ​​the rock mass under different stress states. However, there is currently a lack of mechanical models that specifically reflect the essential properties of layered rock mass in the scope of equivalent continuous media. Traditional mechanical models generally only consider isotropic characteristics or only consider some characteristics of layered rock mass, such as deformation anisotropy, different failure mechanisms of rock blocks and layers, etc. Therefore, there are problems with the effectiveness and accuracy of the simulation results when conducting numerical simulation of layered rock mass engineering. Summary of the invention

[0003] The purpose of the present invention is to provide a method and system for constructing a layered rock mass constitutive model, aiming to solve the above-mentioned problems in the prior art.

[0004] An embodiment of the present invention provides a method for constructing a constitutive model of a layered rock mass, comprising:

[0005] Establish basic assumptions based on the fundamental mechanical characteristics of layered rock masses; and

[0006] Based on the basic assumptions, a constitutive model of layered rock mass is constructed in an equivalent continuous medium framework according to the mechanical properties of the layered rock mass with progressive degradation.

[0007] An embodiment of the present invention provides a system for constructing a layered rock mass constitutive model, comprising:

[0008] A foundation module that establishes basic assumptions based on the fundamental mechanical characteristics of layered rock masses; and

[0009] A construction module is used to construct a constitutive model of layered rock mass according to the mechanical properties of the layered rock mass with progressive degradation within an equivalent continuous medium framework based on the basic assumptions.

[0010] The use of the embodiments of the present invention can include the following beneficial effects: the embodiments of the present invention highlight the contribution of the layers in the deformation and failure process of the rock mass, reflect the influence of the stress state change on the mechanical parameters of the rock mass, realize the progressive deterioration of the mechanical parameters of the layered rock mass with the equivalent plastic strain during the deformation and failure process, and improve the rationality and accuracy of the numerical calculation simulation results carried out during the construction of the layered rock mass engineering. BRIEF DESCRIPTION OF THE DRAWINGS

[0011] In order to more clearly illustrate one or more embodiments of this specification or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments recorded in this specification. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative labor.

[0012] Figure 1 is a flow chart of a method for constructing a constitutive model of a layered rock mass according to an embodiment of the present invention;

[0013] Figure 2 is a schematic diagram of the coordinate system relationship of an embodiment of the present invention;

[0014] Figure 3 is a schematic diagram of joint groups and geometric parameters in the observation area of ​​an embodiment of the present invention;

[0015] Figure 4 is a schematic diagram of shear and tensile yield criteria and their dividing lines in an embodiment of the present invention;

[0016] Figure 5 is a schematic diagram of the positive domain and negative domain of the layer flow law of an embodiment of the present invention;

[0017] Figure 6 Schematic diagram of stress-strain curve and mechanical parameter evolution curve during rock mass damage process of an embodiment of the present invention (a is stress-strain curve, b is mechanical parameter evolution curve);

[0018] Figure 7 is a schematic diagram of a numerical model of a uniaxial compression specimen according to an embodiment of the present invention;

[0019] Figure 8 Schematic diagram of the relationship between the apparent elastic modulus and the layer inclination angle according to an embodiment of the present invention;

[0020] Fig. 9 Schematic diagram of the relationship between the uniaxial compressive strength of the strong layer rock and the layer inclination angle according to an embodiment of the present invention;

[0021] Fig.10 Schematic diagram of the relationship between the uniaxial compressive strength of weak-layer rock and the inclination angle of the layer in an embodiment of the present invention;

[0022] Fig.11 Schematic diagram of the geometric model for numerical calculation of the Brazilian splitting test according to an embodiment of the present invention;

[0023] Fig.12 Schematic diagram of Brazilian splitting failure mode of strong layer rock in an embodiment of the present invention;

[0024] Fig.13 Schematic diagram of the Brazilian splitting failure mode of weak-plane rock in an embodiment of the present invention;

[0025] Fig.14 Schematic diagram of axial stress-strain curves of rock at different confining pressures and different layer inclination angles in an embodiment of the present invention;

[0026] Fig.15 Schematic diagram of axial stress-strain curves of weak-layer rock under different confining pressures and different layer inclination angles according to an embodiment of the present invention;

[0027] Fig.16 Schematic diagram of the relationship between the triaxial compressive strength of an embodiment of the present invention and the inclination angle of the layer and the confining pressure;

[0028] Fig.17 Schematic diagram comparing numerical simulation results and indoor test results under uniaxial and triaxial compression conditions of an embodiment of the present invention (a and d represent stress-strain curves; b and e represent octahedral shear strain rates; c and f represent rock sample failure modes);

[0029] Fig.18 is a schematic diagram comparing the Brazilian splitting numerical test results and the indoor test results of an embodiment of the present invention;

[0030] Fig.19 is a schematic diagram of a unit failure mode under a Brazilian splitting stress state according to an embodiment of the present invention;

[0031] Fig. 20 Schematic diagram of the distribution relationship among the inclination angles LRFD, RFDm and RFDj of different layers under the uniaxial compression test conditions of an embodiment of the present invention;

[0032] Fig.21 is a schematic diagram of a circular tunnel calculation model and boundary conditions according to an embodiment of the present invention;

[0033] Fig. 22 Schematic diagram of maximum displacement values ​​(inclination angle 60°) obtained by calculation and simulation of different constitutive models according to an embodiment of the present invention;

[0034] Fig.23 Schematic diagram of calculating the vault, arch foot and maximum displacement value of a circular tunnel using the LRDM model according to an embodiment of the present invention;

[0035] Fig.24 is a schematic diagram of the displacement field distribution of the surrounding rock of a circular tunnel according to an embodiment of the present invention;

[0036] Fig.25 is a schematic diagram of the maximum principal stress field distribution of the surrounding rock of a circular tunnel according to an embodiment of the present invention;

[0037] Fig.26 Schematic diagram of the distribution of plastic zones of the surrounding rock of a circular tunnel according to an embodiment of the present invention;

[0038] Fig. 27 Schematic diagram of a system for constructing a layered rock mass constitutive model according to an embodiment of the present invention. DETAILED DESCRIPTION

[0039] In order to enable those skilled in the art to better understand the technical solutions in one or more embodiments of this specification, the following will be combined with the drawings in one or more embodiments of this specification to clearly and completely describe the technical solutions in one or more embodiments of this specification. Obviously, the described embodiments are only part of the embodiments of this specification, not all of the embodiments. Based on one or more embodiments of this specification, all other embodiments obtained by ordinary technicians in this field without creative work should fall within the scope of protection of this document.

[0040] Method Embodiment

[0041] According to an embodiment of the present invention, a method for constructing a constitutive model of a layered rock mass is provided. Figure 1 is a flow chart of a method for constructing a constitutive model of a layered rock mass according to an embodiment of the present invention. Figure 1 As shown, the method for constructing a constitutive model of a layered rock mass according to an embodiment of the present invention specifically includes:

[0042] Step S101, establishing basic assumptions based on basic mechanical characteristics of layered rock mass;

[0043] The basic assumptions include: the layered rock mass is composed of rock blocks and layers, the layer is the only weak surface in the rock mass, and the rock layer extends infinitely;

[0044] The layered rock mass is considered as a transversely isotropic elastic medium before and after failure along the layer surface;

[0045] The occurrence of rock formations can uniquely determine the spatial distribution characteristics of layered rock masses;

[0046] Under low and medium stress conditions, rock formation strength has nothing to do with the intermediate principal stress;

[0047] The attenuation of the dilatancy angle of the layered rock mass is not considered, and only the degradation of the elastic modulus, cohesion and internal friction angle is considered;

[0048] The changes in rock physical and mechanical parameters caused by confining pressure are not considered.

[0049] Step S102, based on the basic assumptions, constructing a constitutive model of the layered rock mass in an equivalent continuous medium framework according to the mechanical properties of the layered rock mass with progressive degradation, specifically includes:

[0050] S41: Incorporating the joint deformation characteristics into the Oda fabric tensor to obtain the joint flexibility tensor, and constructing the equivalent elastic flexibility matrix of the layered rock mass in the local coordinate system through the joint flexibility tensor, specifically including:

[0051] The joint deformation characteristics are integrated into the Oda fabric tensor in the form of deformation stiffness coefficients to obtain the joint flexibility tensor as shown in Formula 1, the eigenvalues ​​and eigenvectors of the joint flexibility tensor are calculated, and the equivalent elastic flexibility matrix of the layered rock mass in the local coordinate system as shown in Formula 2 is constructed based on the eigenvalues ​​and eigenvectors; wherein the equivalent elastic flexibility matrix of the layered rock mass after the dominant joint group is weakened in the local coordinate system is shown in Formula 3;

[0052]

[0053] Where JCT represents the joint flexibility tensor, E(n,JF,r) is the joint probability density function of the joint normal vector, joint stiffness and joint geometric dimensions; JF is the joint deformation stiffness coefficient, r is the statistical average diameter of the joint crack, n is the unit direction vector of the joint crack, ρ v is the volume density of the joint fissure, E′1 is the Young's modulus in the transverse isotropic surface, E′3 is the Young's modulus perpendicular to the transverse isotropic surface, and v 12 is the Poisson’s ratio in the transversely isotropic plane, v 13 is the Poisson’s ratio in the perpendicular transverse isotropic plane, G′ 12 is the shear modulus in the transverse isotropic plane, G′ 13 is the shear modulus perpendicular to the transverse isotropic plane, E1 is the equivalent Young's modulus in the transverse isotropic plane, E3 is the equivalent Young's modulus perpendicular to the transverse isotropic plane, G 13 is the equivalent shear modulus perpendicular to the transverse isotropic plane, G 12 is the equivalent shear modulus in the transverse isotropic plane, D' R is the flexibility matrix of the rock mass, D' J is the local flexibility matrix in the main feature direction, K n0 and K s0 are the benchmark joint normal stiffness and joint shear stiffness respectively.

[0054] S42: Establish the transformation relationship between the elastic stiffness matrix in the global coordinate system, the local coordinate system and the principal stress space coordinate system;

[0055] S43: Determine strength criteria, plastic potential functions, and flow laws applicable to rock masses and layers, including:

[0056] The shear failure of rock blocks and layers adopts the Mohr-Coulomb yield failure criterion, and the tensile failure of rock blocks and layers adopts the Rankine maximum tensile stress yield failure criterion. The plastic potential function of shear failure of rock blocks and layers obeys the non-associated flow law, and the plastic potential function of tensile failure of rock blocks and layers obeys the associated flow law.

[0057] When the rock block is sheared, the principal stress of the rock block after shear plasticity correction is as shown in Formula 4; when the rock block is tensilely damaged, the principal stress of the rock block after tensile plasticity correction is as shown in Formula 5; when the layer is sheared, the stress of the layer after shear plasticity correction is as shown in Formula 6, and the stress increment after shear plasticity correction is as shown in Formula 7; when the layer is tensilely damaged, the stress of the layer after tensile plasticity correction is as shown in Formula 8, and the stress increment after tensile plasticity correction is as shown in Formula 9;

[0058]

[0059]

[0060] in, They are respectively the current maximum principal stress, intermediate principal stress and minimum principal stress when calculating according to elastic stress. are the maximum principal stress, intermediate principal stress and minimum principal stress after shear plasticity correction, They are the maximum principal stress, intermediate principal stress and minimum principal stress after tensile plasticity correction respectively; is the plasticity multiplier in different cases, P ij The corresponding component in the elastic stiffness matrix P in the principal stress space coordinate system, ψ m is the rock dilatancy angle, φ is the internal friction angle; τ O are the stress components in the local coordinate system when calculating according to elastic stress; τ N are the stress components in the local coordinate system after the correction of the plane shear plasticity, C 13 , C 23 , C 33 is the corresponding component in the elastic stiffness matrix C in the local coordinate system, ψ j is the dilatancy angle of the layer; Δσ 1′1′ , Δσ 2′2′ , Δσ 3′3′ , Δσ 1′3′ , Δσ 2′3′ is the stress increment after shear failure of the surface in the local coordinate system, is the component corresponding to the normal stress of the plane when calculating according to elastic stress; are the stress components in the local coordinate system after the tensile plasticity correction of the layer, is the stress increment after tensile failure of the layer in the local coordinate system, is the tensile strength of the layer.

[0061] S44: Using the corresponding internal variables and mechanical parameter evolution laws, the elastic modulus, cohesion and internal friction angle are taken as functions of equivalent plastic strain, including:

[0062] The cohesion and internal friction angle of the layer are calculated according to formula 10 and formula 11;

[0063]

[0064] Among them, c j0 ,φ j0 are the reference values ​​of the cohesion and internal friction angle of the layer, c jN ,φ jN are the current values ​​of the layer cohesion and internal friction angle, c jr ,φ jr are the residual values ​​of the interfacial cohesion and internal friction angle, are the equivalent plastic strain values ​​corresponding to the residual values ​​of the interfacial cohesion and internal friction angle, κ j is the current equivalent plastic strain value of the layer.

[0065] S45: Establish evaluation indicators for the progressive degradation process and degree of layered rock mass, including:

[0066] Formula 12 is used to calculate the rock layer damage index, and the evaluation index of the progressive degradation process and degree of the layered rock mass as shown in Formula 13 is established based on the rock layer damage index;

[0067]

[0068] LRFD = max(RFD m , RFD j ) Formula 13;

[0069] Among them, g i (θ) is the plastic potential function of the layer, and are the equivalent plastic strain on the plane and its corresponding limit equivalent plastic strain, F j Represents the yield function of the layer, RFD j Represents the rock layer damage index, RFD mrepresents the surrounding rock damage index, LRFD represents the evaluation index of the progressive degradation process and degree of layered rock mass, p and q represent the average stress and equivalent shear stress respectively, A j , B j , C j All are the layer strength criterion coefficients.

[0070] The method further comprises:

[0071] The constitutive model of layered rock mass is verified by comparing theoretical solutions, analyzing numerical test results, comparing indoor test results, and simulating circular tunnel excavation effects.

[0072] The above technical solution of the embodiment of the present invention is described in detail below in combination with the specific situation of the method for constructing a constitutive model of a layered rock mass in the embodiment of the present invention.

[0073] The layered rock constitutive model of the embodiment of the present invention discloses a numerical calculation method for the mechanical response of layered rock based on the mechanical properties of progressive degradation of the rock. The method includes: proposing 6 basic assumptions, establishing the conversion relationship between the elastic stiffness matrix in the global coordinate system, the local coordinate system and the principal stress space coordinate system; integrating the joint deformation characteristics into the Oda fabric tensor in the form of the deformation stiffness coefficient JF and constructing the joint flexibility tensor; adopting the strength criterion, plastic potential function and flow law applicable to rock blocks and layers; adopting appropriate internal variables and the evolution law of mechanical parameters, and considering the three mechanical parameters of elastic modulus, cohesion and internal friction angle as functions of equivalent plastic strain; establishing the degradation process and degree evaluation index of the layered rock, highlighting the control role of the layer in the process of rock damage and fracture; considering the mechanical properties of the progressive degradation of the layered rock, constructing the mechanical model of the layered rock in the framework of the equivalent continuous medium; verifying the validity of the layered rock constitutive model based on the comparison of theoretical solutions, the analysis of numerical test results, the comparison of indoor test results and the simulation of circular tunnel excavation effect.

[0074] Embodiments of the present invention include:

[0075] 1. By proposing six basic assumptions, the layered rock mass is regarded as a transversely isotropic medium;

[0076] According to the basic mechanical characteristics of layered rock mass, six basic assumptions are proposed as the premise for establishing the mechanical model of layered rock mass. The basic assumptions are as follows: (1) Layered rock mass is composed of rock blocks and layers. The layer is the only weak surface in the rock mass and the rock layer is infinitely extended; (2) Layered rock mass is regarded as a transversely isotropic elastic medium before and after destruction along the layer; (3) The rock layer dip can uniquely determine the spatial distribution characteristics of the layered rock mass; (4) Under medium and low stress conditions, the strength of the rock layer is independent of the intermediate principal stress; (5) The attenuation of the shear dilation angle of the layered rock mass is not considered, and only the deterioration of the elastic modulus E, cohesion c and internal friction angle φ is considered; (6) The changes in the physical and mechanical parameters of the rock caused by the confining pressure are not considered.

[0077] 2. Integrate the joint deformation characteristics into the Oda fabric tensor in the form of deformation stiffness coefficient JF to construct the joint flexibility tensor; specifically, integrate the joint deformation characteristics into the Oda fabric tensor in the form of deformation stiffness coefficient JF and propose the concept of joint flexibility tensor, and construct the equivalent elastic flexibility matrix D′ of the layered rock mass in the local coordinate system. Among them, the deformation stiffness coefficient JF can be calculated by formula (1).

[0078]

[0079] Among them, K n0 and K s0 are the benchmark joint normal stiffness and joint shear stiffness, K n and K s are the normal stiffness and shear stiffness of the joint, respectively.

[0080] It involves the equivalent elastic flexibility matrix of the layered rock mass in the local coordinate system. The matrix adopts the form of the equivalent elastic flexibility matrix of the rock mass cut by three sets of mutually orthogonal structural planes, which is specifically shown in formula (2):

[0081]

[0082] Where E′1 is the Young's modulus in the transverse isotropic surface, E′3 is the Young's modulus perpendicular to the transverse isotropic surface, and v 12 is the Poisson’s ratio in the transversely isotropic plane, v 13 is the Poisson’s ratio in the perpendicular transverse isotropic plane, G′ 12 is the shear modulus in the transverse isotropic plane, G′ 13 is the shear modulus perpendicular to the transverse isotropic plane.

[0083] 3. Establish the conversion relationship between the elastic stiffness matrices in the global coordinate system, local coordinate system and principal stress space coordinate system; specifically, it includes: establishing the conversion relationship between the elastic stiffness matrix A in the global coordinate system, the elastic stiffness matrix C in the local coordinate system and the elastic stiffness matrix P in the principal stress space coordinate system. That is, the deformation parameter matrices involved are the elastic stiffness matrix A in the global coordinate system, the elastic stiffness matrix C in the local coordinate system and the elastic stiffness matrix P in the principal stress space coordinate system. The above three deformation parameter matrices are not the same concept as the joint deformation characteristics and the equivalent elastic flexibility matrix D' of the layered rock mass involved. Involving three coordinate systems, such as Figure 2 As shown. The first is the global coordinate system, which belongs to the right-hand coordinate system, the positive direction of the X-axis is due east, the positive direction of the Y-axis is due north, and the positive direction of the Z-axis points upward; the second is the local coordinate system, the positive direction of the T-axis points to the direction of the layer, the positive direction of the N-axis points to the layer normal, the S-axis and the other two axes form a right-hand coordinate system, and the projection of the S-axis on the horizontal plane points to the layer inclination; the third is the principal stress space coordinate system, in which the positive direction of the P1 axis represents the direction of the maximum principal stress, the positive direction of the P2 axis represents the direction of the intermediate principal stress, and the positive direction of the P3 axis represents the direction of the minimum principal stress, and they together form a right-hand coordinate system. The transformation relationship between the second-order tensor in different coordinate systems obeys the tensor coordinate transformation law, as shown in formula (3):

[0084] a ij = l im l jn a mn (i,j,m,n=1,2,3) (3)

[0085] Among them, a ij represents the transformed coordinate system, a mn represents the coordinate system before transformation, l im , l jn They are all direction cosines of the coordinate system before transformation in the coordinate system after transformation. i, j, m, and n vary from 1 to 3. i and j are free indicators, and m and n are dummy indicators.

[0086] The coordinate transformation matrix R between the global coordinate system and the local coordinate system can be calculated by formula (4):

[0087]

[0088] Where l1, l2, l3, m1, m2, m3, n1, n2, n3 are the direction cosine values ​​of the three coordinate axes of the local coordinate system in the global coordinate system, which can be calculated by equation (5). The subscript numbers 1, 2, and 3 correspond to the coordinates of the T axis, N axis, and S axis in the local coordinate system, respectively.

[0089]

[0090] Where dip is the dip angle of the layer, dd is the inclination of the layer, then the direction cosine vectors l, m, n of the coordinate axes of the local coordinate system in the global coordinate system can be expressed as:

[0091]

[0092] Assume that the direction cosine vectors u, v, and w of the coordinate axes of the principal stress space coordinate system in the global coordinate system are as shown in formula (7):

[0093]

[0094] Where u1, u2, u3, v1, v2, v3, w1, w2, and w3 are the direction cosine components of the three coordinate axes of the principal stress space coordinate system in the global coordinate system.

[0095] The coordinate transformation matrix Q between the principal stress space coordinate system and the local coordinate system can be calculated according to formula (8):

[0096]

[0098] In the formula, l ′ 1. l ′ 2. l ′ 3. m ′ 1. m ′ 2. m ′ 3.n ′ 1.n ′ 2.n ′ 3 is the direction cosine component of the coordinate axis of the local coordinate system in the principal stress space coordinate system, which can be obtained by the following formula (9), and the subscript numbers 1, 2, and 3 correspond to the maximum principal stress direction axis, the intermediate principal stress direction axis, and the minimum principal stress direction axis in the principal stress space coordinate system, respectively.

[0099]

[0100] Wherein, l, m, n are the direction cosine vectors of the coordinate axes of the local coordinate system in the global coordinate system; u, v, w are the direction cosine vectors of the coordinate axes of the principal stress space coordinate system in the global coordinate system; l1, l2, l3, m1, m2, m3, n1, n2, n3 are the direction cosine values ​​of the three coordinate axes of the local coordinate system in the global coordinate system; u1, u2, u3, v1, v2, v3, w1, w2, w3 are the direction cosine components of the three coordinate axes of the principal stress space coordinate system in the global coordinate system respectively.

[0101] 4. Adopt strength criteria, plastic potential function and flow law applicable to rock blocks and layers; adopt appropriate internal variables and evolution laws of mechanical parameters; regard the three mechanical parameters of elastic modulus, cohesion and internal friction angle as functions of equivalent plastic strain; specifically include: the shear failure of rock blocks and layers adopts the Moore-Coulomb yield failure criterion to determine whether shear yield occurs; the tensile failure of rock blocks and layers adopts the Rankine maximum tensile stress yield failure criterion to determine whether tensile yield occurs; the plastic potential function of compressive shear failure of rock blocks and layers obeys the non-associated flow law, which is used to determine the direction of the shear plastic strain increment of rock blocks and layers; and the plastic potential function of tensile failure of rock blocks and layers obeys the associated flow law, which is used to determine the direction of the tensile plastic strain increment of rock blocks and layers. The user-defined equivalent plastic strain in the finite element software FLAC3D is used as the internal variable to describe and evaluate the plastic deformation behavior of layered rock materials. Four hardening parameters are defined to reflect the changes in the internal structure and performance of layered rock materials after plastic deformation. The evolution law of the three mechanical parameters of elastic modulus, cohesion and internal friction angle of layered rock during the damage yield fracture process is regarded as a function of the equivalent plastic strain. The hardening parameters are obtained by integrating the equivalent plastic strain and used to calculate the degradation values ​​of elastic modulus, cohesion and internal friction angle.

[0102] 5. Establish evaluation indicators for the degradation process and degree of layered rock mass, highlighting the controlling role of the layer in the process of rock mass damage and fracture; consider the mechanical properties of the progressive degradation of layered rock mass, construct a layered rock mass mechanical model within the framework of equivalent continuous medium, and conduct verification and analysis on it; specifically, based on the surrounding rock destruction index RFDm, an evaluation index LRFD for the degradation process and degree of layered rock mass is proposed, and the LRFD index value is calculated by equivalent plastic strain. This index incorporates the degradation process and degree of the rock layer into the evaluation of the rock mass fracture degree; using Microsoft Visual Studio 2010 software as the development platform, the constitutive characteristics of the proposed layered rock mass mechanical model are compiled through the supporting development interface provided by the existing software platform, thereby realizing the program development of the mechanical model, which can be used for numerical simulation calculations.

[0103] The surrounding rock damage index RFDm is involved, which can be calculated by formula (10):

[0104]

[0105] in:

[0106]

[0107] Where F represents the strength failure criterion of the surrounding rock. Before the peak failure (F<0), g(θ) is the plastic potential function, θ is the Lode angle, p and q are the average stress and equivalent shear stress, respectively, and A, B, and C are the strength criterion coefficients of the rock block. For the post-peak failure stage (F≥0), and They are the current equivalent plastic strain of the rock material and its corresponding limit equivalent plastic strain. The RFDm index at the peak failure stage actually reflects the damage development induced by the volume and shape changes of the rock material. The value range of the RFDm evaluation index is [0, 2]. When RFDm<1, it reflects the degree of damage development before the surrounding rock peak; when RFDm=1, it indicates the degree of damage development when the surrounding rock is at the peak; when RFDm>1, it reflects the degree of damage development after the surrounding rock peak; when RFDm>2, the surrounding rock stress state has reached the residual stress level.

[0108] 6. Verify the rationality and accuracy of the constitutive model of layered rock mass based on comparison of theoretical solutions, analysis of numerical test results, comparison of indoor test results, and simulation of circular tunnel excavation effects. Specifically, compare the differences between theoretical solutions and numerical solutions, use numerical simulation to reproduce the effects of indoor tests, establish a circular tunnel geometry model, analyze the simulation effects of the proposed model and other mechanical models in tunnel surrounding rock excavation, and prove the rationality and effectiveness of the mechanical model of layered rock mass.

[0109] Specifically, the embodiment of the present invention includes the following steps:

[0110] 1. Establish the joint flexibility tensor

[0111] The joint deformation characteristics are integrated into the Oda fabric tensor as deformation stiffness coefficient, thus proposing the concept of joint flexibility tensor. The expression of Oda fabric tensor (fabric tensor) is shown in equation (12).

[0112]

[0113] In the formula, ρ v is the volume density of the joints and fissures, which can be expressed by m (v) / V is calculated, where m (v) is the number of joints and fissures whose center points are located in the volume V, r is the statistical average diameter of the joints and fissures, n is the unit direction vector of the joints and fissures, E(n,r) is the joint probability density function of n and r, and Ω is the solid angle of the unit sphere surface (Ω=4π).

[0114] If the information of each joint and crack is known, the Oda fabric tensor can be simplified to equation (13).

[0115]

[0116] In the formula, and r(k) are the volume density and equivalent diameter of the kth group of joints, respectively.

[0117] The Oda fabric tensor can well characterize the geometric characteristics of the spatial distribution of joints in engineering rock mass, but its shortcoming is that it fails to consider the deformation characteristics of the joints themselves. The geometric distribution characteristics of joints and the deformation characteristics of joints themselves have a significant impact on the mechanical properties and deformation anisotropy of engineering rock mass. Therefore, the joint stiffness information is introduced into the fabric tensor to construct the joint flexibility tensor, which can be calculated by formula (14).

[0118]

[0119] where E(n,JF,r) is the joint probability density function of the joint normal vector, joint stiffness and joint geometric dimensions; JF is the joint deformation stiffness coefficient, which is defined by equation (1) above.

[0120] Assuming that the joint distribution geometric parameters and the joint deformation stiffness coefficient are independent of each other, the joint flexibility tensor can be expressed by equation (15).

[0121]

[0122] Where E(n) represents the probability density function of the unit normal vector n of the joint crack, E(JF) represents the probability density function of the deformation stiffness coefficient JF, and E(r) represents the probability density function of the statistical average diameter of the joint crack.

[0123] If the information of each joint inside the observation area can be obtained, the joint flexibility tensor JCT can be obtained by formula (16):

[0124]

[0125] In the formula, the superscript k represents the sequence number of the joint, JF (k) and n (k) denote the deformation stiffness coefficient and normal vector of the kth joint, s (k) represents the joint area, which can be expressed by π(r(k)) 2 / 4 is calculated.

[0126] The joints in engineering rock mass usually appear in groups or clusters regularly and have a certain combination pattern. A joint group is generally composed of joints with basically the same occurrence and basically the same mechanical properties. Its geometric distribution characteristics are usually reflected by the joint spacing, connectivity and orientation. If the joints in the study area can be approximately characterized by several joint groups, the joint flexibility tensor can be expressed by equation (17).

[0127]

[0128] In the formula, m represents the total number of dominant joint groups, superscript i represents the sequence number of dominant joint groups, and N i 、s i , JF i and n (i) They represent the number of joints, joint area, joint deformation stiffness coefficient and joint orientation vector in the i-th group of joints respectively.

[0129] Assuming the observation area is cylindrical, L and D * Respectively represent the length and diameter of the observation area, such as Figure 3 As shown, the volume in the study area can be expressed by formula (18).

[0130]

[0131] The average area of ​​the i-th group of joints can be expressed as:

[0132]

[0133] Where, d i represents the average diameter of the i-th group of joints.

[0134] Substituting equations (18) and (19) into equation (17) at the same time, the joint flexibility tensor can be expressed as equation (20):

[0135]

[0136] Assume that the angle between the observation line and the normal vector of the joint surface is α i , then we can get:

[0137]

[0138] Then formula (20) can be expressed as formula (22):

[0139]

[0140] In the formula, the first term is the joint group density, joint spacing The relationship between and joint density is shown by formula (23):

[0141]

[0142] At this time, the joint flexibility tensor can be expressed by formula (24):

[0143]

[0144] In the formula, represents the joint spacing of the i-th group of joints.

[0145] Analogous to the equivalent elastic flexibility matrix of rock mass cut by three sets of mutually orthogonal structural planes established by Amadei and Goodman, the equivalent elastic flexibility matrix of layered rock mass in the local coordinate system can be constructed as shown in the following formula (25). First, the joint flexibility tensor in the study area can be obtained based on the statistical data of the joint information in the study area, and the eigenvalue and eigenvector of the joint flexibility tensor can be calculated; then the joint flexibility tensor can be expressed as the local flexibility matrix D' in the main characteristic direction. J , and then based on the deformation superposition of rock blocks and joints, the equivalent elastic flexibility matrix D' of the layered rock mass in the local coordinate system can be deduced. T1, T2, T3 are the eigenvalues ​​of the joint flexibility tensor, D' R is the flexibility matrix of the rock block, where the specific form of D' is shown in formula (2).

[0146] D'=D' R +D' J (25)

[0147] If there are still several groups of dominant joints in the layered rock mass, we can deal with it as follows: the dominant joint group has a weakening effect on the deformation parameters of the layered rock mass, so the equivalent elastic modulus of the layered rock mass can be calculated by the joint flexibility tensor as a reference value for numerical simulation. Then we can derive equations (26) and (27):

[0148]

[0149] Let E1 = min(E a ,E b ), E3=E c , G 13 =min(G ac ,G ab ), then the equivalent elastic flexibility matrix D of the layered rock mass after the weakening of the dominant joint group in the local coordinate system can be obtained, as shown in formula (28):

[0150]

[0151] Where E1 is the equivalent Young's modulus in the transverse isotropic surface, E2 is the equivalent Young's modulus perpendicular to the transverse isotropic surface, and v 12 is the Poisson’s ratio in the transversely isotropic plane, v 13 is the Poisson’s ratio in the perpendicular transverse isotropic plane, G 13 is the equivalent shear modulus perpendicular to the transverse isotropic plane, G 12 is the equivalent shear modulus in the transverse isotropic plane,

[0152] 2. Establish elastic behavior relationship based on the transverse isotropic characteristics of layered rock mass

[0153] Generally speaking, before the layered rock mass fails, it can be equivalent to a transversely isotropic elastic material, and the relationship between the stress increment and the strain increment in the global coordinate system can be expressed by formula (29).

[0154]

[0155] Where A is the elastic stiffness matrix in the global coordinate system, which can be calculated according to formula (30).

[0156] [A]=[R][C][R] T (30)

[0157] In the formula, the elastic stiffness matrix C in the local coordinate system can be obtained by inverting the elastic flexibility matrix D in the local coordinate system, which can be obtained according to formula (31), and the coordinate transformation matrix R can be calculated according to formula (4).

[0158]

[0159] 3. Constructing Strength Criteria for Layered Rock Mass

[0160] In the field of geotechnical engineering, since the Mohr-Coulomb criterion and the Rankine maximum tensile stress criterion are widely used by many scholars and their parameters are relatively easy to obtain, for the layered rock mechanics model of the embodiment of the present invention, in order to better describe the strength characteristics of the transversely isotropic rock material, the shear failure of the rock block and the layer adopts the Mohr-Coulomb yield failure criterion to determine whether shear yield occurs, and the tensile failure of the rock block and the layer adopts the Rankine maximum tensile stress yield failure criterion to determine whether tensile yield occurs. The yield function of the rock block is composed of the shear failure yield function and the tensile failure yield function The composition is shown in formula (32) and formula (33).

[0161]

[0162] In the formula, σ1 and σ3 are the first and third principal stresses respectively, c m ,φ m are the cohesion and internal friction angle of the rock block respectively. is the tensile strength of the rock block, and its maximum value

[0163] Shear failure yield function and the tensile failure yield function The composition is shown in formula (34) and formula (35).

[0164]

[0165] In the formula, c j ,φ j are the cohesion and internal friction angle of the layer, σ 1′3′ , σ 2′3′ is the shear stress in the direction parallel to the plane, σ 3′3′ is the tensile stress in the direction perpendicular to the plane, is the tensile strength of the layer.

[0166] 4. Use plastic potential functions and flow laws applicable to rock blocks

[0167] The plastic potential function of the rock block is expressed by two functions and Description, shear plastic flow and extensional plastic flow are defined respectively, and They are used to determine the directions of the shear plastic strain increment and the tensile plastic strain increment of the rock block. Function It satisfies the non-associated flow law and can be expressed by formula (36).

[0168]

[0169] In the formula, ψ m is the dilatancy angle of the rock block,

[0170] function It conforms to the law of associated flow and can be expressed by formula (37).

[0171]

[0172] Because the shear and tensile composite yield failure criterion is used, it is necessary to determine the dividing line between shear failure and tensile failure under different stress states. Here, the bisector of the angle between the two strength envelopes of shear failure and tensile failure is used as the dividing line, such as Figure 4 shown.

[0173] The coordinates of the intersection of the Moore-Coulomb criterion and the Rankine criterion can be expressed by formula (38).

[0174]

[0175] Where c is the cohesion, φ is the internal friction angle, σ t is the tensile strength,

[0176] The slope of the angle bisector equation can be calculated using equation (39).

[0177]

[0178] Then the angle bisector can be expressed by equation (40).

[0179]

[0180] ①First, consider the shear failure of rock blocks and transform the plastic potential function Taking partial derivatives of the three principal stresses σ1, σ2, and σ3, we get:

[0181]

[0182] In the principal stress space coordinate system, the Hooke's law increment represented by the generalized stress and stress increment can be expressed by equation (42).

[0183]

[0184] Where P ij The corresponding component in the elastic stiffness matrix P in the principal stress space coordinate system, are the elastic strain increments in the directions of maximum principal stress, intermediate principal stress and minimum principal stress, respectively; Δσ1, Δσ2 and Δσ3 are the stress increments in the directions of maximum principal stress, intermediate principal stress and minimum principal stress, respectively.

[0185] For the convenience of subsequent reference, the above formula (42) is rewritten as formula (43).

[0186]

[0187] In formula (43), Replace with You can get:

[0188]

[0189] Where P 11 , P 21 , P 31 , P 13 , P 23 , P 33 The corresponding component in the elastic stiffness matrix P in the principal stress space coordinate system, ψ m is the dilatancy angle of the rock mass.

[0190] Then the principal stress after shear plasticity correction can be expressed by equation (45).

[0191]

[0192] In the formula, is the plasticity multiplier, They are respectively the current maximum principal stress, intermediate principal stress and minimum principal stress when calculating according to elastic stress; They are the maximum principal stress, intermediate principal stress and minimum principal stress after shear plasticity correction respectively; the principal stress corrected according to the consistency condition satisfies the shear yield function Then the plasticity multiplier can be obtained by equation (46).

[0193]

[0194] ② Considering the case of tensile failure, the plastic potential function Taking partial derivatives of the three principal stresses σ1, σ2, and σ3, we get:

[0195]

[0196] In formula (43), Replace with You can get:

[0197]

[0198] Where P 13 , P 23 , P 33 They respectively represent the corresponding components in the elastic stiffness matrix P in the principal stress space coordinate system.

[0199] Then the principal stress after tensile plasticity correction can be expressed by formula (49).

[0200]

[0201] In the formula, They are respectively the current maximum principal stress, intermediate principal stress and minimum principal stress when calculating according to elastic stress; They are the maximum principal stress, intermediate principal stress and minimum principal stress after tensile plasticity correction respectively; is the plasticity multiplier and can be calculated by formula (50).

[0202]

[0203] 5. Use plastic potential functions and flow laws applicable to the layer

[0204] The plastic potential function of the layer is a composite function that takes into account shear and tension ( and ), which are used to determine the direction of the shear plastic strain increment and the direction of the tensile plastic strain increment of the layer. Potential function Satisfying the non-associated flow law, the potential function Satisfying the associated flow law, it can be expressed by equation (51) and equation (52) respectively.

[0205]

[0206] In the formula, τ j is the shear stress on the surface, σ 3′3′ is the normal stress on the plane, ψ j is the shear dilatancy angle of the layer.

[0207] The following method can be used to determine the unique expression of the layer plastic potential function that satisfies the flow law. Function H(σ 3′3′ ,τ) is defined as (σ 3′3′ ,τ) plane and The angle bisector of Figure 5 As shown in (the principal stress is negative for compression and positive for tension). The function is selected based on the positive and negative domains of the function, as shown in formula (53).

[0208]

[0209] Where τ is the shear stress on the surface, and are two constants and can be calculated by equation (54) and equation (55).

[0210]

[0211] In the formula, φ j is the internal friction angle of the layer.

[0212] Assuming that a certain stress state breaks through the composite yield function, in (σ 3′3′ ,τ) plane indicates that the point is located inside domain 1 and domain 2, which represent the positive and negative domains of H = 0 respectively. If the stress point is located in domain 1, it means that the point is in shear yield state, and the stress point is located on the curve On the plastic potential function The derived flow law; if the stress point is in domain 2, it means that the point is in a tensile yield state, and the new stress point satisfies Applicable to plastic potential function Derived flow law.

[0213] ①First, consider the case of shear failure in the layer, and transform the plastic potential function For the local stress component σ 1′1′ , σ 2′2′ , σ 3′3′ , τ j Taking partial derivatives respectively, we get:

[0214]

[0215] The Hooke's law increment expressed by generalized stress and stress increment in the local coordinate system can be calculated by equation (57).

[0216]

[0217] In the formula, C ij represents the corresponding component in the elastic stiffness matrix C in the local coordinate system of the plane; Δγ e are the elastic strain increments in the local coordinate system; G is the shear modulus. For the convenience of subsequent reference, the above formula (57) is rewritten as formula (58).

[0218]

[0219]

[0220] In formula (58), Δγ e Replace with We can get formula (59).

[0221]

[0222] Then the stress after the correction of the layer shear plasticity can be expressed by formula (60).

[0223]

[0224] In the formula, τ O are the stress components in the local coordinate system when calculating according to elastic stress; τ N are the stress components in the local coordinate system after the correction of the plane shear plasticity; C 13 , C 23 , C 33 is the corresponding component in the elastic stiffness matrix C in the local coordinate system; ψ j is the layer dilatancy angle, is the plasticity multiplier and can be calculated by formula (61).

[0225]

[0226] From τ N and τ O The new shear stress component expression on the plane is derived as follows:

[0227]

[0228] In the formula, is the shear stress component after correction for shear plasticity; It is the component of normal stress in the plane when calculating according to elastic stress.

[0229] Therefore, the correction of plastic stress increment after layer shear failure can be expressed by formula (63).

[0230]

[0231] In the formula, Δσ 1′1′ , Δσ 2′2′ , Δσ 3′3′ , Δσ 1′3′ , Δσ 2′3′ is the stress increment after shear failure of the surface in the local coordinate system; is the layer plasticity multiplier; C 13 , C 23 , C 33 is the corresponding component in the elastic stiffness matrix C in the local coordinate system; ψ j is the layer dilatancy angle; It is the component corresponding to the normal stress of the plane when calculating according to elastic stress.

[0232] ② Now consider the case of tensile failure in the layer, and change the plastic potential function Taking partial derivatives of the local stress components, we get:

[0233]

[0234] In formula (58), Δγ e Replace with You can get:

[0235]

[0236] In the formula, C 13 , C 23 , C 33 is the corresponding component in the elastic stiffness matrix C in the local coordinate system.

[0237] Therefore, the stress after plastic correction of layer tensile failure can be calculated by formula (66).

[0238]

[0239] In the formula, are the stress components in the local coordinate system after the tensile plasticity correction of the layer; C is the stress components in the local coordinate system when calculating according to elastic stress; 13 , C 23 , C 33is the corresponding component in the elastic stiffness matrix C in the local coordinate system; is the tensile plasticity multiplier, which can be calculated by equation (67) according to the consistency condition.

[0240]

[0241] Therefore, the plastic stress correction increment after layer tensile failure can be calculated by formula (68).

[0242]

[0243] In the formula, is the stress increment after tensile failure of the layer in the local coordinate system; is the stress component in the local coordinate system when calculating according to elastic stress; C is the tensile strength of the layer; 13 , C 23 , C 33 is the corresponding component in the elastic stiffness matrix C in the local coordinate system.

[0244] 6. Use appropriate internal variables and mechanical parameter evolution laws

[0245] In the layered rock mass mechanical model of the embodiment of the present invention, the self-defined equivalent plastic strain in the FLAC3D software is used as an internal variable to describe and evaluate the plastic deformation behavior of the layered rock mass material, and four hardening parameters are defined: Reflects the changes in the internal structure and properties of layered rock materials after undergoing plastic deformation. represents the plastic shear strain of the rock mass; It represents the plastic volume tensile strain of the rock mass; represents the plastic shear strain of the layer; Represents the volumetric plastic tensile strain of the plane. These defining parameters are the sum of the plastic strain increments of the element. The element hardening increment is calculated as the volume average of all tetrahedral hardening increments within the element. For the current calculation step, the square root of the second invariant of the plastic shear-strain increment deviator tensor is one tetrahedral shear hardening increment.

[0246] For rocks:

[0247]

[0248] in:

[0249]

[0250] In the formula, are the equivalent plastic shear strain increments in the directions of maximum and minimum principal strain, respectively; is the volumetric plastic shear strain increment.

[0251] The levels are:

[0252]

[0253] In the formula, are the equivalent plastic shear strain increments of the plane in the local coordinate system.

[0254] The increment of tetrahedral tensile hardening is the increment of plastic tensile strain. For rock blocks, we have:

[0255]

[0256] In the formula, is the equivalent plastic tensile strain increment in the direction of minimum principal strain.

[0257] For the levels:

[0258]

[0259] In the formula, is the equivalent plastic tensile strain increment of the layer in the local coordinate system.

[0260] In the embodiment of the present invention, it is considered that the process of rock mass damage and microcrack development is a process of continuous separation of rock mass microstructure, that is, the rock mass cohesion (c) decreases with the equivalent plastic strain and the effective friction (φ) increases with the plastic strain. In the mechanical model of the embodiment of the present invention, the degree of plastic damage after rock mass yield is expressed by the equivalent plastic strain. The elastic modulus (E), cohesion (c) and internal friction angle (φ) in the process of rock mass yield fracture damage are considered as equivalent plastic strain. The sudden brittle yield fracture process of hard rock after unloading under high stress can be divided into four stages: elastic deformation, damage yield, brittle fracture and complete destruction. The typical stress-strain curve is shown in Figure 6 (a) In this process, the mechanical parameters of the rock mass evolve continuously from the damage yield stage, and the initial value and residual value of the mechanical parameters in the elastic deformation stage and the complete destruction stage respectively remain unchanged.

[0261]

[0262] In the formula, the mechanical parameters E0, c0, and φ0 are the initial values ​​of the rock mass before damage yield; and is the degradation value of mechanical parameters under certain equivalent plastic strain; and The equivalent plastic strain is A function of may be used in the form of a linear function, a high-order function, or a power function to reflect the evolution law of the mechanical parameters with the equivalent plastic strain. The embodiment of the present invention intends to use a linear function, as shown in the following equations (75) to (78):

[0263]

[0264]

[0265] In the formula, E m0 、c m0 and φ m0 are the reference values ​​of rock elastic modulus, cohesion and internal friction angle respectively; E mN 、c mN and φ mN are the current elastic modulus, cohesion and internal friction angle of the rock block; E mr 、c mr and φ mr They are the residual elastic modulus, residual cohesion and residual internal friction angle values ​​respectively; and are the equivalent plastic strain values ​​corresponding to the elastic modulus, cohesion and internal friction angle when they reach the residual value; κ m is the current value of the equivalent plastic strain of the rock block.

[0266] For the cohesion of the layer c j and internal friction angle φ j Then we have:

[0267]

[0268] In the formula, c j0 ,φ j0 are the reference values ​​of the cohesion and internal friction angle of the layer, respectively; c jN ,φ jN are the current values ​​of the layer cohesion and internal friction angle, respectively; c jr ,φ jr are the residual values ​​of layer cohesion and internal friction angle, respectively; are the equivalent plastic strain values ​​corresponding to the residual values ​​of the interfacial cohesion and internal friction angle; κ j is the current equivalent plastic strain value of the layer.

[0269] 7. Establish evaluation indexes for the degradation process and degree of layered rock mass

[0270] The surrounding rock damage index RFDm is proposed based on the theory of continuum mechanics and the stress and strain analysis of true triaxial tests of rock. It can be used to evaluate the degree of damage development at a local position in the rock. This index can intuitively and quantitatively characterize the location and degree of surrounding rock rupture. At the same time, it can reasonably unify the typical stage of the rock yield / rupture curve with the support optimization of the engineering surrounding rock, thereby realizing the organic combination of the plastic failure zone of numerical calculation and the depth of the loose zone measured on site and the visualization depth of the surrounding rock. Its calculation expression is shown in the above formula (10):

[0271]

[0272] in:

[0273]

[0274] Where F represents the yield function of the rock block. Before the peak failure (F<0), g(θ) is the plastic potential function, θ is the Lode angle, p and q are the average stress and equivalent shear stress, σ1, σ2, σ3 are the maximum principal stress, the intermediate principal stress and the minimum principal stress, respectively. A, B, and C are the strength criterion coefficients of the rock block. For the post-peak failure stage (F≥0), and They are the current equivalent plastic strain of the rock material and its corresponding limit equivalent plastic strain. The RFDm index at the peak failure stage actually reflects the damage development induced by the volume and shape changes of the rock material. The value range of the RFDm evaluation index is [0, 2]. When RFDm<1, it reflects the degree of damage development before the peak of the surrounding rock; when RFDm=1, it indicates the degree of damage development of the surrounding rock at the peak; when RFDm>1, it reflects the degree of damage development after the peak of the surrounding rock; when RFDm>2, the stress state of the surrounding rock has reached the residual stress level.

[0275] Similar to the definition of surrounding rock destruction index RFDm, for the layer yield failure of layered rock mass, the rock layer destruction index RFDj is defined, and its specific expression is as follows:

[0276]

[0277] In the formula, g i (θ) is the plastic potential function of the layer, and are the equivalent plastic strain on the level and its corresponding limit equivalent plastic strain, Fj represents the yield function of the level, A j , B j , C j All are the layer strength criterion coefficients.

[0278] Rock block failure and layer failure are the two main types of failure of layered rock mass. Because rock mass always fails along or through its weakest part, the failure of layered rock mass depends on rock failure or layer failure or both. Therefore, the larger value between RFDm and RFDj is used as the LRFD index of layered rock mass to evaluate the cracking process and degree of layered rock mass. Its expression is:

[0279] LRFD = max(RFD m , RFD j ) (82)

[0280] 8. Comparison of theoretical solutions

[0281] The uniaxial compression specimen numerical model was established in the FLAC3D software platform to compare and analyze the differences between the numerical calculation results and theoretical analysis of the layered rock mass mechanical model. The height of the calculation model is 100 mm, the diameter is 50 mm, the unit type is mainly tetrahedral unit, the number of nodes is 5341, and the total number of units is 25241. Figure 7 As shown, the cohesion c of the layer j and tensile strength Each of them sets a very large numerical simulation of the situation where the layer bonding is very strong, and this numerical sample is collectively referred to as strong layer rock, and the numerical samples with normal assignment are collectively referred to as weak layer rock. The Fish programming language is used to write a function to control the axial loading rate of the upper and lower end faces of the cylindrical sample. After the simulation is completed, the numerical solution of the apparent elastic modulus and uniaxial compressive strength of the numerical sample at different layer inclination angles is calculated, and its comparison with the analytical solution is as follows Figure 8 , Fig. 9 , Fig.10 shown.

[0282] According to the results of numerical simulation, it can be known that the error range between the numerical solution and the analytical solution of the apparent Young's modulus is 0.8% to 9.1%, with an average error of 5%. The analytical solution is smaller than the numerical solution, and the numerical values ​​of the two are basically consistent with the trend of the change of the layer inclination angle, and both show a U-shaped curve change law as a whole. Among them, the mathematical expression of the analytical solution is:

[0283]

[0284] Where α is the inclination angle of the layer; E1 is the Young's modulus in the transverse isotropic surface; E3 is the Young's modulus perpendicular to the transverse isotropic surface; v 13 is the Poisson’s ratio in the perpendicular transverse isotropic plane; G 13 is the shear modulus perpendicular to the transverse isotropic plane.

[0285] The analytical solution to the uniaxial compressive strength of strong-layer rock adopts the strength formula that describes the change of rock strength with loading angle (also known as TK strength criterion), and its mathematical expression is:

[0286]

[0287] Where c and φ are the cohesion and internal friction angle of the rock block, respectively; α is the angle between the normal vector of the plane and the direction of the maximum principal stress; n and k are two parameters, which are defined as functions of the rock deformation parameters and can be used as fitting parameter values. In this solution, k = 0.97 and n = 2.3 are taken.

[0288] According to the calculation results, the error range between the numerical solution and the analytical solution is between 1.94% and 26.14%, and both show a U-shaped variation trend with the inclination of the plane. The maximum error between the numerical solution and the analytical solution is relatively large: first, the analytical solution reflects the ultimate load stress of a certain unit grid; second, the numerical solution reflects the overall strength of many unit grids after macroscopic equivalence, and the rock material is regarded as a transversely isotropic medium. Therefore, under the action of the stress field, not all unit grids can reach the ultimate bearing stress, and some units cannot reach the peak stress intensity.

[0289] For the theoretical solution of the compressive strength of rock in weak layers, the composite form of the TK strength formula and the Jaeger single weak plane strength formula is used, and the strength of a single weak plane can be calculated by formula (85). When the TK strength formula controls the compressive strength of rock in weak layers, the analytical solution is basically greater than the numerical solution; when the weak layers control the compressive strength of rock (the layer inclination range is approximately 40 to 60°), the analytical solution and the numerical solution are basically the same, but when the layer inclination is 70° and 80°, the theoretical solution is smaller than the numerical solution. The change law of the uniaxial compressive strength of rock in weak layers with the inclination angle is generally characterized by a U-shaped curve, which is basically consistent with many experimental results.

[0290]

[0291] In the formula, σ jc represents the compressive strength of rock with a single weak plane, c j represents the cohesion of the weak surface, φ j represents the internal friction angle of the weak surface, and β represents the inclination angle of the weak surface.

[0292] 9. Brazilian splitting numerical test results

[0293] The computational geometry model for the Brazilian splitting numerical simulation test is as follows Fig.11As shown in the figure, for the convenience of loading, two upper and lower loading plates are set. The diameter of the calculation model is 50mm, the height is 25mm, the unit type is all tetrahedron, the number of nodes is 21606, the number of units is 108573, and the FISH language is used to write functions to control the axial movement of the loading plate and thus control the loading rate. Two cases of strong layer rock and weak layer rock are simulated respectively, and the relationship between their failure mode and layer inclination is shown in the figure below. Fig.12 and Fig.13 As shown in the figure. The failure mode of strong layered rock is rock block tensile failure, and there is no tensile failure of the layer. From 0° to 90°, the inclination angle is basically split along the center line of the disk, which is a typical Brazilian splitting failure mode. For weak layered rock, ideal splitting failure surfaces appear at 0° and 90° inclination angles, but at 0°, rock block tensile failure is dominant, while at 90°, layer tensile failure is dominant. This is completely different from strong layered rock. At 30° and 60° inclination angles, not only the main splitting failure surface deviating from the center line of the disk appears, but also the secondary fracture surface roughly parallel to the layer appears and is dominated by layer tensile failure. The main failure surface at 30° inclination angle is rock block tensile failure, and at 60° inclination angle, the main fracture surface is composed of three composite failures: layer tensile failure, rock block-layer chain failure, and rock block tensile failure.

[0294] 10. Results of conventional triaxial compression numerical tests

[0295] The numerical geometric model of the conventional triaxial compression test is the same as that of the uniaxial compression numerical test, simulating the triaxial compression tests of strong layer rock and weak layer rock under 10MPa, 20MPa, and 30MPa confining pressures, respectively. The axial stress-strain curves of strong layer rock and weak layer rock under different confining pressures and different layer inclinations are shown in Figure 2. Fig.14 and Fig.15 As shown in the figure. When the plane inclination angle is the same, the peak axial stress and residual axial stress increase with the increase of confining pressure, and the apparent elastic modulus of the simulated rock sample does not change significantly. This is because the mechanical model does not consider the influence of confining pressure on the deformation modulus. When the confining pressure is the same, the peak axial stress of the sample with a plane inclination angle of 90° is the largest, while that of the sample with a plane inclination angle of 60° is the smallest. The apparent elastic modulus and compressive strength of the samples with different plane inclination angles are different, indicating that the deformation modulus and compressive strength are affected by the plane inclination angle and have obvious anisotropic characteristics, which is consistent with the previous test results.

[0296] The relationship between the triaxial compressive strength of strong layer rock and weak layer rock and the layer inclination and confining pressure is as follows: Fig.16As shown. The overall variation trend of triaxial compressive strength with the inclination angle of the layer obtained by numerical simulation is basically consistent with the theoretical solution. The compressive strength calculated by the embodiment of the present invention is taken as the average stress of all units. For strong layer rock, the numerical solution of its triaxial compressive strength is smaller than the corresponding theoretical solution, which is due to the size effect (the model unit will not reach the peak strength at the same time and be damaged) and the influence of strain threshold. For weak layer rock, the numerical solution is smaller than the theoretical solution when the inclination angle of the layer is small (0°~40°), and the numerical solution is greater than the theoretical solution when the inclination angle is 50°~80°.

[0297] 11. Comparison of indoor test results

[0298] (I) Compression test

[0299] Uniaxial compression and conventional compression indoor tests were carried out on two rock samples with layer inclinations of 45° and 90°, respectively. At the same time, the corresponding numerical geometric model was established to carry out computational simulation.

[0300] Numerical and laboratory test stress-strain curves, shear strain rates, and failure modes of rock specimens Fig.17 As shown. Fig.17 (a)(d) It can be seen that the layered rock mechanical model constructed by the embodiment of the present invention is able to characterize the stress-strain behavior of layered rock under uniaxial compression and conventional triaxial compression stress states. The indoor test results are basically consistent with the overall change trend of the stress-strain curve obtained by numerical calculation simulation. The shortcoming is that it fails to reflect the mechanical behavior of the initial nonlinear compaction stage. The simulation results of the conventional triaxial test (20MPa confining pressure) are highly consistent with the test curve, especially the pre-peak yield and post-peak strain softening stages are basically consistent. The failure mode simulated by the numerical test generally has current shear and tensile failure along the layers and rock blocks, which is not easy to distinguish. Therefore, the current shear deformation characterized by the octahedral shear strain rate is used to reflect the current failure state instead of directly using the failure unit. By Fig.17 (b)(c) and (e)(f) show that the width range of the shear failure concentration zone obtained by numerical simulation results is larger than the width range of the shear slip surface obtained by indoor tests, while the azimuth inclination of the shear strain concentration is roughly the same as the azimuth inclination of the shear slip surface of the rock sample, which indicates that the restriction effect at the end of the sample is not significant.

[0301] (II) Brazilian splitting test

[0302] The diameter of the Brazilian splitting numerical test calculation model is 50mm and the height is 25mm. The Brazilian splitting numerical test results and the indoor test results are shown in the figure. Fig.18 As shown, the failure mode is Fig.19As shown. It can be seen that: when the layer inclination angle is different, the tensile failure zone starts from the upper and lower ends of the circular sample and then gradually expands to the center until it penetrates the entire sample; the tensile failure zone of the samples with 0° and 90° layer inclination angles is basically formed along the center line of the sample. The difference is that when the layer inclination angle is 0°, the rock block tensile failure is mainly caused, while when the layer inclination angle is 90°, the layer tensile failure is mainly caused. The tensile failure concentration area of ​​the 0° layer inclination sample is slightly larger than that of the 90° layer inclination; the rock block tensile failure zone of the 45° layer inclination sample has the largest range and has a certain small angle with the center line. At the same time, the sample also has a layer tensile failure zone with an angle roughly the same as the layer inclination angle, which shows that the existence of the layer affects the failure mode of the rock material to a certain extent. Compared with the results of the indoor test, the numerical test results at an inclination angle of 90° have the highest consistency, and the tensile failure zone is basically consistent with the splitting surface of the rock sample. At an inclination angle of 45°, the tensile failure zone of the numerical test forms an angle with the center line of the sample and has a certain range, while the fracture surface of the indoor test presents an arc surface. The numerical simulation results at an inclination angle of 0° have a poor consistency with the indoor test results. The numerical test can obtain the typical failure morphology of Brazilian splitting, while the failure surface obtained in this Brazilian test does not have typical splitting characteristics. This may be because the stress field of the layered rock sample deviates from the stress analytical solution of isotropic homogeneous rock under the influence of the directional arrangement of rock minerals. The above results show that the proposed model can better reproduce the central fracture of Brazilian splitting, thereby verifying the effectiveness of the mechanical model under indirect tensile stress state.

[0303] 12. Analysis of evaluation indexes of layered rock mass degradation process and degree

[0304] The numerical simulation results of uniaxial compression test are used to analyze the degradation process and degree evaluation index (LRFD) of layered rock mass. The distribution relationship of LRFD, RFDm and RFDj of samples with different layer inclination angles under uniaxial compression test conditions is shown in the figure. Fig. 20As shown in the figure, it can be seen that when the inclination angle increases from 0° to 90°, the LRFD index can well reflect the transformation process of the rock sample failure mode, that is, the inclination angle is 0°, the rock block is mainly controlled, the inclination angle is 30°, the rock block and the plane are jointly controlled, the inclination angle is 60°, the plane is mainly controlled, and the inclination angle is 90°, the rock block is controlled, but the vertical plane aggravates the failure. The RFDm index can still well characterize the failure area of ​​the rock material when the inclination angle is 0°, but when there are other inclination angles, RFDm cannot well reflect the failure mode and fracture degree of the sample, and RFDj can just make up for the shortcomings of RFDm. When there are certain plane inclination angles (such as 30° and 60°), FRDj can well reflect the failure mode of the rock sample, such as shear slip failure at an inclination angle of 60°, and the failure is mainly controlled by the rock block at an inclination angle of 90°, but the existence of the plane can aggravate the fracture of the sample. In summary, based on the original RFDm index, the LRFD index incorporates the evaluation index RFDj for evaluating the layer damage process and degree, and combines the two into one evaluation index, which can be used to reasonably evaluate the degradation process and degree of layered rock mass.

[0305] 13. Simulation of circular tunnel excavation effect

[0306] In order to analyze the simulation effect of the proposed layered rock mass mechanical model (LRDM) and other mechanical models on tunnel surrounding rock excavation, a numerical calculation model for circular tunnel excavation simulation was established, with a total of 16,000 units, such as Fig.21 As shown. Among them, the model size is 20m (X) × 0.1m (Y) × 20m (Z), the diameter of the circular tunnel is 2m, the distance from the center of the tunnel to the boundary of the calculation model is much greater than 5 times the radius, and the model thickness is 0.1m, so it can be regarded as a plane strain analysis. Stress boundary conditions are used at the top and a vertical stress of 10MPa is applied. Fixed boundary conditions are used in the horizontal direction, the longitudinal axis direction and the bottom of the model to constrain its normal displacement and vertical displacement. The deadweight of the rock mass within the calculation range of the model is not considered during the analysis. The transversely isotropic ubiquitous joint model (AEUbi), the Mohr Coulomb ideal elastoplastic model (MC), the strain softening Mohr Coulomb elastoplastic model (SSMC), the ubiquitous joint model (Ubi), the strain softening ubiquitous joint model (SSUbi) and the model (LRDM) proposed in the embodiment of the present invention built-in in the FLCA3D software are used to simulate the circular tunnel excavation.

[0307] 1. Comparative analysis of surrounding rock deformation field

[0308] The maximum displacements calculated by different constitutive models are as follows: Fig. 22 As shown in the figure, the circular tunnel vault, arch foot and maximum displacement values ​​calculated by LRDM model are as follows Fig.23Comparison of the calculation results of surrounding rock displacement field of layered rock mass mechanics model (LRDM) and strain softening uniform joint model (SSUbi), uniform joint model (Ubi), and transversely isotropic uniform joint model (AEUbi) Fig.24 As shown by Fig. 22 , 23 as well as Fig.24 It can be seen that: ① Using the LRDM model for simulation calculation, when the layer inclination angle gradually increases from 0° to 90°, the displacement field of the surrounding rock of the circular tunnel undergoes a significant shift change, and the area with larger deformation gradually shifts from the top arch to the arch shoulder and arch foot. This shows that the deformation field distribution of the surrounding rock after the excavation of the circular tunnel is significantly affected by the rock formation dip; ② When the layer inclination angle is the same, on the one hand, the deformation field shape distribution of the LRDM model and the AEUbi model has a high similarity. This is because both the LRDM model and the AEUbi model consider the mechanical properties of the transversely isotropic elastic medium before the layered rock mass is destroyed, so the shape distribution of the displacement field has a good similarity. On the other hand, the deformation field shape distribution of the LRDM model is not very similar to that of the SSUbi model and the Ubi model. This is because although the SSUbi model and the Ubi model consider the influence of the layer joints, the layered rock mass is still considered as an isotropic elastic medium. However, in general, the displacement value of the LRDM model is greater than that of the SSUbi model, the Ubi model and the AEUbi model. This is because the LRDM model considers the degradation of the surrounding rock deformation parameters with the equivalent plastic strain. ③ When the inclination angle of the layer is 60°, the displacement value calculated by the LRDM model is the largest, the displacement values ​​of the SSUbi model, Ubi model and SSMC model are similar, while the displacement values ​​of the MC model and AEUbi model are much smaller than those of the LRDM model, indicating that the calculation results of these other mechanical models are more conservative than those of LRDM, and the LRDM simulation results are closer to the actual deformation of the layered rock mass.

[0309] (II) Comparative analysis of surrounding rock stress field

[0310] The calculation results of the maximum principal stress field of the tunnel surrounding rock using the layered rock mass mechanics model (LRDM) and the strain softening uniform joint model (SSUbi), the uniform joint model (Ubi), and the transversely isotropic uniform joint model (AEUbi) are shown in the figure. Fig.25 As shown. Fig.25 It can be seen that: ① When the LRDM model is used for simulation, when the inclination angle of the layer gradually increases from 0° to 90°, the distribution of the maximum principal stress field of the surrounding rock of the circular tunnel shows anisotropic characteristics that are significantly affected by the occurrence of the rock formation; ② The stress relaxation zone of the rock around the tunnel gradually shifts from the arch crown concentration area to the arch shoulder and arch foot, and the concentration range of the stress relaxation zone is larger than that of the SSUbi model and the Ubi model, while the tunnel area of ​​the AEUbi model is basically a stress concentration area, which is inconsistent with the actual situation.

[0311] (III) Comparative analysis of the distribution of surrounding rock plastic zones

[0312] The calculation results of the plastic zone distribution of tunnel surrounding rock using the layered rock mass mechanics model (LRDM) and the strain softening ubiquitous joint model (SSUbi), the ubiquitous joint model (Ubi), and the transversely isotropic ubiquitous joint model (AEUbi) are shown in the figure. Fig.26 As shown. Fig.26 It can be seen that: ① When the LRDM model is used for calculation, the inclination angle of the layer gradually increases from 0° to 90°, and the plastic zone of the surrounding rock of the circular tunnel undergoes obvious displacement and rotation changes, with obvious anisotropic characteristics; ② When the inclination angle of the layer is the same, the failure mode of the plastic zone of the LRDM model is also significantly different from that of other models. Around the hole of the circular tunnel, the calculation results of the LRDM model show that the main shear-stretching composite failure of the rock block and the layer occurs, and some areas are accompanied by rock shear failure. However, the calculation results of the SSUbi model and the Ubi model show that the main shear failure of the rock block and the layer occurs, but the shear failure of the layer occurs in some areas. For the AEUbi model, the main shear-stretching composite failure of the joint surface occurs, and no failure occurs in some areas.

[0313] (IV) Comparative Analysis and Summary

[0314] In summary, the deformation field, stress field and plastic zone distribution of the LRDM model are significantly affected by the occurrence of the rock formation and have obvious anisotropic characteristics. Compared with other models (SSUbi, Ubi, AEUbi, MC, SSMC), the deformation field, stress field and plastic zone calculated by the LRDM model are more consistent with the actual situation of layered rock mass failure. The specific reasons are as follows:

[0315] ① Although the SSUbi model takes into account the changes in the strength parameters of rock blocks and layers with equivalent plastic strain, it ignores the mechanical characteristics that layered rock mass is usually regarded as a transversely isotropic elastic medium before failure, and does not consider the degradation of deformation parameters with equivalent plastic strain after rock mass failure;

[0316] ② The Ubi model takes into account the influence of the layer on the basis of the Mohr-Coulomb elastoplastic model, but it is still only for isotropic elastic media, and it does not consider the changes of the strength parameters of the rock blocks and layers with the equivalent plastic strain during the deformation and failure process;

[0317] ③The AEUbi model takes into account the mechanical property that the layered rock mass is regarded as a transversely isotropic elastic medium before failure. However, in the plastic failure correction, only the plastic correction of the layer is considered without considering the plastic deformation of the rock block, let alone the degradation characteristics of the mechanical parameters of the rock mass.

[0318] ③ Both the SSMC model and the MC model only consider the isotropic situation. Although the SSMC model takes into account the degradation of the rock strength parameters, both the SSMC model and the MC model do not consider the influence of the layered rock surface and ignore the degradation of the rock deformation parameters.

[0319] ④ The LRDM model incorporates the mechanical property that the layered rock mass can usually be regarded as a transversely isotropic elastic medium before failure into the analysis of the deformation and failure process of the layered rock mass. It also considers the influence of stress state changes on the mechanical parameters of the layered rock mass, which is manifested as the change of deformation parameters and strength parameters with the change of equivalent plastic strain. At the same time, the elastic-plastic deformation of the rock block and the layer is considered separately, and the influence of the layer on the deformation is highlighted. Table 1 lists the feature comparisons of the LRDM model with the SSUbi model, Ubi model, AEUbi model, SSMC model, and MC model in the existing software platform.

[0320] Table 1 Comparison of LRDM model characteristics with other models

[0321]

[0322]

[0323] System Example

[0324] According to an embodiment of the present invention, a system for constructing a constitutive model of a layered rock mass is provided. Fig. 27 Schematic diagram of a system for constructing a layered rock mass constitutive model according to an embodiment of the present invention. Fig. 27 As shown, the layered rock mass constitutive model construction system according to an embodiment of the present invention specifically includes:

[0325] Basic module 2700, used to establish basic assumptions based on the basic mechanical characteristics of layered rock masses;

[0326] The construction module 2702 is used to construct a constitutive model of the layered rock mass based on the basic assumptions and in the framework of an equivalent continuous medium according to the mechanical properties of the layered rock mass with progressive degradation.

[0327] The embodiment of the present invention is a system embodiment corresponding to the above-mentioned method embodiment. The specific operations of each module can be understood by referring to the description of the method embodiment, which will not be repeated here.

[0328] In summary, the constitutive model of layered rock mass in the embodiment of the present invention provides a numerical calculation method for the mechanical response of layered rock mass based on the mechanical properties of the progressive degradation of the rock mass. According to the six basic assumptions proposed, the mechanical properties of the progressive degradation of the layered rock mass are considered, and a mechanical model of the layered rock mass is established within the framework of an equivalent continuous medium. The rationality and accuracy of the mechanical model of the layered rock mass are verified based on comparison of theoretical solutions, analysis of numerical test results, comparison of indoor test results, and simulation of circular tunnel excavation effects. The embodiment of the present invention highlights the contribution of the layer in the process of rock mass deformation and destruction, reflects the influence of stress state changes on the physical and mechanical parameters of the rock mass, realizes the progressive degradation of the mechanical parameters of the layered rock mass with the equivalent plastic strain during the deformation and destruction process, and improves the rationality and accuracy of the numerical calculation simulation results carried out during the construction of layered rock mass engineering.

[0329] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit it. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the aforementioned embodiments, or replace some or all of the technical features therein with equivalents. However, these modifications or replacements 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 constructing a constitutive model of a layered rock mass, characterized in that include: Establish basic assumptions based on the basic mechanical characteristics of layered rock masses; as well as Based on the basic assumptions, a constitutive model of layered rock mass is constructed in an equivalent continuous medium framework according to the mechanical properties of the layered rock mass with progressive degradation.

2. The method according to claim 1, characterized in that The method further comprises: The constitutive model of layered rock mass is verified by comparing theoretical solutions, analyzing numerical test results, comparing indoor test results, and simulating circular tunnel excavation effects.

3. The method according to claim 1, characterized in that The basic assumptions include: the layered rock mass is composed of rock blocks and layers, the layer is the only weak surface in the rock mass, and the rock layer extends infinitely; The layered rock mass is considered as a transversely isotropic elastic medium before and after failure along the layer surface; The occurrence of rock formations can uniquely determine the spatial distribution characteristics of layered rock masses; Under low and medium stress conditions, rock formation strength has nothing to do with the intermediate principal stress; The attenuation of the dilatancy angle of the layered rock mass is not considered, and only the degradation of the elastic modulus, cohesion and internal friction angle is considered; The changes in rock physical and mechanical parameters caused by confining pressure are not considered.

4. The method according to claim 1, characterized in that: Based on the basic assumptions, the constitutive model of layered rock mass is constructed according to the mechanical properties of the progressive degradation of layered rock mass within the framework of equivalent continuous medium, specifically including: S41: Incorporating the joint deformation characteristics into the Oda fabric tensor to obtain the joint flexibility tensor, and constructing an equivalent elastic flexibility matrix of the layered rock mass in a local coordinate system through the joint flexibility tensor; S42: Establish the transformation relationship between the elastic stiffness matrix in the global coordinate system, the local coordinate system and the principal stress space coordinate system; S43: Determine strength criteria, plastic potential functions, and flow laws applicable to rock masses and layers; S44: Using the corresponding internal variables and mechanical parameter evolution laws, the elastic modulus, cohesion and internal friction angle are taken as functions of equivalent plastic strain; S45: Establish evaluation indicators for the process and degree of progressive degradation of layered rock masses.

5. The method according to claim 4, characterized in that The S41 specifically includes: The joint deformation characteristics are integrated into the Oda fabric tensor in the form of deformation stiffness coefficients to obtain the joint flexibility tensor as shown in Formula 1, the eigenvalues ​​and eigenvectors of the joint flexibility tensor are calculated, and the equivalent elastic flexibility matrix of the layered rock mass in the local coordinate system as shown in Formula 2 is constructed based on the eigenvalues ​​and eigenvectors; wherein the equivalent elastic flexibility matrix of the layered rock mass after the dominant joint group is weakened in the local coordinate system is shown in Formula 3; Where JCT represents the joint flexibility tensor, E(n,JF,r) is the joint probability density function of the joint normal vector, joint stiffness and joint geometric dimensions; JF is the joint deformation stiffness coefficient, r is the statistical average diameter of the joint crack, n is the unit direction vector of the joint crack, ρ v is the volume density of the joint fissure, E′1 is the Young's modulus in the transverse isotropic surface, E′3 is the Young's modulus perpendicular to the transverse isotropic surface, and v 12 is the Poisson’s ratio in the transversely isotropic plane, v 13 is the Poisson’s ratio in the perpendicular transverse isotropic plane, G′ 12 is the shear modulus in the transverse isotropic plane, G′ 13 is the shear modulus perpendicular to the transverse isotropic plane, E1 is the equivalent Young's modulus in the transverse isotropic plane, E3 is the equivalent Young's modulus perpendicular to the transverse isotropic plane, G 13 is the equivalent shear modulus perpendicular to the transverse isotropic plane, G 12 is the equivalent shear modulus in the transverse isotropic plane, D' R is the flexibility matrix of the rock mass, D' J The local flexibility matrix for the main feature direction, k n0 and K s0 are the benchmark joint normal stiffness and joint shear stiffness respectively.

6. The method according to claim 4, characterized in that The S43 specifically includes: The shear failure of rock blocks and layers adopts the Mohr-Coulomb yield failure criterion, and the tensile failure of rock blocks and layers adopts the Rankine maximum tensile stress yield failure criterion. The plastic potential function of shear failure of rock blocks and layers obeys the non-associated flow law, and the plastic potential function of tensile failure of rock blocks and layers obeys the associated flow law.

7. The method according to claim 6, characterized in that The S43 specifically includes: When the rock block is sheared, the principal stress of the rock block after shear plasticity correction is as shown in Formula 4; when the rock block is tensilely damaged, the principal stress of the rock block after tensile plasticity correction is as shown in Formula 5; when the layer is sheared, the stress of the layer after shear plasticity correction is as shown in Formula 6, and the stress increment after shear plasticity correction is as shown in Formula 7; when the layer is tensilely damaged, the stress of the layer after tensile plasticity correction is as shown in Formula 8, and the stress increment after tensile plasticity correction is as shown in Formula 9; in, They are respectively the current maximum principal stress, intermediate principal stress and minimum principal stress when calculating according to elastic stress. are the maximum principal stress, intermediate principal stress and minimum principal stress after shear plasticity correction, They are the maximum principal stress, intermediate principal stress and minimum principal stress after tensile plasticity correction respectively; is the plasticity multiplier in different cases, P ij The corresponding component in the elastic stiffness matrix P in the principal stress space coordinate system, ψ m is the rock dilatancy angle, φ is the internal friction angle; are the stress components in the local coordinate system when calculating according to elastic stress; are the stress components in the local coordinate system after the correction of the plane shear plasticity, C 13 , C 23 , C 33 is the corresponding component in the elastic stiffness matrix C in the local coordinate system, ψ j is the dilatancy angle of the layer; Δσ 1′1′ , Δσ 2′2′ , Δσ 3′3′ , Δσ 1′3′ , Δσ 2′3′ is the stress increment after shear failure of the surface in the local coordinate system, is the component corresponding to the normal stress of the plane when calculating according to elastic stress; are the stress components in the local coordinate system after the tensile plasticity correction of the layer, is the stress increment after tensile failure of the layer in the local coordinate system, is the tensile strength of the layer.

8. The method according to claim 4, characterized in that The S44 specifically includes: The cohesion and internal friction angle of the layer are calculated according to formula 10 and formula 11; Among them, c j0 ,φ j0 are the reference values ​​of the cohesion and internal friction angle of the layer, c jN ,φ jN are the current values ​​of the layer cohesion and internal friction angle, c jr ,φ jr are the residual values ​​of the interfacial cohesion and internal friction angle, are the equivalent plastic strain values ​​corresponding to the residual values ​​of the interfacial cohesion and internal friction angle, κ j is the current equivalent plastic strain value of the layer.

9. The method according to claim 4, characterized in that The S45 specifically includes: Formula 12 is used to calculate the rock layer damage index, and the evaluation index of the progressive degradation process and degree of the layered rock mass as shown in Formula 13 is established based on the rock layer damage index; LRFD=max(RFD m , RFD j ) Formula 13; Among them, g i (θ) is the plastic potential function of the layer, and are the equivalent plastic strain on the plane and its corresponding limit equivalent plastic strain, F j Represents the yield function of the layer, RFD j Represents the rock layer damage index, RFD m represents the surrounding rock damage index, LRFD represents the evaluation index of the progressive degradation process and degree of layered rock mass, p and q represent the average stress and equivalent shear stress respectively, A j , B j , C j All are the layer strength criterion coefficients.

10. A system for constructing a layered rock mass constitutive model, characterized in that include: Basic module, used to establish basic assumptions based on the basic mechanical characteristics of layered rock masses; as well as A construction module is used to construct a constitutive model of layered rock mass based on the basic assumptions in an equivalent continuous medium framework according to the mechanical properties of the layered rock mass with progressive degradation.

Images