Three-dimensional numerical simulation method for simulating soil-rock mixed shear characteristics based on simulation data

Abstract

This invention discloses a three-dimensional numerical simulation method for the shear properties of soil-rock hybrid materials based on simulation data, belonging to the field of numerical simulation in geotechnical engineering. The method constructs a three-dimensional geometric model containing a randomly polyhedral rock block skeleton embedded with the surrounding soil matrix. An adaptive mesh refinement algorithm is used to discretize the rock block corner regions and the soil-rock interface layer at multiple scales, generating a numerical simulation mesh with contact elements. This mesh is imported into a finite element platform, and dynamic shear loads are applied according to a preset shear rate and normal stress. During the simulation, the normal separation and tangential slip at the contact interface are extracted in real time, and the friction coefficient and cohesion are dynamically adjusted accordingly to form a closed-loop updated numerical simulation result of the soil-rock hybrid shear properties, obtaining the rock block displacement field, soil strain field, and contact interface damage distribution.

Description

Technical Field

[0001] This invention relates to the field of numerical simulation technology in geotechnical engineering, specifically a three-dimensional numerical simulation method for the shear properties of mixed soil and rock based on simulation data. Background Technology

[0002] The shear mechanical behavior of soil-rock mixtures is crucial for slope stability assessment and geological hazard prevention, and numerical simulation is an important means to reveal their failure mechanisms. In existing technologies, simulation methods for the shear characteristics of soil-rock mixtures typically idealize the rock block as a regular geometric shape, discretize it using a uniform mesh, and assign a fixed friction coefficient and cohesion to the soil-rock interface throughout the shearing process. This approach cannot reflect the local stress concentration effect caused by the sharp edges of real irregular rock blocks, nor can it describe the progressive damage caused by the opening and sliding of the contact interface during shearing. Uniform meshes often exhibit mesh distortion and computational divergence when faced with severe deformation gradients near the edges, while fixed contact parameters neglect key mechanical responses such as frictional decay caused by interface separation and loss of cohesion induced by slip accumulation. When the rock block rotates or debonds, the actual bearing characteristics of the contact surface fundamentally change. Simulation results with fixed parameters will overestimate the shear strength and form a distorted one-time brittle fracture failure mode, failing to reproduce the residual strength stage and the gradual evolution of the shear zone. Therefore, how to construct a high-precision and computationally stable discretization scheme for the sharp edges of irregular polyhedral rock blocks, and how to make the friction and cohesion properties of the soil-rock interface continuously evolve with the local normal separation and tangential slip state in shear simulation, have become urgent problems to be solved to improve the realism of the simulation of shear characteristics of soil-rock mixtures. Summary of the Invention

[0003] The purpose of this invention is to provide a three-dimensional numerical simulation method for the shear characteristics of soil-rock mixtures based on simulation data, in order to solve the problems of missing interface damage evolution and insufficient calculation accuracy at sharp corners caused by fixed contact parameters and uniform grids in existing simulation technologies, and to achieve high-fidelity simulation of the entire process of soil-rock mixtures from local interface failure to the formation of overall shear bands.

[0004] To achieve the above objectives, the present invention provides the following technical solution: The present invention provides a three-dimensional numerical simulation method for the shear characteristics of soil-rock mixtures based on simulation data, comprising: constructing a three-dimensional geometric model of a soil-rock mixture by computer-aided design, including a randomly generated polyhedral rock block skeleton embedded with the surrounding soil matrix; based on the three-dimensional geometric model of the soil-rock mixture, performing multi-scale discretization processing on the rock block corner regions and the soil-rock contact interface layer through a mesh adaptive densification algorithm to generate a three-dimensional numerical simulation mesh model with contact elements; importing the three-dimensional numerical simulation mesh model into a finite element simulation platform, and applying dynamic shear loads according to preset shear rates and normal stress boundary conditions to simulate the stress transmission and displacement evolution of the soil-rock mixture during the shearing process; during the simulation, extracting the normal separation and tangential slip at the soil-rock contact interface in real time, and dynamically adjusting the friction coefficient and cohesion of the contact interface according to the normal separation and tangential slip, forming a closed-loop updated numerical simulation result of the shear characteristics of the soil-rock mixture. This method can accurately reproduce the entire process of nonlinear slippage, debonding, and progressive damage of the soil-rock mixture under shear load.

[0005] In a preferred embodiment, when constructing a three-dimensional geometric model of the soil-rock mixture, polyhedral rock blocks are randomly generated in three-dimensional space based on a preset rock block volume fraction and grain size distribution curve. The vertex coordinates and surface normal vectors of each rock block are recorded. The vertex coordinates and surface normal vectors of each rock block are input into a Poisson disk sampling algorithm to generate uniformly distributed contact seed points on the surface of the rock blocks. Using the contact seed points as anchor points, the boundary face of the surrounding soil matrix is ​​offset towards the rock block direction by a preset embedding depth, so that the soil matrix and the surface of the rock block produce overlapping embedding areas, generating a three-dimensional geometric model of soil and rock interlocking. This modeling method realistically reflects the interlocking and interlocking morphology between the rock blocks and the soil, providing a geometric basis for subsequent accurate simulation.

[0006] Furthermore, during the adaptive mesh refinement process, the edges and vertices of each polyhedral rock block in the 3D geometric model are identified, and the dihedral angle of each edge is calculated. When the dihedral angle is less than a preset angle threshold, the edge is marked as an acute edge. A preset layer thickness is then extended outward along the normal direction from the acute edge to generate a local refinement domain containing the edge transition zone. Within this local refinement domain, the side length of the tetrahedral mesh is reduced to one-third of the global mesh side length, and hexahedral interface elements are arranged layer by layer along the normal direction at the soil-rock interface layer. This significantly improves the calculation accuracy of stress concentration and interface contact behavior at the rock block edges without significantly increasing the overall mesh count.

[0007] When applying dynamic shear load, the lower boundary nodes of the three-dimensional numerical simulation mesh model are completely fixed, and a constant normal stress is applied to the upper boundary nodes. A shear displacement that increases linearly with time is applied along the horizontal direction to the upper boundary nodes. The ratio of the shear displacement to the total height of the model in each time step is calculated to obtain the real-time shear strain. When the real-time shear strain reaches a preset peak strain threshold, the shear rate is reduced to one-quarter of the initial shear rate, and shear displacement continues to be applied until the total shear strain reaches a preset termination strain value. This variable-rate loading method can fully capture the structural response characteristics before and after the shear peak softening stage.

[0008] In the stress transfer and displacement evolution simulation, the Gaussian point stress tensor of all tetrahedral mesh elements is calculated at each time step, and the equivalent nodal force of each node is obtained by shape function interpolation based on the Gaussian point stress tensor. The equivalent nodal force of each node is substituted into the global stiffness matrix to solve for the nodal displacement increment, and it is determined whether the normal tensile stress of the soil-rock interface element exceeds the tensile strength of the soil. If the normal tensile stress exceeds the tensile strength of the soil, the interface element is marked as a debonding element, and the tangential stiffness of the debonding element is reduced to zero, thereby simulating the disappearance of tangential bearing capacity after the interface is normally separated.

[0009] When extracting the normal separation and tangential slip at the soil-rock interface in real time, all soil-rock interface elements marked as debonding elements are traversed to obtain the coordinates of the master node on the rock side and the slave node on the soil side of the interface element. The magnitude of the coordinate difference between the master node on the rock side and the slave node on the soil side at the current time step is calculated, and the magnitude of the coordinate difference at the initial time step is subtracted to obtain the normal separation. The change in the coordinate difference between the master node on the rock side and the slave node on the soil side along the shear direction at the current time step is calculated, and the displacement contributed by elastic deformation is subtracted to obtain the tangential slip. Preferably, the displacement contributed by elastic deformation is calculated based on the tangential stress and the initial tangential stiffness at the current time step, thereby separating the irreversible interface slip.

