Freeze-thaw stratum grouting numerical analysis modeling method and system

Abstract

The invention discloses a freezing and thawing stratum grouting numerical analysis modeling method and system, particularly relates to the technical field of hydraulic engineering numerical simulation, and is used for solving the problem of grouting parameter prediction inaccuracy caused by neglecting the dynamic coupling relationship between the thermodynamic behavior of a phase-change material and a seepage field and a damage field in an existing grouting model. A dynamic correlation system of phase change energy storage performance parameters and seepage field-damage evolution is established, a slurry diffusion path priority channel weight is corrected based on a freeze-thaw interface migration track, a bidirectional iteration equation of a seepage field and a diffusion migration field is coupled, and multi-physical field spatial-temporal scale cooperative solution is realized through a separated iteration algorithm; through multi-dimensional coupling of thermodynamic characteristics of the phase-change material, heterogeneous evolution of seepage paths and damage mechanical response, the prediction precision of the mechanical characteristics of the grouting stone body and the slurry diffusion range under the action of freeze-thaw cycle is remarkably improved, and the leakage and structural instability risks are effectively reduced.

Description

Technical Field

[0001] The present invention relates to the technical field of numerical simulation of hydraulic engineering. More specifically, the present invention relates to a numerical analysis and modeling method and system for freeze-thaw stratum grouting. Background Art

[0002] In hydraulic engineering applicable to cold regions and with periodic freeze-thaw action (such as canals, reservoirs), the rock and soil mass undergoes dynamic expansion of seepage channels and cumulative structural damage due to temperature cycling. Traditional grouting models are mostly based on the assumption of a static seepage field, simplify the slurry diffusion as the flow of a homogeneous medium at a constant temperature, and use empirical constitutive relations to describe the mechanical behavior of the grouted stone body, without quantifying the feedback mechanism of the heat storage and release process of the phase change material on the seepage field and structural damage of the surrounding rock.

[0003] In the prior art, since the grouting model does not couple the thermodynamic behavior of the phase change material with the dynamic coordination relationship between the stratum seepage field and damage evolution, it will lead to inaccurate prediction of the non-uniform mechanical properties of the grouted stone body and the slurry diffusion range under low-temperature environments, resulting in a deviation between the grouting design parameters and the actual working conditions of hydraulic engineering with periodic freeze-thaw action, and aggravating the risks of leakage and structural instability. Summary of the Invention

[0004] In order to overcome the above-mentioned defects of the prior art, the embodiments of the present invention provide a numerical analysis and modeling method and system for freeze-thaw stratum grouting to solve the problems raised in the above background art.

[0005] To achieve the above object, the present invention provides the following technical solutions: A numerical analysis and modeling method for freeze-thaw stratum grouting, comprising the following steps: S1. Determine the material parameters of the phase change grouting composite material, where the material parameters include the phase change energy storage performance parameters of the composite slurry material; S2. Generate the correction parameters of the elastoplastic damage constitutive model according to the correlation relationship between the phase change energy storage performance parameters and the damage parameters of the grouted stone body; S3. Based on the correction parameters and the mechanical property parameters of the grouted stone body, construct a three-dimensional mechanical model corresponding to the elastoplastic damage constitutive model; S4. Establish a coupled model of the temperature field - seepage field - stress field of the freeze-thaw stratum based on the phase change energy storage performance parameters and obtain dynamic coupling data; S5. Analyze the influence of the freeze-thaw interface fluctuation on the slurry diffusion path according to the dynamic coupling data, and generate diffusion path optimization parameters; S6. Based on the diffusion path optimization parameters and the dynamic characteristic parameters of the interaction between the stratum and the slurry, construct a diffusion and migration model of the composite slurry; S7. Perform multi - physical field coupling analysis on the three - dimensional mechanical model, the temperature field - seepage field - stress field coupling model, and the diffusion migration model to obtain the numerical analysis model for grouting in frozen - thawed strata.

[0006] In a preferred embodiment, the material parameters further include the composition ratio parameters of the composite grouting material and the micro - structure parameters of the grouting stone body: Based on the mass proportion of each component of the composite grouting material, determine the composition ratio parameters of the composite grouting material; By observing the pore distribution characteristics of the solidified grouting stone body through a scanning electron microscope, extract the porosity, average pore diameter, and connectivity parameters as the micro - structure parameters of the grouting stone body.

[0007] In a preferred embodiment, the material parameters include the phase - change energy storage performance parameters of the composite grouting material. Through the phase - change material thermal cycle experiment, obtain the latent heat of phase - change value and the phase - change temperature range of the composite grouting material under different temperature gradients as the phase - change energy storage performance parameters.

[0008] In a preferred embodiment, generate the correction parameters of the elastoplastic damage constitutive model according to the correlation between the phase - change energy storage performance parameters and the damage parameters of the grouting stone body, including: Through the triaxial shear test, obtain the axial strain - stress curve and shear failure mode of the grouting stone body under different confining pressure conditions, and extract the shear strength, residual strength, and damage softening coefficient as the damage parameters of the grouting stone body; Based on the latent heat of phase - change value and the phase - change temperature range, establish a non - linear regression relationship between the phase - change energy storage performance parameters and the damage softening coefficient, and obtain the correction factor of the damage evolution rate varying with the latent heat of phase - change value through least - squares fitting; According to the porosity and connectivity parameters in the micro - structure parameters of the grouting stone body, set the weight coefficient of the correction factor when the porosity is lower than the preset porosity threshold to generate the correction parameters of the elastoplastic damage constitutive model including the influence of pore structure.

[0009] In a preferred embodiment, based on the correction parameters and the mechanical property parameters of the grouting stone body, construct a three - dimensional mechanical model corresponding to the elastoplastic damage constitutive model, including: Based on the porosity and average pore diameter in the micro - structure parameters of the grouting stone body, define the initial elastic modulus and Poisson's ratio of the grouting stone body as the mechanical property parameters; Input the damage softening coefficient, correction factor, and weight coefficient in the correction parameters into the elastoplastic damage constitutive model, and correlate the latent heat of phase - change value and porosity's attenuation effect on the elastic modulus through the damage evolution equation; The elastoplastic damage constitutive model is discretized into a three-dimensional grid model by the finite element method, and a three-dimensional mechanical model including the coupling effect of phase change energy storage and damage is generated based on the mechanical property parameters and the damage evolution equation.

[0010] In a preferred embodiment, a coupled model of the temperature field - seepage field - stress field of the frozen-thawed stratum is established based on the phase change energy storage performance parameters, and dynamic coupling data is obtained, including: A temperature field model of the frozen-thawed stratum is established based on the latent heat value of phase change and the phase change temperature range, and the non-linear variation relationship between the temperature gradient and the latent heat release rate of phase change is defined. A seepage field model is constructed based on the porosity and connectivity parameters in the microstructure parameters of the grouting stone body. The water conductivity coefficient of the stratum is adjusted by the porosity, and the priority channel weight of the water migration path is corrected by the connectivity parameter. Based on the initial elastic modulus and Poisson's ratio in the mechanical property parameters of the grouting stone body, combined with the temperature gradient change and seepage rate, the stress distribution of the surrounding rock is calculated through the coupled equation of the temperature field - seepage field - stress field. The temperature gradient change, seepage rate and stress distribution of the surrounding rock are output as dynamic coupling data.

