Transversely isotropic rock viscoelastic-plastic creep analysis method

Abstract

A transversely isotropic rock viscoelastoplastic creep analysis method. The method establishes a transversely isotropic Nishihara viscoelastoplastic model, and modifies the Drucker-Prager yield and plastic potential function with Hill-type orthogonal anisotropic equivalent stress; a standard polygonal proportional boundary mother unit library is pre-constructed in the parameter space, the shape function and derivative are obtained in advance through Laplace equation, weighted residual value and Hamiltonian system characteristic decomposition, and are reused in actual physical polygons through isoparametric mapping; rheological parameters are identified by least square method according to creep test, and time increment solving is carried out according to trial, yield discrimination, viscoplasticity correction and state updating by using generalized midpoint format and Perzyna overstress flow rule. The method can reduce repeated characteristic solving in nonlinear iteration, and is suitable for long-term deformation analysis of cracks, heterogeneous interfaces and local dense rock mass domains.

Description

Technical Field

[0001] This invention relates to the fields of numerical analysis in geotechnical engineering, rock rheology and nonlinear computational mechanics, and specifically to a method for analyzing the viscoelastic-plastic creep of transversely isotropic rocks. Background Technology

[0002] Underground caverns, tunnels, slopes, underground storage facilities, and other rock engineering structures will experience creep deformation that develops over time under long-term loads. For layered, bedding, or foliated rock masses, the elastic modulus, shear modulus, Poisson effect, strength, and yield characteristics of the material differ significantly in different directions. This directionality is coupled with the viscoelasticity that develops over time and the viscoplastic flow after yielding, resulting in a significant path correlation between long-term stress redistribution and deformation accumulation.

[0003] In existing rock creep analysis, the classical Nishihara model can describe instantaneous elasticity, initial creep, and steady viscoplastic flow using Hooke elastomers, Kelvin viscoelastics, and Bingham viscoplastics. However, if isotropic parameters are used directly, it is difficult to reflect the differences in stiffness, Poisson coupling, and yield in different directions of layered rocks. If anisotropy is introduced only in the elastic stage, it is difficult to uniformly handle pressure-sensitive yielding and time-dependent flow after yielding.

[0004] The finite element method (FEM) has a mature constitutive integral and incremental solution system in regular domains. However, when the computational domain includes cracks, pores, inclusions, material interfaces, or locally refined regions, a large number of triangular or quadrilateral elements are often required to ensure element shape quality and interface compatibility. Especially near crack tips and heterogeneous interfaces, where stress gradients are high, traditional regular meshes need to be significantly refined and interpolation orders increased, leading to increased preprocessing workload, degrees of freedom, and computational costs.

[0005] The scaled boundary finite element method can decompose the governing equations into radial and circumferential parts using scaling centers, and form a high-precision semi-analytical solution based on boundary discretization. It has a significant advantage for elastic or linear viscoelastic problems; however, in plastic and viscoplastic problems, stiffness, overstress state, and internal variables evolve with time and space within the domain, and the radial response usually no longer satisfies a simple scale separation relationship, requiring domain integration and incremental constitutive updates.

[0006] Furthermore, traditional polygon scale boundary solving methods, which establish eigenvalue problems for each physical unit with different geometries, incur repetitive solution overhead during large-scale nonlinear time iterations. Therefore, there is an urgent need for a computational method that unifies transversely isotropic viscoelastic-plastic constitutive models, reusable polygon scale boundary interpolation, isoparametric mapping, and steady-time integration to improve the geometric adaptability and computational efficiency of long-term creep problems in complex rock masses. Summary of the Invention

[0007] The main objective of this invention is to provide a method for analyzing the viscoelastic-plastic creep of transversely isotropic rocks. This method combines transversely isotropic Nishihara viscoelastic-plastic constitutive model, Drucker-Prager yield description with Hill-type orthotropic equivalent stress correction, a pre-built standard polygonal scale boundary element library, isoparametric mapping, and time-incremental integration based on Perzyna hyperstress theory. This approach addresses the challenges of characterizing directional creep in complex rock mass domains, the high overhead of repetitive feature solving, and the high cost of mesh processing for cracks and heterogeneous interfaces.

