Phase field simulation method for buffer material gas breakthrough path of high-level radioactive waste repository

Abstract

The invention discloses a phase field simulation method for a buffer material gas breakthrough path of a high-level radioactive waste repository, and relates to the technical field of high-level radioactive waste disposal. According to the method, a heat-water-force (THM) coupled phase field cohesion model (PF-CZM) is constructed, deformation, temperature and water-gas two-phase flow of a saturated buffer material are comprehensively considered, a stable fluid source item is introduced to improve the robustness of the model, and formation of dominant flow channels, temperature distribution, pore pressure and crack phase field evolution in the gas breakthrough process can be accurately captured. Through comparison with a classical analytical solution and a numerical solution, the reliability of the model is verified, the model is used for simulating the gas breakthrough process of the saturated buffer material of the high-level radioactive waste repository, and the influence of parameters such as temperature, boundary rigidity and material heterogeneity on gas breakthrough is analyzed; theoretical support is provided for optimization of a buffer material and design of a deep geological repository, and the problems that an existing model neglects the temperature influence, and a dominant flow channel cannot be accurately simulated are solved.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of high-level waste disposal safety assessment, and particularly relates to a phase field simulation method for a gas breakthrough path of a high-level waste disposal repository buffer material. BACKGROUND

[0002] The safe disposal of high-level nuclear waste is a key challenge for the sustainable development of nuclear energy, and a deep geological disposal repository is an internationally recognized feasible solution. As a core component of the engineering barrier, the buffer material needs to have low permeability, high water retention, and self-healing ability, and the commonly used material is GMZ bentonite. During the operation of the disposal repository, a large amount of gas is generated by processes such as metal corrosion, water radiolysis, and microbial degradation, and the gas migrates in the saturated buffer material and forms high pressure, which can easily induce crack propagation and preferential flow channels, threatening the integrity of the engineering barrier.

[0003] The existing gas breakthrough simulation methods have many limitations: the traditional hydraulic fracturing model ignores the influence of temperature on the properties of gas and the mechanical behavior of materials; the conventional phase field method does not consider the quasi-brittle softening characteristics of the buffer material, and the calculation stability is insufficient; most models assume that the material is homogeneous, and it is difficult to reflect the control effect of the heterogeneity of the actual geological material on the crack path. For example, the phase field model proposed by Guo et al. only considers water-force coupling and does not include temperature effects; Liaudat et al. use zero-thickness interface elements to simulate cracks, which requires pre-setting the crack location and cannot predict complex crack branching and expansion. Therefore, it is of great significance to develop a gas breakthrough simulation method that considers multi-field coupling, material heterogeneity, and calculation stability. SUMMARY

[0004] The application aims to provide a phase field simulation method for a gas breakthrough path of a high-level waste disposal repository buffer material, which solves the problems of existing models that ignore temperature effects, have poor calculation stability, and cannot reflect material heterogeneity.

[0005] The technical solution adopted by the application is: a phase field simulation method for a gas breakthrough path of a high-level waste disposal repository buffer material, and the specific operation steps are as follows: By constructing a THM-coupled phase field cohesive zone model PF-CZM, the deformation of the saturated buffer material, the gas temperature, and the water-gas two-phase flow are comprehensively considered, a stable fluid source term is introduced to improve the robustness of the model, and the formation of the preferential flow channel, the temperature distribution, the pore pressure, and the crack phase evolution in the gas breakthrough process are accurately captured; the specific operation steps are as follows: Step 1: Construct a THM-coupled phase field cohesive zone model PF-CZM, which includes control equations for displacement field, phase field, temperature field, and gas-water pressure field; Step two: define the phase field approximation function of the buffer material crack, regularize the sharp crack by using the crack surface density function, and describe the material damage evolution by the energy degradation function; Step three: establish the mass conservation equation of gas-water two-phase flow, describe the fluid seepage velocity based on Darcy's law, and define the water retention curve and relative permeability of the intact buffer material and the crack area by using the van-Genuchten model; Step four: introduce the energy conservation equation under the condition of local thermal equilibrium, calculate the heat conduction flux of the buffer material by Fourier's law, and consider the influence of heat convection and temperature on the fluid density and viscosity; Step five: establish the evolution equation of intrinsic porosity, and then obtain the evolution equation of intrinsic permeability, update the crack permeability by the phase field variable and strain, further identify the crack area and matrix area by the conversion function, and give different seepage parameters; Step six: discretize the control equation by using the finite element method, introduce the polynomial pressure projection stabilization term in the gas-water pressure field equation, use the implicit backward difference format for time integration, and solve the coupled equation group by the separation scheme; Step seven: set the simulation boundary conditions and initial parameters, input the physical and mechanical parameters of the buffer material (length scale parameter, crack geometry function, etc. 、 、 、 、 、 、 、 、 、 、 、 ), gas generation rate and boundary stiffness coefficient, and obtain the crack propagation path, pressure distribution, saturation and temperature evolution results in the gas breakthrough process by numerical simulation.