[0011] In a preferred embodiment, the influence of the fluctuation of the frozen-thawed interface on the slurry diffusion path is analyzed based on the dynamic coupling data, and diffusion path optimization parameters are generated, including: Based on the temperature gradient change data, the dynamic fluctuation range and fluctuation frequency of the frozen-thawed interface are extracted, and the interface migration trajectory of ice-water phase change in the temperature field is identified by the gray threshold segmentation algorithm. Combined with the seepage rate data and the initial elastic modulus in the mechanical property parameters of the grouting stone body, the migration probability of the priority channel and the diffusion radius correction coefficient of the slurry diffusion path caused by the fluctuation of the frozen-thawed interface are calculated. Based on the interface migration trajectory, migration probability and diffusion radius correction coefficient, diffusion path optimization parameters including the priority channel weight and diffusion direction adjustment factor are generated.

[0012] In a preferred embodiment, a diffusion and migration model of the composite slurry is constructed based on the diffusion path optimization parameters and the dynamic characteristic parameters of the interaction between the stratum and the slurry, including: The priority channel weight and diffusion direction adjustment factor in the diffusion path optimization parameters are input into the porous medium seepage equation, and the permeability tensor and flow direction deflection angle of the slurry in the priority channel are defined. Based on the porosity in the microstructure parameters of the grouting stone body and the seepage rate data, combined with the time-varying characteristic equation of the slurry viscosity, the effective viscosity and flow resistance coefficient of the slurry in the heterogeneous pores are calculated. The finite volume method is used to couple and discretize the seepage equation of porous media and the viscosity time-varying characteristic equation to generate a composite slurry diffusion and migration model, and output slurry pressure distribution, velocity vector and diffusion front evolution data.

[0013] In a preferred embodiment, a multi-physical field coupling analysis is performed on the three-dimensional mechanical model, the temperature field-seepage field-stress field coupling action model and the diffusion and migration model to obtain a numerical analysis model for grouting in frozen-thawed strata, including: Convert the surrounding rock stress distribution data output by the three-dimensional mechanical model into node stress field data matching the spatio-temporal resolution of the temperature field-seepage field-stress field coupling action model, and define the stress-temperature feedback coefficient of the mechanical-thermal coupling interface; Based on the seepage rate data in the dynamic coupling data and the slurry pressure distribution data output by the diffusion and migration model, establish a two-way coupling equation for the seepage field-diffusion field, and iteratively update the seepage rate and slurry pressure through the time step synchronization rule; Input the stress-temperature feedback coefficient and the two-way coupling equation of the seepage field-diffusion field into the multi-field coupling solver, and use the split iteration algorithm to synchronously solve the three-dimensional mechanical model, the temperature field-seepage field-stress field coupling action model and the diffusion and migration model to generate a numerical analysis model for grouting in frozen-thawed strata.

[0014] On the other hand, the present invention provides a numerical analysis modeling system for grouting in frozen-thawed strata, including: Phase change material parameter determination module: Determine the material parameters of the phase change grouting composite material, and the material parameters include the phase change energy storage performance parameters of the composite slurry material; Constitutive correction parameter module: Generate correction parameters for the elastoplastic damage constitutive model according to the correlation between the phase change energy storage performance parameters and the damage parameters of the grouting stone body; Three-dimensional mechanical modeling module: Based on the correction parameters and the mechanical property parameters of the grouting stone body, construct a three-dimensional mechanical model corresponding to the elastoplastic damage constitutive model; Multi-field coupling modeling module: Establish a temperature field-seepage field-stress field coupling action model for frozen-thawed strata based on the phase change energy storage performance parameters and obtain dynamic coupling data; Diffusion path optimization module: Analyze the influence of the frozen-thawed interface fluctuation on the slurry diffusion path according to the dynamic coupling data, and generate diffusion path optimization parameters; Diffusion and migration modeling module: Based on the diffusion path optimization parameters and the dynamic characteristic parameters of the interaction between the stratum and the slurry, construct a diffusion and migration model of the composite slurry; Multi-field coupling analysis module: Perform multi-physical field coupling analysis on the three-dimensional mechanical model, the temperature field-seepage field-stress field coupling action model and the diffusion and migration model to obtain a numerical analysis model for grouting in frozen-thawed strata.

[0015] Compared with the prior art, the present invention has the following beneficial effects: 1. By establishing a dynamic correlation system between the phase change energy storage performance parameters, seepage field, and damage evolution, it is possible to accurately depict the dynamic interaction of temperature gradient changes, pore structure deterioration, and slurry diffusion paths during the freeze-thaw cycle. Through the three-dimensional mechanical reconstruction of the elastoplastic damage constitutive model and the real-time feedback of the freeze-thaw interface fluctuation effect, the full-cycle simulation prediction of the non-uniform mechanical properties of the grouting stone body and the slurry diffusion range is realized, significantly improving the model's analytical ability for non-steady seepage and structural damage evolution in freeze-thaw strata, and providing a numerical analysis basis closer to the actual working conditions for grouting design in cold-region water conservancy projects.

[0016] 2. Through the multi-field collaborative iterative calculation of the seepage field - stress field - temperature field, it breaks through the limitation of the simplified assumptions of the static model for the thermodynamic behavior of phase change materials, enabling the optimization parameters of the slurry diffusion path to be dynamically adjusted with the migration of the freeze-thaw interface. This adaptive model correction mechanism not only strengthens the correlation mapping between the microscopic structure and macroscopic mechanical response of the grouting stone body, but also forms a complete technical link from material performance parameter calibration to engineering grouting effect prediction through the deep integration of the diffusion migration model and multi-field coupling data, improving the reliability of grouting process optimization in complex hydro-thermal coupling environments, reducing leakage risks, and enhancing structural stability. BRIEF DESCRIPTION OF THE DRAWINGS

[0017] Figure 1 is a flowchart of a numerical analysis modeling method for grouting in freeze-thaw strata according to the present invention; Figure 2 is a structural schematic diagram of a numerical analysis modeling system for grouting in freeze-thaw strata according to the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0018] The following will clearly and completely describe the technical solutions in the embodiments of the present invention with reference to the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all of the embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0019] Example 1: Figure 1 A numerical analysis modeling method for grouting in freeze-thaw strata according to the present invention is given, which includes the following steps: S1. Determine the material parameters of the phase change grouting composite material, where the material parameters include the phase change energy storage performance parameters of the composite slurry material; S2. Generate the correction parameters of the elastoplastic damage constitutive model according to the correlation relationship between the phase change energy storage performance parameters and the damage parameters of the grouting stone body; S3. Based on the correction parameters and the mechanical property parameters of the grouting stone body, construct a three-dimensional mechanical model corresponding to the elastoplastic damage constitutive model; S4. Based on the phase change energy storage performance parameters, establish a coupled model of temperature field - seepage field - stress field for the frozen-thawed stratum and obtain dynamic coupling data; S5. Analyze the influence of the frozen-thawed interface fluctuation on the slurry diffusion path according to the dynamic coupling data, and generate optimized parameters for the diffusion path; S6. Based on the optimized parameters of the diffusion path and the dynamic characteristic parameters of the interaction between the stratum and the slurry, construct a diffusion and migration model for the composite slurry; S7. Conduct multi-physical field coupling analysis on the three-dimensional mechanical model, the temperature field - seepage field - stress field coupling model and the diffusion and migration model to obtain a numerical analysis model for grouting in the frozen-thawed stratum.