[0008] To achieve the above objectives, this invention provides a method for analyzing the viscoelastic-plastic creep of transversely isotropic rocks, comprising the following steps: S1. The total creep strain increment of rock is decoupled into three series-connected strain increment components: elastic, viscoelastic, and viscoplastic. For plane strain and plane stress conditions, transversely isotropic unit compliance matrices are constructed respectively, and elastic anisotropy is characterized by Young's modulus, in-plane shear modulus, and Poisson coupling parameters. The incremental constitutive relations of the Hooke elastic, Kelvin viscoelastic, and Bingham viscoplastic branches are established in sequence, and finally a unified incremental constitutive equation of the three branches is formed in series. S2. A Drucker-Prager yield function framework is constructed using mean normal stress and equivalent shear stress, and cohesion and internal friction angle are introduced to reflect the pressure-sensitive yielding characteristics of soil and rock. The equivalent shear stress is weighted by direction using Hill-type orthotropic equivalent stress, and the strength directionality is characterized by orthotropic yield and shear yield parameters to obtain a transversely isotropic modified Drucker-Prager yield function. A plastic potential function is established using isotropic structure, and the dilatation angle parameter is introduced to adapt to associated and unassociated flow rules. S3. Construct a multi-topological standard polygon mother element library in the parameter space, and uniformly adopt counterclockwise boundary node numbering; select the geometric center as the scaling center, discretize the boundary into line elements and divide it into triangular sectors, and establish a radial-circular proportional boundary coordinate system; solve the Laplace equation, derive the Euler-Cauchy equation by the weighted residual method, transform the Euler-Cauchy equation into a Hamiltonian system, and obtain the shape function and its derivative matrix of the mother element proportional boundary after eigenvalue decomposition, and pre-store them to form a reusable shape function library; S4. Generate a polygonal mesh for the actual rock mass computational domain, and call the corresponding mother element according to the topological matching of the physical polygonal nodes; adopt isoparametric mapping, and establish the coordinate interpolation relationship of the physical space through pre-stored shape functions; use the Jacobian matrix to realize the spatial transformation of the derivative of the shape function, derive the strain matrix and domain integral expression of the physical element, and complete the spatial discretization of the nonlinear problem. S5. Based on the time-strain curves of uniaxial or triaxial creep tests, the Nishihara creep expression under constant stress is parameterized and reorganized into a fitting form containing exponential decay terms, constant terms, and linear time terms; the least squares method is used to fit the test data, and the Hooke elastic modulus, Kelvin elastic modulus, Kelvin viscosity coefficient, and Bingham viscosity coefficient are inverted and identified; combined with the rock bedding direction, transversely isotropic elastic, Poisson coupling, and yield anisotropic parameters are set. S6. Using the state of the previous time step as a known quantity, establish a consistent tangent matrix; establish an incremental equilibrium equation based on the principle of virtual work, and express the total strain increment as the displacement increment to obtain the global incremental equation; assume no yielding at each time step, calculate the trial stress state, and then obtain the trial yield value. When the trial state exceeds the static yield surface, correct the current stress; the trial state determination and viscoplasticity correction process are executed sequentially in each time step until the entire loading process is completed.

[0009] The beneficial effects of this invention are: (1) This invention places the transverse isotropic compliance relation, the Nishihara three-branch rheological model and the Hill-type equivalent stress correction Drucker-Prager criterion in the same incremental constitutive framework, so that instantaneous elasticity, initial viscoelastic creep, pressure-sensitive yielding and time-dependent viscoplastic flow after yielding can be considered simultaneously, and allows Young's modulus, shear modulus, Poisson coupling and yield parameters in different directions to be adjusted independently.

[0010] (2) In this invention, a limited number of standard polygonal mother units are pre-constructed in the parameter space, and the proportional boundary shape function and derivative matrix of the mother unit are pre-solved and stored. In actual nonlinear time analysis, physical polygons with the same topology can directly reuse the pre-constructed data through isoparametric mapping, thereby avoiding the problem of repeatedly establishing large-scale eigenvalues ​​for each physical unit and each time step.

[0011] (3) The present invention uses polygonal elements with different topologies such as triangles, quadrilaterals and pentagons to nonagons to describe crack tips, heterogeneous inclusions and local mesh transitions. It can increase boundary nodes and perform local refinement in high gradient regions, and maintain coarser discretization in areas far from critical regions, thereby reducing the cost of generating meshes for complex geometries.

[0012] (4) This invention employs a generalized midpoint scheme, Perzyna overstress theory, and an incremental algorithm of “trial calculation-yield criterion-overstress correction-state update”. Under complex geometric conditions such as cracks, the proposed method achieves consistent displacement and stress responses with a significantly fewer element count than traditional finite element convergence models. Attached Figure Description

