A discrete element simulation method of freeze-thaw process considering latent heat and dynamic evolution of thermodynamic properties of water ice

Abstract

The application discloses a discrete element freezing-thawing simulation method considering water ice phase change latent heat and dynamic evolution of thermodynamic properties, and belongs to the cross field of rock mass engineering mechanics and computer numerical simulation technology. The method comprises the following steps: constructing a rock discrete element model, and filling water particles in the fissure; in the phase change temperature interval, the thermodynamic parameters of the water particles are dynamically updated based on the equivalent heat capacity model to simulate the phase change latent heat effect, the mechanical parameter step update and the radius increment update of the water particles are triggered synchronously, and the radius increment update is decoupled from the thermal expansion and contraction calculation in the non-phase change period; and the stiffness attenuation calculation is performed on the ice particles to simulate the irreversible mechanical damage in the limited space. The application solves the problems that the dynamic change of physical and mechanical parameters, the latent heat effect and the ice body pressure failure are ignored in the existing discrete element freezing-thawing simulation, realizes the complete coupling of the phase change-stress field-damage, and significantly improves the accuracy of the simulation of the frost heaving force and the crack propagation evolution process in the rock freezing-thawing process.

Description

Technical Field

[0001] This invention belongs to the interdisciplinary field of rock mass engineering mechanics and computer numerical simulation technology, and specifically relates to a discrete element freeze-thaw simulation method that takes into account both the latent heat of water-ice phase change and the dynamic evolution of thermodynamic properties. Background Technology

[0002] In cold-region rock engineering, repeated freeze-thaw phase changes caused by water trapped in rock fissures generate frost heave stress, leading to internal rock damage and continuous strength deterioration. This is a core contributing factor to rock cracking and structural instability. Research on the freeze-thaw damage mechanism of rocks has significant engineering value for stability control in cold-region rock engineering. Particle flow discrete element method (PFC) can effectively simulate the mechanical behavior of discontinuous media and the dynamic crack propagation process. With its unique micromechanical advantages, it has become the mainstream method for numerical research on rock freeze-thaw damage.

[0003] Currently, the mainstream PFC freeze-thaw simulation method in this field simulates the freeze-thaw effect by driving particle volume deformation through the coefficient of thermal expansion. By statically assigning values ​​to particle thermodynamic parameters, it achieves a simplified simulation of the rock freeze-thaw process and can initially reproduce the macroscopic phenomenon of rock mass freeze-thaw cracking. However, the existing simulation system simplifies the multi-physics coupling mechanism of water-ice phase transition, resulting in several key technical defects and causing the simulation results to deviate significantly from the actual physical process.

[0004] First, existing methods neglect the abrupt change in stiffness during the water-ice phase transition. The core mechanism of actual frost heave is that the stiffness of particles increases significantly after the water phase transition, leading to a surge in contact forces caused by particle overlap within a confined space, which in turn compresses the rock mass and causes damage, rather than simply volume expansion. Traditional models rely solely on the coefficient of thermal expansion to represent particle volume changes, lacking the stiffness abrupt change mechanism, thus distorting the microscopic nature of frost heave damage and causing inaccurate mechanism simulations.

[0005] Secondly, the latent heat effect of phase transition is not coupled. The water-ice phase transition process exhibits isothermal latent heat release / absorption characteristics, which are crucial for controlling the evolution of the freeze-thaw temperature field and the advancement of the freeze-thaw front. Existing simulations ignore the effect of latent heat and directly skip the phase transition temperature plateau stage, which easily leads to instantaneous spurious thermal deformation and thermal stress, resulting in a distortion of the evolution sequence of frost heave force.

[0006] Furthermore, the static assignment of thermodynamic parameters does not conform to actual working conditions. The thermal conductivity, specific heat capacity, and other thermodynamic parameters of water and ice differ significantly, while the parameters of the traditional model remain fixed and cannot be dynamically updated with the phase transition of particulate water and ice. This disrupts the dynamic feedback mechanism between the temperature field and the phase transition field, further exacerbating the quantitative deviation of the temperature field and frost heave force distribution.

