Deep rock hydraulic coupling constitutive model and finite element simulation method

Abstract

The invention provides a deep rock hydraulic coupling constitutive model and finite element simulation method. The method comprises the steps that a first relation is determined according to the equivalent bulk modulus of a multiphase medium composite material; determining a second relationship according to the multiphase medium deformation coordination relationship; constructing a coupling equation based on the first relation and the second relation; obtaining a yield criterion and a flow rule according to a Drucker-Prager criterion; obtaining a plastic multiplier according to a consistency principle and a flow rule, and updating the stiffness matrix according to the plastic multiplier; determining coupling equation model parameters through a deep rock mechanical hydraulic coupling test; and performing language writing on the coupling equation through a programming language to obtain a finite element model. According to the method, the dynamic coupling effect among the porosity, the water pressure and the stress environment is quantified, the programming language is combined, secondary development of a finite element algorithm is achieved, and the safety of engineering design and the reliability of safety evaluation are improved.

Description

Technical Field

[0001] The present invention relates to the field of rock mechanics and engineering technology, and in particular to a deep rock mass hydraulic coupling constitutive model and a finite element simulation method. Background Art

[0002] In geological and geotechnical engineering, the mechanical behavior of rock plays a crucial role in project safety and stability. This is particularly true in engineering fields involving fluid-rock interactions, such as groundwater resource development and oil and gas extraction. Changes in rock porosity and their mutual influence on rock mechanical properties cannot be ignored. Traditional rock mechanics models typically only consider the unidirectional effect of rock mechanical behavior on porosity changes, failing to account for the interrelationships among porosity, water pressure, and the stress environment. The complex interactions among rock porosity, water pressure, and the stress environment significantly impact the overall stability of a project. For example, during groundwater resource development, changes in water pressure can cause changes in rock porosity, which in turn affects the rock's mechanical properties and stress state. Simultaneously, rock deformation can in turn affect water seepage paths and porosity, forming a dynamic coupling relationship. Failure to accurately describe this coupling interaction can lead to significant deviations in engineering design and safety assessments, potentially posing engineering risks. Therefore, it is crucial to design a hydraulically coupled constitutive model for deep rock masses and use finite element simulation methods. Summary of the Invention

[0003] The purpose of this invention is to provide a hydraulic coupling constitutive model of deep rock mass and a finite element simulation method, which quantifies the dynamic coupling between porosity, water pressure and stress environment, and combines it with a programming language to achieve secondary development of the finite element algorithm and improve the safety of engineering design and the reliability of safety assessment.

[0004] To achieve the above object, the present invention provides the following solutions:

[0005] A deep rock mass hydraulic coupling constitutive model and finite element simulation method includes the following steps:

[0006] Determining a first relationship based on an equivalent bulk modulus of a multiphase medium composite material; the first relationship being a relationship between a rock bulk modulus and a rock skeleton bulk modulus;

[0007] The second relationship is determined based on the deformation coordination relationship of the multiphase medium; the second relationship is the relationship between the porosity of the rock under different stress states and the initial porosity and the volume strain of the multiphase medium;

[0008] Based on the first and second relations, the coupling equations are constructed through Biot theory;

[0009] The yield criterion and flow law of the coupled equations are obtained according to the Drucker-Prager criterion;

[0010] Calculate according to the consistency principle and flow law to obtain the plastic multiplier, and update the stiffness matrix of the coupling equation according to the plastic multiplier;

[0011] The model parameters of the updated coupling equations are determined through deep rock mechanics and hydraulic coupling tests;

[0012] The coupled equations after determining the model parameters are written in a programming language to obtain a finite element model.

[0013] Optionally, determining the first relationship according to the equivalent bulk modulus of the multiphase dielectric composite material includes:

[0014] By formula Determine the equivalent bulk modulus of a multiphase dielectric composite material, where K e is the bulk modulus of the composite medium, N is the total number of components, C i is the ratio of the volume of the i-th medium to the total volume of the composite medium, K i is the bulk modulus of the i-th medium;

