Simulation method and system based on tailings-based material parameters, and storage medium

Abstract

This invention relates to the field of simulation, and in particular to a simulation method, system, and storage medium based on tailings-based material parameters. Specifically, it includes: acquiring mechanical test data of tailings-based materials under different confining pressures, water content states, and loading paths; calculating dilatation-related indices; constructing a dilatation response sequence; identifying key feature points and extracting slope and curvature features; establishing the correlation between state parameters, critical state indices, and plastic flow direction based on these features; constructing a multi-factor-related dilatation function; calculating the path memory factor; constructing historical variables to correct critical state indices; determining the structural attenuation coefficient during path switching; correcting the yield surface, hardening, and flow parameters; determining particle interlocking representations based on normalized cumulative dilatation strain; updating shear and bulk modulus based on the attenuation coefficient; updating parameters according to the plastic strain increment correction rule selected based on the calculated stress location; completing the numerical simulation; and outputting the stress field, strain field, and volumetric response results under a given boundary.

Description

Technical Field

[0001] This application belongs to the field of simulation, and in particular relates to a simulation method, system and storage medium based on tailings-based material parameters. Background Technology

[0002] Tailings substrate materials, as special granular materials, exhibit volumetric strain due to interparticle friction, tumbling, and rearrangement under varying confining pressures, water content, and loading paths. Geotechnical constitutive models cannot accurately represent the shear dilatation response of tailings substrate materials induced by microstructural changes under multi-field correlation, leading to significant discrepancies between predicted stress, strain, and volumetric strain responses in engineering simulations and actual conditions. This paper combines critical state soil mechanics with shear dilatation response representation techniques, utilizing state parameters, plastic strain increment ratios, and the ratio of deviatoric stress to mean stress to improve simulation parameters for tailings substrate materials. By detecting the volumetric strain behavior during material shearing, the predictive power of an optimized model based on non-correlated flow rules is constructed. However, due to the lack of identification of the slope and curvature characteristics of the complete dilatation response sequence, it is impossible to construct a dilatation function that is associated with confining pressure, water content, and loading path. When switching loading paths, existing methods cannot calculate path memory factors and water content history variables, and lack a dimensionless structural attenuation representation based on the plastic work density variation law, resulting in distortion of the correction of yield surface scale, hardening parameters, and flow parameters. Moreover, when updating shear and bulk modulus, the particle interlocking effect caused by cumulative dilatation volumetric strain is omitted, and no plastic strain increment correction rule is formulated according to the position of the trial stress relative to a specific dilatation stage in each stress calculation step, causing the parameter changes of the existing model to deviate from the actual physical process, and it is impossible to output simulation response results under given boundary conditions. Summary of the Invention

[0003] The first aspect of this invention provides a simulation method based on tailings-based material parameters, comprising the following steps: The test data of axial strain, volumetric strain, mean stress and deviatoric stress of tailings-based materials under different confining pressures, water content and loading paths were obtained. The ratio of deviatoric stress to mean stress and the ratio of plastic volumetric strain increment to plastic shear strain increment were calculated to construct a shear dilatation response sequence. Based on the shear dilatation response sequence, the shear dilatation inflection point, zero shear dilatation point and the post-peak fallback segment were identified, and the corresponding slope, curvature and segment length features were extracted. Based on the aforementioned characteristic quantities, the correlation between state parameters, critical state indices, and plastic flow direction is established, and a dilatation function that varies with confining pressure, water-bearing state, and loading path is constructed. The path memory factor is calculated based on the changes in the zero dilatation point under different confining pressures and loading paths, and historical variables are constructed in conjunction with the water-bearing state to correct the critical state indices. When the loading path is switched, the structural attenuation coefficient is determined based on the change in the ratio of plastic work density to reference energy density of adjacent path segments, and the yield surface scale parameters, hardening parameters, and flow parameters are corrected. Based on the normalized cumulative shear dilatation strain of the current calculation step and the cumulative volumetric strain characteristics that distinguish between shrinkage and swelling states, the particle interlocking representation is determined. The shear modulus and bulk modulus are updated together with the structural attenuation coefficient and the compaction hardening factor. Within each stress step, the corresponding plastic strain increment correction rule is selected according to the position of the trial stress relative to the shear dilatation inflection point and the zero shear dilatation point, and the strength parameters, deformation parameters, and flow parameters are updated. Based on the updated parameters, numerical simulation of tailings-based materials is performed to obtain the stress field, strain field, and volumetric strain response results under given boundary conditions.

[0004] A second aspect of this invention provides a simulation system based on tailings-based material parameters, comprising the following modules: The extraction module is used to acquire axial strain, volumetric strain, mean stress, and deviatoric stress test data of tailings-based materials under different confining pressures, water-bearing states, and loading paths. It calculates the ratio of deviatoric stress to mean stress and the ratio of plastic volumetric strain increment to plastic shear strain increment, and constructs a shear dilatation response sequence. Based on the shear dilatation response sequence, it identifies the shear dilatation inflection point, zero shear dilatation point, and post-peak fallback segment, and extracts the corresponding slope, curvature, and segment length features. The correction module is used to establish the correlation between state parameters, critical state indices and plastic flow direction based on the characteristic quantities, and to construct a dilatation function that varies with confining pressure, water-bearing state and loading path; to calculate the path memory factor based on the changes in the zero dilatation point under different confining pressures and loading paths, and to construct historical variables in combination with the water-bearing state for correcting the critical state indices; when the loading path is switched, to determine the structural attenuation coefficient based on the change in the ratio of plastic work density to reference energy density of adjacent path segments, and to correct the yield surface scale parameters, hardening parameters and flow parameters. The update module is used to determine the particle interlocking representation based on the normalized cumulative shear dilatation strain of the current calculation step and the cumulative volumetric strain characteristics that distinguish between shrinkage and swelling states. It also updates the shear modulus and bulk modulus in conjunction with the structural attenuation coefficient and the compaction hardening factor. Within each stress step, it selects the corresponding plastic strain increment correction rule based on the position of the trial stress relative to the shear dilatation inflection point and the zero shear dilatation point, and updates the strength parameters, deformation parameters, and flow parameters. Based on the updated parameters, it performs numerical simulation of tailings-based materials to obtain the stress field, strain field, and volumetric strain response results under given boundary conditions.

[0005] This invention constructs a shear dilatation response sequence by acquiring multi-condition test data and extracting stage-specific characteristic quantities to represent the shear dilatation variation law of tailings-based materials. Combining the influence of multiple factors such as confining pressure, water content, and loading path, the critical state indices are corrected by calculating the path memory factor and historical variables. During the loading path switching stage, the degradation process of the material's internal structure is represented by the structural attenuation coefficient, thereby correcting the yield surface, hardening, and key flow parameters. Simultaneously, by combining particle interlocking representation quantities and the correction rule for plastic strain increment within the stress step, adjustments are made to the shear modulus, bulk modulus, and strength parameters. This addresses the shortcomings of conventional models in representing mechanical properties and improves the reliability of numerical simulations of stress field, strain field, and bulk response under specific boundary conditions. Attached Figure Description