[0010] The process of dynamically adjusting the friction coefficient and cohesion of the contact interface based on the normal separation and tangential slip is as follows: Substitute the normal separation into an exponential decay function to calculate the normal separation decay factor, and multiply the initial friction coefficient by the normal separation decay factor to obtain the friction coefficient at the current time step; divide the tangential slip by a preset characteristic slip length to obtain the tangential cumulative damage variable, and multiply the initial cohesion by the difference between the initial and the tangential cumulative damage variable to obtain the cohesion at the current time step; when the tangential cumulative damage variable is greater than a preset failure threshold, the cohesion at the current time step is forcibly set to zero. Through this closed-loop update method, the contact interface parameters evolve in real time following the separation and slip states, accurately reflecting the inherent mechanism of the interface strength gradually weakening with damage accumulation.

[0011] To accurately describe the directional characteristics of the contact interface, the three-dimensional numerical simulation mesh model with contact elements establishes a local coordinate system for each contact element in the following manner: For each contact element on the soil-rock contact interface, the outward normal direction of the rock block side triangular facet in that contact element is taken as the first axis direction of the local contact coordinate system; the projection direction of the horizontal shear loading direction onto the plane containing the rock block side triangular facet is taken as the second axis direction of the local contact coordinate system; the cross product direction of the first axis direction and the second axis direction is taken as the third axis direction of the local contact coordinate system, and the normal separation, first tangential slip, and second tangential slip are projected onto the first, second, and third axis directions for calculation, respectively. This local coordinate system definition method allows the interface normal separation and the slip along the two tangential directions to be extracted and evaluated independently.

[0012] The closed-loop updated numerical simulation results of the soil-rock hybrid shear characteristics include the rock block displacement field, soil strain field, and damage distribution cloud map of the soil-rock interface. To further quantitatively analyze the simulation results, the following post-processing is performed: The equivalent plastic strain values ​​of all tetrahedral mesh elements within the shear zone are extracted, and elements with equivalent plastic strain values ​​exceeding a preset plastic threshold are marked as failed elements; the total volume of all failed elements is calculated, and the ratio of this total volume to the total volume of the shear zone is calculated to obtain the shear failure volume fraction; a curve of the shear failure volume fraction versus shear displacement is plotted, and the displacement point corresponding to the peak shear strength of the soil-rock hybrid is determined based on the peak point of the second derivative of the shear failure volume fraction in the curve. This post-processing method provides a clear and repeatable quantitative indicator for objectively determining the peak shear strength from the perspective of volumetric damage rate.

[0013] As a preferred technique to ensure numerical stability during simulation, the application of shear displacement also includes adaptive adjustment of the loading step size: after each time step, the maximum norm of the global nodal unbalanced force is calculated. When the maximum norm exceeds a preset convergence threshold, the shear displacement increment of the current time step is halved and recalculated. When the maximum norm of the global nodal unbalanced force for three consecutive time steps is less than one-tenth of the convergence threshold, the shear displacement increment of the next time step is doubled. The ratio of the actual shear displacement increment to the preset initial shear displacement increment is recorded as the loading efficiency coefficient for that time step. This adaptive step size adjustment effectively avoids redundant calculations caused by a fixed small step size while ensuring computational convergence accuracy, significantly improving the computational efficiency of the entire dynamic shear simulation process.

[0014] The technical effects and advantages provided by the present invention in the above technical solution are as follows:

[0015] By identifying the dihedral angles of the edges of each polyhedral rock block in the geometric model of the soil-rock mixture, acute edges with angles less than a threshold are marked. A pre-defined layer thickness is then extended outwards along the normal direction from these acute edges to generate a local refinement domain for the edge transition zone. Within this local refinement domain, the side length of the tetrahedral mesh is forcibly reduced to one-third of the global mesh side length. Simultaneously, hexahedral interface elements are layered along the normal direction at the soil-rock interface layer. This multi-scale discretization process allows for refined analysis of high stress gradients and drastic displacement changes near the edges. The layered hexahedral interface elements ensure the numerical stability of stress transmission in the normal and tangential directions at the contact surface, avoiding computational interruptions or result oscillations caused by severe element distortion at sharp edges in traditional uniform meshes. This enables accurate capture of the local failure initiation point triggered by rock block rotation and stress concentration in the early stages of shearing. The coordinate difference between the master node on the rock side and the slave node on the soil side of each contact unit is extracted in real time. The normal separation amount and the tangential slip amount after deducting the contribution of elastic deformation are calculated. The normal separation amount is substituted into the exponential decay function to obtain the normal separation decay factor. The friction coefficient at the current moment is obtained by multiplying the factor by the initial friction coefficient. At the same time, the tangential slip amount is divided by the characteristic slip length to form the tangential cumulative damage variable. The cohesion at the current moment is obtained by multiplying the initial cohesion by a factor minus the damage variable. When the damage variable exceeds the failure threshold, the cohesion is forcibly fixed to zero. By adopting this real-time parameter update mechanism that is directly related to the separation amount and slip amount, the mechanical properties of the contact interface can smoothly and continuously degrade with the local separation and sliding process during the simulation. The process from the initial bonded state through progressive damage to complete frictional sliding is naturally presented. This gives the nonlinear strength response that varies with shear displacement and the actual shear band expansion path, eliminating the erroneous estimation of residual strength after the peak value and the distortion of the fracture mode under fixed parameter settings. Attached Figure Description

[0016] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments recorded in this invention. For those skilled in the art, other drawings can be obtained based on these drawings.

[0017] Figure 1 This is a flowchart of a three-dimensional numerical simulation method for the shear properties of mixed soil and rock;

[0018] Figure 2 This is a flowchart of the construction process for a 3D numerical simulation mesh model of a soil-rock mixture;

[0019] Figure 3 This is a flowchart of the adaptive loading and step size adjustment process for numerical simulation of direct shear of soil-rock mixture;

[0020] Figure 4This is a flowchart of the soil-rock interface debonding determination and tangential stiffness reduction process.

[0021] Figure 5 This is a flowchart of the extraction of separation amount and dynamic adjustment of contact parameters in the debonding unit at the soil-rock contact interface;

[0022] Figure 6 This is a flowchart for calculating the shear failure volume fraction and determining the peak strength of soil-rock mixtures;

[0023] Figure 7 This is a comparison diagram of the equivalent grain size distribution of polyhedral rock blocks based on grain size gradation curves;

[0024] Figure 8 It is the curve of the loading efficiency coefficient changing with time step during shear loading;

[0025] Figure 9 This is a diagram showing the distribution of normal tensile stress in the debonding unit at the soil-rock interface;

[0026] Figure 10 These are curves showing the relationship between the friction coefficient and cohesion at the soil-rock interface as a function of shear displacement.

[0027] Figure 11 It is the curve of shear displacement corresponding to shear failure volume fraction and its second derivative. Detailed Implementation

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

[0029] See Figure 1 This invention provides a three-dimensional numerical simulation method for the shear characteristics of soil-rock mixtures based on simulation data. The method constructs a three-dimensional geometric model containing a randomly multifaceted rock block skeleton and an embedded soil matrix, and adaptively refines the mesh to generate a three-dimensional numerical simulation mesh model with contact elements. The mesh model is then imported into a finite element simulation platform, and dynamic shear loads are applied according to preset shear rates and normal stress boundary conditions to simulate stress transfer and displacement evolution of the soil-rock mixture during shearing. During the simulation, the normal separation and tangential slip at the soil-rock contact interface are extracted in real time, and the friction coefficient and cohesion of the contact interface are dynamically adjusted to form a closed-loop updated numerical simulation result of the soil-rock mixture shear characteristics.