[0006] The characteristics of the present application are also, The crack phase field approximation function and the energy degradation function in step two are as follows: Based on the phase field fracture theory, the sharp crack is regularized and approximated by using the crack surface density function, and the crack surface density function is Satisfies: (1) (2) (3) Wherein, is the crack phase field variable, is the gradient operator, b is the length scale parameter, is the crack geometry function and , satisfy , is a scale parameter, is an integrand function, is an integral variable; The material damage is described by a monotonically decreasing energy degradation function, which satisfies : (4) (5) where, is an energy degradation auxiliary function, p is a power parameter, is a polynomial function; The parameter is obtained by the crack geometry function and the energy degradation function : , , : (6) where, is the elastic modulus, is the critical energy release rate, is the tensile strength, is the initial slope of the softening curve, is the crack limit displacement; The above parameters and are determined by the corresponding cohesive force law, when considering the exponential softening curve and , we can get: (7) The control equations of the displacement field and the phase field in step one are as follows: For damaged solids, the local energy functional can be expressed by the phase field d and the elastic strain tensor : (8) According to the local energy functional, the effective stress tensor is obtained: (9) The thermal strain is expressed as: (10) where, is the fourth-order elastic tensor, thermal expansion coefficient, T is the current temperature, T0 is the reference temperature, F is the second order identity tensor; For quasi-static fracture of a solid under small strain, the mechanical equilibrium equation is: (11) where, σ is the stress tensor, f is the body force vector, Ω is the solid domain, n is the unit normal vector on the boundary of the solid domain, t is the surface traction, N is the Neumann boundary, u is the displacement vector, u0 is the prescribed displacement vector, D is the Dirichlet boundary; The phase field evolution equation is obtained by the Kuhn-Tucker loading and unloading conditions: (12) where, ∂φ / ∂t is the first order derivative of the phase field variable with respect to time, G is the energy dissipation function; According to and the variational derivative of the crack surface density function , the effective energy release rate is defined to consider different mechanical behaviors under tensile and compressive stress states: (13) (14) (15) (16) where, G is the energy release rate, σ1 is the principal stress, Eeff is the effective elastic modulus, λ is the Lame constant, Mac is the Macaulay bracket; The governing equation and boundary conditions of the phase field are derived from the analysis: (17) where, Ωc is the crack band, n is the unit normal vector on the outer boundary of the crack band, ∂Ωc is the outer boundary of the crack band; The phase field must satisfy the irreversibility condition in the phase field fracture theory, which is numerically implemented by the effective energy release rate replaced by its maximum value H : (18) where, H is the maximum value of the effective crack driving force.

[0007] The mass conservation equation for gas-liquid two-phase flow in step three considers the porosity variation, fluid compressibility and thermal expansion effect, and the expression is: (19) where, is the density of water, is the porosity, is the compressibility of water, is the saturation of water, is the capillary pressure, is the gas pressure, is the seepage velocity of water, is the volumetric strain, is the thermal expansion coefficient of water, is the temperature, is the density of gas, is the compressibility of gas, is the saturation of gas, is the seepage velocity of gas, is the water pressure, is the thermal expansion coefficient of gas, is time; The flow rates of water and gas follow Darcy's law: (20) where, is the intrinsic permeability, is the relative permeability of water, is the relative permeability of gas, is the dynamic viscosity of water, is the acceleration of gravity, is the dynamic viscosity of gas; The mass conservation equation applies Dirichlet boundary conditions and Neumann boundary conditions, and the expression is: (21) where, is the water pressure, is the boundary water pressure, is the unit normal vector of the boundary, is the density of water, is the seepage velocity, is the mass flux of water.

[0008] The van-Genuchten model described in Step three defines the water retention curve and relative permeability of the intact buffer material and the fracture region: Water retention curve: (22) where, is the effective saturation, is the characteristic pressure, subscript i is taken as m for the intact buffer material region, and f for the fracture region, n and m are shape parameters that satisfy ; According to equation (22), the effective saturation is a function of the capillary pressure The first derivative of the effective saturation with respect to the capillary pressure is: (23) The effective water saturation is given by: (24) where, is the water residual saturation, is the gas residual saturation; The capillary pressure is given by: (25) The relative permeability of the porous matrix: (26) The relative permeability of the fracture: (27) where, is the relative permeability of water in the fracture, is the relative permeability of gas in the fracture.

[0009] The energy conservation equation described in Step four expresses the heat conduction by Fourier's law: (28) (29) where, subscript eff is the effective value, is the specific heat capacity of water at constant pressure, is the density of the solid matrix, is the specific heat capacity of the solid matrix at constant pressure, is the effective thermal conductivity, is the thermal conductivity of water, is the thermal conductivity of gas, thermal conductivity of the solid matrix, for the heat flux; The energy conservation equation applies Dirichlet and Neumann boundary conditions, expressed as: (30) where, is the boundary temperature, is the surface heat flux.

