A numerical analysis modeling method and system for grouting in freeze-thaw strata

Abstract

The present invention discloses a numerical analysis and modeling method and system for grouting in freeze-thaw strata, specifically relating to the technical field of numerical simulation in water conservancy engineering, and is used to solve the problem of inaccurate prediction of grouting parameters caused by the existing grouting model ignoring the dynamic coupling relationship between the thermodynamic behavior of phase change materials and the seepage field and damage field; by establishing a dynamic correlation system between the phase change energy storage performance parameters and the seepage field-damage evolution, correcting the weight of the priority channel of the slurry diffusion path based on the freeze-thaw interface migration trajectory, coupling the two-way iterative equations of the seepage field and the diffusion migration field, and realizing the collaborative solution of the multi-physical field space-time scale through a split iterative algorithm; through the multi-dimensional coupling of the thermodynamic characteristics of phase change materials, the heterogeneous evolution of seepage paths and the mechanical response of damage mechanics, it significantly improves the prediction accuracy of the mechanical properties of grouting stones and the slurry diffusion range under freeze-thaw cycles, and effectively reduces the risks of leakage and structural instability.

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 grouting in frozen-thawed strata. Background Art

[0002] In hydraulic engineering applicable to cold regions and subjected to periodic freeze-thaw action (such as channels, reservoirs), the rock and soil mass undergoes dynamic expansion of seepage channels and cumulative structural damage due to temperature cycles. 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 phase change materials 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 phase change materials with the dynamic synergy relationship between the formation 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, embodiments of the present invention provide a numerical analysis and modeling method and system for grouting in frozen-thawed strata to solve the problems raised in the above background art.

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

[0006] A numerical analysis and modeling method for grouting in frozen-thawed strata, comprising the following steps:

[0007] S1. 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;

[0008] 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 grouted stone body;

[0009] 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;

[0010] S4. Establish a coupled model of the temperature field - seepage field - stress field in the frozen-thawed formation based on the phase change energy storage performance parameters and obtain dynamic coupling data;

[0011] 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;

[0012] S6. Construct a diffusion and migration model of the composite grout based on the optimized parameters of the diffusion path and the dynamic characteristic parameters of the interaction between the formation and the grout;

[0013] S7. Conduct a 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 frozen-thawed strata.

[0014] In a preferred embodiment, the material parameters further include the composition ratio parameters of the composite grout material and the microscopic structure parameters of the grouting stone body:

[0015] Determine the composition ratio parameters of the composite grout material based on the mass proportion of each component of the composite grout material;

[0016] Observe the pore distribution characteristics of the solidified grouting stone body through a scanning electron microscope, and extract the porosity, average pore diameter, and connectivity parameters as the microscopic structure parameters of the grouting stone body.

[0017] In a preferred embodiment, the material parameters include the phase change energy storage performance parameters of the composite grout 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 grout material under different temperature gradients as the phase change energy storage performance parameters.

[0018] 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:

[0019] 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;

[0020] 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;

[0021] 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 to generate the correction parameters of the elastoplastic damage constitutive model including the influence of pore structure.

[0022] In a preferred embodiment, based on the correction parameters and the mechanical characteristic parameters of the grouting stone body, construct a three-dimensional mechanical model corresponding to the elastoplastic damage constitutive model, including:

[0023] Based on the porosity and average pore diameter in the microscopic structure parameters of the grouting stone body, define the initial elastic modulus and Poisson's ratio of the grouting stone body as the mechanical characteristic parameters;

[0024] 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;

[0025] Discretize the elastoplastic damage constitutive model into a three-dimensional grid model by using the finite element method, and generate a three-dimensional mechanical model including the coupling effect of phase change energy storage - damage based on the mechanical property parameters and the damage evolution equation.

[0026] In a preferred embodiment, establish a coupled model of the temperature field - seepage field - stress field of the frozen-thawed formation based on the phase change energy storage performance parameters and obtain dynamic coupling data, including:

[0027] Establish a temperature field model of the frozen-thawed formation 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 release rate of the latent heat of phase change;

[0028] Construct a seepage field model based on the porosity and connectivity parameters in the microstructure parameters of the grouting stone body, adjust the formation water conductivity coefficient through the porosity, and correct the priority channel weight of the water migration path through the connectivity parameter;

[0029] 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, calculate the surrounding rock stress distribution through the temperature field - seepage field - stress field coupling equation;

[0030] Output the temperature gradient change, seepage rate, and surrounding rock stress distribution as dynamic coupling data.

[0031] In a preferred embodiment, 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, including:

[0032] Based on the temperature gradient change data, extract the dynamic fluctuation range and fluctuation frequency of the frozen-thawed interface, and identify the interface migration trajectory of ice-water phase change in the temperature field through the gray threshold segmentation algorithm;

[0033] 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 priority channel and the diffusion radius correction coefficient of the slurry diffusion path caused by the frozen-thawed interface fluctuation;

[0034] Generate diffusion path optimization parameters including the priority channel weight and diffusion direction adjustment factor based on the interface migration trajectory, migration probability, and diffusion radius correction coefficient.

[0035] In a preferred embodiment, construct a diffusion and migration model of the composite slurry based on the diffusion path optimization parameters and the dynamic characteristic parameters of the interaction between the formation and the slurry, including:

[0036] Input the priority 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 grout in the priority channel;

[0037] Based on the porosity and seepage rate data in the microstructure parameters of the grouting stone body, and combined with the time-varying characteristic equation of the grout viscosity, calculate the effective viscosity and flow resistance coefficient of the grout in the heterogeneous pores;

[0038] 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 composite grout diffusion and migration model, and output the grout pressure distribution, flow velocity vector and diffusion front evolution data.

[0039] In a preferred embodiment, perform 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 frozen-thawed strata, including:

[0040] 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 model, and define the stress-temperature feedback coefficient of the mechanical-thermal coupling interface;

[0041] Based on the seepage rate data in the dynamic coupling data and the grout 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 grout pressure through the time step synchronization rule;

[0042] 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 iterative algorithm to synchronously solve the three-dimensional mechanical model, the temperature field-seepage field-stress field coupling model and the diffusion and migration model to generate a numerical analysis model for grouting in frozen-thawed strata.

[0043] On the other hand, the present invention provides a numerical analysis modeling system for grouting in frozen-thawed strata, including:

[0044] 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 grout material;

[0045] 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;

[0046] 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;

[0047] Multi-field coupling modeling module: Based on the performance parameters of phase change energy storage, establish a coupled model of temperature field - seepage field - stress field for freeze-thaw strata and obtain dynamic coupling data;

[0048] Diffusion path optimization module: Analyze the influence of freeze-thaw interface fluctuations on the slurry diffusion path based on dynamic coupling data, and generate diffusion path optimization parameters;

[0049] 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, construct a diffusion and migration model of the composite slurry;

[0050] Multi-field coupling analysis module: Conduct multi-physical field coupling analysis on the three-dimensional mechanical model, temperature field - seepage field - stress field coupling model, and diffusion and migration model to obtain a numerical analysis model for freeze-thaw strata grouting.

[0051] Compared with the prior art, the present invention has the following beneficial effects:

[0052] 1. By establishing a dynamic correlation system between the performance parameters of phase change energy storage, seepage field, and damage evolution, it is possible to accurately depict the dynamic interaction of temperature gradient changes, pore structure deterioration, and slurry diffusion path during the freeze-thaw cycle. Through the three-dimensional mechanical reconstruction of the elastoplastic damage constitutive model and the real-time feedback of the influence of freeze-thaw interface fluctuations, the full-cycle simulation and prediction of the non-uniform mechanical properties of the grouting stone body and the slurry diffusion range are realized, significantly enhancing 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.