[0030] Example 1: In specific implementation, please refer to Figure 2 The process of constructing a three-dimensional geometric model of a soil-rock mixture, in which a randomly generated polyhedral rock block skeleton is embedded in the surrounding soil matrix, is carried out in stages. Pre-defined rock block volume fraction and particle size distribution curves are obtained. The rock block volume fraction represents the percentage of polyhedral rock blocks in the total volume of the soil-rock mixture. The particle size distribution curve is determined based on the actual soil-rock mixture sieving test results. The horizontal axis of the curve represents the logarithm of the particle size, and the vertical axis represents the cumulative percentage content of particles smaller than a certain size. In three-dimensional space, a random method is used to generate polyhedral rock blocks. The specific steps are as follows: Within the specified cuboid region of the model boundary, the equivalent grain size of the polyhedral rock block to be generated is determined based on the grain size distribution curve. The equivalent grain size is taken as the diameter of a sphere with the same volume as the polyhedral rock block. Then, the total number of polyhedral rock blocks required is calculated based on the volume fraction of the rock blocks. The initial center position of each polyhedral rock block is randomly generated, and the number of vertices of the polyhedral rock block is also randomly generated, with the number of vertices ranging from 6 to 20. Using the initial center position as a reference, the spatial offset of each vertex is randomly generated. The offset is no greater than half of the equivalent grain size, so that the polyhedral rock block presents an irregular convex polyhedral shape. Collision detection is performed on the generated initial polyhedral rock blocks. If there is an overlap with the already generated polyhedral rock blocks, the center position and vertex offset are regenerated until the non-overlap condition is met. The above process is repeated until all polyhedral rock blocks are generated. Record the vertex coordinates of each polyhedral rock block and the face normal vectors of each triangular face in each polyhedral rock block. The face normal vectors of the triangular faces all point to the outside of the polyhedral rock block.

[0031] The vertex coordinates and surface normal vectors of each polyhedral rock block are input into the Poisson disk sampling algorithm to generate uniformly distributed contact seed points on the surface of the polyhedral rock block. The specific implementation process of the Poisson disk sampling algorithm is as follows: Define the minimum distance parameter between contact seed points. The minimum distance parameter is set according to the minimum finite element mesh size in the soil matrix, and is taken as 0.8 times the minimum finite element mesh size. Randomly select an initial point from the surface of the polyhedral rock block as the first contact seed point and add it to the activity list. When the activity list is not empty, randomly select a current seed point from the activity list and generate several candidate points around the current seed point. The candidate points are located in the intersection area of ​​the spherical shell range with the current seed point as the center, the minimum distance parameter as the inner radius, and twice the minimum distance parameter as the outer radius and the surface of the polyhedral rock block. Check the distance between the candidate point and all confirmed contact seed points one by one. If the distance between a candidate point and any confirmed contact seed point is greater than the minimum distance parameter, then add the candidate point as a contact seed point and add it to the activity list. If the candidate points do not meet the conditions, then remove the current seed point from the activity list. Iterate continuously until the activity list is empty, thereby obtaining a uniformly distributed set of contact seed points on the surface of the polyhedral rock block.

[0032] Using the generated contact seed points as anchor points, the boundary of the surrounding soil matrix is ​​offset towards the polyhedral rock block direction by a preset embedding depth, creating an overlapping embedding zone between the soil matrix and the polyhedral rock block surface. In the computer-aided design environment, the initial boundary of the soil matrix is ​​defined by the outer surface of a cuboid or cylinder that completely surrounds all polyhedral rock blocks. Each contact seed point is projected onto the initial boundary of the soil matrix to obtain the corresponding initial soil boundary node. The initial soil boundary node is then moved along the inner normal direction of the polyhedral rock block surface at that node by an embedding depth value. The embedding depth is set to 0.15 times the median particle size of the equivalent particle size fractal dimension distribution of the polyhedral rock block, and not less than 0.1 times the smallest unit size of the soil matrix. The new soil matrix boundary surface formed after the movement constitutes an interpenetrating embedding zone with the polyhedral rock block surface. Boolean operations are used to attribute the volume of the embedding zone to the overlapping portion of both the soil matrix and the rock block, generating a three-dimensional geometric model of the interpenetrating soil and rock.

[0033] After obtaining the aforementioned three-dimensional geometric model of the soil-rock mixture, a multi-scale discretization process was performed on the rock block corner regions and the soil-rock interface layer using an adaptive mesh refinement algorithm, generating a three-dimensional numerical simulation mesh model with contact elements. First, acute-angle edge identification was performed: traversing every edge of all polyhedral rock blocks in the three-dimensional geometric model, for any edge, the face normal vectors of the two triangular faces sharing that edge were obtained, and the angle between the two face normal vectors was calculated and recorded as the dihedral angle of that edge. When the dihedral angle of an edge is less than a preset angle threshold, that edge is marked as an acute-angle edge. The angle threshold was set to 75 degrees, based on the fact that for most rock materials, the local stress concentration coefficient increases sharply when the dihedral angle is less than 75 degrees, requiring mesh refinement to capture high stress gradients.

[0034] Centered on each acute-angled edge, a predetermined layer thickness is extended outward along the outward normal direction of the triangular facets on both sides of the acute-angled edge, generating a local consolidation domain containing the transition zone of the acute angle. The layer thickness consists of two parts: one part is a function of the characteristic length of the polyhedral rock block containing the acute-angled edge, where the characteristic length is taken as the radius of the circumscribed sphere of the polyhedral rock block, and the basic component of the layer thickness is taken as 0.05 times the characteristic length; the other part is related to the total thickness of the hexahedral interface unit layer of the soil-rock interface layer, and is taken as 1.5 times the total thickness of the hexahedral interface unit layer. The final layer thickness is obtained by adding these two components. Centered on each point on the acute-angled edge and with the layer thickness as the extension distance, a ring-shaped prism-shaped region surrounding the acute-angled edge is constructed. The part of this region that intersects with the soil-rock interface layer is defined as the local consolidation domain.

[0035] Within the localized encrypted domain, mesh size control is implemented, reducing the edge length of the tetrahedral mesh to one-third of the global mesh edge length. The global mesh edge length is taken as one-tenth of the total dimension of the 3D geometric model of the soil-rock mixture along the shear direction. Hexahedral interface elements are layered along the normal direction at the soil-rock contact interface layer, with a fixed number of three layers, each with equal thickness, and the total thickness is taken as 0.3 times the minimum edge length of the local tetrahedral mesh. The hexahedral interface elements are in the form of zero-thickness or finite-thickness interface elements, connecting the hexahedral element surfaces of the polyhedral rock block side and the soil matrix side. Through the above processing, the final 3D numerical simulation mesh model with contact elements is generated.

[0036] See Figure 7 In the figure, the horizontal axis represents the particle size in millimeters, using a logarithmic scale, ranging from approximately 1e-2 mm (0.01 mm) to 1e2 mm (100 mm). The vertical axis represents the cumulative percentage content of a certain particle at the corresponding particle size, in percentages, ranging from 0% to 100%. The curve uses a dashed line to represent the target gradation curve, and the scattered points represent the distribution of the equivalent particle size of the polyhedral rock blocks generated in Example 1.