[0015] When the multiphase medium is a dual medium consisting of pores and rock skeleton, the first relationship is determined according to the equivalent bulk modulus of the multiphase medium composite material; the expression of the first relationship is: K = K s (1-φ)+K w φ; where K is the rock bulk modulus, K s is the rock skeleton bulk modulus, φ is the rock porosity, K w is the bulk modulus of water.

[0016] Optionally, determining the second relationship according to the deformation coordination relationship of the multiphase medium includes:

[0017] Determine the change in rock volume based on the change in initial porosity;

[0018] The porosity is calculated based on the volume change of the rock;

[0019] The porosity is updated according to different stress states of the rock to obtain the second relationship.

[0020] Optionally, based on the first relationship and the second relationship, a coupling equation is constructed by Biot theory, including:

[0021] Based on the static equilibrium relationship, the third relationship is determined according to the porosity; the third relationship is the equilibrium relationship between rock skeleton stress and pore water pressure;

[0022] Based on the Lemaitre strain equivalence principle, the damage variable is calculated according to the third relationship;

[0023] The fourth and fifth relationships are determined based on the damage variables; the fourth relationship is the relationship between the net water pressure of the stress state of the rock skeleton and the net water pressure of the stress state of the rock, and the fifth relationship is the relationship between the second invariant of the effective deviator stress and the second invariant of the deviator stress of the stress state of the rock;

[0024] The third, fourth and fifth relations are integrated to obtain the coupling equation.

[0025] Optionally, the yield criterion and flow law of the coupled equations are derived according to the Drucker-Prager criterion, including:

[0026] The yield criterion of the coupled equation is obtained through the yield function of the Drucker-Prager criterion; the expression of the yield criterion is: Where f′ is the yield function, η is the first material constant, k0 is the second material constant, σ is the stress tensor, s is the deviatoric stress tensor, and J2′ is the second invariant of the effective deviatoric stress;

[0027] The flow law of the coupled equation is obtained based on the non-associated flow law; the expression of the flow law is: in, is the plastic strain, λ is the plastic multiplier, and g is the plastic function.

[0028] Optionally, a plastic multiplier is calculated based on the consistency principle and the flow law, and the stiffness matrix of the coupled equation is updated based on the plastic multiplier, including:

[0029] The direction of the return vector is determined by the regression mapping algorithm and the flow law to obtain the updated stress;

[0030] The scaling factor is obtained based on the updated stress and the second invariant of the effective deviatoric stress; the scaling factor is determined by the plastic multiplier;

[0031] Based on the consistency principle, the correlation equation between the plastic multiplier and the updated stress is constructed and solved by the Newton–Raphson iterative method to obtain the plastic multiplier.

[0032] Update the deviatoric stress tensor and net water pressure according to the scaling factor and calculate the damage coefficient;

[0033] The updated deviatoric stress tensor and the differential form of the net water pressure are obtained through the partial derivative components of the Jacobian matrix in the coupling equations;

[0034] Linearize the consistency principle and combine it with the differential form to obtain the elasticity relationship;

[0035] The parameters of the stiffness matrix are updated according to the elastic relationship.

[0036] Optionally, calculations are performed based on the consistency principle and the flow law to obtain a plastic multiplier, and the stiffness matrix of the coupling equation is updated based on the plastic multiplier. This also includes: correcting the direction of the return vector in the vertex area of ​​Drucker-Prager, verifying the consistency principle through the corrected direction, and updating the stress and stiffness matrices of the vertex area through the verified consistency principle.

[0037] Optionally, the model parameters of the updated coupling equation are determined through deep rock mechanics hydraulic coupling tests, specifically: the model parameters of the updated coupling equation are determined through sandstone triaxial compression tests and rock CT tests with a water pressure of 8 MPa and a confining pressure of 14 MPa; the model parameters include: elastic modulus, Poisson's ratio, internal friction angle, cohesion and porosity.

[0038] Optionally, the programming language is Fortran.

[0039] Optionally, the finite element model is an ABAQUS finite element model.