[0006] Figure 1 A flowchart of the first embodiment; Figure 2 This is a schematic diagram for identifying dilatation features. Figure 3 This is a schematic diagram illustrating the modulus decay. Figure 4 This is a schematic diagram comparing the effects of changes in deviatoric stress. Detailed Implementation

[0007] Exemplary embodiments will now be described in detail, examples of which are illustrated in the accompanying drawings. When the following description relates to the drawings, unless otherwise indicated, the same numbers in different drawings represent the same or similar elements. The embodiments described in the following exemplary embodiments do not represent all embodiments consistent with this specification. Rather, they are merely examples of apparatuses and methods consistent with some aspects of this specification as detailed in the appended claims.

[0008] It should be understood that the terms “comprising” and “having”, and any variations thereof, in the embodiments of this specification are intended to cover but not exclude inclusion. For example, a product or device that includes a series of components is not necessarily limited to those components that are explicitly listed, but may include other components that are not explicitly listed or that are inherent to such product or device.

[0009] Example 1 In Embodiment 1 of the present invention, as Figure 1 As shown, a simulation method based on tailings-based material parameters includes: S1: Obtain test data of axial strain, volumetric strain, mean stress and deviatoric stress of tailings base material under different confining pressures, water content, and loading paths; calculate the ratio of deviatoric stress to mean stress and the ratio of plastic volumetric strain increment to plastic shear strain increment; and construct a shear dilatation response sequence.

[0010] Mechanical tests were conducted on tailings-based materials under different confining pressures, initial moisture contents, and loading paths including conventional triaxial compression, proportional loading, and cyclic loading, using a dynamic and static triaxial testing system. Real-time time-series data of axial strain, volumetric strain, mean stress, and deviatoric stress were acquired at different time steps using sensors. The NumPy library was used in Python to read the experimental data files, and the raw time-series data was input into the signal processing module of the SciPy library, where a low-pass filtering algorithm was used to remove high-frequency environmental noise interference. Hooke's law, combined with the material's elastic modulus and Poisson's ratio, was used to isolate the elastic strain component, obtaining the increments of plastic volumetric strain and plastic shear strain. Based on array calculations, the stress ratio was calculated by dividing the deviatoric stress at each data point by the mean stress, and the shear dilatation ratio was calculated by dividing the increment of plastic volumetric strain by the increment of plastic shear strain. The stress ratios and shear dilatation ratios at different loading steps were concatenated and combined to generate a multidimensional time-series matrix, forming a complete shear dilatation response sequence.

[0011] In an optional embodiment, the acquisition of axial strain, volumetric strain, mean stress, and deviatoric stress test data of tailings-based materials under different confining pressures, water-bearing states, and loading paths, and the calculation of the ratio of deviatoric stress to mean stress, the ratio of plastic volumetric strain increment to plastic shear strain increment, and the construction of a shear dilatation response sequence, includes: Within each loading step, the triaxial principal stress data and corresponding principal strain data of the tailings-based material are read; The mean stress and deviatoric stress are calculated based on the triaxial principal stress data, and the stress ratio is obtained by dividing the deviatoric stress by the mean stress. The plastic volumetric strain increment and the plastic shear strain increment are calculated based on the principal strain data, and the shear dilatation ratio is obtained by dividing the plastic volumetric strain increment by the plastic shear strain increment. The stress ratio and shear dilatation ratio corresponding to each loading step are arranged in the loading order to construct a shear dilatation response sequence that represents the relationship between stress state and volumetric deformation.

[0012] In acquiring triaxial principal stress data and corresponding principal strain data, taking a conventional consolidated drained triaxial compression test (CD) as an example, the confining pressure range is typically set to 100 kPa to 800 kPa, for example, three test levels: 100 kPa, 200 kPa, and 400 kPa. The sampling frequency is set to 1 Hz, and the major principal stresses recorded by the data acquisition device are read in real time. and minor principal stress In conventional triaxial intermediate principal stress , and axial strain measured by extensometer and radial strain Within each tiny loading step, the mean principal stress of the current state is calculated. and generalized deviatoric stress Calculate the stress ratio For example, when a certain incremental step reads... , When, it can be calculated , Thus, the stress ratio can be obtained. .

[0013] The elastic and plastic components in the total strain are separated based on the generalized Hooke's law. The initial bulk modulus K is determined based on the elastic unloading test results of the tailings-based material, typically in the range of 20 MPa to 50 MPa, with a preferred example of 30 MPa, along with the initial shear modulus. Typically, the pressure is in the range of 10 MPa to 30 MPa, with a preferred example of 15 MPa, thereby calculating the elastic volumetric strain increment of this loading step. and elastic shear strain increment Based on the total strain increment, i.e. , Subtracting the elastic component obtained above, we obtain the plastic volumetric strain increment representing irreversible deformation. With plastic shear strain increment Divide the two values ​​to calculate the shear dilatation ratio during the instantaneous loading step. ,when A value greater than 0 indicates that the material is undergoing volume shrinkage, while a value less than 0 indicates that it is undergoing volume expansion. The extracted discrete state parameter coordinates are calculated over N consecutive loading increments, for example, N=2000. By storing the data in an array according to the absolute time sequence of the laboratory loading, a high-resolution one-dimensional time series of shear dilatation response containing the complete stress ratio change history and the corresponding volumetric rate change was successfully constructed.

[0014] S2, based on the shear dilatation response sequence, identify the shear dilatation inflection point, zero shear dilatation point, and post-peak fallback segment, and extract the corresponding slope, curvature, and segment length features.

[0015] The shear dilatation response sequence is input into the Pandas library to construct a data frame structure. The peak finding algorithm from the SciPy library is used to perform first-order differentiation on the shear dilatation ratio sequence and detect the locations where the derivative sign changes. The location where the derivative changes from positive to negative is identified as the shear dilatation inflection point, and the point where the shear dilatation ratio is equal to zero is extracted as the zero shear dilatation point. The post-peak decline segment is obtained from the data segment after the stress ratio reaches its peak in the shear dilatation response sequence. The polynomial feature algorithm from the Scikit-learn library is used to perform quadratic polynomial regression fitting on the neighborhood of the shear dilatation inflection point, the neighborhood of the zero shear dilatation point, and the post-peak decline segment. The first-order derivative of each fitted term is calculated as the slope at the corresponding feature location, and the curvature is obtained by calculating the related terms of the second-order derivative. Simultaneously, the Euclidean distance between the first and last data points of the corresponding segment in the state space is calculated using the linear algebra module of the NumPy library, and this Euclidean distance is used as the segment length feature.

[0016] In an optional embodiment, the step of identifying the dilatation inflection point, zero dilatation point, and post-peak pullback segment based on the dilatation response sequence, and extracting the corresponding slope, curvature, and segment length features, includes: The derivative of the curve of shear dilatation ratio versus shear strain in the shear dilatation response sequence is taken, and the point where the derivative turns from negative to positive and the shear dilatation ratio reaches a minimum value is identified as the shear dilatation inflection point; Traverse the shear dilatation response sequence and identify the stress state point corresponding to the shear dilatation ratio crossing zero as the zero shear dilatation point; Identify the stress-strain softening interval after the deviatoric stress reaches its peak value, and define the interval as the post-peak decline segment; The tangent slopes before and after the dilatation inflection point, the curvature near the zero dilatation point, and the shear strain span corresponding to the post-peak fall-off segment are extracted as segment length features.