[0010] The evolution equation of the intrinsic porosity in Step 5 is: (31) Considering the cubic relationship between intrinsic permeability and porosity: (32) The generated crack can be regarded as two parallel plates, and the anisotropic crack permeability is: (33) (34) (35) where, is the Biot coefficient, is the volumetric compression coefficient of the particles, is the average pressure, the solid thermal expansion coefficient, h e is the unit size; The conversion function in Step 5 is: (36) (37) where, is the damage interpolation function, is the hyperbolic tangent function, is the physical quantity of the crack phase, is the physical quantity of the matrix phase, is the transition zone width parameter, is the damage critical value; The finite element weak form of the displacement field, phase field, gas, liquid, and heat transfer in Step 6 is: Displacement field: (38) where, is the reference pore pressure, is the linear strain tensor of the displacement variation, is the body force vector, is the volume,​ is the surface force vector, is the area; is the phase field: (39) is the water pressure field: (40) is the gas pressure field: (41) is the temperature field: (42) A polynomial pressure projection based stabilization term is added to the water and gas mass continuity equations, (43) where, G is the shear modulus, is the element volume, is the trial function of the pore pressure, is the first order derivative of the pore pressure with respect to time.

[0011] The numerical solving module in step seven is integrated in the COMSOL software platform, the solid mechanics module is used to solve the displacement field, the Darcy law module is used to solve the gas and water pressure field, the porous medium heat transfer module is used to solve the temperature field, the Poisson equation module is used to solve the phase field evolution, and the historical maximum value of the effective driving of the phase field is determined by using the state variable H .

[0012] The spatial correlation non-uniform field generation method of the rotating belt method is adopted, the variance factor and the reference variance coefficient are introduced, the discreteness of the random field is quantified, and the heterogeneity of the buffer material is accurately simulated; the phase field method can clearly present the details of crack propagation: as the variance factor increases, the crack propagation path gradually changes, and even bifurcation phenomenon occurs, which reflects the high-precision simulation capability of the crack dynamic evolution process in the non-uniform medium; combined with the gas pressure response, the coupling simulation of the crack phase field and the gas seepage field is realized, and the synergistic effect of heterogeneity on the multi-process of "crack propagation-gas breakthrough" can be systematically studied.

[0013] The present application has the following advantages: Theoretical innovation: THM coupling and phase field cohesive zone model PF-CZM are combined for the first time, the influence of temperature on gas properties and material damage is considered, and the defects of the existing model ignoring multi-field coupling are solved; Numerical advantage: the stabilization term is introduced, low-order elements are allowed to be used, the calculation efficiency is improved, and the numerical robustness is ensured at the same time; Engineering value: It can quantitatively analyze the impact of multiple parameters on gas breakthrough, providing direct basis for the optimization of buffer material composition (such as adjusting the compaction degree of bentonite to change the permeability) and the design of disposal repository structure (such as optimizing the barrier contact stiffness), reducing the risk of radioactive leakage and ensuring the safety of nuclear waste disposal. Attached Figure Description

[0014] To more clearly illustrate the embodiments of the present invention, the accompanying drawings used in the embodiments of the present invention will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0015] Figure 1 This is a flowchart of the method of the present invention; Figure 2 This is the calculation model of the present invention; Figure 3 It is a distribution diagram of multiphysics variables during the gas breakthrough process; Figure 4 It is a diagram showing the energy evolution of the gas breakthrough process in the buffer material; Figure 5 This is a diagram showing the evolution of water pressure at various points along the boundary; Figure 6 This is a spatial distribution map of the average pressure at the model boundary I at different times; Figure 7 This is a graph showing the evolution of gas pressure at different temperatures; Figure 8 This is a diagram showing the effect of the variance factor on the damaged area. Detailed Implementation

[0016] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without creative effort are within the scope of protection of the present invention.

[0017] Example 1 This embodiment provides a phase-field simulation method for the gas breakthrough path of the buffer material in a high-level radioactive waste disposal repository. The following is a detailed description of the above scheme in conjunction with simulation results: This invention uses saturated buffer materials as the research object, focusing on the thermal-thermal-mechanical coupling effect during gas breakthrough. A numerical model is constructed using the COMSOL platform to verify the effectiveness of the phase-field cohesion model (PF-CZM) and analyze the influence of key parameters on gas breakthrough. In the simulation, the crack internal pressure is 5 × 10⁻⁶. 3Pa / s, with constant heat flux Q = -3 x 10 3 W / m 2 , fixed pore pressure of 0, fixed temperature T = 20 °C, due to the symmetry of the system, only a quarter of the sample was simulated. The computational domain was discretized using bilinear quadrilateral elements with the maximum element size of (b / 2), where b is specified as 10 mm, the total number of elements was 14350, the implicit backward difference format was used for time stepping, the initial time step was 0.001 day, and the maximum time step was 5 days. A stopping condition was added, i.e., max (u d > 0.4) on boundary IV, to determine whether the gas had completely broken through. The key parameters are listed in Table 1. Here, to reduce the amount of calculation, the permeability of the developing fracture was set to be 1 x 10 8 times the permeability of the porous matrix. See Table 2 for detailed information on the boundary conditions of the model. The inner boundary V of the model was specified as 1 in the phase field to simulate the pre-existing fracture, and in the gas flow field, this boundary was given a constant mass flux.