[0053] 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 diffusion path optimization parameters 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 and migration model and multi-field coupling data, improving the reliability of grouting process optimization in complex hydro-thermal-mechanical coupling environments, reducing leakage risks, and enhancing structural stability. BRIEF DESCRIPTION OF THE DRAWINGS

[0054] Figure 1 It is a flowchart of a numerical analysis modeling method for freeze-thaw strata grouting according to the present invention;

[0055] Figure 2 It is a structural schematic diagram of a numerical analysis modeling system for freeze-thaw strata grouting according to the present invention. DETAILED DESCRIPTION OF THE INVENTION

[0056] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with 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.

[0057] Embodiment 1: Figure 1 A numerical analysis and modeling method for grouting in freeze-thaw strata of the present invention is given, which includes the following steps:

[0058] S1. 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 grout material;

[0059] 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;

[0060] 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;

[0061] S4. Establish a coupled model of the temperature field - seepage field - stress field in the freeze-thaw strata based on the phase change energy storage performance parameters and obtain dynamic coupling data;

[0062] S5. Analyze the influence of the freeze-thaw interface fluctuation on the grout diffusion path according to the dynamic coupling data, and generate diffusion path optimization parameters;

[0063] S6. Based on the diffusion path optimization parameters and the dynamic characteristic parameters of the interaction between the stratum and the grout, construct a diffusion and migration model of the composite grout;

[0064] 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 and migration model to obtain a numerical analysis model for grouting in freeze-thaw strata.

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

[0066] Based on the mass ratio of each component of the composite grout material, determine the composition ratio parameters of the composite grout material;

[0067] By observing the pore distribution characteristics of the cured grouting stone body through a scanning electron microscope, extract the porosity, average pore diameter and connectivity parameters as the microscopic structure parameters of the grouting stone body.

[0068] S1. 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 grout material, including:

[0069] Through the thermal cycling experiment of the phase change material, the latent heat of phase change 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;

[0070] The phase change energy storage performance parameters of the composite slurry material are obtained through the thermal cycling experiment of the phase change material. Specifically, the differential scanning calorimetry is used to perform the heating and cooling cycle test on the composite slurry material. The temperature change range is set from -30°C to 25°C, and the heating and cooling rate is 5°C / min. The heat flow curves of the heat absorption and heat release processes of the material in each cycle are recorded. The latent heat of phase change value is calculated according to the integral areas of the endothermic peak and the exothermic peak in the heat flow curve, and the phase change start temperature, peak temperature and end temperature are extracted 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 from -10°C to 5°C, it adapts to the typical temperature fluctuation range of the freeze-thaw formation.

[0071] Based on the mass ratio of each component of the composite slurry material, the composition ratio parameters are determined. Specifically, the cement-based material, phase change microcapsules and water reducer are mixed in a preset ratio. The composition ratio parameters are obtained by weighing the masses of each component and calculating their percentages in the total mass. For example, when the cement-based material accounts for 60%, the phase change microcapsules account for 30%, and the water reducer accounts for 10%, this ratio parameter is recorded 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.

[0072] The pore distribution characteristics of the grouted stone body after curing are observed by a scanning electron microscope, and the porosity, average pore diameter and connectivity parameters are extracted as the microstructural parameters of the grouted stone body. Specifically, the cured grouted stone body sample is cut into thin slices of 5 mm×5 mm×2 mm, and after sputter coating, it is placed under the scanning electron microscope. Microstructural images are taken at magnification multiples from 500 to 2000 times. The pores in the image are binarized and segmented using image processing software. The porosity is obtained by statistically calculating the proportion of the pore area, the average value of the equivalent diameters of all pores is calculated as the average pore diameter, and the connectivity parameter is generated by judging the connectivity between pores through the pore skeleton extraction algorithm. For example, when the cement-based material accounts for 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.

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

[0074] 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:

[0075] Obtain the axial strain-stress curve and shear failure mode of the grouting stone body under different confining pressure conditions through triaxial shear tests, and extract the shear strength, residual strength and damage softening coefficient as the damage parameters of the grouting stone body;

[0076] Based on the latent heat value of phase change and the phase change temperature range, establish the 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 value of phase change by least squares fitting;

[0077] 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 pore structure.

[0078] 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, apply three groups of confining pressures of 0.5 MPa, 1.0 MPa, and 1.5 MPa to the grouting stone body specimen using a rock mechanics testing machine, 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, calculate the ratio of the peak strength to the residual strength as the damage softening coefficient, and 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.

[0079] Based on the latent heat of phase change value and the phase change temperature range in step S1, a non-linear regression relationship between the phase change energy storage performance parameters and the damage softening coefficient is established. Specifically, the test data of the damage softening coefficient of the grouting stone body corresponding to different latent heat of phase change values are fitted, and the least squares method is used 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. The correction factor of the damage evolution rate with the change of the latent heat of phase change value is determined 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.

[0080] According to the porosity and connectivity parameters in the microscopic structure parameters of the grouting stone body in step S1, the weight coefficient of the correction factor when the porosity is lower than the preset porosity threshold is set. Specifically, by statistically analyzing the correlation data between the damage softening coefficient and the correction factor of specimens with different porosities, the porosity threshold range is determined. When the porosity is lower than this threshold, the weight coefficient is adjusted according to the negative correlation between the connectivity parameter and the correction factor to generate the correction parameters of the elastoplastic damage constitutive model including the influence of the 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 the grouting stone body specimens. For example, by performing a normal distribution analysis on the porosity data of 30 groups of specimens, the mean value minus twice the standard deviation is taken as the threshold lower limit.

[0081] 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, it corresponds to the in-situ stress level of the typical burial depth of 5 - 15 meters in the cold region canal lining project. 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 parameters, while in this embodiment, the dynamic influence of the heat storage and release scale on the damage evolution is quantified by introducing the latent heat of phase change value. In the setting of the weight coefficient, the joint calibration method of the porosity and connectivity parameters corrects the deviation that only depends on the porosity 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%.

[0082] 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:

[0083] Based on the porosity and average pore size 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;

[0084] The damage softening coefficient, correction factor, and weight coefficient in the correction parameters are input into the elastoplastic damage constitutive model, and the attenuation effect of the latent heat of phase change value and porosity on the elastic modulus is correlated through the damage evolution equation;

[0085] The finite element method is used to discretize the elastoplastic damage constitutive model into a three-dimensional grid model, and 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.

[0086] When constructing the 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 size 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, the empirical formula of 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 size 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 combining the Poisson's ratio recommended value table in the material handbook. The equivalent sphere model is a geometric simplification method that equates the pores to spherical cavities, and the pore shape factor is calculated by the ratio of the average pore size to the sphere diameter.

[0087] The damage softening coefficient, correction factor, and weight coefficient in the correction parameters are input into the elastoplastic damage constitutive model, and the attenuation effect of the latent heat of phase change value and porosity on the elastic modulus is correlated through the damage evolution equation. Specifically, the damage evolution equation takes the damage softening coefficient as the reference elastic modulus attenuation rate, the correction factor as the nonlinear 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 of the elastic modulus attenuation rate and the correction parameters of specimens with different latent heat of phase change values under the same confining pressure needs to be greater than 0.9.

[0088] The elastoplastic damage constitutive model is discretized into a three-dimensional grid model by 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 hexahedral element meshes are divided with an element size ranging from 5 mm to 10 mm. The grid 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. The initial elastic modulus, Poisson's ratio, and parameters of the damage evolution equation are input into the material property module of the finite element software, and the nonlinear equations are solved by the Newton-Raphson iterative algorithm to generate a three-dimensional mechanical model that includes the coupling effect of phase change energy storage - damage. The division of hexahedral element meshes and the Newton-Raphson iterative method are conventional techniques in finite element analysis. In this embodiment, the heterogeneous characteristics of the grouting stone body in the freeze-thaw stratum are adapted by limiting the element size range and the type of input parameters.