[0020] The material parameters also include the composition ratio parameters of the composite slurry material and the microscopic structure parameters of the grouting stone body.

[0021] Based on the mass ratio of each component of the composite slurry material, determine the composition ratio parameters of the composite slurry material; By observing the pore distribution characteristics of the solidified grouting stone body through a scanning electron microscope, extract the porosity, average pore size and connectivity parameters as the microscopic structure parameters of the grouting stone body.

[0022] S1. Determine the material parameters of the phase change grouting composite material. The material parameters include the phase change energy storage performance parameters of the composite slurry material, including: Through the phase change material thermal cycle experiment, obtain the latent heat of phase change value and the phase change temperature range of the composite slurry material under different temperature gradients as the phase change energy storage performance parameters; Obtain the phase change energy storage performance parameters of the composite slurry material through the phase change material thermal cycle experiment. Specifically, use differential scanning calorimetry to conduct heating and cooling cycle tests on the composite slurry material. Set the temperature change range to -30°C to 25°C, and the heating and cooling rate to 5°C / min. Record the heat flow curves of the heat absorption and heat release processes of the material in each cycle. Calculate the latent heat of phase change value according to the integral areas of the endothermic peak and exothermic peak in the heat flow curve, and extract the starting temperature, peak temperature and ending temperature of the phase change to form the phase change temperature range. The latent heat of phase change value and the phase change temperature range are used as the phase change energy storage performance parameters. For example, when the latent heat of phase change value is higher than 200 J / g, it is determined that the material has significant heat storage and release capacity, and when the phase change temperature range is -10°C to 5°C, it is adapted to the typical temperature fluctuation range of the frozen-thawed stratum.

[0023] Determine the composition ratio parameters based on the mass proportion of each component of the composite slurry material. Specifically, mix the cement-based material, phase change microcapsules, and water reducer in a preset ratio, and obtain the composition ratio parameters by weighing the mass of each component and calculating its percentage in the total mass. For example, when the proportion of the cement-based material is 60%, the phase change microcapsules are 30%, and the water reducer is 10%, record the ratio parameter as 60:30:10. The composition ratio parameters control the slurry curing rate through the proportion of the cement-based material, adjust the heat storage performance through the proportion of the phase change microcapsules, and affect the initial fluidity of the slurry through the proportion of the water reducer.

[0024] Observe the pore distribution characteristics of the grouting stone body after curing through a scanning electron microscope, and extract the porosity, average pore diameter, and connectivity parameters as the microstructural parameters of the grouting stone body. Specifically, cut the cured grouting stone body sample into thin slices of 5 mm×5 mm×2 mm, place it under the scanning electron microscope after sputtering, take microstructural images at a magnification of 500 to 2000 times, perform binary segmentation on the pores in the image using image processing software, calculate the proportion of the pore area to obtain the porosity, calculate the average value of the equivalent diameters of all pores as the average pore diameter, and judge the connectivity between pores through a pore skeleton extraction algorithm to generate the connectivity parameter. For example, when the proportion of the cement-based material is 60%, the porosity is about 12% to 15%, the average pore diameter is 20 to 50 μm, and the connectivity parameter is 0.3 to 0.6.

[0025] The temperature range of -30°C to 25°C and the heating and cooling rate of 5°C / min in the above differential scanning calorimetry correspond to the actual freeze-thaw cycle conditions in cold region projects. The multi-scale image acquisition at a magnification of 500 to 2000 times in the scanning electron microscope observation is used to adapt to the heterogeneous distribution characteristics of the micro-pores of the grouting stone body. The pore skeleton extraction algorithm is an optimization of the existing technology based on image gray threshold segmentation. Specifically, set the gray threshold range to 120 to 180 to distinguish pores from solid matrices.

[0026] S2. Generate the correction parameters of the elastoplastic damage constitutive model according to the correlation between the phase change energy storage performance parameters and the damage parameters of the grouting stone body, including: Obtain the axial strain-stress curve and shear failure mode of the grouting stone body under different confining pressure conditions through a triaxial shear test, and extract the shear strength, residual strength, and damage softening coefficient as the damage parameters of the grouting stone body; Based on the latent heat of phase change value and the phase change temperature range, establish a non-linear regression relationship between the phase change energy storage performance parameters and the damage softening coefficient, and obtain the correction factor of the damage evolution rate varying with the latent heat of phase change value through least squares fitting; According to the porosity and connectivity parameters in the microscopic structure parameters of the grouting stone body, set the weight coefficient of the correction factor when the porosity is lower than the preset porosity threshold, and generate the correction parameters of the elastoplastic damage constitutive model including the influence of the pore structure.

[0027] When generating the correction parameters of the elastoplastic damage constitutive model, it is realized through the following steps: Obtain the axial strain-stress curve and shear failure mode of the grouting stone body under different confining pressure conditions through triaxial shear tests. For example, specifically, use a rock mechanics testing machine to apply three groups of confining pressures of 0.5 MPa, 1.0 MPa, and 1.5 MPa to the grouting stone body specimen, apply axial pressure at a loading rate of 0.1 mm / min until the specimen fails, record the corresponding relationship curve of axial strain and stress, determine the shear strength according to the peak point of the curve, take the stress value when the stress drops to the stable stage as the residual strength, and calculate the ratio of the peak strength to the residual strength as the damage softening coefficient. The shear strength, residual strength, and damage softening coefficient are used as the damage parameters of the grouting stone body. For example, when the confining pressure is 1.0 MPa, the shear strength is 8.5 MPa, the residual strength is 3.2 MPa, and the damage softening coefficient is 0.38.

[0028] Based on the latent heat of phase change value and phase change temperature range in step S1, establish a non-linear regression relationship between the phase change energy storage performance parameter and the damage softening coefficient. Specifically, fit the test data of the damage softening coefficient of the grouting stone body corresponding to different latent heat of phase change values, and use the least squares method to construct a quadratic polynomial regression equation with the latent heat of phase change value as the independent variable and the damage softening coefficient as the dependent variable. Determine the correction factor of the damage evolution rate with the change of the latent heat of phase change value according to the change of the regression equation slope. For example, when the latent heat of phase change value is 200 J / g, the correction factor is 1.2, indicating that the improvement of the phase change energy storage performance increases the damage evolution rate by 20%, while when the latent heat of phase change value is 150 J / g, the correction factor drops to 0.9.