[0017] The characteristics of the variation of shear dilatation ratio with shear strain in tailings-based materials and the identification results of key feature points are as follows: Figure 2 As shown. When processing the constructed dilatation response sequence, the central difference numerical algorithm is used to calculate the dilatation ratio D sequence relative to the plastic shear strain. discrete first derivative sequence After smoothing the derivative curve by applying a moving average filter function with a window size of 5 points, the system searches point by point within the curve for a result satisfying the first derivative. The sign of the number changes from negative across zero to positive, and the corresponding second derivative is calculated. Data feature nodes greater than zero. Points satisfying the second derivative test criterion are identified as shear dilatation inflection points. This physical state typically corresponds to the initial extreme state of the maximum volume expansion rate of tailings material; for example, the minimum shear dilatation ratio D at this point may reach -0.25. The program scans the entire shear dilatation ratio D sequence array sequentially from the starting point, locating the value of D as it first decreases from the positive range, representing the macroscopic contraction stage, and then drops and crosses the absolute zero line into the negative range, representing the critical physical moment of the macroscopic expansion stage. The program then records the corresponding instantaneous stress state parameters at this moment, such as the stress ratio read at this time. It is usually defined as the phase transformation stress ratio and marked as the critical zero shear dilatation point.

[0018] Based on this identification, the coordinates of the global maximum point of the historical data sequence of deviatoric stress q are retrieved, for example, the peak deviatoric stress. The corresponding shear strain The algorithm then determines the post-peak decline segment based on whether the deviatoric stress value shows a continuous monotonic decrease or fluctuating decay after crossing the peak point until it approaches a constant residual stress. For example, a stress-strain interval is defined as one where the cumulative stress decrease exceeds 5% of the peak value. The local feature extraction algorithm module is then activated: short window subsets of shear strain spanning 0.5% are extracted before and after the locked shear dilatation inflection point. Least squares regression is used for local linear regression fitting to extract the tangent slopes before and after the inflection point. For example, the slope of the steep drop segment before the turning point ,and For example, the slope of the gradually rising segment after the turning point Extract the coordinates of 10 consecutive discrete data points before and after the zero dilatation point, and use cubic polynomial interpolation to calculate the local geometric radius of curvature at the zero point of the extreme value. For example, curvature calculation results The shear strain at the residual steady-state initiation point, for example... Subtracting the peak shear strain from the 8% value yields the absolute difference. The shear strain span length is used as a characteristic quantity to represent the energy dissipation characteristics of material ductility softening. All parameters with explicit physical meaning are passed to the equations as floating-point variables, serving as strong constraint boundary conditions to modify the shape of the constitutive model's response.

[0019] S3. Based on the aforementioned characteristic quantities, establish the correlation between state parameters, critical state indices, and plastic flow direction, and construct a shear dilatation function that varies with confining pressure, water content, and loading path.

[0020] The difference between the current porosity and the critical porosity under the corresponding average stress is defined as the state parameter. The extracted slope, curvature, and segment length features, along with the corresponding tensors of confining pressure, water cut, and loading path, are used as input feature vectors. A multilayer perceptron neural network is constructed using the PyTorch framework. Experimental observation shear dilatation ratios under different conditions are used as supervisory labels to train the multilayer perceptron neural network using backpropagation gradient descent. By iteratively updating the weight matrix, the network output layer establishes a mapping relationship between the state parameters, the critical state index representing the relationship between the critical porosity and stress, and the plastic flow direction representing the direction of plastic strain increment. Based on the trained multilayer perceptron network structure and weights, the network black-box model is transformed into an explicit multivariable function equation using a Taylor series expansion algorithm, thereby outputting a shear dilatation function that changes in real time with the current confining pressure, current water cut, and current loading path.

[0021] S4 calculates the path memory factor based on the changes in the zero shear dilatation point under different confining pressures and loading paths, and constructs historical variables in conjunction with the water-bearing state to correct the critical state index.

[0022] The stress ratio and cumulative plastic strain corresponding to the zero-shear dilatation point under different confining pressures and loading paths are read. The deviation between the stress ratio at the zero-shear dilatation point under the current loading path and that under the initial reference path is calculated. The path memory factor is obtained by performing an exponentially weighted integral on the deviation during path turning using the integral function of the SciPy library. The path memory factor is multiplied by the saturation parameter corresponding to the water-bearing state, and the plastic deformation history scalar represented by the cumulative plastic work is added. After processing with a normalization algorithm, a historical variable containing the correlation between the current state and the historical path is generated. Based on the original mathematical expression of the critical state index, the logarithmic term decay function of the historical variable with the natural logarithm as the product correction term is used to perform real-time multiplicative numerical correction on the critical state index, thereby obtaining the corrected critical state index.

[0023] S5, when the loading path is switched, the structural attenuation coefficient is determined based on the change in the ratio of plastic work density to reference energy density of adjacent path segments, and the yield surface scale parameters, hardening parameters and flow parameters are corrected.

[0024] In a subroutine of the finite element method (FEM) program, the strain increment direction tensor is monitored in real time. When the deflection angle of the strain increment tensor in phase space exceeds a set threshold, a loading path switch is determined to have occurred. A numerical integration algorithm is called to calculate the inner product of the deviatoric stress and plastic shear strain increment in two adjacent path segments before and after the path switch, and the entire path is integrated to obtain the plastic work density of the adjacent path segments. The calculated plastic work density is divided by the reference energy density constant corresponding to material failure, and the difference between the two ratios is calculated using an absolute value function as the energy dissipation increment. A negative exponential distribution algorithm is used to convert this energy dissipation increment into a structural attenuation coefficient between 0 and 1. Linear interpolation function equations are established for the initial yield surface scale parameters, initial hardening parameters, and initial flow parameters, respectively, with respect to the structural attenuation coefficient. Multiplication operations are used to perform real-time reduction calculations on the initial values ​​of the relevant parameters, completing the attenuation correction for the yield surface scale parameters, hardening parameters, and flow parameters.

[0025] S6. Based on the normalized cumulative shear dilatation strain of the current calculation step and the cumulative volumetric strain characteristics that distinguish between volumetric contraction and volumetric expansion states, determine the particle interlocking representation quantity, and update the shear modulus and bulk modulus together with the structural attenuation coefficient and the compaction hardening factor.

[0026] The plastic volumetric strain tensor and plastic shear strain tensor of the current computation time step in the finite element solver are extracted. The cumulative shear dilatation strain and the cumulative volumetric strain value distinguishing the deformation direction from the initial time step to the current step are calculated using the Euler backward integration algorithm. The cumulative volumetric contraction strain and cumulative shear dilatation strain are extracted and divided by the corresponding critical failure thresholds to achieve data normalization. An interlocking contact model between particles at the material's mesoscale is constructed using the normalized cumulative volumetric contraction strain to obtain a particle interlocking hardening representation for the degree of particle contact tightness in tailings. Simultaneously, the normalized cumulative shear dilatation strain is used as an independent parameter to characterize the structural loosening damage caused by particle rearrangement. The initial shear modulus and initial bulk modulus of the material in its undamaged state are extracted. The stiffness strengthening factor is constructed by combining the particle interlocking hardening expression with an exponential function. The initial shear modulus and initial bulk modulus are multiplied by the stiffness strengthening factor, and then multiplied by a weakening coefficient formula that includes the structural attenuation coefficient and the normalized cumulative shear dilatation strain. The updated shear modulus and bulk modulus at the current time step are calculated by comprehensively considering the competition mechanism between compaction hardening and shear dilatation softening.