[0007] Finally, the mechanism of ice particle stress-induced damage degradation is missing. Ice bodies within confined fracture spaces undergo irreversible stiffness decay under high-pressure compression. However, existing models simplify ice bodies as ideal elastic or elastoplastic bodies, lacking damage evolution characteristics. This leads to spurious stiffness recovery during the frost heave force decline phase, resulting in severe distortion in the calculation of cumulative rock mass damage under multiple freeze-thaw cycles.

[0008] In summary, existing PFC freeze-thaw simulation technology staticizes and simplifies the dynamically coupled phase change freeze-thaw process, failing to accurately reproduce the microscopic mechanism and multi-physics field evolution law of water-ice phase change in rock fractures in cold regions. The numerical simulation accuracy is insufficient and cannot meet the needs of accurate long-term stability analysis of rock mass engineering in cold regions. Summary of the Invention

[0009] This invention addresses the problem in existing discrete element freeze-thaw simulations of PFC particles that fail to comprehensively consider abrupt changes in particle stiffness, release of latent heat of phase change, dynamic coupling evolution of particle thermodynamic properties, and stress-induced damage mechanisms during the water-ice phase transition, leading to significant deviations between simulated frost heave curves, temperature fields, and crack propagation modes and real-world conditions. It proposes a discrete element freeze-thaw simulation method that takes into account both the latent heat of water-ice phase change and the dynamic evolution of thermodynamic properties.

[0010] To achieve the above-mentioned technical improvements, this invention provides a discrete element freeze-thaw simulation method that takes into account both the latent heat of water-ice phase change and the dynamic evolution of thermodynamic properties, comprising:

[0011] A discrete element model of rock containing prefabricated fracture structures was constructed, and water particles were filled into the fracture structures.

[0012] Within the phase change temperature range, the thermodynamic parameters of water particles are dynamically updated based on the equivalent heat capacity model to simulate the latent heat effect of phase change.

[0013] Within the phase change temperature range, the mechanical parameters of water particles are simultaneously updated by step and the radius is updated by increment. The radius update is decoupled from the thermal expansion and contraction calculation during the non-phase change period.

[0014] Stiffness decay calculations are performed on frozen water particles based on strain accumulation to simulate irreversible mechanical damage in a confined space.

[0015] As one implementation method, the construction of a rock discrete element model containing a prefabricated fracture structure and the filling of water particles within the fracture structure includes: generating rock particles within a defined model area and establishing interparticle bonding; deleting rock particles at specific locations to generate the prefabricated fracture structure; and filling the fracture structure space with water particles based on the geometric volume and porosity of the prefabricated fracture and establishing contact with the rock particles on the fracture wall.

[0016] As one implementation method, after generating rock particles and water particles, the method further includes: obtaining the macroscopic thermal conductivity of rock, water and ice; converting the macroscopic thermal conductivity into microscopic particle contact thermal resistance based on the relationship between macroscopic thermal conductivity and microscopic particle contact thermal resistance; and assigning the microscopic particle contact thermal resistance to the corresponding rock particles and water particles for heat conduction calculation.

[0017] As one implementation method, the method of dynamically updating the thermodynamic parameters of water particles based on the equivalent heat capacity model to simulate the latent heat effect of phase change includes: treating water particles as independent thermodynamic units, and dynamically adjusting the specific heat capacity and thermal conductivity of water particles according to the real-time temperature and phase change state; simulating the absorption or release of latent heat during the phase change process through the dynamic adjustment of specific heat capacity and thermal conductivity, so as to generate a temperature plateau period in the temperature field evolution.