[0029] According to the porosity and connectivity parameters in the microscopic structure parameters of the grouting stone body in step S1, set the weight coefficient of the correction factor when the porosity is lower than the preset porosity threshold. Specifically, by statistically analyzing the correlation data between the damage softening coefficient of specimens with different porosities and the correction factor, determine the porosity threshold range. When the porosity is lower than this threshold, adjust the weight coefficient according to the negative correlation between the connectivity parameter and the correction factor to generate the correction parameters of the elastoplastic damage constitutive model considering the influence of pore structure. For example, the preset porosity threshold is 15%. When the porosity is 10% and the connectivity parameter is 0.4, the weight coefficient is set to 0.8, indicating that the effect of the correction factor is weakened by 20% under the conditions of low porosity and low connectivity. If the connectivity parameter is 0.6, the weight coefficient is adjusted to 1.1. The preset porosity threshold is determined by statistically analyzing the porosity distribution of grouting stone body specimens. For example, perform a normal distribution analysis on the porosity data of 30 groups of specimens, and take the mean minus twice the standard deviation as the lower limit of the threshold.

[0030] The confining pressure range of 0.5 MPa to 1.5 MPa in the triaxial shear test corresponds to the stress conditions of the surrounding rock from the shallow layer to the middle and deep layers of the frozen-thawed stratum. For example, the typical buried depth of 5 - 15 meters in the cold region channel lining project corresponds to the in-situ stress level. The least squares fitting is a well-known regression analysis method. In this embodiment, by limiting the quadratic polynomial form and taking the latent heat of phase change value as the core variable, the problem that the traditional linear regression cannot characterize the non-linear correlation between phase change energy storage and damage is solved. For example, the traditional method only correlates the phase change temperature and the strength parameter, while in this embodiment, the dynamic influence of the heat storage and release scale on damage evolution is quantified by introducing the latent heat of phase change value. In the method for jointly calibrating the porosity and connectivity parameters in setting the weight coefficient, the deviation that solely depends on the porosity is corrected by introducing the connectivity parameter. For example, when the porosity is the same but the difference in the connectivity parameter exceeds 0.2, the adjustment range of the weight coefficient is ±15%.

[0031] S3. Based on the correction parameters and the mechanical property parameters of the grouting stone body, construct a three-dimensional mechanical model corresponding to the elastoplastic damage constitutive model, including: Define the initial elastic modulus and Poisson's ratio of the grouting stone body as mechanical property parameters based on the porosity and average pore diameter in the microscopic structure parameters of the grouting stone body; Input the damage softening coefficient, correction factor, and weight coefficient in the correction parameters into the elastoplastic damage constitutive model, and correlate the latent heat of phase change value and the attenuation effect of porosity on the elastic modulus through the damage evolution equation; Use the finite element method to discretize the elastoplastic damage constitutive model into a three-dimensional grid model, and generate a three-dimensional mechanical model including the phase change energy storage-damage coupling effect based on the mechanical property parameters and the damage evolution equation.

[0032] When constructing a three-dimensional mechanical model corresponding to the elastoplastic damage constitutive model, it is achieved through the following steps: Based on the porosity and average pore diameter in the microscopic structure parameters of the grouting stone body, the initial elastic modulus and Poisson's ratio of the grouting stone body are defined as mechanical property parameters. Specifically, according to the negative correlation between porosity and elastic modulus, an empirical formula for elastic modulus is obtained by fitting the laboratory uniaxial compression test data. For example, when the porosity is 10%, the elastic modulus is 8.5 GPa, and for every 5% increase in porosity, the elastic modulus decreases by 15%. The average pore diameter is converted into a pore shape factor through the equivalent sphere model, and the initial Poisson's ratio is determined to be 0.25 by referring to the recommended value table of Poisson's ratio in the material handbook. The equivalent sphere model is a geometric simplification method that equates pores to spherical cavities, and the pore shape factor is calculated by the ratio of the average pore diameter to the sphere diameter.

[0033] Input the damage softening coefficient, correction factor, and weight coefficient in the correction parameters into the elastoplastic damage constitutive model, and relate the latent heat of phase change value and the attenuation effect of porosity on the elastic modulus through the damage evolution equation. Specifically, the damage evolution equation uses the damage softening coefficient as the reference elastic modulus attenuation rate, the correction factor as the non-linear gain coefficient of the latent heat of phase change value, and the weight coefficient adjusts the attenuation amplitude according to the connectivity parameter when the porosity is lower than the preset threshold. For example, when the damage softening coefficient is 0.38, the correction factor is 1.2, and the weight coefficient is 0.8, the elastic modulus attenuation rate is calculated as 0.38×1.2×0.8 = 0.365, indicating that the elastic modulus decays to 36.5% of the initial value under the current conditions. The damage evolution equation is calibrated through multiple sets of triaxial test data. For example, the goodness of fit of the linear regression between the elastic modulus attenuation rate of specimens with different latent heat of phase change values and the correction parameters under the same confining pressure needs to be greater than 0.9.

[0034] The elastoplastic damage constitutive model is discretized into a three-dimensional grid model using the finite element method. Specifically, a three-dimensional solid model is established according to the actual geometric size of the grouting stone body, and a hexahedral element mesh is divided with an element size ranging from 5 mm to 10 mm. The mesh density is determined based on the aggregate particle size distribution of the grouting stone body. For example, when the maximum aggregate particle size is 3 mm, the element size is set to 5 mm to ensure that at least 3 elements cover the aggregate. Input the initial elastic modulus, Poisson's ratio, and damage evolution equation parameters into the material property module of the finite element software, and solve the non-linear equations through the Newton-Raphson iterative algorithm to generate a three-dimensional mechanical model including the coupling effect of phase change energy storage-damage. The hexahedral element mesh division and the Newton-Raphson iterative method are conventional techniques in finite element analysis. In this embodiment, the non-homogeneous characteristics of the grouting stone body in the freeze-thaw stratum are adapted by limiting the element size range and the input parameter type.

[0035] The empirical relationship between porosity and elastic modulus is fitted based on laboratory compression test data. For example, uniaxial compression tests are conducted on 20 groups of specimens with different porosities, and an exponential function relationship between elastic modulus and porosity is obtained by fitting. The equivalent sphere model is a geometric model used in the prior art to simplify the pore shape. In this embodiment, the pore shape factor is calculated by substituting the average pore diameter into the model, and the initial parameters are determined by combining with the recommended Poisson's ratio table in the material handbook. The hexahedral element mesh generation and iterative calculation in the finite element method are conventional numerical methods. In this embodiment, by limiting the element size range (5 mm to 10 mm) and the input parameter type (parameters of the damage evolution equation), it is adapted to the heterogeneous mechanical properties of the grouted stone body in the freeze-thaw stratum.