[0037] The scattered dots in the figure show a trend of increasing cumulative percentage content with increasing particle size. The overall curve shape is highly consistent with the target gradation curve, indicating that the method of randomly generating polyhedral rock blocks described in Example 1 can accurately simulate the preset particle size gradation characteristics. Specifically, in the particle size range of approximately 0.01 mm to 0.1 mm, the cumulative percentage content increases slowly, reflecting a small proportion of fine particles; in the particle size range of approximately 0.1 mm to 10 mm, the cumulative percentage content rises rapidly, indicating that particles in this particle size range account for the main proportion of the total volume; when the particle size exceeds 10 mm, the cumulative percentage content approaches 100%, corresponding to the maximum particle size of the particle volume distribution.

[0038] This figure compares the equivalent grain size distribution of the generated rock blocks with the target gradation curve, demonstrating the accuracy and effectiveness of the random generation algorithm in Example 1, which calculates the number and size of polyhedral rock blocks based on the grain size gradation curve and the volume fraction of the rock blocks. This distribution forms the basis for subsequent 3D geometric model construction and finite element mesh generation, ensuring the realism and representativeness of the grain size distribution in the numerical simulation model, and providing reliable geometric input data for simulating the shear mechanical behavior of soil-rock mixtures.

[0039] Example 2: In specific implementation, please refer to Figure 3The three-dimensional numerical simulation mesh model with contact elements, generated using an adaptive mesh refinement algorithm, is imported into the finite element simulation platform. The finite element simulation platform uses general-purpose finite element software with explicit dynamics and implicit statics solving capabilities. During the import process, the continuity of node numbers for tetrahedral and hexahedral interface elements, as well as the completeness of element type definitions, are checked. The lower boundary nodes of the three-dimensional numerical simulation mesh model are completely fixed by constraining the translational degrees of freedom of all nodes along the three orthogonal coordinate axes, but not constraining the rotational degrees of freedom. A constant normal stress is applied to the upper boundary nodes of the three-dimensional numerical simulation mesh model. The magnitude of the normal stress is set based on the actual ground stress value under the simulated working condition; the direction of the applied normal stress is perpendicular to the upper boundary plane and remains unchanged throughout the simulation, while its magnitude remains constant.

[0040] A shear displacement, increasing linearly with time, is applied horizontally at the upper boundary node, perpendicular to the normal stress direction. This shear displacement is achieved by applying horizontal displacement boundary conditions at the upper boundary node. The ratio of the shear displacement to the total model height is calculated at each time step to obtain the real-time shear strain. The total model height is the vertical distance from the lower boundary plane to the upper boundary plane of the 3D numerical simulation mesh model, and the shear displacement is the displacement component of the upper boundary node in the shear direction at the current time step. The real-time shear strain is calculated as follows:

[0041]

[0042] in, Indicates the first Real-time shear strain at the end of each time step, dimensionless quantity; Indicates the first Shear displacement at the end of each time step, in millimeters; This represents the total height of the model, in millimeters. Shear displacement. This represents the cumulative displacement of the boundary nodes in the shear direction at the current time step, increasing from zero with loading time; the total height of the model. The original dimensions of the three-dimensional geometric model of the soil-rock mixture along the normal direction.

[0043] When real-time shear strain When the preset peak strain threshold is reached, the shear rate is reduced to one-quarter of the initial shear rate. The peak strain threshold is set to 0.15, based on the median of the shear failure strain range obtained from numerous direct shear tests and numerical experiments on soil-rock mixtures. A value of 0.15 ensures that the peak strength remains within this strain range. The initial shear rate is determined according to the quasi-static loading principle, using one ten-thousandth of the total model height as the shear displacement increment at each time step, making the inertial effect of the entire shear process negligible. The shear displacement increment at each time step corresponding to the reduced shear rate is one-quarter of the initial shear displacement increment. Shear displacement continues to be applied until the total shear strain reaches the preset termination strain value, set to 0.30. This value is based on the fact that the residual strength of soil-rock mixtures typically stabilizes in the strain range of 0.25 to 0.35. A value of 0.30 ensures complete capture of the mechanical behavior of softening to the residual stage after the peak strain.

[0044] The entire process of applying shear displacement also includes adjusting the loading step size. After each time step, the maximum norm of the global nodal unbalanced force is extracted. The global nodal unbalanced force refers to the residual force vector formed by the difference between the internal stress equivalent nodal force and the externally applied nodal force. The maximum norm is the maximum absolute value of all components in the residual force vector. When the maximum norm of the global nodal unbalanced force exceeds a preset convergence threshold, the shear displacement increment of the current time step is halved and the current time step is recalculated. The convergence threshold is set to one-thousandth of the total externally applied load, which is the product of the normal stress and the upper boundary area plus the sum of the equivalent nodal forces corresponding to the current shear displacement. After halving the shear displacement increment, the iteration solution of the current time step is performed again, and the maximum norm of the global nodal unbalanced force is checked again. If it still exceeds the convergence threshold, it is halved again, with a maximum of five halvings per time step. If convergence is still not achieved after five halvings, the calculation is terminated and an error message is output.

[0045] When the maximum norm of the global nodal unbalanced forces is less than one-tenth of the convergence threshold for three consecutive time steps, the shear displacement increment for the next time step is doubled. The doubled shear displacement increment shall not exceed twice the initial shear displacement increment to avoid excessively large shear displacement increments in a single time step, which could lead to a decrease in solution accuracy. If the maximum norm of the global nodal unbalanced forces exceeds the convergence threshold again in any time step, the doubling operation is immediately stopped, and the shear displacement increment is restored to that of the previous time step.

[0046] During the simulation, the ratio of the actual shear displacement increment to the preset initial shear displacement increment at each time step is recorded and used as the loading efficiency coefficient for that time step. The loading efficiency coefficient reflects the magnitude of the actual load at each time step relative to the reference load. The loading efficiency coefficient for each time step is denoted as . , The expression is ,in, Indicates the first The actual shear displacement increment at each time step, in millimeters; This represents the initial shear displacement increment, in millimeters, and is taken as one ten-thousandth of the total model height. Loading efficiency coefficient. The loading efficiency coefficients of all time steps are adjusted and arranged in the order of time steps to form a loading efficiency coefficient curve, which is used for subsequent analysis of the stability of the solution process and identification of periods of drastic changes in contact state.

[0047] See Figure 8 In the graph, the horizontal axis represents the time step number, and the vertical axis represents the loading efficiency coefficient. The curve represents the change in the ratio of the actual shear displacement increment to the initial shear displacement increment at each time step. According to Example 2, the loading efficiency coefficient reflects the dynamic adjustment state of the loading step size of the shear displacement boundary conditions during the finite element simulation.

[0048] The loading efficiency coefficient remained around 1.0 for most of the time, indicating that the loading step size remained unchanged and the simulation ran stably. At a few specific time steps, the loading efficiency coefficient showed significant jumps, reaching approximately 2.0 or decreasing to approximately 0.5, reflecting the doubling or halving of the loading step size. Specifically: a sudden increase in the loading efficiency coefficient from 1.0 to 2.0 indicates that the maximum norm of the global node unbalanced force was below one-tenth of the convergence threshold for several consecutive time steps, triggering a doubling of the loading step size, and the doubled step size did not exceed the initially set maximum limit of twice; a sharp drop in the loading efficiency coefficient from 1.0 to 0.5 indicates that the maximum norm of the global node unbalanced force exceeded the convergence threshold in the current time step, leading to a halving of the loading step size and a re-iteration of the solution.

[0049] Throughout the loading process, the loading efficiency coefficient fluctuated repeatedly around 1.0, accompanied by sporadic doubling and halving operations. This reflects the adaptive adjustment mechanism of the loading step size based on the global node unbalanced force residual, effectively balancing computational efficiency and solution accuracy. The figure clearly illustrates the dynamic adjustment strategy of the loading step size during the shear loading process in Example 2, verifying that the finite element simulation platform can reasonably adjust the loading step size under the preset convergence threshold, ensuring the convergence and computational efficiency of the shear simulation process.