[0018] As one implementation method, the phase transition temperature range is defined as [ , The specific heat capacity of water particles is updated according to the following formula:

[0019] ;

[0020] In the formula, and These are the specific heat capacities of water and ice, respectively. To account for the equivalent volumetric heat capacity of the latent heat of phase change, The latent heat per unit mass of water; The density of water, The real-time temperature of the particles; for thermal conductivity, when the temperature of the water particles decreases to... At that time, its thermal conductivity was switched from that of water phase to that of ice phase.

[0021] As one implementation method, the step update of the mechanical parameters and the radius increment update of water particles are triggered simultaneously within the phase change temperature range, and the radius increment update is decoupled from the thermal expansion and contraction calculation during the non-phase change period. This includes: when water particles enter the phase change temperature range, the water particles are transformed into ice particles, causing a step increase in the stiffness of the ice particles, while a fixed increment update is assigned to the radius of the ice particles; outside the phase change temperature range, the radius change caused by temperature change is calculated, and the radius change is linearly superimposed with the fixed increment of the radius to achieve physical separation of phase change expansion and thermal expansion and contraction.

[0022] As one implementation method, the stiffness update formula is as follows:

[0023] ;

[0024] In the formula, The initial simulated normal stiffness when water particles freeze into ice particles; The simulated normal stiffness of ice particles; This is the stiffness amplification factor, and ≥4.

[0025] As one implementation method, when the particle is outside the phase transition temperature range, the particle radius drives the volume change based on the coefficient of thermal expansion; when the particle is inside the phase transition temperature range, the particle radius is directly assigned a fixed radius increment proportional to the radius when entering the phase transition temperature range to simulate the volume expansion caused by the water-ice phase transition; when the particle is outside the phase transition temperature range, the radius change caused by the temperature change is calculated and linearly superimposed with the fixed radius increment within the phase transition temperature range to obtain the current radius of the particle.

[0026] As one implementation method, the stiffness attenuation calculation for frozen water particles based on strain accumulation includes: obtaining the real-time strain generated by the compression of frozen water particles in a confined space; calculating the effective strain accumulation that damages stiffness when the real-time strain exceeds the strain threshold that does not damage stiffness; and calculating the real-time normal stiffness of the ice particles based on the effective strain accumulation using a step-type stiffness descent function, so that the real-time normal stiffness decreases monotonically with the increase of the effective strain accumulation.

[0027] As one implementation method, the step stiffness reduction function is as follows:

[0028] ;

[0029] In the formula, The simulated normal stiffness of ice particles. is the minimum stiffness of the ice particles; h is the stiffness reduction factor; For strain accumulation that damages stiffness, where, To respond in real time, Let be the critical strain value, which satisfies: ,in, This is the critical stress for frost heave damage.

[0030] Compared with the prior art, the present invention has the following advantages:

[0031] 1. More realistic simulation of physical mechanisms: The equivalent volumetric heat capacity method is introduced to accurately realize the temperature plateau period caused by the release of latent heat of phase change, and correct the deviation between the timing of frost heave force and the real physical process.

[0032] 2. Eliminate physical disconnect of parameters: Decouple the volume expansion caused by phase change from the rate of temperature change, correct the error in traditional simulation that the expansion only continues during the cooling process, ensure the synchronization of stiffness change and volume expansion, and improve the accuracy of simulating the principle of fracture water freeze-thaw.

[0033] 3. An irreversible damage mechanism was constructed: a stiffness reduction function based on strain accumulation was introduced to effectively simulate the irreversible mechanical damage of ice under constrained compression, avoiding the problem of false recovery of ice stiffness after the frost heave force subsides.

[0034] 4. Constructing a multi-physics coupling channel: Real-time dynamic synchronous coupling of temperature field, phase transition field, and mechanical field was realized, providing a high-precision numerical simulation tool for the study of freeze-thaw damage mechanism of rock fractures in cold regions and the evaluation of engineering stability. Attached Figure Description

[0035] Figure 1 This is a schematic flowchart of a discrete element freeze-thaw simulation method that takes into account both the latent heat of water-ice phase change and the dynamic evolution of thermodynamic properties, according to Embodiment 1 of the present invention. Detailed Implementation

[0036] This embodiment provides a discrete element freeze-thaw simulation method that takes into account both the latent heat of water ice phase change and the dynamic evolution of thermodynamic properties, including the following steps:

[0037] S1. Construct a discrete element model of rock containing prefabricated fracture structures and fill the fracture structures with water particles.

[0038] Rock samples are taken from the mining area to be tested, rock mass specimens are made, and the thermodynamic and mechanical parameters of the rock are determined. The thermodynamic parameters include thermal conductivity and specific heat capacity, while the mechanical parameters include elastic modulus and the normal stiffness, tangential stiffness, and bond strength of the rock structural surfaces.

[0039] A discrete element model of rock is constructed. The model area is set according to the geometric outline of the rock specimen. Rock particles are generated in the model area according to the rock material type and particle distribution of the mining area to be tested. The bonding between particles is established by using a parallel bonding model. The location and conditions of the fractures in the field are obtained. Rock particles are deleted at the corresponding positions in the rock discrete element model to generate fracture structures. Water particles are filled into the fracture structure space and establish contact with the rock particles on the fracture wall.

[0040] Based on the measured thermodynamic and mechanical parameters of the rock, combined with the thermodynamic and mechanical parameters of water and ice, parameters are set for the discrete element model of the rock. The set parameters include at least the initial thermal resistance values ​​of different material particles and the initial radius of water particles. Based on the relationship between macroscopic thermal conductivity and microscopic particle contact thermal resistance, the thermal conductivity is converted into microscopic particle contact thermal resistance. The microscopic particle contact thermal resistance is assigned to the corresponding rock and water particles for heat conduction calculations to solve the mapping problem from macroscopic thermal properties to microscopic contact attributes. The relationship between macroscopic thermal conductivity and microscopic particle contact thermal resistance is as follows:

[0041] ;

[0042] In the formula, For particle contact thermal resistance, and Let be the thermal conductivity of the two contacting particles. The contact area is defined as the area between the two components. By transforming parameters from macroscopic to microscopic, it is possible to ensure that the heat conduction behavior of the discrete element model is consistent with the actual working conditions, eliminate the thermal field calculation deviation caused by the scale effect, and avoid the temperature conduction time sequence deviating from the actual process due to directly assigning empirical thermal resistance values.

[0043] S2. Within the phase change temperature range, the thermodynamic parameters of water particles are dynamically updated based on the equivalent heat capacity model to simulate the latent heat effect of phase change.

[0044] Treating water particles as independent thermodynamic units, within the phase transition temperature range, the specific heat capacity and thermal conductivity of the water particles are dynamically adjusted based on real-time temperature and phase transition state. Through this dynamic adjustment of specific heat capacity and thermal conductivity, the absorption or release of latent heat during the phase transition process is simulated, thereby generating a temperature plateau in the temperature field evolution. Specifically:

[0045] Boundary conditions with temperature variations over time are applied to the particles at the boundary of the discrete element model. Using a heat conduction module, the temperature of each particle is updated in real-time at each computational step. Within each computational step, each water particle is traversed, and its thermodynamic properties are independently updated based on its current temperature and phase state to accurately simulate the energy transfer and absorption during the water-ice phase transition. This generates a temperature plateau period during the phase transition, representing the absorption and release of latent heat. The thermodynamic properties include specific heat capacity and thermal conductivity.

[0046] Define the phase transition temperature range as [ , Within the phase transition temperature range, for specific heat capacity, according to the equivalent heat capacity calculation method, the general formula for calculating latent heat capacity is:

[0047] ;

[0048] In the formula, and These are the specific heat capacities of water and ice, respectively. To account for the equivalent volumetric heat capacity of the latent heat of phase change, The latent heat per unit mass of water; The density of water, This represents the real-time temperature of the particles. The liquid phase ratio is defined as the volume fraction of unfrozen water in a single water particle to make this macroscopic general formula applicable to microscopic water particles. Replace the liquid phase fraction in the general formula for calculating latent heat capacity. ;in, The definition is as follows:

[0049] ;

[0050] Will Differentiating with respect to T, we get Substituting into the general formula for calculating latent heat capacity, the updated formula for the specific heat capacity of water particles within the phase change temperature range is as follows:

[0051] ;

[0052] The significance of this formula lies in the fact that when the temperature of a water particle enters the phase transition temperature range, its specific heat capacity is set to increase significantly to characterize the inhibitory effect of latent heat release on temperature change and to achieve a temperature plateau during the phase transition process. Outside the phase transition temperature range, when the water particle temperature is above 273 K, the specific heat capacity of the aqueous phase is used; when the temperature is below 272.8 K, the specific heat capacity of the ice phase is used. The specific heat capacity of the water particle is then updated according to the following formula:

[0053] .

[0054] Regarding thermal conductivity, when the temperature of the water particles is greater than... When the water particle temperature drops to 273 K, the thermal conductivity of the water phase is used, which is 0.58 W / (m·K); when the water particle temperature drops to 273 K, the thermal conductivity is switched to that of the ice phase, which is 2.22 W / (m·K). This phase-state-based step switching ensures that the heat conduction calculation strictly corresponds to the actual physical process, thus ensuring the accurate simulation of the freeze-thaw front advance velocity.

[0055] Accurate simulation of the temperature plateau period is crucial for determining the advance rate of the freeze-thaw front. If this effect is ignored, the simulated cooling rate will be much faster than the actual rate, causing the subsequent frost heave force to be generated significantly earlier in time, which in turn affects the judgment of the rock mass damage evolution law.

[0056] S3. Within the phase transition temperature range, the mechanical parameters of water particles are simultaneously updated by step and the radius increment is updated, and the radius increment update is decoupled from the thermal expansion and contraction calculation during the non-phase transition period.

[0057] When water particles enter the phase transition temperature range, the particle phase transition conversion module is invoked to transform the water particles into ice particles. This simultaneously triggers a step increase in mechanical parameters and a fixed increment update in radius, simulating the sudden change in stiffness and volume expansion during water freezing. This achieves complete coupling of the "temperature field-phase transition field-stress field" during the phase transition process. Outside the phase transition temperature range, an independent thermal expansion coefficient module calculates the radius change caused by temperature variations and linearly superimposes this radius change with the fixed radius increment to achieve physical separation of phase transition expansion and thermal expansion and contraction. Specifically:

[0058] Regarding particle stiffness, when water particles freeze into ice particles, their mechanical parameters, including stiffness, increase dramatically. The update formula is as follows:

[0059] ;

[0060] In the formula, The initial simulated normal stiffness (in N / m) when water particles freeze into ice particles. The simulated normal stiffness of the ice particles (in N / m); This is the stiffness amplification factor. The range of values ​​is determined based on the stiffness ratio of ice particles to water particles. ≥4, the specific value is determined by optimization based on simulation results.

[0061] For the dynamic updating of particle radius, in each calculation step, the particle phase transition module iterates through each water particle, quantifies the phase state of the water or ice particle based on the particle's real-time temperature, and iterates through the current real-time radius of each water particle. When the particle is outside the phase transition temperature range, the particle radius depends on the coefficient of thermal expansion driving the volume change; when the particle is within the phase transition temperature range, the particle radius is directly assigned a fixed radius increment proportional to the radius when entering the phase transition temperature range, to simulate the volume expansion caused by the water-ice phase transition; when the temperature is higher than the upper limit of the phase transition temperature range, the radius change caused by the temperature change is calculated to obtain the particle's current radius; when the temperature is lower than the lower limit of the phase transition temperature range, the radius change caused by the temperature change is linearly superimposed with the fixed radius increment within the phase transition temperature range to obtain the particle's current radius. The radius update formula within the phase transition temperature range is as follows:

[0062] ;

[0063] In the formula, This is the real-time radius (in meters) of the water particles when they enter the phase change temperature range, i.e., the radius when the particle temperature just drops to 273K. This is the preset growth rate coefficient.

[0064] In traditional discrete element freeze-thaw simulations, the coefficient of thermal expansion is typically used to drive particle volume changes to simulate freeze-thaw expansion. However, when the latent heat of phase change is introduced, the rate of temperature change of water particles during the phase change approaches zero. If the traditional thermal expansion formula, which relies on the rate of temperature change, is still used, the volume expansion of the particles will also approach zero, which clearly contradicts the physical fact that water volume inevitably increases dramatically when it freezes. This embodiment decouples the phase change expansion mechanism from the thermal expansion and contraction calculations based on the coefficient of thermal expansion during the non-phase change period, eliminating the temporal disorder of volume expansion caused by confusing the two in traditional methods. Simultaneously, this embodiment achieves synchronous coupling between the abrupt change in particle stiffness during the phase change and volume expansion, simulating the microscopic nature of the sudden change in stiffness and volume expansion during the instant of water freezing.

[0065] S4. Perform stiffness attenuation calculations on frozen water particles based on strain accumulation to simulate irreversible mechanical damage in a confined space.

[0066] The system acquires the real-time strain of frozen water particles (i.e., ice particles) under compression within a confined space. When the real-time strain exceeds a critical strain value that does not damage stiffness, the effective cumulative strain that damages stiffness is calculated. Based on the effective cumulative strain, the real-time normal stiffness of the ice particles is calculated using a step-type stiffness descent function, ensuring that the real-time normal stiffness decreases monotonically with increasing effective cumulative strain. Specifically:

[0067] Real-time strain of ice particles is monitored at each computation time step. When real-time response Exceeding the critical strain value without damaging stiffness At this time, the ice particle damage constitutive module is activated, and the mechanical response behavior of the ice particles is quantified through a step-type stiffness descent function. The step-type stiffness descent function is as follows:

[0068] ;

[0069] In the formula, is the minimum stiffness of ice particles, expressed in N / m; h is the stiffness reduction factor. For strain accumulation that damages stiffness, where, Critical stress caused by frost heave With initial stiffness jointly determined, that is It should be noted that the stiffness reduction coefficient h was manually calibrated based on the mechanical characteristics of rocks and ice, combined with measured data of frost heave force, through parameter sensitivity analysis. Calculations show that when h is set to 50, the simulated frost heave force attenuation curve and ice damage evolution process have the highest degree of matching with the actual situation.

[0070] Traditional discrete element method (DEM) simulations often simplify ice bodies as ideal elastic or elasto-plastic bodies, leading to spurious stiffness recovery. In actual physical processes, when fissure water freezes and expands under strong constraint from surrounding rock, irreversible structural damage occurs within the ice body, such as microcracks, lattice misalignment, or creep. If a conventional elastic model is used, the stiffness of ice particles instantly returns to its initial undamaged state when the temperature rises or fissures expand, causing the frost heave force to be unloaded. This clearly violates the objective law of material damage accumulation. This embodiment introduces accumulated strain as a damage-driving variable and constructs a step-type stiffness decrease function. During the simulation, it ensures that after ice particle stiffness is damaged, the stiffness only monotonically decreases with increasing accumulated strain, rather than increasing as expected. This accurately simulates the process of ice bodies undergoing compression within confined fissure spaces, resulting in microcracks and creep damage. This has significant engineering guidance for predicting cumulative fatigue damage in rocks under multiple freeze-thaw cycles.

[0071] The key features that distinguish this invention from existing technologies are as follows: First, it abandons the erroneous practice of using the coefficient of thermal expansion to drive phase change volume expansion, as is common in traditional methods. Instead, it strictly distinguishes between phase change (fixed radius increment) and pure thermal expansion and contraction (coefficient of thermal expansion), greatly improving the accuracy of the simulation results in terms of physical principles. Second, it achieves synchronous coupling between the abrupt change in particle stiffness and volume expansion, overcoming the limitation of existing technologies that rely solely on volume expansion to calculate frost heave force, realizing the synchronous occurrence of particle stiffness abrupt change and volume expansion in the discrete element method. Third, addressing the deficiency of statically assigning thermodynamic parameters (specific heat capacity, thermal resistance) in traditional simulations, it establishes a particle thermodynamic... The real-time dynamic update mechanism of parameters precisely controls the advance speed of the freeze-thaw front, corrects the deviations in the timing of frost heave force generation and temperature field distribution caused by neglecting the latent heat of phase change, and realizes the accurate evolution of the temperature field-phase change field. By introducing strain accumulation as a driving variable, a step-type stiffness descent function is constructed. In the simulation process, it is realized that after the stiffness of ice particles is damaged, the stiffness only decreases monotonically with the increase of accumulated strain, rather than increasing in a way that contradicts the objective phenomenon. It accurately simulates the process of ice body being squeezed in a confined fracture space to generate microcracks and creep damage, which has high engineering guiding significance for predicting the cumulative fatigue damage of rocks under multiple freeze-thaw cycles.

[0072] The above description represents a preferred embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any equivalent substitutions or modifications made by those skilled in the art within the scope of the technology disclosed in the present invention, based on the technical solution and inventive concept of the present invention, should be covered within the scope of protection of the present invention.

Claims

1. A discrete-element freeze-thaw simulation method that takes into account both the latent heat of water-ice phase change and the dynamic evolution of thermodynamic properties, characterized in that, include: A discrete element model of rock containing prefabricated fracture structures was constructed, and water particles were filled into the fracture structures. Within the phase change temperature range, the thermodynamic parameters of water particles are dynamically updated based on the equivalent heat capacity model to simulate the latent heat effect of phase change. Specifically, water particles are regarded as independent thermodynamic units, and their specific heat capacity and thermal conductivity are dynamically adjusted according to the real-time temperature and phase change state. Through the dynamic adjustment of specific heat capacity and thermal conductivity, the absorption or release of latent heat during the phase change process is simulated to generate a temperature plateau period in the temperature field evolution. Within the phase change temperature range, the mechanical parameters of water particles are simultaneously updated by step and the radius is updated by increment. The radius update is decoupled from the thermal expansion and contraction calculation during the non-phase change period. Stiffness decay calculations are performed on frozen water particles based on strain accumulation to simulate irreversible mechanical damage in a confined space.

2. The discrete element freeze-thaw simulation method according to claim 1, characterized in that, The construction of a rock discrete element model containing a prefabricated fracture structure and the filling of the fracture structure with water particles includes: generating rock particles within a defined model area and establishing interparticle bonding; deleting rock particles at specific locations to generate the prefabricated fracture structure; and filling the fracture structure space with water particles based on the geometric volume and porosity of the prefabricated fracture and establishing contact with the rock particles on the fracture wall.

3. The discrete element freeze-thaw simulation method according to claim 2, characterized in that, After generating rock and water particles, the process also includes: obtaining the macroscopic thermal conductivity of the rock, water, and ice; converting the macroscopic thermal conductivity into microscopic particle contact thermal resistance based on the relationship between the macroscopic thermal conductivity and the microscopic particle contact thermal resistance; and assigning the microscopic particle contact thermal resistance to the corresponding rock and water particles for heat conduction calculations.

4. The discrete element freeze-thaw simulation method according to claim 1, characterized in that, Define the phase transition temperature range as [ , The specific heat capacity of water particles is updated according to the following formula: ； In the formula, and These are the specific heat capacities of water and ice, respectively. To account for the equivalent volumetric heat capacity of the latent heat of phase change, The latent heat per unit mass of water; The density of water, The real-time temperature of the particles; for thermal conductivity, when the temperature of the water particles decreases to... At that time, its thermal conductivity was switched from that of water phase to that of ice phase.

5. The discrete element freeze-thaw simulation method according to claim 1, characterized in that, Within the phase change temperature range, the simultaneous triggering of step updates and radius increment updates of the mechanical parameters of water particles, and the decoupling of the radius increment update from the thermal expansion and contraction calculation during the non-phase change period, includes: when water particles enter the phase change temperature range, the water particles are transformed into ice particles, causing a step increase in the stiffness of the ice particles, while a fixed increment update is assigned to the radius of the ice particles; outside the phase change temperature range, the radius change caused by temperature changes is calculated, and the radius change is linearly superimposed with the fixed increment of the radius to achieve physical separation of phase change expansion and thermal expansion and contraction.

6. The discrete element freeze-thaw simulation method according to claim 5, characterized in that, The stiffness update formula is as follows: ； In the formula, The initial simulated normal stiffness when water particles freeze into ice particles; The simulated normal stiffness of ice particles; This is the stiffness amplification factor, and ≥4.

7. The discrete element freeze-thaw simulation method according to claim 5, characterized in that, When the particle is outside the phase transition temperature range, the particle radius depends on the coefficient of thermal expansion to drive the volume change; when the particle is inside the phase transition temperature range, the particle radius is directly assigned a fixed radius increment proportional to the radius when entering the phase transition temperature range to simulate the volume expansion caused by the water-ice phase transition; when outside the phase transition temperature range, the radius change caused by temperature change is calculated and linearly superimposed with the fixed radius increment inside the phase transition temperature range to obtain the current radius of the particle.

8. The discrete element freeze-thaw simulation method according to claim 1, characterized in that, The calculation of stiffness decay of frozen water particles based on strain accumulation includes: obtaining the real-time strain generated by the compression of frozen water particles in a confined space; calculating the effective strain accumulation that damages stiffness when the real-time strain exceeds the strain threshold that does not damage stiffness; and calculating the real-time normal stiffness of the ice particles based on the effective strain accumulation using a step-type stiffness descent function, so that the real-time normal stiffness decreases monotonically with the increase of the effective strain accumulation.

9. The discrete element freeze-thaw simulation method according to claim 8, characterized in that, The step-type stiffness descent function is as follows: ； In the formula, The simulated normal stiffness of ice particles. is the minimum stiffness of the ice particles; h is the stiffness reduction factor; For strain accumulation that damages stiffness, where, To respond in real time, Let be the critical strain value, which satisfies: ,in, This is the critical stress for frost heave damage.

Images