[0027] In an optional embodiment, the step of determining the particle interlocking representation based on the normalized cumulative shear dilatation strain of the current calculation step and the cumulative volumetric strain characteristic that distinguishes between shrinkage and swelling states, and updating the shear modulus and bulk modulus in conjunction with the structural attenuation coefficient and the compaction hardening factor, includes: Calculate the normalized cumulative shear dilatation strain of the current calculation step, and extract the attenuation function value constructed by accumulating the plastic volumetric strain increment when the current calculation step is in the dilatation state as the cumulative volumetric strain feature. At the same time, construct the compaction hardening factor based on the volumetric strain increment when the current calculation step is in the contraction state. The normalized cumulative shear dilatation strain and the cumulative volumetric strain feature are weighted and summed to obtain a particle interlocking representation that represents the degree of damage to the internal particle skeleton reconstruction. The initial shear modulus and initial bulk modulus are multiplied by a non-negative attenuation coefficient constructed from the particle interlocking representation and the structural attenuation coefficient, and further multiplied by the compaction hardening factor to obtain the updated shear modulus and bulk modulus under the current stress step.

[0028] In each stress iteration step of the finite element analysis or self-written numerical calculation program, all plastic volumetric strain increment data from the initial loading state to the current calculation step are extracted. The absolute values ​​of the strain increments in the dilatation stage (i.e., the expansion state with negative increment values) are integrated and accumulated, and then divided by the reference maximum dilatation strain value of the ultimate failure state measured by material testing. For example, the reference value is set to 0.05. This yields the normalized cumulative dilatation strain truncated in the interval [0,1]. If the current calculation step yields... The sign condition of the total volumetric strain increment at each discrete increment step throughout the entire mechanical loading history is determined, and the integral of the volumetric strain increment during the shearing stage (compaction state) is extracted to construct a compaction hardening factor reflecting the densification of the skeleton. On the other hand, only the absolute values ​​of the plastic volumetric strain increments during the expansion stage (expansion state) are extracted and purely numerically accumulated, i.e. The cumulative amount of body expansion is then substituted into a preset monotonically increasing negative exponential decay function, with the preferred function being... The fitting parameter c used to control the convergence of the material damage evolution rate can preferably be 15.0. After calculation by this function mapping, a floating-point value between 0 and 1 is output as a characteristic variable of volumetric swelling damage. For example, by substituting the current cumulative volume of body expansion to obtain... .

[0029] To represent the tumbling and rearrangement behavior and the degree of interlocking damage of the micro-particle skeleton of tailings-based materials during shearing, a set of normalized weighting coefficients is used to perform a linear weighted combination operation on the two damage-related parameters mentioned above. For example, key weighting coefficients for normalized cumulative shear dilatation strain are set through indoor test calibration. The auxiliary weighting coefficient for cumulative body variation characteristics is 0.6. A value of 0.4 ensures compliance. Using the conservation conditions, calculate the scalar result of the weighted summation. =0.322 is used as the particle intercalation representation. The dimensionless structural attenuation coefficient stored in the previous analysis step is retrieved. For example, if the current value is 0.80, construct a nonnegative decay multiplier formula for controlling the degradation of the elastic modulus. The stiffness degradation rate control constant 'a' set in the model parameters can be taken as an empirical value of 1.2, so the current reduction factor R can be calculated to be approximately 0.69.

[0030] The user-defined initial macroscopic shear modulus G (e.g., 20 MPa) and initial macroscopic bulk modulus K (e.g., 40 MPa) are multiplied by the generated non-negative attenuation coefficient R, and then multiplied by the compaction hardening factor reflecting the volumetric shrinkage effect. For example, if a significant compaction and consolidation process occurred during the initial loading phase, the current... The updated and effective true shear modulus G = 15.87 MPa and bulk modulus K = 31.74 MPa were obtained under the current stress step. The evolution of bulk modulus and shear modulus of the tailings-based material during the entire loading process, alternating between axial strain and material bulk deformation state, is as follows: Figure 3 As shown.

[0031] Those skilled in the art should know that the current time step corresponds to either a compaction state or an expansion state, and the corresponding process can be used for calculation based on the current time step. In another embodiment, the normalized cumulative shear dilatation strain of the current calculation step is calculated, and a determination is made based on the volumetric strain state of the current calculation step: if the current calculation step is in a volumetric expansion state, the current plastic volumetric strain increment is used to update the attenuation function value as the cumulative volumetric strain characteristic, and the compaction hardening factor is kept at the value of the previous calculation step; if the current calculation step is in a volumetric contraction state, the current volumetric strain increment is used to update the compaction hardening factor, and the cumulative volumetric strain characteristic is kept at the value of the previous calculation step; the normalized cumulative shear dilatation strain and the cumulative volumetric strain characteristic are weighted and summed to obtain a particle interlocking representation quantity representing the degree of internal particle skeleton reconstruction damage; the initial shear modulus and the initial volumetric modulus are multiplied by a non-negative attenuation coefficient constructed from the particle interlocking representation quantity and the structural attenuation coefficient, and further multiplied by the compaction hardening factor to obtain the updated shear modulus and volumetric modulus under the current stress step.

[0032] S7. Within each stress step, select the corresponding plastic strain increment correction rule based on the position of the trial stress relative dilatation inflection point and the zero dilatation point, and update the strength parameters, deformation parameters, and flow parameters.

[0033] In the material subroutine of the Abaqus finite element software, the elastic prediction module is called to calculate the elastic test stress state point and extract the shear dilatation inflection point stress threshold and the zero shear dilatation point stress threshold of the current state. A logical judgment algorithm is used to compare the elastic test stress state point with each threshold. When the elastic test stress is less than the shear dilatation inflection point threshold, the first set of strain correction rules including the compressive hardening function is selected for backward Euler mapping, iteratively solving for the plastic strain increment. When the elastic test stress is between the shear dilatation inflection point and the zero shear dilatation point, the tangent modulus correction algorithm is called to match the plastic strain increment calculation before the phase transition. When the elastic test stress crosses the zero shear dilatation point, the associated flow rule is used to match the plastic strain increment calculation corresponding to the fully plastic stage. The actual plastic strain increment after convergence of the iterative calculation and the current cumulative plastic work state quantity are substituted into the Mohr-Coulomb failure criterion and the elastic constitutive model, and Jacobian matrix updates and tensor state variable updates are performed on the strength parameters, deformation parameters, and flow parameters.