[0013] Figure 1This is a flowchart illustrating the method for analyzing the viscoelastic-plastic creep of transversely isotropic rocks based on polygonal scale boundary finite element method of the present invention. Figure 2 A schematic diagram of the Nishihara model; Figure 3 The standard polygon mother cell library and representative topology diagrams are pre-constructed in the parameter space, where: (a) a set of square multi-nodes, (b) a set of regular polygons; Figure 4 The diagram shows the parameter space discretization, where: (a) is a schematic diagram of line element discretization; (b) is a schematic diagram of triangular sector discretization. Figure 5 The first embodiment is a schematic diagram of the geometric boundary conditions and polygonal mesh of a square plate, wherein: (a) schematic diagram of the geometric boundary conditions of the square plate; (b) schematic diagram of the polygonal mesh of the square plate; Figure 6 The first embodiment - a comparison diagram of creep displacement history of monitoring points of square plate, wherein: (a) creep comparison diagram of monitoring point A; (b) creep comparison diagram of monitoring point B; (c) creep comparison diagram of monitoring point C; Figure 7 The second embodiment is a schematic diagram of the geometric boundary conditions and polygonal mesh of the cracked rectangular plate, wherein: (a) schematic diagram of the geometric boundary conditions of the cracked rectangular plate; (b) schematic diagram of the polygonal mesh of the cracked rectangular plate; Figure 8 The second embodiment is a comparison diagram of the creep displacement history of the monitoring points of the cracked rectangular plate, wherein: (a) creep comparison diagram of monitoring point A; (b) creep comparison diagram of monitoring point B. Detailed Implementation

[0014] It should be understood that the specific embodiments described herein are for illustrative purposes only and are not intended to limit the scope of the invention. The realization of the invention's objectives, functional characteristics, and advantages will be further explained in conjunction with the embodiments and with reference to the accompanying drawings.

[0015] As shown in the attached diagram Figure 1 This is a flowchart illustrating the method for analyzing the viscoelastic-plastic creep of transversely isotropic rocks based on polygonal scale boundary finite element method of the present invention. Figure 2 A schematic diagram of the Nishihara model; Figure 3 This is a pre-constructed library of standard polygonal parent elements and a representative topology diagram in the parameter space; Figure 4 This is a schematic diagram of the scaled boundary coordinate system and the discretization of the triangular sector in the parameter space. Figure 5 This is a schematic diagram of the geometric boundary conditions and polygonal mesh of the first embodiment - a square plate; Figure 6 This is a comparison chart of the creep displacement history of the monitoring points of the square plate in the first embodiment; Figure 7This is a schematic diagram of the geometric boundary conditions and polygonal mesh of the second embodiment - a cracked rectangular plate; Figure 8 This is a comparison chart of the creep displacement history of the monitoring points on the cracked rectangular plate in the second embodiment. The detailed method includes the following steps: S1. The total creep strain increment of rock is decoupled into three series-connected strain increment components: elastic, viscoelastic, and viscoplastic. For plane strain and plane stress conditions, transversely isotropic unit compliance matrices are constructed respectively, and elastic anisotropy is characterized by Young's modulus, in-plane shear modulus, and Poisson coupling parameters. The incremental constitutive relations of the Hooke elastic, Kelvin viscoelastic, and Bingham viscoplastic branches are established in sequence, and finally a unified incremental constitutive equation of the three branches is formed in series. S2. A Drucker-Prager yield function framework is constructed using mean normal stress and equivalent shear stress, and cohesion and internal friction angle are introduced to reflect the pressure-sensitive yielding characteristics of soil and rock. The equivalent shear stress is weighted by direction using Hill-type orthotropic equivalent stress, and the strength directionality is characterized by orthotropic yield and shear yield parameters to obtain a transversely isotropic modified Drucker-Prager yield function. A plastic potential function is established using isotropic structure, and the dilatation angle parameter is introduced to adapt to associated and unassociated flow rules. S3. Construct a multi-topological standard polygon mother element library in the parameter space, and uniformly adopt counterclockwise boundary node numbering; select the geometric center as the scaling center, discretize the boundary into line elements and divide it into triangular sectors, and establish a radial-circular proportional boundary coordinate system; solve the Laplace equation, derive the Euler-Cauchy equation by the weighted residual method, transform the Euler-Cauchy equation into a Hamiltonian system, and obtain the shape function and its derivative matrix of the mother element proportional boundary after eigenvalue decomposition, and pre-store them to form a reusable shape function library; S4. Generate a polygonal mesh for the actual rock mass computational domain, and call the corresponding mother element according to the topological matching of the physical polygonal nodes; adopt isoparametric mapping, and establish the coordinate interpolation relationship of the physical space through pre-stored shape functions; use the Jacobian matrix to realize the spatial transformation of the derivative of the shape function, derive the strain matrix and domain integral expression of the physical element, and complete the spatial discretization of the nonlinear problem. S5. Based on the time-strain curves of uniaxial or triaxial creep tests, the Nishihara creep expression under constant stress is parameterized and reorganized into a fitting form containing exponential decay terms, constant terms, and linear time terms; the least squares method is used to fit the test data, and the Hooke elastic modulus, Kelvin elastic modulus, Kelvin viscosity coefficient, and Bingham viscosity coefficient are inverted and identified; combined with the rock bedding direction, transversely isotropic elastic, Poisson coupling, and yield anisotropic parameters are set. S6. Using the state of the previous time step as a known quantity, establish a consistent tangent matrix; establish an incremental equilibrium equation based on the principle of virtual work, and express the total strain increment as the displacement increment to obtain the global incremental equation; assume no yielding at each time step, calculate the trial stress state, and then obtain the trial yield value. When the trial state exceeds the static yield surface, correct the current stress; the trial state determination and viscoplasticity correction process are executed sequentially in each time step until the entire loading process is completed.