[0089] The empirical relationship between porosity and elastic modulus is fitted based on laboratory compression test data. For example, uniaxial compression tests are performed on 20 groups of specimens with different porosities, and an exponential function relationship between the 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 the recommended Poisson's ratio value table in the material handbook. The division of hexahedral element meshes and iterative calculations in the finite element method are conventional numerical methods. In this embodiment, the heterogeneous mechanical characteristics of the grouting stone body in the freeze-thaw stratum are adapted by limiting the element size range (5 mm to 10 mm) and the type of input parameters (parameters of the damage evolution equation).

[0090] 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, the elastic modulus decay rate and the damage softening coefficient show a benchmark linear relationship; when the correction factor is increased to 1.2 due to the latent heat of phase change value, and the weight coefficient is decreased 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 benchmark value, reflecting the dynamic coupling effect of phase change energy storage and pore structure. The element size range is determined according to the characteristic size of the grouting stone body (such as the aggregate particle size of 1 - 3 mm) to ensure a balance between grid accuracy and calculation efficiency.

[0091] S4. Establish a coupled model of the temperature field - seepage field - stress field in the freeze-thaw stratum based on the phase change energy storage performance parameters and obtain dynamic coupling data, including:

[0092] Establish a temperature field model of the freeze-thaw stratum based on the latent heat of phase change value and the phase change temperature range, and define the nonlinear variation relationship of the temperature gradient with the latent heat release rate of phase change;

[0093] 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 through the porosity, and the priority channel weight of the moisture migration path is corrected through the connectivity parameter;

[0094] 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 surrounding rock stress distribution is calculated through the temperature field-seepage field-stress field coupling equation;

[0095] The temperature gradient change, seepage rate, and surrounding rock stress distribution are output as dynamic coupling data.

[0096] When establishing the temperature field-seepage field-stress field coupling action model of the frozen-thawed formation and obtaining the dynamic coupling data, it is realized through the following steps: Based on the latent heat of phase change value and phase change temperature range in step S1, a temperature field model of the frozen-thawed formation is established. Specifically, the latent heat of phase change value is input as the heat source term into the heat conduction equation, and the latent heat release rate is set 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 of temperature. The growth coefficient is determined by fitting the test data of differential scanning calorimetry.

[0097] Based on the porosity and connectivity parameters in the microstructural parameters of the grouting stone body in step S1, a seepage field model is constructed. Specifically, the formation water conductivity is calculated through the porosity, and the calculation formula is the modified Darcy's law, where the water conductivity is proportional to the square of the porosity. For example, when the porosity is 10%, the water conductivity is 1×10 -5 m / s, and for every 5% increase in porosity, the water conductivity is increased to 1.25×10 -5 m / s; the priority channel weight of the seepage path is corrected through the connectivity parameter, and the weight coefficient calculation formula 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.

[0098] 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, the stress field distribution is calculated. Specifically, a temperature field-seepage field-stress field coupling equation is established, and the thermal expansion strain caused by the temperature gradient, the change in pore water pressure caused by the seepage rate are related to the initial elastic modulus and Poisson's ratio, and the coupling equation is solved by the finite element method. For example, when the initial elastic modulus is 8.5 GPa, the corresponding thermal expansion coefficient is 1.2×10 -5 / ℃, the thermal stress generated per 1℃ / m increase in the temperature gradient is 8.5 GPa × 1.2 × 10 -5 × 1 = 0.102 MPa; the pore water pressure drops by 5 kPa corresponding to a seepage rate of 0.1 m / s, and the net stress change of the surrounding rock is calculated by the effective stress principle.