[0018] Table 1. Related parameters of the bentonite

[0019] Table 2. Boundary conditions of the model

[0020] It is assumed that both hydraulic and mechanical properties can obey the lognormal distribution. This simulation uses MATLAB to generate a two-dimensional spatial autocorrelation random seed that satisfies the normal distribution, and then it is imported into COMSOL as a random seed through the interpolation function int(x). It is converted into a lognormal distribution of heterogeneous field with parameters such as permeability, elastic modulus, tensile strength, and inlet capillary pressure by the following formula.

[0021]

[0022] where, is the mean of the normal distribution of the random field, is the standard deviation of the normal distribution of the random field; Combining the above simulation, we can get Figure 3 the distribution of multi-physical field variables during the gas breakthrough process. As time increases, gas pressure gradually accumulates, and when it reaches a certain value, cracks begin to appear. Due to the existence of the heterogeneous field, the phase field cracks are not smooth and do not extend along a straight line. The water saturation in the crack propagation path is close to 0, because the gas that breaks through along the high-permeability crack channel flows rapidly and occupies the entire crack area. The high-temperature gas flows along the developed cracks, resulting in a temperature distribution closely related to the crack morphology. This process also changes the thermal stress distribution, which in turn affects the crack propagation path.

[0023] Figure 4 It can be seen that the elastic strain energy , the fracture energy and the pressure-dependent energy all show a growing trend over time, but the growth rate varies. The fracture energy continues to increase as the crack expands, and its value reaches about 1 J / m when the gas completely breaks through. The elastic strain energy first increases and then remains almost constant during the rapid expansion of the crack. There is a slight decrease until it approaches complete breakthrough. Its value is almost equal to the fracture energy. The evolution of energy related to pressure can be divided into three stages: slow increase, rapid rise and slow rise. Its final value is more than four times the fracture energy. Thermal energy is three orders of magnitude larger than the above three energies, and its evolution is similar to that of elastic energy.

[0024] (44) Figure 6 The spatial distribution of average pressure at the model boundary I at different times is shown. The average pressure distribution has no obvious rule. However, there is a significant difference in pressure distribution between the pre-crack area and the porous matrix. When the time is from 100 days to 300 days, the average pressure in the crack area increases significantly, while the average pressure in the matrix has no obvious change. When the time is from 300 days to 600 days, the average pressure in the crack area gradually decreases, but the average pressure in the entire matrix area increases significantly. Figure 5 The water pressure changes of the four observation points at the upper left and lower right boundaries of the model are analyzed. It can be seen from the figure that the crack is along the upper left boundary, so the peak pressure of point #1 is higher than that of point #3 in the figure, and the peak pressure of point #2 is higher than that of point #4.

[0025] With the increase of temperature, there are significant differences in the crack phase field. In particular, when the temperature rises from 100 ℃ to 110 ℃, the main crack changes from upward deflection to downward deflection. At the same time, a circle of cracks is formed along the inner boundary I. When the gas temperature is 120 ℃, the crack splits at the initial stage and finally expands upward.

[0026] As shown in Figure 7 , when the gas mass flux remains unchanged, the higher the boundary stiffness, the higher the peak gas pressure, but the time to reach the peak is gradually delayed. However, Figure 8It can be seen from the evolution of the damage zone that the earliest cracking occurs at the lowest pressure peak, i.e. the lowest stiffness. Over time, although the gas pressure is also rising, it is not enough to drive the crack to continue to expand. Until close to the 9000th day, the calculation is terminated. When both the gas mass flux and the boundary stiffness change, the greater the boundary stiffness and the gas mass flux, the greater the corresponding peak, but the time to reach the peak is advanced. It can be observed from the damage evolution curve that under high boundary stiffness and high gas mass flux, the crack first initiates, but lacks sufficient energy to drive the crack to expand until about the 3600th day, when the gas completely breaks through. There is a delay of 2550 days compared to the earliest breakthrough.