[0016] In step S1, the total creep strain increment of the rock is decoupled into three series strain increment components: elastic, viscoelastic, and viscoplastic. in, This represents the total strain increment. , and These represent the elastic strain increment, viscoelastic strain increment, and viscoplastic strain increment, respectively. Indicates increment, subscript n This represents the nth time step.

[0017] For plane strain and plane stress conditions, transversely isotropic unit compliance matrices are constructed respectively, and elastic anisotropy is characterized by directional Young's modulus, in-plane shear modulus, and Poisson coupling parameters: in, To view the transverse isotropic unit compliance matrix, , , and These are Young's moduli in different directions. It is the in-plane shear modulus. and The parameters are used to describe Poisson coupling in different directions, and the reverse Poisson coupling parameters are determined based on the reciprocal work relationship.

[0018] Establish the incremental constitutive relation for the Hooke elastic branch: in, Let be the Young's modulus in the y-direction of the Hooke elastic branch. This represents the stress increment.

[0019] Establish the incremental constitutive relation for the Kelvin viscoelastic branch: in, , , It is the viscoelastic strain rate. It is the stress increment. Let be the Young's modulus in the y-direction of the Kelvin spring. Kelvin viscosity coefficient The length of the nth time step. For the generalized midpoint integral parameter.

[0020] Establish the incremental constitutive relation for the Bingham viscoplastic branch: in It is viscoplastic strain rate. It is the tangent matrix corresponding to the viscoplastic strain rate, and its specific expression is: in, Bingham viscosity coefficient, For Macaulay brackets, It is the overstress function. F Let be the yield function. G Let be the plastic potential function. , and Represent , F and G stress Find the partial derivative, where the superscript T represents the matrix transpose.

[0021] The final result is a unified incremental constitutive equation with three branches connected in series: In step S2, the Drucker-Prager yield function framework is constructed using mean normal stress and equivalent shear stress: in, q For equivalent shear stress, p For average normal stress, a This is the pressure coefficient.

[0022] Introducing cohesion and internal friction angle to reflect the pressure-sensitive yielding characteristics of soil and rock: in, The initial yield strength, It is cohesion. It is the internal friction angle.

[0023] Hill-type orthotropic equivalent stress is used to weight the equivalent shear stress by direction, and the strength directionality is characterized by orthogonal yield and shear yield parameters. in, To view the isotropic equivalent shear stress laterally, It is the transverse isotropic pressure coefficient. The orthogonal yield parameter, This is the shear yield parameter. , and These are the normal stress components in the x, y, and z directions, respectively. , and These are the shear stress components in the yz plane, zx plane, and xy plane, respectively.

[0024] The transversely isotropic modified Drucker-Prager yield function is obtained as follows: A plastic potential function is established using a homogeneous structure, and a dilatation angle parameter is introduced to adapt to both associated and unassociated flow rules: in, It is the shear expansion angle, when When the correlation flow rule is applied, when The non-associated flow rule is used at this time.

[0025] In step S3, a multi-topology standard polygon mother unit library is constructed in the parameter space, and the counterclockwise boundary node numbering is uniformly adopted.

[0026] By selecting the geometric center as the scaling center, discretizing the boundary into line elements and dividing it into triangular sectors, a radial-circumferential scaling boundary coordinate system is established: in, , It is the coordinate value of any point in the parameter space. It is a radial coordinate. It is a circumferential coordinate system. and These are the coordinates of the scaling center. and These are the boundary coordinates, obtained by interpolation using a one-dimensional linear element shape function: in, It is a shape function matrix. , It is the node coordinate vector.