[0040] According to a specific embodiment provided by the present invention, the present invention discloses the following technical effects: a deep rock mass hydraulic coupling constitutive model and finite element simulation method provided by the present invention, comprising: determining a first relationship based on the equivalent bulk modulus of a multiphase medium composite material; determining a second relationship based on the deformation coordination relationship of the multiphase medium; constructing a coupling equation based on the first and second relationships using Biot theory; obtaining the yield criterion and flow law of the coupling equation according to the Drucker-Prager criterion; calculating based on the consistency principle and the flow law to obtain a plastic multiplier, and updating the stiffness matrix of the coupling equation based on the plastic multiplier; determining the model parameters of the updated coupling equation through deep rock mechanics hydraulic coupling testing; and programming the coupling equation after the model parameters are determined using a programming language to obtain a finite element model. This method quantifies the dynamic coupling between porosity, water pressure, and stress environment, and combines it with a programming language to achieve secondary development of the finite element algorithm, thereby improving the safety of engineering design and the reliability of safety assessment. BRIEF DESCRIPTION OF THE DRAWINGS

[0041] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0042] Figure 1This is a flow chart of the coupled constitutive model finite element simulation method of the present invention;

[0043] Figure 2 A flow chart for updating the stiffness matrix of the coupling equation of the present invention;

[0044] Figure 3 This is a comparison diagram of the full stress-strain curve of the sandstone triaxial compression test and the numerical simulation curve according to an embodiment of the present invention.

[0045] Figure 4 Writing a logic diagram for the Fortran language of the present invention;

[0046] Figure 5 Schematic diagram of parameter variation characteristics of the constitutive model of the present invention in three-dimensional principal stress space. DETAILED DESCRIPTION

[0047] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of the present invention.

[0048] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, the present invention is further described in detail below with reference to the accompanying drawings and specific embodiments.

[0049] like Figure 1 As shown, the present invention provides a deep rock mass hydraulic coupling constitutive model and finite element simulation method, comprising the following steps:

[0050] Step 100: determining a first relationship based on the equivalent bulk modulus of the multiphase medium composite material; the first relationship is the relationship between the rock bulk modulus and the rock skeleton bulk modulus;

[0051] Specifically, through the formula Determine the equivalent bulk modulus of a composite material of a multiphase medium, where K e is the bulk modulus of the composite medium, N is the total number of components, C i is the ratio of the volume of the i-th medium to the total volume of the composite medium, K i is the bulk modulus of the i-th medium. When the rock is a dual medium consisting of pores and rock skeleton, the relationship between the bulk modulus of the rock and the rock skeleton is determined according to the equivalent bulk modulus of the multiphase medium composite material, and the expression is:

[0052] K=K s (1-φ)+K w φ;

[0053] Where K is the rock bulk modulus, K s is the rock skeleton bulk modulus, φ is the rock porosity, K w is the bulk modulus of water.

[0054] Step 200: Determine a second relationship based on the deformation coordination relationship of the multiphase medium; the second relationship is the relationship between the porosity of the rock under different stress states, the initial porosity, and the volume strain of the multiphase medium;

[0055] Specifically, from the definition of rock porosity:

[0056]

[0057] Where V is the volume of the rock, V s is the volume of the rock skeleton. When the porosity changes from the initial state φ0 to φ, the rock and the rock skeleton produce a certain volume strain, which are:

[0058] ΔV=V0ε v ;

[0059] ΔV s =V s0 ε vs ;

[0060] Among them, ΔV and ΔV s is the volume of rock and rock skeleton changes, V0 and V s0 is the initial volume of rock and rock skeleton, ε v and ε vs is the volume strain of the rock and rock skeleton. When the rock is deformed, the porosity is calculated as follows:

[0061]

[0062] If the solid is incompressible, then ε vs = 0. When the rock is in the elastic stage, the relationship between the volume strain of the rock and the net water pressure is: Where p is the net water pressure of the rock. Integrating the relationship between the incompressible state of the solid and the elastic stage, the porosity is updated to: When the rock is in the plastic stage, update the porosity to: in is the elastic volume strain, is the plastic volume strain. After integrating the updated porosity, the relationship between the porosity of rock under different stress states, the initial porosity and the volume strain of multiphase medium is obtained, and the expression is:

[0063] Step 300: constructing a coupling equation based on the first relationship and the second relationship through Biot theory;

[0064] Specifically, the overburden pressure acts on the cross section of the rock, the skeleton stress acts on the skeleton cross section, and the pore water pressure acts on the pore area. Based on the static equilibrium relationship, the expression is:

[0065] σ ij A z =σ s A s +P w A w ;

[0066] where σ ij , σ s 、P w are rock stress, rock skeleton stress, and pore water pressure, respectively. z 、A s 、A w are the rock cross-sectional area, rock skeleton cross-sectional area, and pore area, respectively. Since the area affected by the rock skeleton stress is almost equal to the rock cross-sectional area, the porosity and equilibrium relationship are combined to obtain the equilibrium relationship between the rock skeleton stress and the pore water pressure, which can be expressed as follows:

[0067] σ′ ij =σ s =σ ij -BP w ;

[0068] in According to the theory of continuum damage mechanics, the damage process is described by the damage variable D. The damage variable is calculated based on the Lemaitre strain equivalence principle. The calculation formula is:

[0069]

[0070] Then, based on the damage variable, we can obtain the relationship between the net water pressure of the stress state of the rock skeleton and the net water pressure of the stress state of the rock, and the relationship between the second invariant of the effective deviator stress and the second invariant of the deviator stress of the stress state of the rock. The calculation formula is as follows:

[0071]

[0072] in is the net water pressure of the stress state of the rock skeleton, p is the net water pressure of the stress state of the rock, is the second deviatoric stress invariant of the stress state of the skeleton, also known as the effective second deviatoric stress invariant, and J2 is the second deviatoric stress invariant of the stress state of the rock. Finally, all the relationships obtained in step 300 are integrated into a coupled equation.

[0073] Step 400: Obtain the yield criterion and flow law of the coupled equation according to the Drucker-Prager criterion;

[0074] Specifically, the yield behavior of rock under complex stress state is measured by the yield function of the Drucker-Prager criterion, that is, the yield criterion, which is expressed as:

[0075]

[0076] Where f′ is the yield function, η is the first material constant, which affects the contribution of the mean stress term in the yield criterion to the yield condition and reflects the sensitivity of the material yield to the mean stress. k0 is the second material constant, which comprehensively reflects the inherent strength characteristics of the material itself and is independent of the specific form of the stress state. It is a strength-related parameter inherent to the material itself. σ is the stress tensor, and s is the deviatoric stress tensor. The volume expansion generated during the development of pores and microcracks causes damage evolution in the rock. As pores and microcracks develop, damage gradually accumulates, resulting in stress reduction, which manifests as a softening phenomenon. This embodiment also constructs a softening equation based on damage mechanics and volumetric strain increments. The expression is:

[0077]

[0078] Where a1 is the material parameter, which is related to the water pressure. is the total strain at time t, is the elastic strain at time t. Then according to the non-associative flow law Construct the flow law of the coupled equations. Assuming that the direction of the plastic strain increment at a point in stress space is unique, that is, it is only related to the stress state of the point and has nothing to do with the direction of the applied stress increment, the flow law expression of the coupled equations is obtained:

[0079]

[0080] in, is the plastic strain, λ is the plastic multiplier, and g is the plastic function. The flow vector of the smooth part of the yield surface is also calculated according to the non-associative flow law. The calculation formula is:

[0081]

[0082] Where I is the unit tensor. The flow vector of the smooth part of the vertebral body is calculated by the flow law of the coupled equation, and the corresponding plastic strain will increase. The calculation formula of its increment is:

[0083]

[0084] Step 500: Calculate according to the consistency principle and flow law to obtain the plastic multiplier, and update the stiffness matrix of the coupling equation according to the plastic multiplier; the specific steps are as follows: Figure 2 Shown, including:

[0085] Step 501: Determine the direction of the return vector using a regression mapping algorithm and flow rules to obtain updated stress;

[0086] Specifically, the update formula of the regression mapping algorithm is:

[0087]