[0027] Introducing a variance factor C v to characterize the dispersion of the random field under the same random seed. The standard deviation of the lognormal distribution of the random field is related to the mean value, and for the permeability, C v is set to 10 and other parameters are set to 0.2. The greater the variance factor, the greater the heterogeneity. As the variance factor increases, the crack propagation path gradually changes, especially when C v = 4, the crack bifurcates halfway and eventually realizes complete breakthrough along one side. As shown in Figure 8 , as C v decreases, the peak of the gas pressure is higher, and the time to reach the peak is gradually delayed. Greater C v corresponds to earlier cracking, and when C v = 0.4, the crack propagation rate is very slow between 1075 days and 1750 days, which belongs to the pressure accumulation stage. When the pressure accumulates to a certain extent, the crack bifurcates.

[0028] (45) wherein, the mean value of the lognormal distribution of the random field; In summary, the present application proposes a novel THM coupled phase field cohesive zone model PF-CZM framework, which overcomes the key limitations in existing models. The proposed method contains three innovations: (1) a new phase field-cohesive zone model PF-CZM specifically suitable for multiphase transport in buffer materials is developed; (2) an advanced stabilization technique is implemented through a bilinear equal-order quadrilateral element to ensure the robustness of the calculation; (3) a comprehensive THM coupling scheme is built within the platform, successfully linking multiple physical processes. This integrated computational framework provides new capabilities for simulating complex gas migration behavior, which was previously unattainable by traditional methods.

[0029] Example 2 The phase field simulation method for the gas breakthrough path of the high-level waste disposal repository buffer material of the application, by constructing the THM coupled phase field cohesive zone model PF-CZM, comprehensively considers the deformation of the saturated buffer material, the gas temperature and the water-gas two-phase flow, introduces a stable fluid source term to improve the robustness of the model, and accurately captures the formation of the dominant flow channel, the temperature distribution, the pore pressure and the crack phase evolution in the gas breakthrough process.

[0030] As Figure 2 shown is a simulation model of the buffer material, a grid, the buffer material is wrapped around the high-level waste tank, the waste tank is filled with glass solidification, and the buffer material is surrounded by granite, and the wrapping thickness of the buffer material is about 350mm.

[0031] Example 3 The phase field simulation method for the gas breakthrough path of the high-level waste disposal repository buffer material of the application, as Figure 1 shown, the specific operation steps are as follows: Step 1: Construct the THM coupled phase field cohesive zone model PF-CZM, which includes the control equations of the displacement field, the phase field, the temperature field and the gas-water pressure field; Step 2: Define the phase field approximation function of the buffer material crack, regularize the sharp crack by using the crack surface density function, and describe the material damage evolution by using the energy degradation function; Step 3: Establish the mass conservation equation of the gas-water two-phase flow, describe the fluid seepage velocity based on Darcy's law, and define the water retention curve and the relative permeability of the intact buffer material and the crack region by using the van-Genuchten model; Step 4: Introduce the energy conservation equation under the local thermal equilibrium condition, calculate the heat conduction flux of the buffer material by using the Fourier law, and consider the influence of heat convection and temperature on the fluid density and viscosity; Step 5: By establishing the evolution equation of the intrinsic porosity, the evolution equation of the intrinsic permeability is obtained, the crack permeability is updated by the phase variable and the strain, and further the crack region and the matrix region are identified by using the conversion function and different seepage parameters are given; Step 6: Discretize the control equation by using the finite element method, introduce the polynomial pressure projection stabilization term in the gas-water pressure field equation, use the implicit backward difference format for time integration, and solve the coupled equation groups by using the separation scheme; Step 7: Set the simulation boundary conditions and initial parameters, input the physical and mechanical parameters of the buffer material, the gas generation rate and the boundary stiffness coefficient, and obtain the crack propagation path, the pressure distribution, the saturation and the temperature evolution results in the gas breakthrough process by numerical simulation.