[0027] Solving the Laplace equation, the Euler-Cauchy equation is derived using the weighted residual method: in, It is the radial analytic function of the scalar potential function, subscript Representative to Find the first partial derivative, index Representative to Find the second partial derivative. , and Here is the boundary integral coefficient matrix: in, It is a boundary Jacobian determinant. and Defined as: Among them, subscript Representative to Find the first partial derivative.

[0028] Transforming the Euler-Cauchy equations into Hamiltonian systems: in, It is the internal flux auxiliary function, Hamiltonian matrix. for: After eigenvalue decomposition, the proportional boundary shape functions of the parent element and their derivative matrices are obtained, and a reusable shape function library is pre-stored. in, It is a proportional boundary shape function matrix. It is a vector of nodal values ​​of the scalar potential function. It is the eigenvector matrix. It is an eigenvalue matrix. Represents the eigenvector matrix The extracted local submatrix contains only rows related to the current line cell (or triangular sector) node. , They represent respectively to and Find the partial derivative. It is an identity matrix.

[0029] In step S4, a polygonal mesh is generated for the actual rock mass calculation domain, and the corresponding parent element is called according to the physical polygon node topology matching.

[0030] Using isoparametric mapping, coordinate interpolation relationships in physical space are established through pre-stored shape functions: in, These are physical space coordinates. , It is the coordinate vector of the physical boundary node.

[0031] Spatial transformation of the derivatives of shape functions is achieved using the Jacobian matrix: in, , Representative to , Find the partial derivative; the superscript -1 indicates finding the inverse matrix. The Jacobian matrix has the following specific form: Derivation of the strain matrix and domain integral expression for the physical unit, thus completing the spatial discretization of the nonlinear problem: in, It is the strain matrix of a polygonal element. It is the number of polygon nodes. It is the differential area of ​​the physical domain. and They represent and The differential, It is the Jacobian determinant, the strain matrix of the polygonal element. block matrix The specific form is as follows: in, , These represent the derivative matrices of the proportional boundary shape functions, respectively. , The i-th component.

[0032] In step S5, the time-strain curves are based on uniaxial or triaxial creep tests: in, It is axial strain. Represents the plastic potential function G right Find the partial derivative. t It is creep time. , , and These represent the first, second, and third principal stresses, respectively.

[0033] The Nishihara creep expression under constant stress is parameterized and reorganized into a fitted form containing exponential decay terms, constant terms, and linear time terms: in, , , and represents the fitting coefficient.

[0034] The least squares method was used to fit the experimental data, and the Hooke elastic modulus, Kelvin elastic modulus, Kelvin viscosity coefficient and Bingham viscosity coefficient were identified by inversion. Combined with the rock bedding direction, transversely isotropic elastic, Poisson coupling and yield anisotropy parameters were set.

[0035] In step S6, the state of the previous time step is used as a known quantity to establish a consistent tangent matrix: in, It is a uniform tangent matrix.

[0036] Incremental equilibrium equations are established based on the principle of virtual work: in, Indicates a time step The internal external force increment vector.

[0037] The total strain increment is expressed as the displacement increment: in, It is the displacement increment.

[0038] The global incremental equation is obtained as follows: Among them, the global tangent stiffness matrix Additional viscoelastic force and additional viscoplastic forces The results are obtained by integrating each physical polygon domain separately and then assembling them.

[0039] Assuming no yielding at each time step, calculate the trial stress state: in, , , These represent the experimental stress increment, the experimental total strain increment, and the experimental viscoelastic strain rate at the nth time step, respectively. This represents the trial stress state at the (n+1)th time step. The superscript "trial" indicates trial, and the subscript "[...]" indicates trial. n +1 represents the (n+1)th time step. It is the elastic and viscoelastic tangent stiffness matrix, and its specific expression is: The yield value was then calculated: in, , and These represent the experimental transverse isotropic equivalent shear stress, experimental mean normal stress, and experimental yield function at the (n+1)th time step, respectively.

[0040] When the calculated state exceeds the static yield surface, the current stress is corrected: in, This represents the stress state at the (n+1)th time step.

[0041] The trial state determination and viscoplasticity correction processes are executed sequentially at each time step until the entire loading process is completed.

[0042] Furthermore, to verify the accuracy of the present invention, a first embodiment—a square viscoelastic plate—is provided for verification.