[0036] The synergistic effect of the correction factor and the weight coefficient in the damage evolution equation is calibrated through multiple groups of control tests. For example, when the correction factor is 1.0 and the weight coefficient is 1.0, there is a baseline linear relationship between the elastic modulus decay rate and the damage softening coefficient; when the correction factor is increased to 1.2 due to the latent heat of phase change value, and the weight coefficient is reduced to 0.8 due to the porosity being lower than the threshold and the connectivity parameter being 0.4, the decay rate is adjusted to 1.2×0.8 = 0.96 times the baseline value, reflecting the dynamic coupling effect of phase change energy storage and pore structure. The unit size range is determined according to the characteristic size of the grouted stone body (such as the aggregate particle size of 1 - 3 mm) to ensure the balance between grid accuracy and calculation efficiency.

[0037] S4. Establish a coupled temperature field - seepage field - stress field model for the freeze - thaw stratum based on the phase change energy storage performance parameters and obtain dynamic coupling data, including: Establish a temperature field model for the freeze - thaw stratum based on the latent heat of phase change value and the phase change temperature range, and define the non - linear variation relationship of the temperature gradient with the latent heat release rate of phase change; Construct a seepage field model based on the porosity and connectivity parameters in the microscopic structure parameters of the grouted stone body. Adjust the formation water conductivity through the porosity and correct the priority channel weight of the water migration path through the connectivity parameter; Based on the initial elastic modulus and Poisson's ratio in the mechanical property parameters of the grouted stone body, combined with the temperature gradient change and seepage rate, calculate the surrounding rock stress distribution through the temperature field - seepage field - stress field coupling equation; Output the temperature gradient change, seepage rate, and surrounding rock stress distribution as dynamic coupling data.

[0038] When establishing a coupled model of the temperature field, seepage field, and stress field in freeze-thaw strata and obtaining dynamic coupling data, it is achieved through the following steps: Based on the latent heat of phase change value and the phase change temperature range in step S1, establish a temperature field model for the freeze-thaw strata. Specifically, input the latent heat of phase change value as a heat source term into the heat conduction equation, and set the latent heat release rate as a piecewise function in combination with the phase change temperature range. A non-linear release mode is adopted within the phase change temperature range, and a constant heat conduction coefficient is adopted outside the range. For example, when the phase change temperature range is from -10°C to 5°C, the latent heat release rate starts at -10°C and terminates at 5°C, and the release rate increases exponentially with the increase in temperature. The growth coefficient is determined by fitting the test data of differential scanning calorimetry.

[0039] Based on the porosity and connectivity parameters in the microscopic structure parameters of the grouting stone body in step S1, construct a seepage field model. Specifically, calculate the formation hydraulic conductivity through the porosity. The calculation formula is the modified Darcy's law, where the hydraulic conductivity is proportional to the square of the porosity. For example, when the porosity is 10%, the hydraulic conductivity is 1×10 -5 m / s. For every 5% increase in porosity, the hydraulic conductivity increases to 1.25×10 -5 m / s; correct the priority channel weight of the seepage path through the connectivity parameter. The calculation formula for the weight coefficient is the ratio of the connectivity parameter to the maximum connectivity. For example, when the maximum connectivity parameter is 0.8, the weight coefficient corresponding to the connectivity parameter of 0.6 is 0.6 / 0.8 = 0.75. The modified Darcy's law and the weight calculation formula are calibrated through the seepage test data of 10 groups of specimens with different porosities and connectivities.

[0040] Based on the initial elastic modulus and Poisson's ratio in the mechanical property parameters of the grouting stone body in step S3, combined with the temperature gradient change and seepage rate, calculate the stress field distribution. Specifically, establish a coupled equation of the temperature field, seepage field, and stress field, and correlate the thermal expansion strain caused by the temperature gradient, the change in pore water pressure caused by the seepage rate with the initial elastic modulus and Poisson's ratio, and solve the coupled equation through the finite element method. For example, when the initial elastic modulus is 8.5 GPa, the thermal expansion coefficient is 1.2×10 -5 / °C, and the thermal stress generated by every 1°C / m increase in the temperature gradient is 8.5 GPa×1.2×10 -5 ×1 = 0.102 MPa; the pore water pressure corresponding to the seepage rate of 0.1 m / s drops by 5 kPa, and the change in the net stress of the surrounding rock is calculated through the effective stress principle.

[0041] Output the temperature gradient change data, seepage rate data, and surrounding rock stress distribution data as dynamically coupled data, with the output format being three-dimensional field data in time series. For example, the temperature gradient data is updated to a matrix of grid node temperature values at each hourly time step, the seepage rate data is marked with the flow velocity direction and magnitude of each unit in vector field format, and the surrounding rock stress data is output as the von Mises stress value of each unit. The time step setting rule is that the temperature field and seepage field are updated synchronously every hour, and the stress field is iteratively calculated every 10 minutes to match the migration rate of the freeze-thaw interface.

[0042] S5. Analyze the influence of the freeze-thaw interface fluctuation on the slurry diffusion path based on the dynamically coupled data, and generate optimization parameters for the diffusion path, including: Based on the temperature gradient change data, extract the dynamic fluctuation range and frequency of the freeze-thaw interface, and identify the interface migration trajectory of ice-water phase change in the temperature field through the gray threshold segmentation algorithm. Combine the seepage rate data and the initial elastic modulus in the mechanical property parameters of the grouting stone body to calculate the migration probability of the preferential channel and the diffusion radius correction coefficient of the slurry diffusion path affected by the freeze-thaw interface fluctuation. Based on the interface migration trajectory, migration probability, and diffusion radius correction coefficient, generate optimization parameters for the diffusion path including the preferential channel weight and diffusion direction adjustment factor.

[0043] When analyzing the influence of the freeze-thaw interface fluctuation on the slurry diffusion path based on the dynamically coupled data, it is achieved through the following steps: Extract the dynamic fluctuation range and frequency of the freeze-thaw interface based on the temperature gradient change data in step S4. Specifically, convert the spatio-temporal distribution matrix of the temperature field output in step S4 into a grayscale image, and the gray value mapping rule is that the gray value increases by 10 for every 1°C / m increase in the temperature gradient. Identify the continuous region with gray values in the range of 50 - 200 as the interface migration trajectory of ice-water phase change through the gray threshold segmentation algorithm.

[0044] The gray threshold segmentation algorithm uses the adaptiveThreshold function of the OpenCV library, with the block size set to 11×11 pixels and the constant offset C = 2 to adapt to the gradient smoothing characteristics of the freeze-thaw interface. For example, when the temperature gradient is -10°C / m, the corresponding gray value is 50, and when it is 5°C / m, the corresponding gray value is 200. After segmentation, the interface trajectory is stored in the form of a polygon contour line.