[0050] Example 3: In specific implementation, please refer to Figure 4The finite element simulation platform performs stress calculations and nodal force solutions at each time step. For all tetrahedral mesh elements in the 3D numerical simulation mesh model, it iterates through all Gaussian integration points within each tetrahedral mesh element. A complete integration scheme is used for the number of Gaussian integration points: one Gaussian integration point for a four-node linear tetrahedral element and five Gaussian integration points for a ten-node quadratic tetrahedral element. At each Gaussian integration point, the strain increment is calculated based on the nodal displacements of the current time step. This strain increment is then converted into a stress increment using the constitutive model and superimposed with the stress tensor from the previous convergence time step to obtain the Gaussian stress tensor for the current time step. The Gaussian stress tensor is stored as a second-order tensor containing six independent stress components.

[0051] After obtaining the stress tensors at each Gaussian point within all tetrahedral mesh elements, the equivalent nodal forces at each node are obtained through shape function interpolation based on the Gaussian point stress tensors. For any tetrahedral mesh element, the equivalent nodal force component generated by the tetrahedral mesh element for one of its nodes is calculated according to the principle of virtual work. Specifically, the stress tensors at all Gaussian integration points within the tetrahedral mesh element are multiplied by the transpose of the strain-displacement matrix at that Gaussian integration point, then multiplied by the integration weight coefficient corresponding to that Gaussian integration point and the volume of the tetrahedral mesh element, and finally summed over all Gaussian integration points. The equivalent nodal force contributions of all tetrahedral mesh elements in the 3D numerical simulation mesh model to the shared node are accumulated to obtain the complete equivalent nodal force vector of the shared node.

[0052] The equivalent nodal forces at each node are substituted into the global stiffness matrix to solve for the nodal displacement increments. The global stiffness matrix is ​​assembled from the element stiffness matrices of all tetrahedral mesh elements and hexahedral interface elements in the 3D numerical simulation mesh model. The element stiffness matrix is ​​constructed using the tangent modulus matrix of the element material at the current time step; the soil matrix uses the elasto-plastic tangent modulus matrix, and the polyhedral rock block uses the linear elastic tangent modulus matrix. The overall equilibrium equations after assembly are:

[0053]

[0054] in, Indicates the first The global stiffness matrix after assembly at each time step is integrated from the tangent modulus matrix of each element, and the unit is Newtons per millimeter. Indicates the first The nodal displacement increment vector to be solved at each time step, in millimeters; Indicates the first The external nodal force vector applied to the node at each time step is obtained by converting the normal stress and shear displacement boundary conditions, and the unit is Newton; Indicates the first The internal nodal force vector at each time step is derived from the stress field equivalent, i.e., the vector assembled from the equivalent nodal forces of all the aforementioned nodes, in Newtons. The overall equilibrium equations are solved using a sparse direct solver or a preconditional conjugate gradient iterative solver to obtain the converged nodal displacement increment vector at the current time step.

[0055] After obtaining the nodal displacement increments and updating the nodal coordinates, it is determined whether the normal tensile stress of the soil-rock interface element exceeds the tensile strength of the soil. The soil-rock interface element is a hexahedral interface element, which has two faces, top and bottom, connecting the soil matrix element and the polyhedral rock block element, respectively. The stress vector at the Gaussian integration point of each hexahedral interface element is extracted, and the stress vector is projected into a normal stress component and two tangential stress components using the local coordinate system of the hexahedral interface element. The normal stress component is positive with tension. When the normal stress component value at any Gaussian integration point of a hexahedral interface element exceeds the tensile strength of the soil, it is determined that the normal tensile stress of that hexahedral interface element exceeds the tensile strength of the soil. The tensile strength of the soil is taken as one-tenth of the zero point of the soil matrix cohesion, and the tensile strength decreases synchronously as the soil matrix cohesion softens with plastic strain.

[0056] If the normal tensile stress of a hexahedral interface element exceeds the tensile strength of the soil, the hexahedral interface element is marked as a debonded element. The debonded element is identified using an integer state variable, initially set to zero to indicate no debonding; after being marked as a debonded element, the state variable is set to 1. Simultaneously with marking the hexahedral interface element as a debonded element, its tangential stiffness is reduced to zero. Tangential stiffness refers to the stiffness coefficients in the two tangential directions within the local coordinate system of the hexahedral interface element. Reducing it to zero involves setting the stiffness coefficients of all associated tangential displacement degrees of freedom in the element stiffness matrix of the hexahedral interface element to zero, retaining only the normal stiffness to prevent normal embedding. In subsequent time steps, hexahedral interface elements marked as debonded elements with their tangential stiffness reduced to zero no longer transmit tangential stress, but the normal direction still participates in the contact calculation to handle possible normal closure cases.

[0057] See Figure 9 In the figure, the horizontal axis represents the contact element number, and the vertical axis represents the normal stress minus the tensile strength of the debonded soil, in kilopascals. Black dots represent the distribution of normal stress values ​​in undebonded elements, black crosses represent the distribution of normal stress values ​​in debonded elements, and the dashed horizontal line represents the reference value of the soil tensile strength (approximately 2 kPa), i.e., the debonding threshold.

[0058] As can be seen from the figure, the normal stress values ​​of most unbonded elements are concentrated in the range below or close to zero, and none exceed the reference value of soil tensile strength, indicating that these elements are in a normal contact state and cannot fail tensilely. The normal stress values ​​of a few unbonded elements are close to but do not exceed the threshold, showing the fluctuation of normal stress at the contact interface.

[0059] The normal stress values ​​of the debonded elements are all significantly higher than the tensile strength threshold of the soil, and the distribution range is large, with the maximum exceeding 4 kPa, indicating that these elements have undergone normal tensile stress failure and meet the debonding criterion. After the state variable of the debonded element changes from 0 to 1, its corresponding normal tensile stress exceeds the tensile strength of the soil, and the tangential stiffness is simultaneously reduced to zero, which is used to simulate the interface debonding phenomenon.

[0060] The numerical distribution of normal stress minus the tensile strength of the soil in the figure intuitively reflects the effective capture of the failure state of the soil-rock contact interface by the three-dimensional numerical simulation method described in Example 3. This figure shows the contact element state determination results after convergence at each time step, verifying the accuracy of the stress calculation and contact element debonding determination algorithm. It achieves the determination of debonding elements and the closed-loop update of tangential stiffness reduction when the normal tensile stress at the contact interface exceeds the tensile strength of the soil, meeting the technical requirements for debonding element identification and processing in Example 3.

[0061] Example 4: In specific implementation, please refer to Figure 5 The process of extracting the normal separation and tangential slip at the soil-rock interface in real time is executed after each convergence time step. All soil-rock interface elements marked as debonding elements are traversed, and the coordinates of the rock block-side master node and the soil-side slave node are obtained for each debonding element. The rock block-side master node refers to the corner node of the hexahedral interface element connected to the polyhedral rock block surface mesh in the debonding element, and the soil-side slave node refers to the corresponding corner node connected to the soil matrix mesh in the debonding element. The two nodes have the same spatial coordinates in the initial undeformed state, forming a contact node pair. When the debonding element separates due to shear or tension, the spatial coordinates of the rock block-side master node and the soil-side slave node no longer coincide. The degree of separation is quantified by comparing the differences between the rock block-side master node coordinates and the soil-side slave node coordinates.