[0043] Example 1: Verification with a square viscoelastic plate. After closing the Bingham viscoplastic branch, the model degenerates into a linear three-parameter solid. Figure 5 The side length of the square plate shown in (a) is l1 = 1m, a uniformly distributed constant tension is applied to the right boundary. P =1 Pa, the left boundary restricts horizontal displacement, the lower boundary restricts vertical displacement, the upper boundary has no external load, plane stress condition is adopted, and the material parameters are: , , Monitoring points A, B, and C are used to examine the displacement response in different directions and at different spatial locations, respectively. Figure 5 (b) in the figure gives the polygon discretization used in this invention.

[0044] Under uniform uniaxial stress, this three-parameter solid has a closed-form displacement solution. Displacements in the x and y directions... and They are respectively: Figure 6 A comparison of the calculation results, closed-form analytical solution, and adaptive proportional boundary finite element results of this invention is presented. The three monitored quantities basically overlap throughout the creep process, which can correctly reflect the transition from initial creep to long-term asymptotic response. At the same time, it reproduces the negative lateral displacement generated by Poisson contraction and the spatial proportional relationship between different monitoring points, indicating that the viscoelastic time integral and polygon interpolation implementations are consistent.

[0045] Furthermore, to verify the applicability of the present invention under complex geometric conditions, a second embodiment is given - viscoelastic-plastic creep of a cracked rectangular plate.

[0046] Example 2: Viscoelastic-plastic creep of a cracked rectangular plate. Figure 7 The dimensions of the rectangular plate shown in (a) are: l 2 = 100mm l 3=50mm, with a horizontal crack of 50mm length inside, and an application applied to the upper boundary. P The uniform strain stress is 2 = 12 MPa, the lower boundary constrains the displacement in the y direction, the right boundary constrains the displacement in the x direction, the creep time is 20 h, and the plane strain condition is adopted; the material parameters are as shown in Table 1.

[0047] Table 1: Representative constitutive parameters of Example 2

[0048] First, conventional finite element mesh convergence analysis was performed using Abaqus, employing an eight-node quadratic plane strain element CPE8. The final creep displacement at the plate edge monitoring point and the y-direction normal stress in the crack tip neighborhood were used as convergence indices. The converged Abaqus model contained 48,890 elements; this invention uses 502 polygonal elements and locally refines the mesh at the crack tip, such as... Figure 7 As shown in (b) of the diagram.

[0049] Figure 8 (a) and Figure 8 Figure (b) shows the creep displacement history of monitoring points A and B. The field distribution and time history obtained by the two methods are similar. On the same Intel Core i9-14900K workstation, the computation time of this invention is 26s, and the convergence time of the Abaqus model is 620s, which is equivalent to about 24 times speedup and 95.8% reduction in computation time. Therefore, this invention reduces computational cost while accurately calculating creep response.

Claims

1. A method for analyzing the viscoelastic-plastic creep of transversely isotropic rocks, characterized in that, The method includes the following steps: S1. The total creep strain increment of rock is decoupled into three series-connected strain increment components: elastic, viscoelastic, and viscoplastic. For plane strain and plane stress conditions, transversely isotropic unit compliance matrices are constructed respectively, and elastic anisotropy is characterized by Young's modulus, in-plane shear modulus, and Poisson coupling parameters. The incremental constitutive relations of the Hooke elastic, Kelvin viscoelastic, and Bingham viscoplastic branches are established in sequence, and finally a unified incremental constitutive equation of the three branches is formed in series. S2. A Drucker-Prager yield function framework is constructed using mean normal stress and equivalent shear stress, and cohesion and internal friction angle are introduced to reflect the pressure-sensitive yielding characteristics of soil and rock. The equivalent shear stress is weighted by direction using Hill-type orthotropic equivalent stress, and the strength directionality is characterized by orthotropic yield and shear yield parameters to obtain a transversely isotropic modified Drucker-Prager yield function. A plastic potential function is established using isotropic structure, and the dilatation angle parameter is introduced to adapt to associated and unassociated flow rules. S3. Construct a multi-topological standard polygon mother element library in the parameter space, and uniformly adopt counterclockwise boundary node numbering; select the geometric center as the scaling center, discretize the boundary into line elements and divide it into triangular sectors, and establish a radial-circular proportional boundary coordinate system; solve the Laplace equation, derive the Euler-Cauchy equation by the weighted residual method, transform the Euler-Cauchy equation into a Hamiltonian system, and obtain the shape function and its derivative matrix of the mother element proportional boundary after eigenvalue decomposition, and pre-store them to form a reusable shape function library; S4. Generate a polygonal mesh for the actual rock mass computational domain, and call the corresponding mother element according to the topological matching of the physical polygon nodes; adopt isoparametric mapping, and establish the coordinate interpolation relationship of the physical space through pre-stored shape functions; use the Jacobian matrix to realize the spatial transformation of the derivative of the shape function, derive the strain matrix and domain integral expression of the physical element, and complete the spatial discretization of the nonlinear problem. S5. Based on the time-strain curves of uniaxial or triaxial creep tests, the Nishihara creep expression under constant stress is parameterized and reorganized into a fitting form containing exponential decay terms, constant terms, and linear time terms; the least squares method is used to fit the test data, and the Hooke elastic modulus, Kelvin elastic modulus, Kelvin viscosity coefficient, and Bingham viscosity coefficient are inverted and identified; combined with the rock bedding direction, transversely isotropic elastic, Poisson coupling, and yield anisotropic parameters are set. S6. Using the state of the previous time step as a known quantity, establish a consistent tangent matrix; establish an incremental equilibrium equation based on the principle of virtual work, and express the total strain increment as the displacement increment to obtain the global incremental equation; assume no yielding at each time step, calculate the trial stress state, and then obtain the trial yield value. When the trial state exceeds the static yield surface, correct the current stress; the trial state determination and viscoplasticity correction process are executed sequentially in each time step until the entire loading process is completed.