[0088] Among them, σ n+1 is the updated stress at step n+1, is the test stress at step n+1, ΔγD e :Ν n+1 is the return vector, and Δγ is the plastic multiplier. The regression mapping algorithm is integrated with the plastic strain increment of the flow vector of the smooth part of the vertebral body to obtain the stress update formula:

[0089]

[0090] Step 502: Obtain a scaling factor based on the updated stress and the second invariant of the effective deviatoric stress; the scaling factor is determined by the plasticity multiplier;

[0091] Specifically, the second invariant of the effective deviatoric stress has the identity relation: Substituting this identity into the stress update formula, we get The updated stress is obtained by scaling the test deviation stress by a factor that depends on Δγ.

[0092] Step 503: constructing a correlation equation between the plastic multiplier and the updated stress based on the consistency principle, and solving it using the Newton–Raphson iterative method to obtain the plastic multiplier;

[0093] Specifically, the consistency principle is expressed as: Then the Newton–Raphson iterative method is used to solve the independent variables Δγ and φ in the consistency principle. The expression of the solution process is:

[0094] Step 504: Update the deviatoric stress tensor and net water pressure according to the scaling factor and calculate the damage coefficient;

[0095] Specifically, the values ​​of the deviatoric stress tensor and net water pressure are updated according to the stress update formula. The expression of the update process is as follows:

[0096]

[0097] And calculate the updated damage coefficient D according to the consistency principle n+1 , the calculation formula is:

[0098]

[0099] Step 505: Obtain the updated deviatoric stress tensor and the differential form of the net water pressure by using the partial derivative components of the Jacobian matrix in the coupling equation;

[0100] Specifically, the calculation formula of the partial derivative component of the Jacobian matrix is ​​as follows:

[0101]

[0102] where μ is Poisson's ratio. The relationship between the deviatoric stress tensor and the net water pressure under a given stress state is expressed as:

[0103]

[0104] Combining this relational expression with the updated deviatoric stress tensor, we can obtain the deviatoric stress tensor under the current stress state. The calculation formula is:

[0105]

[0106] The deviatoric stress tensor is then converted into differential form, expressed as:

[0107]

[0108] in is parallel to The unit tensor of . The conversion process is the same as that of the deviatoric stress tensor, and the differential form expression of the net water pressure is obtained as follows:

[0109] Step 506: Linearize the consistency principle and combine it with the differential form to obtain the elasticity relationship;

[0110] Specifically, the elasticity relationship is expressed as:

[0111]

[0112] This expression can be organized as:

[0113]

[0114] Step 507: Update the parameters of the stiffness matrix according to the elastic relationship.

[0115] Specifically, the stiffness matrix is ​​expressed as:

[0116]

[0117] The expression after updating the stiffness matrix parameters through the elastic relationship is:

[0118]

[0119] Furthermore, the process of updating the coupling equation stiffness matrix also includes: correcting the direction of the return vector in the vertex area of ​​Drucker-Prager, verifying the consistency principle through the corrected direction, and updating the stress and stiffness matrices of the vertex area through the verified consistency principle.

[0120] Furthermore, the constraints within the vertex region of Drucker-Prager are:

[0121]

[0122] J 2,n+1 =0;

[0123] The direction of the corrected return vector is corrected by this constraint condition to obtain the updated stress, which is expressed as: Then, according to the updated stress, the consistency principle is solved to obtain The calculation formula is: And convert it into differential form:

[0124] Setting this differential expression equal to 0 gives:

[0125]

[0126] Combined with the expression of the elastic stiffness matrix Perform comprehensive calculations to obtain the updated stiffness matrix expression:

[0127] Step 600: Determine the model parameters of the updated coupling equation through deep rock mechanics and hydraulic coupling tests;

[0128] Specifically, this embodiment determines the model parameters of the coupling equation based on the results of the rock triaxial seepage test and the rock physical properties test. First, a triaxial compression test of sandstone with a water pressure of 8MPa and a confining pressure of 14MPa is carried out. The full stress-strain curve of the test is compared with the numerical simulation curve. Figure 3 The initial porosity of the sandstone was then determined to be 0.15 through CT testing of the rock. Triaxial compression tests were then conducted on the sandstone under different confining pressures. Based on the plane trajectory of the mean pressure and deviatoric stress, the model parameters were determined to be: elastic modulus 30.60 GPa, Poisson's ratio 0.3069, internal friction angle 30.90°, cohesion 24.90 MPa, and porosity 0.15.