[0099] Output the temperature gradient change data, seepage rate data, and surrounding rock stress distribution data as dynamically coupled data, and the output format is three-dimensional field data in time series. For example, the temperature gradient data is updated to the grid node temperature value matrix at each hourly time step, the seepage rate data is marked with the flow velocity direction and magnitude of each unit in the 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 synchronously updated every hour, and the stress field is iteratively calculated every 10 minutes to match the migration rate of the freeze-thaw interface.

[0100] S5. Analyze the influence of the freeze-thaw interface fluctuation on the slurry diffusion path based on the dynamically coupled data, and generate the diffusion path optimization parameters, including:

[0101] Based on the temperature gradient change data, extract the dynamic fluctuation range and fluctuation 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;

[0102] Combined with the seepage rate data and the initial elastic modulus in the mechanical property parameters of the grouting stone body, 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;

[0103] Based on the interface migration trajectory, migration probability, and diffusion radius correction coefficient, generate the diffusion path optimization parameters including the preferential channel weight and diffusion direction adjustment factor.

[0104] When analyzing the influence of the freeze-thaw interface fluctuation on the slurry diffusion path based on the dynamically coupled data, it is realized through the following steps: Extract the dynamic fluctuation range and fluctuation 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℃ / m increase in the temperature gradient. Identify the continuous area 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.

[0105] The gray threshold segmentation algorithm uses the adaptiveThreshold function of the OpenCV library, the block size is 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℃ / m, the corresponding gray value is 50, and when it is 5℃ / m, the corresponding gray value is 200. After segmentation, the interface trajectory is stored in the form of a polygon contour line.

[0106] Calculate the migration probability of the preferential channel of the slurry diffusion path and the diffusion radius correction coefficient 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 calculation formula for the migration probability 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 calculation formula for the diffusion radius correction coefficient is 1 / (1 + migration probability). For example, when the migration probability is 0.085, the corresponding correction coefficient is 1 / 1.085≈0.922, indicating that the diffusion radius is retained at 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).

[0107] 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, and calculate the spatial intersection area between the two. The weight of the preferential channel within 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 tangent direction of the interface migration trajectory. 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 cell.

[0108] 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:

[0109] 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 the flow direction deflection angle of the slurry in the preferential channel;

[0110] Based on the porosity in the microscopic structure parameters of the grouting stone body and the seepage rate data, and 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;

[0111] 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.

[0112] 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%.

[0113] 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.

[0114] 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.

[0115] 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:

[0116] 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 model, and define the stress-temperature feedback coefficient of the mechanical-thermal coupling interface;

[0117] Based on the seepage rate data in the dynamic coupling data 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;

[0118] 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.

[0119] When performing 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, 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 matching 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.

[0120] 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 of 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, set the seepage rate to be updated every 1 hour in the seepage field and the slurry pressure to be updated every 10 minutes in the diffusion field, and the two synchronize data exchange 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%.

[0121] 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 action 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 is performed 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 adopts the Gauss-Seidel iterative format, and the relaxation factor is set to 0.8 to balance the convergence rate and stability.

[0122] 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:

[0123] 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;

[0124] 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;

[0125] Three-dimensional mechanical modeling module: Based on the modification parameters and the mechanical characteristic parameters of the grouting stone body, construct a three-dimensional mechanical model corresponding to the elastoplastic damage constitutive model;

[0126] Multi-field coupling modeling module: Establish a coupling action model of the temperature field-seepage field-stress field in the freeze-thaw stratum based on the phase change energy storage performance parameters and obtain dynamic coupling data;

[0127] 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;

[0128] 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;

[0129] 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 migration model to obtain a numerical analysis model for freeze-thaw stratum grouting.

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

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

[0132] In addition, the functional modules in the various embodiments of the present 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.

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

[0134] 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 should all be included in the protection scope of the present invention.

Claims