[0045] Calculate the migration probability of the preferential channel and the diffusion radius correction coefficient of the slurry diffusion path due to the freeze-thaw interface fluctuation by combining the seepage rate data in step S4 and the initial elastic modulus in the mechanical property parameters of the grouting stone body in step S3. Specifically, the migration probability calculation formula is the product of the seepage rate (unit: m / s) and the initial elastic modulus (unit: GPa) divided by 10. For example, when the seepage rate is 0.1 m / s and the initial elastic modulus is 8.5 GPa, the migration probability = (0.1×8.5) / 10 = 0.085; the diffusion radius correction coefficient calculation formula is 1 / (1 + migration probability). For example, the migration probability of 0.085 corresponds to a correction coefficient of 1 / 1.085 ≈ 0.922, indicating that the diffusion radius remains 92.2%. The formula coefficient (divided by 10) is determined by the regression analysis of 20 groups of grouting test data, and the regression equation is y = 0.1x (R² = 0.91).

[0046] Generate the diffusion path optimization parameters based on the interface migration trajectory, migration probability, and diffusion radius correction coefficient. Specifically, superimpose the interface migration trajectory polygon on the seepage preferential channel vector field in step S4, calculate the spatial intersection area between the two, and the preferential channel weight in the intersection area = migration probability × seepage preferential channel weight. For example, when the migration probability is 0.085 and the seepage weight is 0.8, the preferential channel weight = 0.085×0.8 = 0.068; the diffusion direction adjustment factor = correction coefficient × the angle of the interface migration trajectory tangent direction. For example, when the correction coefficient is 0.922 and the tangent direction angle is 30°, the adjusted direction = 30°×0.922 ≈ 27.7°. The optimization parameters are stored in JSON format, including the weight value and direction angle value of each finite element grid unit.

[0047] S6. Based on the diffusion path optimization parameters and the dynamic characteristic parameters of the interaction between the formation and the slurry, construct a diffusion and migration model of the composite slurry, including: Input the preferential channel weight and diffusion direction adjustment factor in the diffusion path optimization parameters into the porous medium seepage equation, and define the permeability tensor and flow direction deflection angle of the slurry in the preferential channel; Based on the porosity in the microscopic structure parameters of the grouting stone body and the seepage rate data, combined with the time-varying characteristic equation of the slurry viscosity, calculate the effective viscosity and flow resistance coefficient of the slurry in the heterogeneous pores; Use the finite volume method to couple and discretize the porous medium seepage equation and the time-varying characteristic equation of viscosity to generate a diffusion and migration model of the composite slurry, and output the slurry pressure distribution, velocity vector, and diffusion front evolution data.

[0048] When constructing the diffusion and migration model of the composite slurry, it is achieved through the following steps: the priority channel weight and diffusion direction adjustment factor in the diffusion path optimization parameters of step S5 are input into the porous media seepage equation, specifically, the anisotropy coefficient of the permeability tensor is adjusted according to the priority channel weight, and the diffusion direction adjustment factor is converted into the flow direction deflection angle. For example, when the priority channel weight is 0.8, the main direction component of the permeability tensor is increased to 1.8 times the baseline value, and the secondary direction component remains at the baseline value; when the diffusion direction adjustment factor is 30°, the flow direction deflection angle is set to rotate 30° clockwise along the main permeability direction. The permeability tensor adjustment rule is determined by fitting 10 groups of seepage test data with different priority channel weights. For example, for every increase of 0.1 in the weight, the main direction permeability increases by 20%.

[0049] Based on the porosity in the microstructural parameters of the grouting stone body in step S1 and the seepage rate data in step S4, the effective viscosity and flow resistance coefficient of the slurry in the heterogeneous pores are calculated in combination with the time-varying characteristic equation of the slurry viscosity. Specifically, according to the negative correlation between porosity and effective viscosity, an exponential decay function is established in which the effective viscosity decreases with the porosity, and the seepage rate is used as a linear correction term for the flow resistance coefficient. For example, when the porosity is 10%, the effective viscosity is 500 mPa·s, and the effective viscosity decreases to 70% of the previous value for every 5% increase in porosity; when the seepage rate is 0.1 m / s, the flow resistance coefficient = the benchmark resistance coefficient × (1 + 0.1 × seepage rate), if the benchmark coefficient is 1.2, then the corrected resistance coefficient = 1.2 × 1.01 = 1.212. The exponential decay function and the linear correction term coefficient are determined by regression analysis of 20 sets of slurry flow test data, and the goodness of fit is required to be greater than 0.85.

[0050] The finite volume method is used to couple the porous media seepage equation with the viscosity time-varying characteristic equation for discretization. Specifically, the computational domain is divided into hexahedral control volumes. The unit size is consistent with the three-dimensional mechanical model grid of step S3 (5 mm to 10 mm). The seepage equation and the viscosity equation are solved simultaneously in each control volume, and the slurry pressure distribution, velocity vector and diffusion front evolution data are output. For example, the pressure distribution data is stored as a three-dimensional scalar field in Pascals, the velocity vector field records the velocity magnitude and direction of each unit center point, and the diffusion front evolution data is generated by tracking the unit boundaries where the slurry volume fraction is greater than 0.5. The SIMPLE algorithm is used to handle the pressure-velocity coupling in the finite volume method discretization process, and the convergence residual is set to 1×10 -4 , the time step is 1 second, and the total simulation time is set to 60 minutes according to the actual grouting process.

[0051] S7. Perform multi-physics field coupling analysis on the three-dimensional mechanical model, the temperature field-seepage field-stress field coupling model and the diffusion migration model to obtain a numerical analysis model for grouting in frozen-thaw formations, including: Convert the surrounding rock stress distribution data output by the three-dimensional mechanical model into node stress field data that matches the spatio-temporal resolution of the temperature field-seepage field-stress field coupling model, and define the stress-temperature feedback coefficient of the mechanical-thermal coupling interface; Based on the seepage rate data in the dynamic coupling data of step S4 and the slurry pressure distribution data output by the diffusion migration model, establish a two-way coupling equation for the seepage field-diffusion field, and iteratively update the seepage rate and slurry pressure through the time step synchronization rule; Input the stress-temperature feedback coefficient and the two-way coupling equation of the seepage field-diffusion field into the multi-field coupling solver, and use the segregated iterative algorithm to synchronously solve the three-dimensional mechanical model, the temperature field-seepage field-stress field coupling model, and the diffusion migration model to generate a numerical analysis model for grouting in frozen-thawed strata.

[0052] When performing multi-physical field coupling analysis on the three-dimensional mechanical model, the temperature field-seepage field-stress field coupling model, and the diffusion migration model, it is achieved through the following steps: Convert the surrounding rock stress distribution data output by the three-dimensional mechanical model in step S3 into node stress field data that matches the spatio-temporal resolution of the temperature field-seepage field-stress field coupling model in step S4. Specifically, perform bilinear interpolation on the hexahedral mesh element stress data of the three-dimensional mechanical model and map it to the tetrahedral mesh nodes of the temperature field-seepage field model. Define the stress-temperature feedback coefficient of the mechanical-thermal coupling interface as the temperature change corresponding to every 1 MPa increase in the surrounding rock stress. For example, when the surrounding rock stress increases from 5 MPa to 6 MPa and the temperature rises by 0.2 °C, the stress-temperature feedback coefficient is set to 0.2 °C / MPa. The interpolation method is calculated based on the finite element shape function to ensure that the spatial position deviation between the stress field and the temperature field mesh nodes is less than 0.1 mm.