[0034] In an optional embodiment, the step of selecting the corresponding plastic strain increment correction rule based on the position of the trial stress relative dilatation inflection point and the zero dilatation point within each stress step, and updating the strength parameters, deformation parameters, and flow parameters, includes: Compare the calculated stress ratio of the current stress step with the critical stress ratio corresponding to the dilatation inflection point and the zero dilatation point; If the calculated stress ratio is less than the critical stress ratio corresponding to the dilatation inflection point, it is determined that the body is in the shrinkage stage, and the initial plastic strain increment is calculated using the associated flow rule. If the calculated stress ratio is between the critical stress ratio corresponding to the dilatation inflection point and the zero dilatation point, it is determined that the body is in the transition stage from contraction to expansion, and the anisotropic tensor is used to perform non-associated flow compensation on the initial plastic strain increment. If the calculated stress ratio is not less than the critical stress ratio corresponding to the zero dilatation point, it is determined that it is in the strong dilatation stage. The plastic flow direction is updated according to the corrected critical state index and the dilatation function, and the plastic strain increment is recalculated.

[0035] In the stress integral return mapping algorithm of the underlying code of elastoplastic finite element analysis, after receiving a new strain increment, the constitutive model calculates the predicted trial deviatoric stress under the assumption that the current step is completely elastic loading. and trial mean stress Then, by dividing, the calculated stress ratio is obtained. At this point, the lower limit critical stress ratio of the shear dilatation inflection point, corrected based on a specific porosity ratio and pre-stored in the preceding module, is retrieved at the current Gaussian integration point. For example, set it to an empirical value of 0.85, and the ratio of the upper limit critical stress at the zero shear dilatation point. For example, let's set it to 1.20 and establish a strict three-branch conditional statement system. When the calculated value is... Strictly less than For example, when the calculated value is 0.65, it is determined that the material domain is still in the dense pure volume shrinkage stage. At this time, the plastic potential function g representing the strain direction is forcibly set to be coplanar and equivalent with the yield function f representing the yield boundary. That is, the associated flow rule g=f is activated. The orthogonal mapping matrix is ​​generated by using the normal partial derivative vector of the yield surface in the three-dimensional stress space. The plastic volume and shear strain increment in the initial state are calculated to ensure that the local tangent stiffness matrix in the early stage of loading is strictly symmetric and convergence is achieved.

[0036] If the program determines that the calculated stress ratio falls within the transition range, that is... For example, the result calculated at this time With a value of 1.05, the material particles have entered a critical transition state from volumetric contraction to volumetric expansion, characterized by frictional tumbling, mutual rearrangement, and crossover. To reasonably represent the asymmetric deformation characteristics induced by the microparticle arrangement, the termination correlation flow criterion is replaced by a second-order anisotropic tensor representing the directional characteristics of the micropore distribution. For example, by assigning partial differential values ​​related to the tilt angle of the anisotropic principal axes and the angle between them and the current major principal stress axes, and fusing them with the isotropic plastic potential function, a non-associated potential function surface g≠f that produces tilt deviation can be constructed. Based on this new non-associated partial derivative... An additional anisotropic shear compensation component is added to the original initial plastic strain increment. The flow direction angle is corrected by smooth torsional interpolation. And when detected... When the trigger condition is greater than or equal to 1.20, the material element has crossed the phase transition surface and fully erupted into structural expansion failure, entering a strong dilatation stage; the previous transition period tensor compensation is discarded, and the improved dilatation function formula is fully activated. ,in To correct the current critical stress ratio by incorporating effective historical variables such as real-time porosity, for example, a value of 1.30 is obtained, with the constant A set at 0.6. This equation defines the ratio of the normal gradients of the mean stress and deviatoric stress planes, i.e., determined by the relationship... The direction of the flow vector is forcibly defined as a non-associated mapping path. By assembling a plastically consistent tangent stiffness equation that includes this updated vector, the return point of the current trial stress on the new plastic yield surface is obtained through iterative solution. The corrected and physically compatible real plastic strain increment is output, thereby detecting the local dilatation zone phenomenon in deep surrounding rock.

[0037] S8 performs numerical simulation of tailings-based materials based on the updated parameters, and obtains the stress field, strain field and volumetric response results under given boundary conditions.

[0038] The updated strength, deformation, flow, and state variables are assembled into a material-consistent tangent stiffness matrix and passed to the global finite element solver for global stiffness matrix assembly. Externally given displacement and load boundary conditions are applied to the nodes of the global equations. The Newton-Raphson iterative algorithm is used to solve the equations globally. When the residual force norm in the current calculation increment is less than the tolerance threshold, global convergence is considered achieved, and the stress and strain data of all Gaussian integral points are updated. The iteration continues to advance according to the preset time increment until the set total analysis time step is reached. The finite element post-processing module reads the calculation result file to generate principal stress distribution cloud maps on the nodes and elements to obtain the stress field, generates equivalent plastic strain contour maps to obtain the strain field, and extracts the porosity distribution tensor of the whole model to obtain the volumetric response results, thus realizing the three-dimensional full-process numerical simulation of tailings-based materials.

[0039] In an optional embodiment, the numerical simulation of tailings-based materials based on the updated parameters to obtain the stress field, strain field, and bulk strain response results under given boundary conditions includes: A three-dimensional mesh geometric model of the tailings-based material is established, and initial void ratio, initial stress field and boundary conditions are assigned. Input the initial shear modulus, initial bulk modulus, strength parameters, flow parameters, and variation rules as constitutive model parameters into the finite element calculation program; Within each computational increment step, the constitutive model parameters are called, the shear modulus and bulk modulus are updated according to the current increment step state, the stress increment at the Gaussian integral point is calculated, and the overall stiffness matrix is ​​updated. Solve the global equilibrium equations and output the displacement, stress, strain and volumetric strain response results under each loading stage to complete the numerical simulation of tailings-based materials.

[0040] In the actual execution phase of the numerical simulation analysis software, engineers use a preprocessor to establish the physical geometry and perform mesh discretization based on the in-situ dimensions of real mine tailings ponds or underground backfill bodies. For example, a three-dimensional mesh geometric model containing 100,000 eight-node linear reduced integral hexahedral solid elements is generated. In the attribute assignment stage of the initial analysis step, the global integration points of the model are uniformly assigned or the initial porosity obtained from the test is assigned according to the depth gradient formula. The variable's value range is typically preferred to be between 0.75 and 0.90, such as 0.85 for the surface layer and 0.75 for the bottom layer according to consolidation principles. A vertically downward gravitational acceleration constant and the corresponding geostress coefficient are applied to balance and generate an initial self-weight stress field that conforms to in-situ statics. For example, the calculated initial average principal stress of the bottom constraint surface is approximately 200 kPa. Environmental force boundary conditions are set on the node set, including applying environmental boundary conditions to the bottom nodes. Completely fixed constraints are applied, and roller support constraints with zero normal displacement are applied to the nodes on the four outer surfaces. After preparation, a user-defined material subroutine, such as the UMAT subroutine module written in Fortran, is mounted to load the previously calibrated parameters including initial bulk modulus K=30MPa, initial shear modulus G=15MPa, and critical friction angle. Conventional strength parameters and various flow parameters in the control matrix of the variation equation are imported as static variable arrays into the core library of the finite element solution engine.