[0032] Example 4 Based on Example 3, The crack phase field approximation function and energy degradation function in step two are as follows: Based on phase-field fracture theory, a regularized approximation of sharp cracks is performed using the crack surface density function. satisfy: (1) (2) (3) in, For the crack phase field variables, For gradient operators, b For length scale parameters, The crack geometry function and ,satisfy , For scale parameters, ; Material damage is described using a monotonically decreasing energy degradation function. satisfy: (4) (5) in, For energy degradation auxiliary function, p For power-order parameters, It is a polynomial function; Through crack geometry function and have Energy degradation function Get parameters , , : (6) in, For elastic modulus, The critical energy release rate. For tensile strength, The initial slope of the softening curve is given. This represents the ultimate displacement of the crack. The above parameters and This is determined by the corresponding cohesive force rule, when considering the exponential softening curve and... We can obtain: (7) The governing equations for the displacement field and phase field in step one are as follows: For damaged solids, the local energy functional can be expressed as d and the elastic strain tensor : (8) The effective stress tensor : (9) The thermal strain is expressed as (10) where is the fourth-order elastic tensor, the thermal expansion coefficient, the current temperature, the reference temperature, is the second-order unit tensor; For quasi-static fracture of solids under small strain, the mechanical equilibrium equation is (11) where is the stress tensor, is the body force vector, is the solid domain, is the unit normal vector on the boundary of the solid domain, is the surface force, is the Neumann boundary, is the displacement vector, is the prescribed displacement vector, is the Dirichlet boundary; The phase field evolution equation is obtained by the Kuhn-Tucker loading and unloading conditions: (12) where is the first-order derivative of the phase field variable with respect to time, is the energy failure function; According to and the variational derivative of the crack surface density function , the effective energy release rate is defined considering different mechanical behaviors under tensile and compressive stress states: (13) (14) (15) (16) where is the energy release rate, principal stress, effective elastic modulus, Lame constant, Macaulay brackets; The phase field governing equation and boundary conditions are derived from the analysis: (17) where, fracture zone, unit normal vector of the outer boundary of the fracture zone, outer boundary of the fracture zone; The phase field must satisfy the irreversibility condition, in numerically the effective energy release rate is replaced by its maximum value H : (18) where, H maximum value of the effective crack driving force.

[0033] Example 5 Based on Example 4, The mass conservation equation of gas-liquid two-phase flow in Step Three considers the change of porosity, fluid compressibility and thermal expansion effect, and the expression is: (19) where, density of water, porosity, compressibility coefficient of water, saturation of water, capillary pressure, gas pressure, seepage velocity of water, volume strain, thermal expansion coefficient of water, temperature, density of gas, compressibility coefficient of gas, saturation of gas, seepage velocity of gas, water pressure, thermal expansion coefficient of gas, time; The flow rates of water and gas follow Darcy's law: (20) where, intrinsic permeability, relative permeability of water, relative permeability for gas, dynamic viscosity of water, acceleration of gravity, dynamic viscosity of gas; The mass conservation equation applies Dirichlet boundary conditions and Neumann boundary conditions, expressed as: (21) where, water pressure, boundary water pressure, unit normal vector of the boundary, density of water, seepage velocity, mass flux of water.

[0034] The van-Genuchten model described in Step Three defines the water retention curve and relative permeability of the intact buffer material and the fracture region: Water retention curve: (22) where, effective saturation, characteristic pressure, subscript i take m intact buffer material region, take f then represents the fracture region, n and m shape parameters, satisfying ; According to equation (22), the effective saturation is the first derivative of the capillary pressure : (23) Effective water saturation : (24) where, residual saturation of water, residual saturation of gas; Capillary pressure : (25) Relative permeability of porous matrix: (26) Relative permeability of fracture: (27) where, relative permeability of water in the fracture, Relative permeability of gas in the fracture.

[0035] The energy conservation equation in Step Four expresses the heat conduction by Fourier's law: (28) (29) where subscript eff is the effective value, is the specific heat capacity of water at constant pressure, is the density of the solid matrix, is the specific heat capacity of the solid matrix at constant pressure, is the effective thermal conductivity, is the thermal conductivity of water, is the thermal conductivity of gas, is the thermal conductivity of the solid matrix, is the heat flux; The energy conservation equation applies Dirichlet boundary conditions and Neumann boundary conditions, expressed as: (30) where is the boundary temperature, is the surface heat flux.

[0036] The evolution equation for the intrinsic porosity in Step Five is: (31) The intrinsic permeability is considered to have a cubic relationship with the porosity: (32) The resulting fracture can be viewed as two parallel plates, and the anisotropic fracture permeability is: (33) (34) (35) where is the Biot coefficient, is the bulk volume compressibility of the grains, is the average pressure, is the coefficient of thermal expansion of the solid, h e is the cell size; The conversion function in Step Five is: (36) (37) where is the damage interpolation function, is the hyperbolic tangent function, is the physical quantity of the fracture phase, is the physical quantity of the matrix phase, is the transition zone width parameter, is the damage critical value.

[0037] Example 6 Based on the example 5, The finite element weak forms of the displacement field, phase field, gas, liquid and heat transfer in step six are as follows: Displacement field: (38) where, is the reference pore pressure, is the linear strain tensor of the displacement variation, is the volume force vector, is the volume, is the surface force vector, is the area; Phase field: (39) Water pressure field: (40) Gas pressure field: (41) Temperature field: (42) A polynomial pressure projection based stabilization term is added to the water and gas mass continuity equations, (43) where, G is the shear modulus, is the element volume, is the test function of the pore pressure, is the first order derivative of the pore pressure with respect to time.

[0038] The numerical solution module in step seven is integrated into the COMSOL software platform, the solid mechanics module is used to solve the displacement field, the Darcy law module is used to solve the gas and water pressure field, the porous medium heat transfer module is used to solve the temperature field, the Poisson equation module is used to solve the phase field evolution, and the historical maximum value of the effective driving of the phase field is determined by the state variable H .