2. The method for analyzing the viscoelastic-plastic creep of transversely isotropic rocks according to claim 1, characterized in that, In step S1, the total creep strain increment of the rock is decoupled into three series strain increment components: elastic, viscoelastic, and viscoplastic. in, This represents the total strain increment. , and These represent the elastic strain increment, viscoelastic strain increment, and viscoplastic strain increment, respectively. Indicates increment, subscript n Represents the nth time step; For plane strain and plane stress conditions, transversely isotropic unit compliance matrices are constructed respectively, and elastic anisotropy is characterized by directional Young's modulus, in-plane shear modulus, and Poisson coupling parameters: in, To view the transverse isotropic unit compliance matrix, , , and These are Young's moduli in different directions. It is the in-plane shear modulus. and The parameters are used to describe Poisson coupling in different directions, and the reverse Poisson coupling parameters are determined based on the reciprocal work relationship; Establish the incremental constitutive relation for the Hooke elastic branch: in, Let be the Young's modulus in the y-direction of the Hooke elastic branch. This represents the stress increment; Establish the incremental constitutive relation for the Kelvin viscoelastic branch: in, , , It is the viscoelastic strain rate. It is the stress increment. Let be the Young's modulus in the y-direction of the Kelvin spring. Kelvin viscosity coefficient The length of the nth time step. For generalized midpoint integral parameters; Establish the incremental constitutive relation for the Bingham viscoplastic branch: in It is viscoplastic strain rate. It is the tangent matrix corresponding to the viscoplastic strain rate, and its specific expression is: in, Bingham viscosity coefficient, For Macaulay brackets, It is the overstress function. F Let be the yield function. G Let be the plastic potential function. , and Represent , F and G stress Find the partial derivatives, where the superscript T represents the matrix transpose; The final result is a unified incremental constitutive equation with three branches connected in series: 。 3. The method for analyzing the viscoelastic-plastic creep of transversely isotropic rocks according to claim 2, characterized in that, In step S2, the Drucker-Prager yield function framework is constructed using mean normal stress and equivalent shear stress: in, q For equivalent shear stress, p For average normal stress, a Pressure coefficient; Introducing cohesion and internal friction angle to reflect the pressure-sensitive yielding characteristics of soil and rock: in, The initial yield strength, It is cohesion. It is the internal friction angle; Hill-type orthotropic equivalent stress is used to weight the equivalent shear stress by direction, and the strength directionality is characterized by orthogonal yield and shear yield parameters. in, To view the isotropic equivalent shear stress laterally, It is the transverse isotropic pressure coefficient. The orthogonal yield parameter, This is the shear yield parameter. , and These are the normal stress components in the x, y, and z directions, respectively. , and These are the shear stress components in the yz plane, zx plane, and xy plane, respectively. The transversely isotropic modified Drucker-Prager yield function is obtained as follows: A plastic potential function is established using a homogeneous structure, and a dilatation angle parameter is introduced to adapt to both associated and unassociated flow rules: in, It is the shear expansion angle, when When the correlation flow rule is applied, when The non-associated flow rule is used at this time.