[0041] Once the static implicit or explicit dynamics calculation process is officially triggered, the total loading time is divided into thousands of time increments. Within any calculation increment step, the above subroutine frequently calls the Gaussian integration point of each element within the geometric model. During the call, the principal strain increment matrix passed down from the global equations to the current integration point is extracted. The algorithm extracts historical parameters stored in memory variables after convergence at the previous time step. It then initiates the constitutive inner loop algorithm, simultaneously calculating the attenuated latest shear modulus and bulk modulus based on the material's real-time stress-strain position and discrimination criteria using the aforementioned modulus reduction equation. This leads to the formation of a uniform tangent stiffness operator for the local continuum, and the calculation of the stress increment matrix at the current integration point. The energy dissipation caused by projection onto the plastic yield surface is considered. The main solver integrates the local stiffness of each microscopic Gaussian point based on the shape function integration rule, completing the assembly and replacement of the macroscopic overall stiffness matrix. At the system level, the unbalanced residual force equation of global displacement is eliminated based on the classical Newton-Raphson iterative solution strategy. When a single incremental step converges, the relevant tensors are written to the results database, such as the .odb file. After all loading steps are completed, visualization post-processing software can be used to extract and render the triaxial settlement displacement field of the nodes, the Mises equivalent shear stress distribution map, the equivalent plastic strain zone trajectory line, and the volume expansion region at key locations of the model at each stage, thereby reproducing and predicting the ultimate bearing capacity state and catastrophic sliding surface characteristics of tailings-based building materials within a given period.

[0042] This control and ablation experiment used conventional consolidated drained triaxial compression tests on tailings-based materials as the verification object. The experimental environment was set at a confining pressure of 400 kPa, an initial void ratio of 0.85, and an eight-node linear reduced integral hexahedral solid element mesh. The loading step was strictly controlled at 0.01% axial strain increment. Three simulation groups were used: a conventional control group using an ideal elastoplastic constitutive model without considering stiffness degradation; a constitutive ablation group using only particle interlocking characteristic modulus attenuation but without three-branch plastic flow correction; and a comprehensive verification group using a complete scheme with modulus attenuation and three-branch plastic strain increment correction based on dilatation response sequence. All models were continuously loaded to 20% axial strain, and stress and strain change data were extracted simultaneously. The effects of different simulation schemes and physical tests on the change of deviatoric stress with axial strain in tailings-based materials are compared as follows: Figure 4 As shown.

[0043] Actual indoor data from conventional consolidated drained triaxial compression tests showed that the peak deviatoric stress of the tailings-based sample was 950 kPa, the residual deviatoric stress after the peak was 720 kPa, and the maximum volumetric expansion strain reached 1.50%. Numerical simulation data indicated that the traditional control group calculated a peak deviatoric stress of 1050 kPa without post-peak softening characteristics, with the residual deviatoric stress maintained at 1020 kPa and the maximum volumetric expansion strain only 0.30%. The constitutive ablation group calculated a peak deviatoric stress reduced to 980 kPa, a residual deviatoric stress reduced to 760 kPa, and a maximum volumetric expansion strain increased to 0.80%. The comprehensive verification group calculated a peak deviatoric stress of 955 kPa, a residual deviatoric stress of 725 kPa, and a maximum volumetric expansion strain of 1.45%, with prediction errors for all key mechanical indicators controlled within 5%.

[0044] Traditional control groups suffer from inaccurate stress-strain results due to their inability to represent the failure process of microstructures. The constitutive ablation group, by utilizing modulus decay, reproduces the weakening of stiffness and stress reduction after macroscopic stress on materials, but still exhibits significant bias in predicting volumetric strain rates. The comprehensive validation group, by integrating a three-branch plastic strain increment correction rule, activates anisotropic compensation during the transition from volume contraction to volume expansion and utilizes a non-correlated correction function during the strong shear dilatation period. This aligns with the physical rearrangement mechanism when the stress ratio exceeds the critical point, enabling the constitutive model to detect the deep shear dilatation characteristics of tailings building materials.

[0045] Example 2 Embodiment 2 of the present invention proposes a simulation system based on tailings-based material parameters, comprising the following modules: The extraction module is used to acquire axial strain, volumetric strain, mean stress, and deviatoric stress test data of tailings-based materials under different confining pressures, water-bearing states, and loading paths. It calculates the ratio of deviatoric stress to mean stress and the ratio of plastic volumetric strain increment to plastic shear strain increment, and constructs a shear dilatation response sequence. Based on the shear dilatation response sequence, it identifies the shear dilatation inflection point, zero shear dilatation point, and post-peak fallback segment, and extracts the corresponding slope, curvature, and segment length features. The correction module is used to establish the correlation between state parameters, critical state indices and plastic flow direction based on the characteristic quantities, and to construct a dilatation function that varies with confining pressure, water-bearing state and loading path; to calculate the path memory factor based on the changes in the zero dilatation point under different confining pressures and loading paths, and to construct historical variables in combination with the water-bearing state for correcting the critical state indices; when the loading path is switched, to determine the structural attenuation coefficient based on the change in the ratio of plastic work density to reference energy density of adjacent path segments, and to correct the yield surface scale parameters, hardening parameters and flow parameters. The update module is used to determine the particle interlocking representation based on the normalized cumulative shear dilatation strain of the current calculation step and the cumulative volumetric strain characteristics that distinguish between shrinkage and swelling states. It also updates the shear modulus and bulk modulus in conjunction with the structural attenuation coefficient and the compaction hardening factor. Within each stress step, it selects the corresponding plastic strain increment correction rule based on the position of the trial stress relative to the shear dilatation inflection point and the zero shear dilatation point, and updates the strength parameters, deformation parameters, and flow parameters. Based on the updated parameters, it performs numerical simulation of tailings-based materials to obtain the stress field, strain field, and volumetric strain response results under given boundary conditions.

[0046] In an optional embodiment, the acquisition of axial strain, volumetric strain, mean stress, and deviatoric stress test data of tailings-based materials under different confining pressures, water-bearing states, and loading paths, and the calculation of the ratio of deviatoric stress to mean stress, the ratio of plastic volumetric strain increment to plastic shear strain increment, and the construction of a shear dilatation response sequence, includes: Within each loading step, the triaxial principal stress data and corresponding principal strain data of the tailings-based material are read; The mean stress and deviatoric stress are calculated based on the triaxial principal stress data, and the stress ratio is obtained by dividing the deviatoric stress by the mean stress. The plastic volumetric strain increment and the plastic shear strain increment are calculated based on the principal strain data, and the shear dilatation ratio is obtained by dividing the plastic volumetric strain increment by the plastic shear strain increment. The stress ratio and shear dilatation ratio corresponding to each loading step are arranged in the loading order to construct a shear dilatation response sequence that represents the relationship between stress state and volumetric deformation.