[0039] The above merely describes preferred embodiments of the present application, and is not used to limit the present application. Any modification, equivalent replacement, improvement, etc. made by those skilled in the art within the spirit and principle of the present application shall be included in the protection scope of the present application.

Claims

1. A phase field simulation method of gas breakthrough paths in a buffer material of a high level waste repository, characterized in that, By constructing the THM-coupled phase field cohesive zone model PF-CZM, the deformation of the saturated buffer material, the temperature of the gas, and the water-gas two-phase flow are comprehensively considered, the stable fluid source term is introduced to improve the robustness of the model, and the formation of the dominant flow channel, the temperature distribution, the pore pressure, and the crack phase field evolution in the gas breakthrough process are accurately captured.

2. The phase field simulation method of gas breakthrough paths in a high level waste repository buffer material according to claim 1, characterized in that, The specific operation steps are as follows: Step one: constructing the THM-coupled phase field cohesive zone model PF-CZM, the phase field cohesive zone model PF-CZM includes the control equations of the displacement field, the phase field, the temperature field, and the gas-water pressure field; Step two: defining the phase field approximation function of the buffer material crack, using the crack surface density function to regularize the sharp crack, and describing the material damage evolution through the energy degradation function; Step three: establishing the mass conservation equation of the water-gas two-phase flow, describing the fluid seepage velocity based on Darcy's law, and using the van-Genuchten model to define the water retention curve and the relative permeability of the intact buffer material and the crack region respectively; Step four: introducing the energy conservation equation under the local thermal equilibrium condition, calculating the heat conduction flux of the buffer material through the Fourier law, and considering the influence of temperature on the fluid density and viscosity; Step five: by establishing the evolution equation of the intrinsic porosity, the evolution equation of the intrinsic permeability is obtained, the crack permeability is updated by the phase field variable and the strain, and further through the conversion function to identify the crack region and the matrix region and give different seepage parameters; Step six: using the finite element method to discretize the control equation, introducing the polynomial pressure projection stabilization term in the gas-water pressure field equation, using the implicit backward difference format for time integration, and solving the coupled equation group through the separation scheme; Step seven: setting the simulation boundary conditions and initial parameters, inputting the physical and mechanical parameters of the buffer material, the gas generation rate and the boundary stiffness coefficient, and obtaining the crack propagation path, pressure distribution, saturation and temperature evolution results in the gas breakthrough process through numerical simulation.

3. The phase field simulation method of gas breakthrough paths in a high level waste repository buffer material according to claim 2, characterized in that, The crack phase field approximation function and the energy degradation function in step two are as follows: Based on the phase-field fracture theory, the sharp crack is regularized and approximated by a crack surface density function, the crack surface density function satisfies: (1) (2) (3) wherein is a crack phase variable, is a gradient operator, b is a length scale parameter, is a crack geometry function and satisfies , is a scale parameter, is an integrand function, is an integration variable; The material damage is described by a monotonically decreasing energy degradation function, which is denoted as satisfies: (4) (5) wherein, is an energy degradation auxiliary function, p is a power parameter, is a polynomial function; By the crack geometry function and having an energy degradation function yielding parameters , , : (6) wherein, E is the modulus of elasticity, G is the critical energy release rate, σ is the tensile strength, m is the initial slope of the softening curve, δc is the crack limit displacement; The above parameters and By the corresponding cohesion law, when considering an exponential softening curve and , one obtains: (7)。 4. The phase field simulation method of gas breakthrough paths in a high level waste repository buffer material according to claim 3, characterized in that, The control equations of the displacement field and the phase field in step one are as follows: For damaged solids, the local energy functional The phase field d and the elastic strain tensor can be expressed as: (8) According to the local energy functional, the effective stress tensor is obtained : (9) thermal strain is represented by: (10) wherein, is the fourth order elasticity tensor, thermal expansion coefficient, is the current temperature, is the reference temperature, is the second order unit tensor; For the quasi-static fracture of a solid under small strain, the mechanical equilibrium equation is: (11) wherein, is a stress tensor, is a volume force vector, is a solid domain, is a unit normal vector to the solid domain boundary, is a surface force, is a Neumann boundary, is a displacement vector, is a prescribed displacement vector, is a Dirichlet boundary; The phase field evolution equation is obtained through the Kuhn-Tucker loading and unloading condition: (12) wherein is the first derivative of the phase field variable with respect to time, is the energy destruction function; According to and the variational derivative of the crack surface density function , considering the different mechanical behavior in tension and compression stress states, the effective energy release rate is defined as (13) (14) (15) (16) wherein, is the energy release rate, is the principal stress, is the effective elastic modulus, is the Lame constant, is the Macaulay bracket; The control equation of the phase field and the boundary condition are derived by analysis and deduction: (17) wherein is a crack band, is a unit normal vector to the outer boundary of the crack band, is an outer boundary of the crack band; In the phase field fracture theory the phase field The irreversible condition must be satisfied, in numerical terms the effective energy release rate is replaced by its maximum value H : (18) wherein, H is the maximum value of the effective crack driving force.