[0062] Calculate the magnitude of the coordinate difference between the rock block-side master node and the soil-side slave node at the current time step, which is equivalent to calculating the Euclidean length of the difference between the three-dimensional spatial position vectors of the two nodes at the current time step. Obtain the initial time step coordinate difference, which corresponds to the initial configuration before the application of shear displacement. In this initial configuration, the rock block-side master node and the soil-side slave node are completely coincident, and the magnitude of the initial time step coordinate difference is zero. Subtract the magnitude of the initial time step coordinate difference from the magnitude of the current time step coordinate difference to obtain the normal separation amount. The normal separation amount represents the opening displacement of the debonded element in the normal direction of the contact interface due to tension.

[0063] Calculate the change in coordinate difference between the rock block-side master node and the soil-side slave node along the shear direction at the current time step. The shear direction is the direction of the horizontal shear displacement applied in the boundary conditions. Project the displacement vectors of the rock block-side master node and the soil-side slave node onto the unit vector of the shear direction, respectively, to obtain two scalar displacement values. Subtract these two scalar displacement values ​​to obtain the change in coordinate difference along the shear direction at the current time step. Then, subtract the displacement contributed by elastic deformation from the change in coordinate difference along the shear direction at the current time step to obtain the tangential slip. The displacement contributed by elastic deformation is calculated based on the tangential stress and initial tangential stiffness at the current time step, using the following calculation method:

[0064]

[0065] in, Indicates the first The displacement contributed by elastic deformation at each time step, in millimeters; Indicates the first The tangential stress along the shear direction in the local coordinate system of the debonding unit at each time step is directly extracted from the stress calculation output of the current time step, and the unit is megapascal; This represents the initial tangential stiffness, expressed in megapascals per millimeter. Take the elastic tangential stiffness value of the debonding unit before damage. This elastic tangential stiffness value is equal to the tangential stiffness modulus defined in the material parameters of the hexahedral interface element of the debonding unit divided by the characteristic thickness of the debonding unit. The tangential stiffness modulus is taken as ten times the elastic shear modulus of the soil matrix, and the characteristic thickness is taken as the side length corresponding to the normal direction in the local coordinate system of the debonding unit.

[0066] After extracting the normal separation and tangential slip, the friction coefficient and cohesion of the contact interface are dynamically adjusted based on these values. This dynamic adjustment process is performed before the stress update phase of each time step to ensure that the friction coefficient and cohesion used in the current time step reflect the latest contact state.

[0067] Substituting the normal separation amount into the exponential decay function, the normal separation decay factor is calculated. The exponential decay function is in the form that the normal separation decay factor is equal to the ratio of the negative normal separation amount to the reference separation amount, where the reference separation amount is taken as 0.05 times the characteristic particle size of the soil matrix, and the characteristic particle size of the soil matrix is ​​taken as the median particle size on the soil matrix particle size distribution curve. Multiplying the initial friction coefficient by the normal separation decay factor yields the friction coefficient at the current time step. The initial friction coefficient is obtained from a laboratory direct shear test between the soil matrix and the polyhedral rock block under non-separation conditions, and the friction coefficient ranges from 0.3 to 0.8.

[0068] The tangential slip is divided by the preset characteristic slip length to obtain the tangential cumulative damage variable. The characteristic slip length is set to 0.1 times the minimum equivalent particle size in the polyhedral rock block. The minimum equivalent particle size is extracted from the particle size distribution curve of the polyhedral rock block, taking the equivalent diameter corresponding to the minimum sieve particle size. The initial cohesion is multiplied by the difference between the initial and the tangential cumulative damage variable to obtain the cohesion at the current time step. The initial cohesion is the laboratory-measured cohesion value at the interface between the soil matrix and the polyhedral rock block in the non-slip state, and the unit of cohesion is kilopascals.

[0069] When the cumulative tangential damage variable exceeds a preset failure threshold, the cohesion at the current time step is forcibly set to zero. The failure threshold is set to 0.95, based on the condition that the tangential slip amount corresponding to the complete loss of cohesion reaches 0.95 times the characteristic slip length, indicating that sufficiently large plastic sliding has occurred at the contact interface, causing the bonding effect to completely disappear. After the cohesion is forcibly set to zero, the cohesion of this debonding unit remains zero in subsequent time steps and does not recover.

[0070] See Figure 10 In the figure, the horizontal axis represents shear displacement in millimeters, indicating the cumulative displacement of the upper boundary node of the soil-rock mixture in the shear direction; the left vertical axis corresponds to the solid curve, representing the friction coefficient, ranging from 0 to 0.8; the right vertical axis corresponds to the dashed curve, representing cohesion in kilopascals, ranging from 0 to 60. This figure reflects the dynamic evolution of the friction coefficient and cohesion at the soil-rock interface as the shear displacement increases, closely following the specific method in Example 4 for closed-loop adjustment of the contact interface parameters based on normal separation and tangential slip.

[0071] The friction coefficient curve in the figure shows a monotonically decreasing trend. The initial friction coefficient is approximately 0.75, and as the shear displacement gradually increases from 0 to 30 mm, the friction coefficient gradually decreases to approximately 0.27. This change is consistent with the adjustment mechanism in Example 4, which calculates the friction coefficient by substituting the normal separation amount into the exponential decay function. That is, the normal separation amount increases with the increase of shear displacement, leading to an exponential decay of the friction coefficient, reflecting the physical process of the gradual weakening of the friction capacity of the soil-rock interface due to tensile cracking and shear slip.

[0072] The cohesion curve is represented by a dashed line. The initial cohesion is approximately 50 kPa, and it decreases non-linearly with increasing shear displacement, dropping to zero at approximately 18 mm shear displacement and remaining at zero. This trend is consistent with the process in Example 4 where the tangential cumulative damage variable was calculated based on the tangential slip amount, and the cohesion was gradually reduced using the cumulative shear slip damage as the damage index. The cohesion dropping to zero corresponds to sufficiently large plastic slip at the contact interface. When the preset failure threshold of 0.95 is reached, it indicates that the bonding effect at the contact interface is completely lost, and the slip damage has reached the failure stage.

[0073] Combining the trends of the two curves, it can be seen that when the shear displacement is small, the contact interface maintains a high coefficient of friction and cohesion, and the interface is in a strong bonding state. As the shear displacement increases, the interface gradually undergoes bond failure and slippage, the cohesion decreases rapidly, and the decrease in the coefficient of friction is relatively gradual, reflecting the characteristic that friction and cohesion synergistically control the shear strength of the soil-rock contact interface. After the cohesion reaches zero, the interface is dominated by the friction mechanism for shear behavior, and the continued slow decrease in the coefficient of friction indicates that the normal opening of the contact interface continuously affects the frictional characteristics.

[0074] This figure accurately reflects the technical effect of extracting normal separation and tangential slip in real time in Example 4, dynamically adjusting the friction coefficient and cohesion of the contact interface, and forming a closed-loop update of the numerical simulation results of soil-rock mixed shear characteristics. It verifies the rationality and effectiveness of the evolution of the mechanical parameters of the contact interface corresponding to the exponential decay and cumulative damage model.

[0075] Example 5: In specific implementation, please refer to Figure 6 The process of establishing the local coordinate system for the contact element is completed during the generation of the 3D numerical simulation mesh model with contact elements, providing a unified local coordinate reference for subsequent calculations. For each contact element on the soil-rock interface, the coordinates of three nodes on the rock block side triangle are obtained, and the outward normal direction of the rock block side triangle is determined according to the right-hand rule. The method for calculating the outward normal direction is as follows: take two edge vectors of the rock block side triangle, calculate the cross product of the two edge vectors, obtain the cross product vector, and normalize the cross product vector to obtain the outward normal direction. In a single contact element, if the rock block side triangle is composed of nodes A, B, and C, then the two edge vectors are the edge vector pointing from node A to node B and the edge vector pointing from node A to node C, respectively. The direction of the cross product vector is determined by the right-hand rule order of these two edge vectors, ensuring that the outward normal direction points to the outside of the polyhedral rock block. This outward normal direction is taken as the first axis direction of the contact local coordinate system, and the first axis direction is denoted as axis. ,axis The reference axis is the normal direction of the contact element.