[0129] Step 700: Compile the coupling equations after determining the model parameters using a programming language to obtain a finite element model.

[0130] Specifically, based on the updated model theory formula of the coupling equation, the Fortran language is written in UMAT, and the logic is written as follows Figure 4 As shown, the compiled UMAT subroutine is embedded into the ABAQUS finite element model through the interface of the ABAQUS main program, and then the model is submitted and the calculation is completed. The model parameters obtained in step 600 are then brought into the Edit Material-General-UserMaterial and load modules of the Property module in ABAQUS to form a constitutive model that takes into account water damage and water pressure adjustment, thereby realizing finite element simulation. The parameter change characteristics of this model in the three-dimensional principal stress space are shown as follows Figure 5 shown.

[0131] The beneficial effects of the present invention are as follows:

[0132] 1) By introducing a coupled equation for the interaction between water pressure changes in porous media and solid deformation, a hydraulic coupled constitutive model that considers changes in rock porosity was established. This improved the accuracy of characterizing the effects of rock stress environment, water pressure conditions, and porosity on rock plastic behavior and rock deformation characteristics.

[0133] 2) It overcomes the problem that traditional constitutive models cannot take into account the continuous changes in porosity and pore water pressure with stress state, thereby improving the safety of engineering design and the reliability of safety assessment.

[0134] The various embodiments in this specification are described in a progressive manner, and each embodiment focuses on the differences from other embodiments. The same or similar parts between the various embodiments can be referenced to each other.

[0135] The present invention uses specific examples to illustrate the principles and implementation methods of the present invention. The above examples are only intended to help understand the method and core concept of the present invention. At the same time, those skilled in the art will find that the specific implementation methods and application scopes may vary based on the concept of the present invention. In summary, the contents of this specification should not be construed as limiting the present invention.

Claims

1. A hydraulic coupled constitutive model of deep rock mass and a finite element simulation method, characterized in that: The steps include: Determining a first relationship based on an equivalent bulk modulus of the multiphase medium composite material; the first relationship being a relationship between the rock bulk modulus and the rock skeleton bulk modulus; Determining a second relationship based on the deformation coordination relationship of the multiphase medium; The second relationship is the relationship between the porosity of rock under different stress states and the initial porosity and volume strain of multiphase medium; Based on the first relationship and the second relationship, constructing a coupling equation through Biot theory; The yield criterion and flow law of the coupled equation are obtained according to the Drucker-Prager criterion; Calculating according to the consistency principle and the flow law to obtain a plastic multiplier, and updating the stiffness matrix of the coupling equation according to the plastic multiplier; Determining updated model parameters of the coupling equation through deep rock mechanics and hydraulic coupling tests; The coupling equations after determining the model parameters are written in a programming language to obtain a finite element model.

2. The deep rock mass hydraulic coupling constitutive model and finite element simulation method according to claim 1 is characterized in that: The first relationship is determined according to the equivalent bulk modulus of the multiphase dielectric composite material, including: By formula Determine the equivalent bulk modulus of the multiphase dielectric composite material, where K e is the bulk modulus of the composite medium, N is the total number of components, C i is the ratio of the volume of the i-th medium to the total volume of the composite medium, K i is the bulk modulus of the i-th medium; When the multiphase medium is a dual medium consisting of pores and rock skeletons, the first relationship is determined according to the equivalent bulk modulus of the multiphase medium composite material; the expression of the first relationship is: K = K s (1-φ)+K w φ; where K is the rock bulk modulus, K s is the rock skeleton bulk modulus, φ is the rock porosity, K w is the bulk modulus of water.