[0053] Establish a two-way coupling equation for the seepage field-diffusion field based on the seepage rate data in the dynamic coupling data of step S4 and the slurry pressure distribution data output by the diffusion migration model in step S6. Specifically, input the seepage rate as the boundary condition of the diffusion field model, and at the same time, feedback the slurry pressure distribution data as the source term of the seepage field model, and iteratively update the seepage rate and slurry pressure through the time step synchronization rule. For example, it is set that the seepage rate of the seepage field is updated every 1 hour, and the slurry pressure of the diffusion field is updated every 10 minutes, and the data exchange is synchronized at the whole hour. The two-way coupling equation is constructed based on the law of conservation of mass. For example, the amount of water flowing out in the seepage field is equal to the increase in the pore volume occupied by the slurry in the diffusion field, and the mass conservation error tolerance is set to 5%.

[0054] The stress-temperature feedback coefficient and the bidirectional coupling equation of the seepage field-diffusion field are input into the multi-field coupling solver. The separated iterative algorithm is used to synchronously solve the three-dimensional mechanical model, the temperature field-seepage field-stress field coupling model, and the diffusion migration model. Specifically, within each time step, the equations of the temperature field, seepage field, stress field, and diffusion field are solved in sequence, and the calculation results of the previous step are used as the initial conditions for the next step. The loop iteration continues until the residuals of all physical fields are less than the set threshold. For example, the residual threshold of the temperature field is set to 0.1 °C, the residual threshold of the seepage field flow velocity is set to 0.01 m / s, the residual threshold of the stress field displacement is set to 0.01 mm, and the residual threshold of the diffusion field concentration is set to 0.5%. The separated iterative algorithm uses the Gauss-Seidel iterative format, and the relaxation factor is set to 0.8 to balance the convergence rate and stability.

[0055] Embodiment 2: Figure 2 The structural schematic diagram of a numerical analysis and modeling system for freeze-thaw stratum grouting according to the present invention is given. A numerical analysis and modeling system for freeze-thaw stratum grouting includes: Phase change material parameter determination module: Determine the material parameters of the phase change grouting composite material, and the material parameters include the phase change energy storage performance parameters of the composite slurry material; Constitutive modification parameter module: Generate the modification parameters of the elastoplastic damage constitutive model according to the correlation between the phase change energy storage performance parameters and the damage parameters of the grouting stone body; Three-dimensional mechanics modeling module: Based on the modification parameters and the mechanical property parameters of the grouting stone body, construct a three-dimensional mechanical model corresponding to the elastoplastic damage constitutive model; Multi-field coupling modeling module: Based on the phase change energy storage performance parameters, establish a coupling model of the temperature field-seepage field-stress field in the freeze-thaw stratum and obtain dynamic coupling data; Diffusion path optimization module: Analyze the influence of the freeze-thaw interface fluctuation on the slurry diffusion path according to the dynamic coupling data, and generate diffusion path optimization parameters; Diffusion migration modeling module: Based on the diffusion path optimization parameters and the dynamic characteristic parameters of the interaction between the stratum and the slurry, construct a diffusion migration model of the composite slurry; Multi-field coupling analysis module: Perform multi-physical field coupling analysis on the three-dimensional mechanical model, the temperature field-seepage field-stress field coupling model, and the diffusion migration model to obtain a numerical analysis model for freeze-thaw stratum grouting.

[0056] The above embodiments can be implemented in whole or in part by software, hardware, firmware, or any other combination. When implemented using software, the above embodiments can be implemented in whole or in part in the form of a computer program product.

[0057] Those of ordinary skill in the art will realize that the modules and algorithm steps of each example described in combination with the embodiments disclosed herein can be implemented by electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are executed in a hardware or software manner depends on the specific application of the technical solution and the inventive constraints. Professional technicians can use different methods to implement the described functions for each specific application, but such implementation should not be considered to exceed the scope of this application.

[0058] In addition, the functional modules in each embodiment of this application can be integrated into one processing module, or each module can exist physically alone, or two or more modules can be integrated into one module.

[0059] As described above, this is only the specific implementation manner of this application, but the protection scope of this application is not limited thereto. Any person skilled in the art within the technical scope disclosed in this application can easily think of changes or substitutions, which should all be covered by the protection scope of this application. Therefore, the protection scope of this application should be subject to the protection scope of the claims.

[0060] Finally: The above are only the preferred embodiments of the present invention and are not used to limit the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention shall be included within the protection scope of the present invention.

Claims

1. A numerical analysis modeling method for freeze-thaw stratum grouting, characterized in that: The steps include: S1. Determine the material parameters of the phase change grouting composite material, the material parameters including the phase change energy storage performance parameters of the composite slurry material; S2, generating correction parameters of the elastic-plastic damage constitutive model according to the correlation between the phase change energy storage performance parameters and the damage parameters of the grouting stone body; S3, constructing a three-dimensional mechanical model corresponding to the elastic-plastic damage constitutive model based on the correction parameters and the mechanical characteristic parameters of the grouting stone body; S4. Establish a coupling model of temperature field-seepage field-stress field in freeze-thaw formations based on phase change energy storage performance parameters and obtain dynamic coupling data; S5. Analyze the influence of freeze-thaw interface fluctuation on the slurry diffusion path based on dynamic coupling data, and generate diffusion path optimization parameters; S6. Based on the diffusion path optimization parameters and the dynamic characteristic parameters of the interaction between the formation and the slurry, a diffusion and migration model of the composite slurry is constructed; S7. Perform multi-physics field coupling analysis on the three-dimensional mechanical model, the temperature field-seepage field-stress field coupling model and the diffusion migration model to obtain a numerical analysis model for grouting in frozen-thaw formations.

2. A freeze-thaw stratum grouting numerical analysis modeling method according to claim 1, characterized in that: The material parameters also include the composition ratio parameters of the composite slurry material and the microstructure parameters of the grouting stone body: Determine the composition ratio parameters of the composite slurry material based on the mass proportion of each component of the composite slurry material; The pore distribution characteristics of the grouting stone body after solidification were observed by scanning electron microscopy, and the porosity, average pore size and connectivity parameters were extracted as the microstructural parameters of the grouting stone body.

3. A freeze-thaw stratum grouting numerical analysis modeling method according to claim 2, characterized in that: The material parameters include the phase change energy storage performance parameters of the composite slurry material. Through the phase change material thermal cycle experiment, the phase change latent heat value and the phase change temperature range of the composite slurry material under different temperature gradients are obtained as the phase change energy storage performance parameters.