[0047] In an optional embodiment, the step of identifying the dilatation inflection point, zero dilatation point, and post-peak pullback segment based on the dilatation response sequence, and extracting the corresponding slope, curvature, and segment length features, includes: The derivative of the curve of shear dilatation ratio versus shear strain in the shear dilatation response sequence is taken, and the point where the derivative turns from negative to positive and the shear dilatation ratio reaches a minimum value is identified as the shear dilatation inflection point; Traverse the shear dilatation response sequence and identify the stress state point corresponding to the shear dilatation ratio crossing zero as the zero shear dilatation point; Identify the stress-strain softening interval after the deviatoric stress reaches its peak value, and define the interval as the post-peak decline segment; The tangent slopes before and after the dilatation inflection point, the curvature near the zero dilatation point, and the shear strain span corresponding to the post-peak fall-off segment are extracted as segment length features.

[0048] In an optional embodiment, the step of determining the particle interlocking representation based on the normalized cumulative shear dilatation strain of the current calculation step and the cumulative volumetric strain characteristic that distinguishes between shrinkage and swelling states, and updating the shear modulus and bulk modulus in conjunction with the structural attenuation coefficient and the compaction hardening factor, includes: Calculate the normalized cumulative shear dilatation strain of the current calculation step, and extract the attenuation function value constructed by accumulating the plastic volumetric strain increment when the current calculation step is in the dilatation state as the cumulative volumetric strain feature. At the same time, construct the compaction hardening factor based on the volumetric strain increment when the current calculation step is in the contraction state. The normalized cumulative shear dilatation strain and the cumulative volumetric strain feature are weighted and summed to obtain a particle interlocking representation that represents the degree of damage to the internal particle skeleton reconstruction. The initial shear modulus and initial bulk modulus are multiplied by a non-negative attenuation coefficient constructed from the particle interlocking representation and the structural attenuation coefficient, and further multiplied by the compaction hardening factor to obtain the updated shear modulus and bulk modulus under the current stress step.

[0049] In an optional embodiment, the step of selecting the corresponding plastic strain increment correction rule based on the position of the trial stress relative dilatation inflection point and the zero dilatation point within each stress step, and updating the strength parameters, deformation parameters, and flow parameters, includes: Compare the calculated stress ratio of the current stress step with the critical stress ratio corresponding to the dilatation inflection point and the zero dilatation point; If the calculated stress ratio is less than the critical stress ratio corresponding to the dilatation inflection point, it is determined that the body is in the shrinkage stage, and the initial plastic strain increment is calculated using the associated flow rule. If the calculated stress ratio is between the critical stress ratio corresponding to the dilatation inflection point and the zero dilatation point, it is determined that the body is in the transition stage from contraction to expansion, and the anisotropic tensor is used to perform non-associated flow compensation on the initial plastic strain increment. If the calculated stress ratio is not less than the critical stress ratio corresponding to the zero dilatation point, it is determined that it is in the strong dilatation stage. The plastic flow direction is updated according to the corrected critical state index and the dilatation function, and the plastic strain increment is recalculated.

[0050] In an optional embodiment, the numerical simulation of tailings-based materials based on the updated parameters to obtain the stress field, strain field, and bulk strain response results under given boundary conditions includes: A three-dimensional mesh geometric model of the tailings-based material is established, and initial void ratio, initial stress field and boundary conditions are assigned. Input the initial shear modulus, initial bulk modulus, strength parameters, flow parameters, and variation rules as constitutive model parameters into the finite element calculation program; Within each computational increment step, the constitutive model parameters are called, the shear modulus and bulk modulus are updated according to the current increment step state, the stress increment at the Gaussian integral point is calculated, and the overall stiffness matrix is ​​updated. Solve the global equilibrium equations and output the displacement, stress, strain and volumetric strain response results under each loading stage to complete the numerical simulation of tailings-based materials.

[0051] It should be understood that in the foregoing description of the embodiments in this specification, various features are combined in a single embodiment, drawing, or description for the purpose of simplifying the description and to aid in understanding a feature. However, this does not mean that the combination of these features is necessary, and those skilled in the art, upon reading this specification, may readily identify some of the devices as separate embodiments. That is, the embodiments in this specification can also be understood as an integration of multiple secondary embodiments. And the content of each secondary embodiment is valid even if it contains fewer than all the features of a single foregoing disclosed embodiment.

Claims

1. A simulation method based on tailings-based material parameters, characterized in that, include: The test data of axial strain, volumetric strain, mean stress and deviatoric stress of tailings-based materials under different confining pressures, water content and loading paths were obtained. The ratio of deviatoric stress to mean stress and the ratio of plastic volumetric strain increment to plastic shear strain increment were calculated to construct a shear dilatation response sequence. Based on the shear dilatation response sequence, the shear dilatation inflection point, zero shear dilatation point and the post-peak fallback segment were identified, and the corresponding slope, curvature and segment length features were extracted. Based on the aforementioned characteristic quantities, the correlation between state parameters, critical state indices, and plastic flow direction is established, and a dilatation function that varies with confining pressure, water-bearing state, and loading path is constructed. The path memory factor is calculated based on the changes in the zero dilatation point under different confining pressures and loading paths, and historical variables are constructed in conjunction with the water-bearing state to correct the critical state indices. When the loading path is switched, the structural attenuation coefficient is determined based on the change in the ratio of plastic work density to reference energy density of adjacent path segments, and the yield surface scale parameters, hardening parameters, and flow parameters are corrected. Based on the normalized cumulative shear dilatation strain of the current calculation step and the cumulative volumetric strain characteristics that distinguish between shrinkage and swelling states, the particle interlocking representation is determined. The shear modulus and bulk modulus are updated together with the structural attenuation coefficient and the compaction hardening factor. Within each stress step, the corresponding plastic strain increment correction rule is selected according to the position of the trial stress relative to the shear dilatation inflection point and the zero shear dilatation point, and the strength parameters, deformation parameters, and flow parameters are updated. Based on the updated parameters, numerical simulation of tailings-based materials is performed to obtain the stress field, strain field, and volumetric strain response results under given boundary conditions.

2. The method according to claim 1, characterized in that, The process involves acquiring experimental data on axial strain, volumetric strain, mean stress, and deviatoric stress of tailings-based materials under different confining pressures, water content, and loading paths; calculating the ratio of deviatoric stress to mean stress and the ratio of plastic volumetric strain increment to plastic shear strain increment; and constructing a shear dilatation response sequence, including: Within each loading step, the triaxial principal stress data and corresponding principal strain data of the tailings-based material are read; The mean stress and deviatoric stress are calculated based on the triaxial principal stress data, and the stress ratio is obtained by dividing the deviatoric stress by the mean stress. The plastic volumetric strain increment and the plastic shear strain increment are calculated based on the principal strain data, and the shear dilatation ratio is obtained by dividing the plastic volumetric strain increment by the plastic shear strain increment. The stress ratio and shear dilatation ratio corresponding to each loading step are arranged in the loading order to construct a shear dilatation response sequence that represents the relationship between stress state and volumetric deformation.

3. The method according to claim 1, characterized in that, The step of identifying the dilatation inflection point, zero dilatation point, and post-peak fallback segment based on the dilatation response sequence, and extracting the corresponding slope, curvature, and segment length features, includes: The derivative of the curve of shear dilatation ratio versus shear strain in the shear dilatation response sequence is taken, and the point where the derivative turns from negative to positive and the shear dilatation ratio reaches a minimum value is identified as the shear dilatation inflection point; Traverse the shear dilatation response sequence and identify the stress state point corresponding to the shear dilatation ratio crossing zero as the zero shear dilatation point; Identify the stress-strain softening interval after the deviatoric stress reaches its peak value, and define the interval as the post-peak decline segment; The tangent slopes before and after the dilatation inflection point, the curvature near the zero dilatation point, and the shear strain span corresponding to the post-peak fallback segment are extracted as segment length features.