1. A numerical analysis modeling method for grouting in freeze-thaw strata, characterized in that, It 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 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 triaxial shear tests, 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 value of phase change 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 value of phase change by 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 pore structure; 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. Establish a coupled model of temperature field-seepage field-stress field for the freeze-thaw formation based on the phase change energy storage performance parameters and obtain dynamic coupling data, including: Establish a temperature field model for the freeze-thaw formation based on the latent heat value of phase change 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 grouting 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 grouting 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; 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 formation and the slurry, construct a diffusion and migration model of 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 freeze-thaw formation.

2. A numerical analysis and modeling method for grouting in freeze-thaw strata according to claim 1, characterized in that, The material parameters also include the composition ratio parameters of the composite slurry material and the microscopic structure parameters of the grouting stone body: Based on the mass ratio of each component of the composite slurry material, determine the composition ratio parameters of the composite slurry material; 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 microscopic structure parameters of the grouting stone body.

3. A numerical analysis modeling method for grouting in freeze-thaw strata 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, obtain the latent heat value of phase change and the phase change temperature range of the composite slurry material under different temperature gradients as the phase change energy storage performance parameters.

4. A numerical analysis and modeling method for grouting in freeze-thaw strata according to claim 3, 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 elastoplastic damage constitutive model is constructed, including: 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; The damage softening coefficient, correction factor and weight coefficient in the correction parameters are input into the elastoplastic damage constitutive model, and the attenuation effect of the latent heat value of phase change and porosity on the elastic modulus is correlated 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 phase change energy storage-damage coupling effect is generated based on the mechanical property parameters and the damage evolution equation.

5. A numerical analysis and modeling method for grouting in freeze-thaw strata according to claim 4, characterized in that, According to the dynamic coupling data analysis of the influence of the freeze-thaw interface fluctuation on the slurry diffusion path, 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 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 preferential channel and the diffusion radius correction coefficient of the slurry diffusion path affected by the freeze-thaw interface fluctuation are calculated; Based on the interface migration trajectory, migration probability and diffusion radius correction coefficient, diffusion path optimization parameters including the preferential channel weight and diffusion direction adjustment factor are generated.

6. A numerical analysis modeling method for grouting in freeze-thaw strata according to claim 5, 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 preferential 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 preferential channel are defined; 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, the effective viscosity and flow resistance coefficient of the slurry in the heterogeneous pores are calculated; The porous medium seepage equation and the time-varying characteristic equation of viscosity are coupled and discretized by the finite volume method to generate a diffusion and migration model of the composite slurry, and the slurry pressure distribution, velocity vector and diffusion front evolution data are output.

7. A numerical analysis and modeling method for grouting in freeze-thaw strata according to claim 6, characterized in that, The three-dimensional mechanical model, the temperature field-seepage field-stress field coupling action model and the diffusion and migration model are subjected to multi-physical field coupling analysis to obtain a numerical analysis model for grouting in frozen strata, including: The surrounding rock stress distribution data output by the three-dimensional mechanical model is converted into node stress field data matching the spatio-temporal resolution of the temperature field-seepage field-stress field coupling action 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 and migration model, a two-way coupling equation of the seepage field-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 two-way coupling equation of the seepage field-diffusion field are input into the multi-field coupling solver, and the three-dimensional mechanical model, the temperature field-seepage field-stress field coupling action model and the diffusion and migration model are synchronously solved by the split iteration algorithm to generate a numerical analysis model for grouting in frozen strata.

8. A numerical analysis and modeling system for freeze-thaw stratum grouting, which is used to implement the numerical analysis and modeling method for freeze-thaw stratum grouting according to any one of claims 1-7, characterized in that, Including: Phase change material parameter determination module: 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; 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 mechanical modeling module: Based on the modification parameters and the mechanical characteristic 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 coupled model of the temperature field - seepage field - stress field in the freeze-thaw formation based on the phase change energy storage performance parameters 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 formation 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 coupled model and the diffusion migration model to obtain a numerical analysis model for grouting in the freeze-thaw formation.

Images