4. A freeze-thaw stratum grouting numerical analysis modeling method according to claim 3, characterized in that: According to the correlation between the phase change energy storage performance parameters and the damage parameters of the grouting stone body, the correction parameters of the elastic-plastic damage constitutive model are generated, including: The axial strain-stress curve and shear failure mode of the grouting stone body under different confining pressure conditions were obtained through triaxial shear test, and the shear strength, residual strength and damage softening coefficient were extracted as the damage parameters of the grouting stone body. Based on the phase change latent heat value and phase change temperature range, the nonlinear regression relationship between the phase change energy storage performance parameters and the damage softening coefficient is established, and the correction factor of the damage evolution rate changing with the phase change latent heat value is obtained by least squares fitting. According to the porosity and connectivity parameters in the microstructural parameters of the grouting stone body, the weight coefficient of the correction factor is set when the porosity is lower than the preset porosity threshold, and the correction parameters of the elastic-plastic damage constitutive model including the influence of the pore structure are generated.

5. A freeze-thaw stratum grouting numerical analysis modeling method according to claim 4, characterized in that: Based on the correction parameters and the mechanical property parameters of the grouting stone body, a three-dimensional mechanical model corresponding to the elastic-plastic damage constitutive model is constructed, including: Based on the porosity and average pore size in the microstructural parameters of the grouting stone body, the initial elastic modulus and Poisson's ratio of the grouting stone body are defined as mechanical characteristic parameters. The damage softening coefficient, correction factor and weight coefficient in the correction parameters are input into the elastic-plastic damage constitutive model, and the attenuation effect of the phase change latent heat value and the porosity on the elastic modulus is related through the damage evolution equation; The finite element method is used to discretize the elastic-plastic damage constitutive model into a three-dimensional mesh model. Based on the mechanical property parameters and the damage evolution equation, a three-dimensional mechanical model including the phase change energy storage-damage coupling effect is generated.

6. A freeze-thaw stratum grouting numerical analysis modeling method according to claim 5, characterized in that: Based on the phase change energy storage performance parameters, a coupling model of temperature field, seepage field and stress field of frozen-thaw formations is established and dynamic coupling data is obtained, including: Based on the phase change latent heat value and phase change temperature range, a freeze-thaw formation temperature field model is established to define the nonlinear relationship between the temperature gradient and the phase change latent heat release rate. Based on the porosity and connectivity parameters in the microstructural parameters of the grouting stone body, a seepage field model is constructed, the formation water conductivity is adjusted by the porosity, and the priority channel weight of the water migration path is corrected by the connectivity parameter. Based on the initial elastic modulus and Poisson's ratio of the mechanical characteristic parameters of the grouting stone body, combined with the temperature gradient change and seepage rate, the surrounding rock stress distribution is calculated through the temperature field-seepage field-stress field coupling equation; The temperature gradient change, seepage rate and surrounding rock stress distribution are output as dynamic coupling data.

7. A freeze-thaw stratum grouting numerical analysis modeling method according to claim 6, characterized in that: The influence of freeze-thaw interface fluctuations on the slurry diffusion path is analyzed based on the dynamic coupling data, and the diffusion path optimization parameters are generated, including: Based on the temperature gradient change data, the dynamic fluctuation range and fluctuation frequency of the freeze-thaw interface are extracted, and the interface migration trajectory of the ice-water phase transition in the temperature field is identified through the grayscale threshold segmentation algorithm; Combined with the seepage rate data and the initial elastic modulus in the mechanical characteristic parameters of the grouting stone body, the probability of preferential channel migration and the diffusion radius correction coefficient of the slurry diffusion path due to the freeze-thaw interface fluctuation were calculated. Based on the interface migration trajectory, migration probability and diffusion radius correction coefficient, the diffusion path optimization parameters including the priority channel weight and diffusion direction adjustment factor are generated.

8. A freeze-thaw stratum grouting numerical analysis modeling method according to claim 7, characterized in that: Based on the diffusion path optimization parameters and the dynamic characteristic parameters of the interaction between the formation and the slurry, a diffusion and migration model of the composite slurry is constructed, including: The priority channel weight and diffusion direction adjustment factor in the diffusion path optimization parameters are input into the porous media seepage equation to define the permeability tensor and flow direction deflection angle of the slurry in the priority channel; Based on the porosity and seepage rate data of the microstructural parameters of the grouting stone body and the time-varying characteristic equation of the slurry viscosity, the effective viscosity and flow resistance coefficient of the slurry in the heterogeneous pores are calculated. The finite volume method is used to couple the porous media seepage equation with the viscosity time-varying characteristic equation and discretize it to generate a composite slurry diffusion and migration model, and output the slurry pressure distribution, velocity vector and diffusion front evolution data.

9. A freeze-thaw stratum grouting numerical analysis modeling method according to claim 8, characterized in that: The three-dimensional mechanical model, the temperature field-seepage field-stress field coupling model and the diffusion migration model are subjected to multi-physics field coupling analysis to obtain the frozen-thaw stratum grouting numerical analysis model, including: The surrounding rock stress distribution data output by the three-dimensional mechanical model is converted into node stress field data that matches the temporal and spatial resolution of the temperature field-seepage field-stress field coupling model, and the stress-temperature feedback coefficient of the mechanical-thermal coupling interface is defined; Based on the seepage rate data in the dynamic coupling data and the slurry pressure distribution data output by the diffusion migration model, a bidirectional coupling equation of seepage field and diffusion field is established, and the seepage rate and slurry pressure are iteratively updated through the time step synchronization rule; The stress-temperature feedback coefficient and the seepage field-diffusion field bidirectional coupling equation are input into the multi-field coupling solver, and the separated iterative algorithm is used to simultaneously solve the three-dimensional mechanical model, the temperature field-seepage field-stress field coupling model and the diffusion migration model to generate a numerical analysis model for grouting in freeze-thaw formations.

10. A freeze-thaw stratum grouting numerical analysis modeling system, used to implement a freeze-thaw stratum grouting numerical analysis modeling method according to any one of claims 1 to 9, characterized in that: include: Phase change material parameter determination module: determine the material parameters of the phase change grouting composite material, including the phase change energy storage performance parameters of the composite slurry material; Constitutive correction parameter module: Generates correction parameters of the elastic-plastic damage constitutive model based on the correlation between phase change energy storage performance parameters and grouting stone body damage parameters; 3D mechanical modeling module: Based on the correction parameters and the mechanical characteristic parameters of the grouting stone body, a 3D mechanical model corresponding to the elastic-plastic damage constitutive model is constructed; Multi-field coupling modeling module: establish the coupling model of temperature field-seepage field-stress field of frozen-thaw formation based on phase change energy storage performance parameters and obtain dynamic coupling data; Diffusion path optimization module: Analyze the influence of freeze-thaw interface fluctuation on slurry diffusion path based on dynamic coupling data, and generate diffusion path optimization parameters; Diffusion and migration modeling module: Based on the diffusion path optimization parameters and the dynamic characteristic parameters of the interaction between the formation and the slurry, the diffusion and migration model of the composite slurry is constructed; Multi-field coupling analysis module: Multi-field coupling analysis is performed on the three-dimensional mechanical model, the temperature field-seepage field-stress field coupling model and the diffusion migration model to obtain a numerical analysis model for grouting in frozen-thaw formations.

Images