[0076] The projection of the horizontal shear loading direction onto the plane containing the triangular facet of the rock block is taken as the second axis direction of the local contact coordinate system. The horizontal shear loading direction is the horizontal unit vector of the shear displacement applied in the boundary conditions. Projecting the horizontal shear loading direction unit vector onto the plane containing the triangular facet of the rock block involves calculating the horizontal shear loading direction unit vector and its relationship to the first axis direction. The dot product of the first axis and the horizontal shear loading direction is then subtracted from the unit vector of that dot product. The product of these two values ​​yields the projection vector of the horizontal shear loading direction onto the plane containing the triangular facet of the rock block. This projection vector, after being normalized, is used as the second axis direction of the local contact coordinate system, denoted as axis. ,axis This is the reference axis in the first tangential direction of the contact unit. When the horizontal shear loading direction is perpendicular to the first tangential direction axis... When parallelism causes the projection vector magnitude to be less than the preset minimum tolerance value, the projection of the principal axis direction of the polyhedral rock block onto the side triangular face of the rock block is used as the second axis direction. The preset minimum tolerance value is taken as .

[0077] The cross product of the first and second axis directions is taken as the third axis direction of the contact local coordinate system. The first axis direction is then... The unit vector and the second axis direction axis Perform a cross product operation on the unit vectors, and the result yields the direction of the third axis, which is denoted as the axis. ,axis This is the reference axis in the second tangential direction of the contact unit. At the same time perpendicular to the axis and shaft This forms a right-handed orthogonal local coordinate system. At this point, the local coordinate system of the contact element is established, with three coordinate axes. ,axis ,axis All are unit vectors and are pairwise orthogonal.

[0078] The normal separation, first tangential slip, and second tangential slip are projected onto the first, second, and third axes for calculation. The normal separation is a scalar, and its absolute value represents the distance between the master node on the rock block side and the slave node on the soil side along the axis. The separation distance in the direction, when the rock block side master node is relative to the soil side from the node to the axis The normal separation is positive during positive movement. The first tangential slip is the distance from the rock block side principal node relative to the soil side from the node along the axis. The tangential slip component in the direction is obtained by interpolating the relative displacement vectors of the two nodes with the axis. The unit vector dot product is used. The second tangential slip is the distance from the rock block side master node relative to the soil side from the node along the axis. The tangential slip component in the direction is obtained by interpolating the relative displacement vectors of the two nodes with the axis. The resultant is obtained by the dot product of unit vectors. The combined tangential slip of the first and second tangential slips is calculated by taking the square root of the sum of squares.

[0079] Based on the dynamic adjustment of the normal separation and tangential slip, the friction coefficient and cohesion of the contact elements are updated at each convergence time step, resulting in a closed-loop updated numerical simulation result of the soil-rock hybrid shear characteristics. This closed-loop updated numerical simulation result includes the rock block displacement field, the soil strain field, and the damage distribution cloud map of the soil-rock interface. The rock block displacement field consists of the displacement components of all polyhedral rock block nodes in three directions; the soil strain field consists of the equivalent plastic strain values ​​of all tetrahedral mesh elements of the soil matrix; and the damage distribution cloud map of the soil-rock interface consists of the spatial distribution of the tangential cumulative damage variables of all debonding elements on the soil-rock interface.

[0080] After generating closed-loop updated numerical simulation results of the soil-rock hybrid shear properties, post-processing is performed on the simulation results. The equivalent plastic strain values ​​of all tetrahedral mesh elements within the shear zone are extracted. The shear zone is defined as follows: a horizontal strip with a height of one-tenth of the total model height is taken along the shear direction in the center of the three-dimensional geometric model of the soil-rock hybrid, and all tetrahedral mesh elements within this horizontal strip are included in the shear zone. The equivalent plastic strain value is the accumulated plastic strain invariant of the tetrahedral mesh elements at the current time step, directly extracted from the element state variables output by the finite element simulation platform.

[0081] Tetrahedral mesh elements whose equivalent plastic strain exceeds a preset plasticity threshold are marked as failed elements. The preset plasticity threshold is set to 0.05, based on the fact that soil material exhibits significant plastic flow when the equivalent plastic strain reaches 0.05, indicating localized failure. The total volume of all failed elements is calculated, taking the actual volume of the tetrahedral mesh element after deformation at the current time step, obtained by updating the node coordinates using node displacements. The ratio of the total volume of failed elements to the total volume of the shear zone region is calculated to obtain the shear failure volume fraction. The shear failure volume fraction is calculated as follows:

[0082]

[0083] in, It represents the shear failure volume fraction, a dimensionless quantity, with a value ranging from zero to one; Indicates the total number of destroyed units; Indicates the first The volume of a single destruction unit, expressed in cubic millimeters; This represents the total number of all tetrahedral mesh elements within the shear band region; Indicates the first [unit] within the shear band region The volume of each tetrahedral mesh element, in cubic millimeters. The summation operation iterates through all relevant tetrahedral mesh elements within the shear band region.

[0084] A curve was plotted to show the shear failure volume fraction as a function of shear displacement, where the shear displacement was the cumulative shear displacement at the upper boundary at the end of each time step. The x-axis represents shear displacement in millimeters, and the y-axis represents the shear failure volume fraction, a dimensionless quantity. The second derivative of the shear failure volume fraction curve was numerically calculated using a central difference scheme. Local peak points were searched on the second derivative curve before and after the zero-crossing point from positive to negative values. The shear displacement points corresponding to these peak points were identified as the displacement points corresponding to the peak shear strength of the soil-rock mixture. The displacement points corresponding to the peak shear strength were used to read the corresponding peak shear stress values ​​from the shear stress-shear displacement curve.

[0085] See Figure 11 In the figure, the horizontal axis represents shear displacement in millimeters, ranging from 0 to 30 mm; the left vertical axis represents shear failure volume fraction, ranging from 0 to 1, dimensionless, reflecting the proportion of the volumetric unit within the shear band; the right vertical axis represents the second derivative of the volume fraction, in units of 1 / mm, used to analyze the trend of failure rate variation. The solid curve represents the trend of shear failure volume fraction with shear displacement, the dashed curve represents the change of the second derivative of the shear failure volume fraction, and the triangle marker indicates the shear displacement corresponding to the peak shear strength determined by the second derivative curve.

[0086] As shown in the figure, the shear failure volume fraction gradually increases with the increase of shear displacement. In the initial stage (0-10 mm), the failure volume fraction is low and increases slowly, indicating that local plastic failure in the soil-rock mixture is just beginning to form in the early stage of shearing. In the middle and late stages (10-20 mm), the failure volume fraction increases rapidly, and the curve becomes steeper, indicating that a large number of shear failure units are generated, and significant plastic flow and damage accumulation occur in the soil-rock mixture. At a shear displacement of approximately 15 mm, the second derivative curve of the volume fraction changes from positive to negative and shows a clear peak. The corresponding triangle marker indicates the peak shear strength displacement point. At this point, the failure rate of the soil-rock mixture reaches its maximum, and shear failure enters the peak stage. Subsequently, the rate of increase of the failure volume fraction curve slows down and tends to saturate, eventually approaching 1, indicating that most units in the shear zone region have undergone plastic failure, and the soil-rock mixture enters the residual failure stage.