4. The method according to claim 1, characterized in that, The process of determining particle interlocking representations based on the normalized cumulative shear dilatation strain of the current calculation step and the cumulative volumetric strain characteristics that distinguish between shrinkage and swelling states, and updating the shear modulus and bulk modulus in conjunction with the structural attenuation coefficient and compaction hardening factor, includes: Calculate the normalized cumulative shear dilatation strain of the current calculation step, and extract the attenuation function value constructed by accumulating the plastic volumetric strain increment when the current calculation step is in the dilatation state as the cumulative volumetric strain feature. At the same time, construct the compaction hardening factor based on the volumetric strain increment when the current calculation step is in the contraction state. The normalized cumulative shear dilatation strain and the cumulative volumetric strain feature are weighted and summed to obtain a particle interlocking representation that represents the degree of damage to the internal particle skeleton reconstruction. The initial shear modulus and initial bulk modulus are multiplied by a non-negative attenuation coefficient constructed from the particle interlocking representation and the structural attenuation coefficient, and further multiplied by the compaction hardening factor to obtain the updated shear modulus and bulk modulus under the current stress step.

5. The method according to claim 1, characterized in that, The step of selecting the corresponding plastic strain increment correction rule based on the position of the trial stress relative dilatation inflection point and the zero dilatation point within each stress step, and updating the strength parameters, deformation parameters, and flow parameters, includes: Compare the calculated stress ratio of the current stress step with the critical stress ratio corresponding to the dilatation inflection point and the zero dilatation point; If the calculated stress ratio is less than the critical stress ratio corresponding to the dilatation inflection point, it is determined that the body is in the shrinkage stage, and the initial plastic strain increment is calculated using the associated flow rule. If the calculated stress ratio is between the critical stress ratio corresponding to the dilatation inflection point and the zero dilatation point, it is determined that the body is in the transition stage from contraction to expansion, and the anisotropic tensor is used to perform non-associated flow compensation on the initial plastic strain increment. If the calculated stress ratio is not less than the critical stress ratio corresponding to the zero dilatation point, it is determined that it is in the strong dilatation stage. The plastic flow direction is updated according to the corrected critical state index and the dilatation function, and the plastic strain increment is recalculated.

6. The method according to claim 1, characterized in that, The numerical simulation of tailings-based materials based on the updated parameters yields the stress field, strain field, and bulk strain response results under given boundary conditions, including: A three-dimensional mesh geometric model of the tailings-based material is established, and initial void ratio, initial stress field and boundary conditions are assigned. Input the initial shear modulus, initial bulk modulus, strength parameters, flow parameters, and variation rules as constitutive model parameters into the finite element calculation program; Within each computational increment step, the constitutive model parameters are called, the shear modulus and bulk modulus are updated according to the current increment step state, the stress increment at the Gaussian integral point is calculated, and the overall stiffness matrix is ​​updated. Solve the global equilibrium equations and output the displacement, stress, strain and volumetric strain response results under each loading stage to complete the numerical simulation of tailings-based materials.

7. A simulation system based on tailings-based material parameters, characterized in that, Includes the following modules: The extraction module is used to acquire axial strain, volumetric strain, mean stress, and deviatoric stress test data of tailings-based materials under different confining pressures, water-bearing states, and loading paths. It calculates the ratio of deviatoric stress to mean stress and the ratio of plastic volumetric strain increment to plastic shear strain increment, and constructs a shear dilatation response sequence. Based on the shear dilatation response sequence, it identifies the shear dilatation inflection point, zero shear dilatation point, and post-peak fallback segment, and extracts the corresponding slope, curvature, and segment length features. The correction module is used to establish the correlation between state parameters, critical state indices and plastic flow direction based on the characteristic quantities, and to construct a dilatation function that varies with confining pressure, water-bearing state and loading path; to calculate the path memory factor based on the changes in the zero dilatation point under different confining pressures and loading paths, and to construct historical variables in combination with the water-bearing state for correcting the critical state indices; when the loading path is switched, to determine the structural attenuation coefficient based on the change in the ratio of plastic work density to reference energy density of adjacent path segments, and to correct the yield surface scale parameters, hardening parameters and flow parameters. The update module is used to determine the particle interlocking representation based on the normalized cumulative shear dilatation strain of the current calculation step and the cumulative volumetric strain characteristics that distinguish between shrinkage and swelling states. It also updates the shear modulus and bulk modulus in conjunction with the structural attenuation coefficient and the compaction hardening factor. Within each stress step, it selects the corresponding plastic strain increment correction rule based on the position of the trial stress relative to the shear dilatation inflection point and the zero shear dilatation point, and updates the strength parameters, deformation parameters, and flow parameters. Based on the updated parameters, it performs numerical simulation of tailings-based materials to obtain the stress field, strain field, and volumetric strain response results under given boundary conditions.

8. The system according to claim 7, characterized in that, The process involves acquiring experimental data on axial strain, volumetric strain, mean stress, and deviatoric stress of tailings-based materials under different confining pressures, water content, and loading paths; calculating the ratio of deviatoric stress to mean stress and the ratio of plastic volumetric strain increment to plastic shear strain increment; and constructing a shear dilatation response sequence, including: Within each loading step, the triaxial principal stress data and corresponding principal strain data of the tailings-based material are read; The mean stress and deviatoric stress are calculated based on the triaxial principal stress data, and the stress ratio is obtained by dividing the deviatoric stress by the mean stress. The plastic volumetric strain increment and the plastic shear strain increment are calculated based on the principal strain data, and the shear dilatation ratio is obtained by dividing the plastic volumetric strain increment by the plastic shear strain increment. The stress ratio and shear dilatation ratio corresponding to each loading step are arranged in the loading order to construct a shear dilatation response sequence that represents the relationship between stress state and volumetric deformation.

9. The system according to claim 7, characterized in that, The step of identifying the dilatation inflection point, zero dilatation point, and post-peak fallback segment based on the dilatation response sequence, and extracting the corresponding slope, curvature, and segment length features, includes: The derivative of the curve of shear dilatation ratio versus shear strain in the shear dilatation response sequence is taken, and the point where the derivative turns from negative to positive and the shear dilatation ratio reaches a minimum value is identified as the shear dilatation inflection point; Traverse the shear dilatation response sequence and identify the stress state point corresponding to the shear dilatation ratio crossing zero as the zero shear dilatation point; Identify the stress-strain softening interval after the deviatoric stress reaches its peak value, and define the interval as the post-peak decline segment; The tangent slopes before and after the dilatation inflection point, the curvature near the zero dilatation point, and the shear strain span corresponding to the post-peak fallback segment are extracted as segment length features.

10. A computer-readable storage medium storing a computer program thereon, characterized in that, The computer program, when executed by a processor, implements the method as described in any one of claims 1-6.

Images