3. The deep rock mass hydraulic coupling constitutive model and finite element simulation method according to claim 1 is characterized in that: The second relationship is determined based on the deformation coordination relationship of the multiphase medium, including: determining a change in rock volume according to a change in the initial porosity; Calculating the porosity according to the rock volume change; The porosity is updated according to different stress states of the rock to obtain the second relationship.

4. The deep rock mass hydraulic coupling constitutive model and finite element simulation method according to claim 1, characterized in that: Based on the first relationship and the second relationship, a coupling equation is constructed using Biot theory, including: Based on the static equilibrium relationship, a third relationship is determined according to the porosity; the third relationship is the equilibrium relationship between rock skeleton stress and pore water pressure; Based on the Lemaitre strain equivalence principle, the damage variable is calculated according to the third relationship; Determining a fourth relationship and a fifth relationship based on the damage variable; the fourth relationship is the relationship between the net water pressure of the stress state of the rock skeleton and the net water pressure of the stress state of the rock; and the fifth relationship is the relationship between the second invariant of the effective deviator stress and the second invariant of the deviator stress of the stress state of the rock; The third relationship, the fourth relationship, and the fifth relationship are integrated to obtain the coupling equation.

5. The deep rock mass hydraulic coupling constitutive model and finite element simulation method according to claim 1 is characterized in that: The yield criterion and flow law of the coupled equations are obtained according to the Drucker-Prager criterion, including: The yield criterion of the coupled equation is obtained by the yield function of the Drucker-Prager criterion; the expression of the yield criterion is: Where f′ is the yield function, η is the first material constant, k0 is the second material constant, σ is the stress tensor, s is the deviatoric stress tensor, and J2′ is the second invariant of the effective deviatoric stress; The flow law of the coupling equation is obtained according to the non-associated flow law; the expression of the flow law is: in, is the plastic strain, λ is the plastic multiplier, and g is the plastic function.

6. The deep rock mass hydraulic coupling constitutive model and finite element simulation method according to claim 4, characterized in that: Calculation is performed based on the consistency principle and the flow law to obtain a plastic multiplier, and the stiffness matrix of the coupling equation is updated based on the plastic multiplier, including: Determine the direction of the return vector by using the regression mapping algorithm and the flow law to obtain the updated stress; Obtaining a scaling factor according to the updated stress and the second invariant of the effective deviatoric stress; wherein the scaling factor is determined by the plastic multiplier; Constructing a correlation equation between the plastic multiplier and the updated stress based on the consistency principle, and solving it by a Newton–Raphson iterative method to obtain the plastic multiplier; updating the deviatoric stress tensor and the net water pressure according to the scaling factor and calculating the damage coefficient; Obtaining updated differential forms of the deviatoric stress tensor and the net water pressure through partial derivative components of the Jacobian matrix in the coupling equation; Linearizing the consistency principle and combining it with the differential form to obtain the elasticity relationship; The parameters of the stiffness matrix are updated according to the elastic relationship.

7. The deep rock mass hydraulic coupling constitutive model and finite element simulation method according to claim 6, characterized in that: Calculation is performed according to the consistency principle and the flow law to obtain a plastic multiplier, and the stiffness matrix of the coupling equation is updated according to the plastic multiplier. The method also includes: correcting the direction of the return vector in the vertex area of ​​the Drucker-Prager, verifying the consistency principle through the corrected direction, and updating the stress of the vertex area and the stiffness matrix through the verified consistency principle.

8. The deep rock mass hydraulic coupling constitutive model and finite element simulation method according to claim 1, characterized in that: The updated model parameters of the coupling equation are determined through deep rock mechanics hydraulic coupling tests. Specifically, the updated model parameters of the coupling equation are determined through sandstone triaxial compression tests and rock CT tests with a water pressure of 8 MPa and a confining pressure of 14 MPa. The model parameters include elastic modulus, Poisson's ratio, internal friction angle, cohesion, and porosity.

9. The deep rock mass hydraulic coupling constitutive model and finite element simulation method according to claim 1, characterized in that: The programming language is Fortran.

10. The deep rock mass hydraulic coupling constitutive model and finite element simulation method according to claim 1, characterized in that: The finite element model is an ABAQUS finite element model.

Images