5. The phase field simulation method of gas breakthrough paths in a high level waste repository buffer material according to claim 4, characterized in that, The mass conservation equation of the gas-liquid two-phase flow considers the porosity variation, fluid compressibility and thermal expansion effect, and the expression is: (19) wherein, is the density of water, is the porosity, is the compressibility of water, is the saturation of water, is the capillary pressure, is the gas pressure, is the seepage velocity of water, is the volumetric strain, is the thermal expansion coefficient of water, is the temperature, is the density of gas, is the compressibility of gas, is the saturation of gas, is the seepage velocity of gas, is the water pressure, is the thermal expansion coefficient of gas, is time; The flow velocity of water and gas follows Darcy's law: (20) wherein, is the intrinsic permeability, is the relative permeability of water, is the relative permeability of gas, is the dynamic viscosity of water, is the acceleration of gravity, is the dynamic viscosity of gas; The mass conservation equation applies the Dirichlet boundary condition and the Neumann boundary condition, and the expression is: (21) wherein is the water pressure, is the boundary water pressure, is the unit normal vector of the boundary, is the density of water, is the seepage velocity, is the mass flux of water.

6. The phase field simulation method of gas breakthrough paths in a high level waste repository buffer material according to claim 5, characterized in that, The van-Genuchten model in step three defines the water retention curve and the relative permeability of the intact buffer material and the crack region: The water retention curve: (22) wherein is the effective saturation, is the characteristic pressure, subscript i take m is the complete cushioning material zone, take f then represents the crack zone, n and m is the shape parameter, satisfying ; According to equation (22), the effective saturation The first derivative of the capillary pressure is: (23) Effective water saturation : (24) wherein, is the water residual saturation, is the gas residual saturation; Capillary pressure : (25) The relative permeability of the porous matrix: (26) The relative permeability of the crack: (27) wherein, Krw is the relative permeability of water in the fracture, Krg is the relative permeability of gas in the fracture.

7. A phase field simulation method of gas breakthrough paths in a high level waste repository buffer material according to claim 6, characterized in that, The energy conservation equation in step four is expressed by the Fourier law to represent the heat conduction amount, and the expression is: (28) (29) wherein the subscript eff is the effective value, is the specific heat capacity of water at constant pressure, is the density of the solid matrix, is the specific heat capacity of the solid matrix at constant pressure, is the effective thermal conductivity, is the thermal conductivity of water, is the thermal conductivity of the gas, is the thermal conductivity of the solid matrix, is the heat flux; The energy conservation equation applies the Dirichlet boundary condition and the Neumann boundary condition, and the expression is: (30) where, is the boundary temperature, is the surface heat flux.

8. A phase field simulation method of gas breakthrough paths in a high level waste repository buffer material according to claim 7, characterized in that, The evolution equation of the intrinsic porosity in step five is: (31) Considering the cubic relationship between the intrinsic permeability and the porosity: (32) The resulting fracture can be visualized as two parallel plates, with anisotropic fracture permeability is: (33) (34) (35) in, Biot coefficient, The particle volume compressibility coefficient, For average pressure, Coefficient of thermal expansion of solids h e Unit size; The conversion function described in step five is: (36) (37) wherein, is a damage interpolation function, is a hyperbolic tangent function, is a physical quantity of the fracture phase, is a physical quantity of the matrix phase, is a transition zone width parameter, is a damage threshold value.

9. The phase field simulation method of gas breakthrough paths in a high level waste repository buffer material according to claim 8, characterized in that, The finite element weak forms of the displacement field, phase field, gas, liquid and heat transfer in step six are: Displacement field: (38) wherein, is the reference pore pressure, is the linear strain tensor of the displacement variation, is the volume force vector, is the volume, is the surface force vector, is the area; Phase field: (39) Water pressure field: (40) Gas pressure field: (41) Temperature field: (42) A polynomial pressure projection based stabilization term is added to the water and gas mass continuity equations, (43) wherein, G G is the shear modulus, V is the unit volume, p is a trial function for pore pressure, is the first derivative of the pore pressure with respect to time.

10. The phase field simulation method of gas breakthrough paths in a high level waste repository buffer material according to claim 2, characterized in that, The numerical solution module in step seven is integrated in the COMSOL software platform, a solid mechanics module is used to solve the displacement field, a Darcy law module is used to solve the gas and water pressure field, a porous medium heat transfer module is used to solve the temperature field, and a Poisson equation module is used to solve the phase field evolution, and the historical maximum value of the effective driving of the phase field is determined by using the state variable H .

Images