4. The method for analyzing the viscoelastic-plastic creep of transversely isotropic rocks according to claim 3, characterized in that, In step S3, a multi-topology standard polygon mother unit library is constructed in the parameter space, and counterclockwise boundary node numbering is uniformly adopted. By selecting the geometric center as the scaling center, discretizing the boundary into line elements and dividing it into triangular sectors, a radial-circumferential scaling boundary coordinate system is established: in, , It is the coordinate value of any point in the parameter space. It is a radial coordinate. It is a circumferential coordinate system. and These are the coordinates of the scaling center. and These are the boundary coordinates, obtained by interpolation using a one-dimensional linear element shape function: in, It is a shape function matrix. , It is a vector of node coordinates; Solving the Laplace equation, the Euler-Cauchy equation is derived using the weighted residual method: in, It is the radial analytic function of the scalar potential function, subscript Representative to Find the first partial derivative, index Representative to Find the second partial derivative. , and Here is the boundary integral coefficient matrix: in, It is a boundary Jacobian determinant. and Defined as: Among them, subscript Representative to Find the first partial derivative; Transforming the Euler-Cauchy equations into Hamiltonian systems: in, It is the internal flux auxiliary function, Hamiltonian matrix. for: After eigenvalue decomposition, the proportional boundary shape functions of the parent element and their derivative matrices are obtained, and a reusable shape function library is pre-stored. in, It is a proportional boundary shape function matrix. It is a vector of scalar potential function node values. It is the eigenvector matrix. It is an eigenvalue matrix. Represents the eigenvector matrix The extracted local submatrix contains only rows related to the current line cell node. , They represent respectively to and Find the partial derivative. It is an identity matrix.

5. The method for analyzing the viscoelastic-plastic creep of transversely isotropic rocks according to claim 4, characterized in that, In step S4, a polygonal mesh is generated for the actual rock mass computational domain, and the corresponding parent element is called according to the physical polygon node topology matching. Using isoparametric mapping, coordinate interpolation relationships in physical space are established through pre-stored shape functions: in, These are physical space coordinates. , It is the coordinate vector of the physical boundary node; Spatial transformation of the derivatives of shape functions is achieved using the Jacobian matrix: in, , Representative to , Find the partial derivative; the superscript -1 indicates finding the inverse of the matrix. The Jacobian matrix has the following specific form: Derivation of the strain matrix and domain integral expression for the physical unit, thus completing the spatial discretization of the nonlinear problem: in, It is the strain matrix of a polygonal element. It is the number of polygon nodes. It is the differential area of ​​the physical domain. and They represent and The differential, It is the Jacobian determinant, the strain matrix of the polygonal element. block matrix The specific form is as follows: in, , These represent the derivative matrices of the proportional boundary shape functions, respectively. , The i-th component.

6. The method for analyzing the viscoelastic-plastic creep of transversely isotropic rocks according to claim 5, characterized in that, In step S5, the time-strain curves are based on uniaxial or triaxial creep tests: in, It is axial strain. Represents the plastic potential function G right Find the partial derivative. t It is creep time. , , and These represent the first, second, and third principal stresses, respectively. The Nishihara creep expression under constant stress is parameterized and reorganized into a fitted form containing exponential decay terms, constant terms, and linear time terms: in, , , and These are the fitting coefficients; The least squares method was used to fit the experimental data, and the Hooke elastic modulus, Kelvin elastic modulus, Kelvin viscosity coefficient and Bingham viscosity coefficient were identified by inversion. Combined with the rock bedding direction, transversely isotropic elastic, Poisson coupling and yield anisotropy parameters were set.

7. The method for analyzing the viscoelastic-plastic creep of transversely isotropic rocks according to claim 6, characterized in that, In step S6, the state of the previous time step is used as a known quantity to establish a consistent tangent matrix: in, It is a uniform tangent matrix; Incremental equilibrium equations are established based on the principle of virtual work: in, Indicates a time step Internal external force increment vector; The total strain increment is expressed as the displacement increment: in, It is the displacement increment; The global incremental equation is obtained: Among them, the global tangent stiffness matrix Additional viscoelastic force and additional viscoplastic forces The results are obtained by integrating and assembling each physical polygon domain separately; Assuming no yielding at each time step, calculate the trial stress state: in, , , These represent the experimental stress increment, the experimental total strain increment, and the experimental viscoelastic strain rate at the nth time step, respectively. This represents the trial stress state at the (n+1)th time step. The superscript "trial" indicates trial, and the subscript "[...]" indicates trial. n+1 Represents the (n+1)th time step. It is the elastic and viscoelastic tangent stiffness matrix, and its specific expression is: The yield value was then calculated: in, , and These represent the experimental transverse isotropic equivalent shear stress, experimental mean normal stress, and experimental yield function at the (n+1)th time step, respectively. When the calculated state exceeds the static yield surface, the current stress is corrected: in, This represents the stress state at the (n+1)th time step; The trial state determination and viscoplasticity correction processes are executed sequentially at each time step until the entire loading process is completed.

Images