[0087] This figure reflects the dynamic evolution process of the shear failure volume ratio calculated by extracting the equivalent plastic strain value of the tetrahedral mesh element in the shear zone region and determining the volume of the failed element in Example 5. It intuitively shows the whole process of the soil-rock mixture transforming from an elastoplastic state to a failed state under shear loading, as well as the key shear displacement parameters corresponding to the peak shear strength, providing a quantitative basis for the determination and analysis of the numerical simulation results of closed-loop updates.

[0088] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application.

Claims

1. A three-dimensional numerical simulation method for the shear properties of mixed soil and rock based on simulation data, characterized in that, include: A three-dimensional geometric model of a soil-rock mixture was constructed using computer-aided design, which included a randomly generated polyhedral rock block skeleton and the surrounding soil matrix interpenetrating with each other. Based on the three-dimensional geometric model of the soil-rock mixture, the rock block corner area and the soil-rock interface layer are discretized at multiple scales using a mesh adaptive densification algorithm to generate a three-dimensional numerical simulation mesh model with contact elements. The three-dimensional numerical simulation mesh model is imported into the finite element simulation platform, and dynamic shear load is applied according to the preset shear rate and normal stress boundary conditions to simulate the stress transfer and displacement evolution of the soil-rock mixture during the shearing process. During the simulation, the normal separation and tangential slip at the soil-rock interface are extracted in real time, and the friction coefficient and cohesion of the interface are dynamically adjusted based on the normal separation and tangential slip to form a closed-loop updated numerical simulation result of the soil-rock mixed shear characteristics.

2. The three-dimensional numerical simulation method for soil-rock mixed shear properties based on simulation data as described in claim 1, characterized in that, The construction of a three-dimensional geometric model of a soil-rock mixture, using computer-aided design, includes the interweaving of a randomly generated polyhedral rock block skeleton with the surrounding soil matrix. Specifically, this includes: Based on the preset rock block volume fraction and grain size distribution curve, polyhedral rock blocks are randomly generated in three-dimensional space and the vertex coordinates and surface normal vectors of each rock block are recorded. The vertex coordinates and surface normal vectors of each rock block are input into the Poisson disk sampling algorithm to generate uniformly distributed contact seed points on the surface of the rock block. Using the contact seed point as the anchor point, the boundary of the surrounding soil matrix is ​​offset towards the rock block by a preset embedding depth, so that the soil matrix and the surface of the rock block create an overlapping embedding area, generating a three-dimensional geometric model of soil and rock interembedded.

3. The three-dimensional numerical simulation method for soil-rock mixed shear properties based on simulation data as described in claim 2, characterized in that, Based on the three-dimensional geometric model of the soil-rock mixture, a multi-scale discretization process is performed on the rock block corner region and the soil-rock interface layer using an adaptive mesh refinement algorithm, specifically including: Identify the edges and vertices of each polyhedral rock block in the three-dimensional geometric model, calculate the dihedral angle of each edge, and mark the edge as an acute edge when the dihedral angle is less than a preset angle threshold. Centered on the acute-angled edge, a preset layer thickness is extended outward along the normal direction to generate a local encryption domain containing the edge transition area; Within the local encryption domain, the side length of the tetrahedral mesh is reduced to one-third of the side length of the global mesh, and hexahedral interface elements are arranged in layers along the normal direction at the soil-rock interface layer.

4. The three-dimensional numerical simulation method for soil-rock mixed shear properties based on simulation data as described in claim 1, characterized in that, The three-dimensional numerical simulation mesh model is imported into the finite element simulation platform, and dynamic shear loads are applied according to preset shear rate and normal stress boundary conditions, specifically including: The lower boundary nodes of the three-dimensional numerical simulation mesh model are completely fixed, and a constant normal stress is applied to the upper boundary nodes. A shear displacement that grows linearly with time is applied to the upper boundary node along the horizontal direction. The ratio of the shear displacement to the total height of the model in each time step is calculated to obtain the real-time shear strain. When the real-time shear strain reaches the preset peak strain threshold, the shear rate is reduced to one-quarter of the initial shear rate, and shear displacement is continued to be applied until the total shear strain reaches the preset termination strain value.

5. The three-dimensional numerical simulation method for soil-rock mixed shear properties based on simulation data according to claim 4, characterized in that, The stress transfer and displacement evolution of the simulated soil-rock mixture during shearing specifically include: Within each time step, the Gaussian point stress tensor of all tetrahedral mesh elements is calculated, and the equivalent nodal force of each node is obtained by shape function interpolation based on the Gaussian point stress tensor. Substitute the equivalent nodal forces of each node into the global stiffness matrix to solve for the nodal displacement increment, and determine whether the normal tensile stress of the soil-rock interface element exceeds the tensile strength of the soil. If the normal tensile stress exceeds the tensile strength of the soil, the contact interface unit is marked as a debonding unit, and the tangential stiffness of the debonding unit is reduced to zero.

6. The three-dimensional numerical simulation method for soil-rock mixed shear characteristics based on simulation data according to claim 5, characterized in that, The real-time extraction of the normal separation and tangential slip at the soil-rock interface specifically includes: Traverse all soil-rock contact interface elements marked as debonding elements, and obtain the coordinates of the master node on the rock block side and the slave node on the soil side of the interface element. Calculate the magnitude of the coordinate difference between the master node on the rock block side and the slave node on the soil side at the current time step, and subtract the magnitude of the coordinate difference at the initial time step to obtain the normal separation amount; Calculate the change in coordinate difference between the master node on the rock block side and the slave node on the soil side along the shear direction at the current time step, subtract the displacement contributed by elastic deformation, and obtain the tangential slip.

7. The three-dimensional numerical simulation method for soil-rock mixed shear properties based on simulation data according to claim 6, characterized in that, When calculating the tangential slip, the displacement contributed by the elastic deformation is calculated based on the tangential stress and initial tangential stiffness at the current time step.

8. The three-dimensional numerical simulation method for soil-rock mixed shear properties based on simulation data according to claim 6, characterized in that, The friction coefficient and cohesion of the contact interface are dynamically adjusted based on the normal separation amount and the tangential slip amount, specifically including: Substitute the normal separation amount into the exponential decay function to calculate the normal separation decay factor, and multiply the initial friction coefficient by the normal separation decay factor to obtain the friction coefficient at the current time step. Divide the tangential slip by the preset characteristic slip length to obtain the tangential cumulative damage variable, and multiply the initial cohesion by the difference between the tangential cumulative damage variable and the initial cohesion to obtain the cohesion at the current time step. When the tangential cumulative damage variable is greater than the preset failure threshold, the cohesion at the current time step is forcibly set to zero.

9. The three-dimensional numerical simulation method for soil-rock mixed shear properties based on simulation data according to claim 1, characterized in that, The three-dimensional numerical simulation mesh model with contact elements establishes the local coordinate system of the contact elements in the following manner: For each contact unit on the soil-rock interface, the outward normal direction of the side triangle of the rock block in the contact unit is taken as the first axis direction of the local contact coordinate system. The projection direction of the horizontal shear loading direction onto the plane containing the triangular facet of the rock block is taken as the second axis direction of the contact local coordinate system; The cross product direction of the first axis direction and the second axis direction is taken as the third axis direction of the contact local coordinate system, and the normal separation amount, the first tangential slip amount and the second tangential slip amount are projected onto the first axis, the second axis and the third axis direction for calculation.

10. The three-dimensional numerical simulation method for soil-rock mixed shear properties based on simulation data according to claim 1, characterized in that, The closed-loop updated numerical simulation results of soil-rock hybrid shear properties include rock block displacement field, soil strain field, and soil-rock contact interface damage distribution cloud map.

Images