A method for simulating crack propagation in triso particle coating layers

Abstract

This invention discloses a method for simulating crack propagation in TRISO particle coatings. The main steps are as follows: 1. Establish a geometric model of the particle coating layer; 2. Set the physical properties of each layer material; 3. Determine the service environment parameters of the fuel; 4. Establish irradiation behavior and fission gas release models; 5. Establish a solid mechanics-phase field model and verify the accuracy of the model through classical fracture mechanics problems; 6. Introduce the geometric defects to be studied into each coating layer and calculate them; 7. After calculation, output the calculation results. This invention considers the influence of geometric defects, material irradiation behavior, and fission gas pressure on crack propagation in particle coating layers, and can accurately predict the initiation and propagation of cracks in the coating layer, providing numerical basis for TRISO particle safety assessment and structural optimization.

Description

Technical Field

[0001] This invention relates to the field of nuclear fuel technology, and more specifically to a TRISO particle coating crack propagation simulation method. Background Technology

[0002] As an accident-resistant fuel, TRISO particles, due to their excellent high-temperature resistance and fission product containment capacity, can be used in high-temperature gas-cooled reactors and pressurized water reactor systems. A typical TRISO particle consists of a fuel core and four cladding layers: a buffer layer, an inner pyrolytic carbon layer (IPyC), a silicon carbide layer (SiC), and an outer pyrolytic carbon layer (OPyC). The SiC layer, as the primary pressure-bearing layer and safety barrier, plays a crucial role in maintaining the structural integrity of the fuel particle. During reactor operation, the materials of each TRISO particle layer are subjected to the coupled effects of various irradiation and mechanical factors, such as irradiation effects, fission gas pressure, thermal stress, and residual stress, which can lead to cracks or even complete fracture of the cladding layer. The formation and propagation of cracks significantly weaken the mechanical properties of the cladding layer, increase the risk of radioactive product leakage, and thus affect the service safety of the fuel elements. Due to the high-temperature and high-irradiation environment, the crack evolution process is difficult to observe directly. Therefore, applying the phase-field method to simulate crack propagation in the TRISO particle cladding layer is of great significance for revealing the crack evolution mechanism and ensuring fuel safety. Summary of the Invention

[0003] The technical problem this invention aims to solve is the difficulty in accurately predicting the initiation and propagation of cracks in the coating layer of TRISO particles under high temperature, strong irradiation, and the coupling effect of fission gases. It provides a simulation method for TRISO particle coating layer crack propagation based on a coupling of solid mechanics and fracture phase-field methods, comprehensively considering the influence of geometric defects, material irradiation behavior, and fission gas pressure. This method establishes a geometric model of the multi-layered coating structure and combines it with the temperature field, irradiation effect, and fission gas release behavior during fuel service to realistically depict the complex service environment of the coating layer. Simultaneously, it leverages the advantages of the fracture phase-field method in describing continuous crack evolution, combined with the accuracy of the solid mechanics model in stress-strain analysis, to uniformly handle the crack initiation, propagation, and interaction processes, providing a reliable theoretical basis and numerical tool for the structural integrity assessment and optimization design of TRISO particles.

[0004] To achieve the above objectives, the present invention adopts the following technical solution: A method for simulating crack propagation in TRISO particle-coated layers, comprising the following steps: Step 1: Establishment of a two-dimensional rotational axisymmetric geometric model of the three outer coating layers of TRISO particles; Step 2: Setting the physical properties of the three-layer coating materials, including the density ρ of pyrolytic carbon and silicon carbide, Young's modulus E, Poisson's ratio υ, crack size parameter l0, and critical energy release rate G. c ; Step 3: Determine the service environment parameters of TRISO particles. The service environment parameters include core thermal power, service life, fuel temperature, fission rate, and fast neutron flux. Calculate the burnup depth and fast neutron flux of the fuel at a certain moment based on the service environment parameters to obtain the material properties under the corresponding environmental conditions and calculate the fission gas release. Step 4: Establish irradiation behavior and fission gas release models to simulate the possible irradiation deformation, irradiation creep, and irradiation swelling of nuclear fuel in the reactor due to neutron irradiation, as well as the release of fission gas generated by uranium atom fission and its effect on the cladding layer. Step 5: Establishment and verification of the fracture phase field method. Establish the governing equations and solve them using the finite element method. The governing equations include the solid mechanics governing equations and the energy conservation equations. Verify the reliability of the method using typical fracture mechanics problems. Step 6: Introduce the geometric imperfections to be studied into the two-dimensional rotational axisymmetric geometric model and solve the governing equations; Step 7: Calculation complete, output the calculation results, and evaluate the crack propagation.

[0005] In one possible embodiment, the two-dimensional rotationally axisymmetric geometric model in step 2 includes a buffer layer of a certain thickness, an inner dense pyrolytic carbon layer, a silicon carbide layer, and an outer dense pyrolytic carbon layer. To facilitate the introduction of anisotropic material properties and the calculation of radial and tangential stress and strain, while reducing unnecessary computational costs, local spherical coordinates (r, θ, φ) and global two-dimensional axisymmetric coordinates (ρ) are established. s The transformation relationship of z, φ) is shown in equation (1): (1) In the formula, and These are the radial and tangential basis vectors in spherical coordinates. and These are the basis vectors for the polar radius and height directions in a two-dimensional axisymmetric coordinate system.

[0006] In one possible embodiment, the service environment parameters of the fuel particles in step 3 include core thermal power, service life, fuel temperature, fission rate, and fast neutron flux rate, thereby calculating the burn-up depth and fast neutron flux of the fuel at a certain moment. This is used to obtain the material properties under the corresponding environmental conditions and calculate the fission gas release. The burn-up depth is calculated by the following formula. (2) In the formula: Fuel consumption depth For fission rate, Initial uranium atom density; The fission rate and the initial uranium atom density are determined by equations (3)-(4): (3) (4) In equations (3)-(4): q represents the fuel enrichment degree. For effective fission cross section, For uranium atom density, Where is the fast neutron flux, D is the theoretical density ratio of the fuel (i.e., the ratio of actual density to theoretical density), and TD is the initial theoretical density of the fuel. Let Avogadro's constant be 1. is the molar mass of uranium dioxide; Fast neutron flux is calculated by the following formula: (5) In the formula: Let be the fast neutron fluence rate, and c be the conversion factor between line power and neutron fluence rate. The linear power density of the fuel.

[0007] In one possible embodiment, the irradiation behavior and fission gas release model in step 4 includes the irradiation deformation, irradiation creep, and irradiation swelling that may occur in nuclear fuel within the reactor due to neutron irradiation, as well as the effect of fission gas generated by uranium atom fission on the cladding layer. The irradiation deformation is calculated using the following formula: (6) (7) In equations (6)-(7), and 1 / s represents the radial and tangential irradiation strain rates of dense pyrolytic carbon, respectively. This parameter describes the rate of volume change of pyrolytic carbon materials in the radial and tangential directions under neutron irradiation. Since pyrolytic carbon is anisotropic, its irradiation deformation behavior varies in different directions. This formula is obtained by fitting experimental data of neutron irradiation and can be expressed as a function of fast neutron flux, reflecting the evolution rate of the material's microstructure with cumulative irradiation. The irradiation creep of a material is calculated by the following formula: (8) (9) (10) In equations (8)-(10): The irradiation creep rate of dense pyrolytic carbon represents the rate at which pyrolytic carbon materials undergo irradiation-induced creep under external stress in an irradiated environment. In other words, it is the strain per unit time under constant stress and neutron irradiation. This parameter reflects the viscoplastic response characteristics of the material under irradiation. K is the creep constant, and σ... i Let υ be the i-th principal stress (i = 1, 2, 3) in the local coordinate system. c denoted as Poisson's ratio of pyrolyzed carbon under irradiation conditions, and T as temperature; The irradiation swelling of silicon carbide is expressed by the following formula. (11) (12) In equations (11)-(12): σ represents the irradiation creep rate of silicon carbide, characterizing the creep deformation rate of silicon carbide materials under neutron irradiation. e This is the effective stress, i.e., the Mies stress; Φ0 represents the volumetric strain of silicon carbide caused by irradiation swelling. This parameter is used to characterize the volume change of silicon carbide materials under neutron irradiation due to microstructural evolution such as point defect accumulation and lattice expansion. Φ0 is the neutron flux parameter. Due to nuclear fission and the special microstructure of uranium dioxide, the fission gas molecules and atoms are bound within the UO2 grains after being generated. Only after the bubbles aggregate, expand, and connect to a certain extent can they break free from the constraints and be released into free space, thereby generating gas pressure on the inner surface of the dense pyrolytic carbon layer. Before and after the concentration reaches saturation, the grain boundary gas concentration is calculated by the model proposed by Cui. This model is based on the diffusion equation of equivalent spherical grains, the incomplete sink boundary conditions of grain boundaries, and the gas mass conservation relationship. The analytical solutions of the grain boundary gas concentration before and after the concentration reaches saturation are obtained by Laplace transform and inverse transform, as shown in equations (13)-(17): (13) (14) (15) (16) (17) In equations (13)-(17): N BS The concentration of the grain boundary gas before saturation is given, where 'a' is the average grain radius and 'Y' is the average number of fission gas atoms and molecules produced in a single fission event. λ is the fission rate, b' is the resolubility of gas at the grain boundary, and λ is the resolubility of gas at the grain boundary. b Let N be the full width at half maximum (FWHM) of the resoluble gas distribution function, D' be the effective diffusion coefficient, and N' be the effective diffusion coefficient. AS The concentration of the grain boundary gas after the concentration reaches saturation, ts The time it takes for the gas concentration to reach saturation; After calculating the fission gas release, the gap width between the buffer layer and the inner dense pyrolytic carbon layer, and the average temperature of the air gap, the fission gas pressure acting on the inner surface of the inner dense pyrolytic carbon layer can be obtained by averaging the pore size of the buffer layer material and the temperature and volume of the gas in the gap between the buffer layer and the inner dense pyrolytic carbon layer using the ideal gas law, as shown in equation (18). (18) In the formula: and This refers to the amount of substance of the fission gas and CO. and This serves as a buffer layer and a space within the air gap to hold the gas. and This refers to the temperature of the gas in the buffer layer and air gap.

[0008] In one possible embodiment, the solid mechanics governing equation in step 5 is shown below. (19) In the formula: For Cauchy stress tensor, It is a volume force; The volumetric strain is calculated based on the small strain assumption, as shown in the following equation: (20) In the formula: Let be the initial volume of the domain. The volume of the domain after strain occurs. , , These are the strain tensors in the x, y, and z directions, respectively; Elastic strain is calculated using equations (21)-(23). (twenty one) (twenty two) (twenty three) In equations (21)-(23): For total strain, and For elastic and inelastic strain, It is a displacement vector. The initial strain is determined by the initial conditions. External strain, which is related to boundary conditions. and For thermal strain and creep strain; In one possible embodiment, step 5 involves constructing an energy conservation equation. First, an expression for the total potential energy is considered. After considering energy dissipation caused by fracture, when a crack Γ appears in a domain Ω made of a certain material, the total potential energy ψ... t As shown in the following formula: (twenty four) In the formula: ψ ε Elastic strain energy density, ε strain tensor, G c Critical energy release rate, S crack surface area, b volume force tensor, u displacement tensor, f surface force tensor; To address the difficulty of integrating the fracture energy term, a phase field p is defined, and a crack surface density function γ(p, p) to approximate the fracture energy term in equation (24) by representing the crack region Г: (25) In the formula: l0 is the crack size parameter; The field variable p=1 indicates that the material at this point is in a completely fractured state, while p=0 indicates that the point is in an undamaged state. Since crack propagation is mainly driven by tensile load, in order to avoid the non-physical contribution of the compressive state to crack propagation, it is necessary to decompose the elastic strain energy density into tensile and compressive parts, and only the tensile part participates in the phase field evolution, while the compressive part does not contribute. Based on the positive and negative decomposition of the strain tensor, the elastic strain energy density decomposition is shown in equations (26)-(28): (26) (27) (28) In equations (26)-(28): ψ ε0 For elastic strain energy density without considering fracture, ψ ε0+ and ψ ε0- Let λ be the tensile and compressive components of the elastic strain energy density, λ and μ be the first and second Lamé constants, and tr[A] be the trace of matrix A. <x> ± This represents the operation (x ± |x|) / 2 on x, where ε + and ε - ψ represents the positive and negative components of the strain tensor. ε (ε,p) is the elastic strain energy density function considering fracture; The total Lagrange potential energy is the sum of the fracture energy term, the elastic strain energy term, the volumetric internal energy term, and the energy term caused by the external load: (29) In the formula: ρ is the material density, and t is time; By variational method, the strong formal equations of solid mechanics and phase-field models are obtained: (30) In the formula: σ ij Stress components, x-position components; If the above phase field model is used directly, when the load decreases or is removed, the phase field value at the location where the crack has formed will decrease directly, indicating that the crack has recovered or even disappeared, which contradicts the irreversibility of cracks in reality. To ensure the monotonically increasing phase field value, a historical strain energy field H(x,t) is defined to record the tensile component of the maximum elastic strain energy density in the historical solution time step: (31) Replacing the tensile component of the elastic strain energy in equation (30) with the historical strain energy H, the strong form equation of the phase field model is shown in equation (32): (32) In the solid mechanics equations, as shown in equation (33), the Cauchy stress σ ij It is also decomposed accordingly, and the Cauchy stress σ tensor is expressed as equation (34): (33) (34) In equation (34): I is the identity matrix; The weak form equations of solid mechanics and the phase field model are as shown in equations (35)-(36): (35) (36) In equations (35)-(36): δx represents a trial function of the variable x; By solving classical fracture mechanics problems and comparing the results with those of other experimental and numerical methods for crack propagation, the reliability of the model and the differences are analyzed, and the reliability of the proposed method is verified.

[0009] In one possible embodiment, step 6 involves introducing and calculating the geometric defects to be studied in each coating layer, including first introducing defects common to particles during production and transportation, such as gaps in dense pyrolytic carbon layers and residual pores in silicon carbide layers.

[0010] In one possible embodiment, after the calculation in step 7 is completed, the calculation results are output, stress distribution and phase field evolution images are plotted, and crack propagation is evaluated.

[0011] Compared with the prior art, the present invention has at least the following beneficial effects: This invention provides a method for simulating crack propagation in TRISO particle coatings. It combines the advantages of solid mechanics models in accurately analyzing stress-strain response with the advantages of fracture phase-field methods in continuously characterizing crack initiation and propagation processes. Addressing the cracking problem of TRISO particle coatings under high temperature, strong irradiation, and the coupled effects of fission gas, the method couples material irradiation behavior, fission gas pressure, and crack propagation response. Irradiation deformation and fission gas release models, solid mechanics control equations, and fracture phase-field control equations are established separately. The coupled model is solved using the finite element method, which largely ensures the completeness of the crack propagation analysis, significantly improves model refinement, and guarantees the accuracy of the calculation results.

[0012] Furthermore, an irradiation behavior and fission gas release model is established. By considering the irradiation deformation, irradiation creep, and irradiation swelling behavior of pyrolytic carbon materials and silicon carbide layers, and combining the generation, diffusion, accumulation, and release processes of fission gas inside fuel particles, an accurate mapping of the service environment to the cladding load response is achieved. This ensures the accuracy of the calculation of internal pressure load and material inelastic evolution in crack propagation analysis.

[0013] Furthermore, using solid mechanics theory, by establishing stress balance equations, geometric relationships, and constitutive relationships, the stress-strain response of each coating layer of TRISO particles under complex service conditions is calculated, ensuring the accuracy of the analysis of the mechanical behavior of the coating layers.

[0014] Furthermore, using the fracture phase field method to handle the crack initiation and propagation process in the coating layer can avoid the limitations of pre-setting crack paths in traditional fracture mechanics methods and fully preserve the continuous evolution characteristics of crack propagation. At the same time, it combines the historical strain energy field to handle the irreversibility of crack evolution and considers the combined effects of geometric defects, irradiation effects and fission gas pressure to ensure the reliability and rationality of TRISO particle coating layer failure assessment.

[0015] Furthermore, this invention can fully leverage the advantages of solid mechanics theory and fracture phase-field method to ensure the completeness of the analysis of crack propagation in TRISO particle coatings. By comprehensively considering factors such as geometric defects, material irradiation behavior, and fission gas pressure, it can achieve a refined simulation of crack initiation, propagation, and failure processes, resulting in more accurate calculation results and ensuring the accuracy of the calculation results. Attached Figure Description

[0016] Figure 1 This is a flowchart of the present invention; Figure 2 This is the geometric model of the present invention; Figure 3 This is a diagram showing the reliability verification results of the present invention; Detailed Implementation To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and specific embodiments. However, it should be understood that the following embodiments are for illustrative purposes only and are not intended to limit the scope of protection of this invention. Furthermore, in the following description, well-known structures, conventional calculation methods, and general numerical solution processes will not be described in detail to avoid unnecessary confusion regarding the core technical content of this invention.

[0017] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments.

[0018] Example: Simulation of crack propagation in TRISO particle-coated layers This embodiment uses a TRISO particle as the research object and applies the technical solution of this invention to simulate the crack propagation of the cladding layer. This fuel particle consists of a fuel core, a buffer layer, an inner dense pyrolytic carbon layer, a silicon carbide layer, and an outer dense pyrolytic carbon layer. During reactor operation, it is subjected to a long-term service environment characterized by high temperature, strong neutron irradiation, and continuous generation of fission products. Due to the combined effects of material irradiation behavior, fission gas release, and local geometric defects in the cladding layer, stress concentration may occur within the particle's cladding layer, inducing crack initiation and propagation. This calculation simulates the stress evolution and crack propagation process of the TRISO particle's cladding layer under typical service conditions, studying its damage distribution characteristics and structural integrity.

[0019] Step 1: Establishment of the geometric model of the particle outer coating layer According to the mainstream TRISO particle model, in this embodiment, the UO2 core diameter is 425 μm, the buffer layer thickness is 100 μm, and the inner surface radius of the inner dense pyrolytic carbon layer is 312.5 μm. The thicknesses of the inner dense pyrolytic carbon layer, silicon carbide layer, and outer dense pyrolytic carbon layer are 35 μm, 40 μm, and 35 μm, respectively. Figure 2 As shown. Due to the requirements of the phase-field method for mesh refinement, in order to reduce unnecessary computational costs and facilitate the introduction of anisotropic material properties and the calculation of radial and tangential stress and strain, a two-dimensional rotational axisymmetric geometric model of the outer triple cladding layer is established. This requires establishing local spherical coordinates (r, θ, φ) and global two-dimensional axisymmetric coordinates (ρ). s The transformation relationship of z, φ) is as shown in equation (1): (1) Step 2: Setting the physical properties of each layer of materials To fully account for the impact of the service environment on the material properties, in this embodiment, the density of the pyrolytic carbon material is 1900 kg / m³. 3 With a Poisson's ratio of 0.21 and Young's modulus being anisotropic, it can be calculated using equations (2)-(3): (2) (3) In equations (2)-(3): and The radial and tangential Young's moduli of pyrolytic carbon. The Pa is 25.5 GPa. Where is density, and BAF is the anisotropy orientation factor of the grains, i.e., sphericity, taken as 1.03. The grain diameter is 3 nm. Where T is the fast neutron flux and T is the temperature; Silicon carbide has a density of 3210 kg / m³ 3 Poisson's ratio is 0.21, and Young's modulus is modeled using the Snead model. (4) For pyrolytic carbon and silicon carbide, the crack length scale l0 is taken as 3 μm and 1 μm respectively, and the critical energy release rate G is... c It can be determined from l0 by the following formula: (5) In the formula: The critical stresses for the materials are 200 MPa for pyrolytic carbon and 667 MPa for silicon carbide.

[0020] Step 3: Determine the service environment parameters of fuel pellets To simulate the behavior of fuel under real-world in-reactor conditions, it is necessary to accurately calculate the core thermal power, service life, fuel temperature, fission rate, and fast neutron flux. Using the burn-up depth and fast neutron flux at a given moment, the material properties under corresponding environmental conditions are obtained, and the release of fission gases is calculated. The burn-up depth is calculated using the following formula: (6) In the formula: Fuel consumption depth For fission rate, Initial uranium atom density; The fission rate and the initial uranium atom density are determined by equations (7)-(8): (7) (8) In equations (7)-(8): q represents the fuel enrichment degree of 5%. For effective fission cross section, For uranium atom density, Where is the fast neutron flux, D is the theoretical density ratio of the fuel (i.e., the ratio of actual density to theoretical density), and TD is the initial theoretical density of the fuel (10970 kg / m³). 3 , Let Avogadro's constant be 1. is the molar mass of uranium dioxide; Fast neutron flux is calculated by the following formula: (9) In the formula: Let be the fast neutron fluence rate, and c be the conversion factor between line power and neutron fluence rate, taken as 3 × 10⁻⁶. 13 (1 / (m 2 •s) / (W / m)), The linear power density of the fuel.

[0021] To fully simulate the real behavior of each layer of material under reactor service conditions, the irradiation behavior and fission gas release model in this embodiment includes the irradiation deformation, irradiation creep, and irradiation swelling that nuclear fuel may experience due to neutron irradiation within the reactor, as well as the effect of fission gas generated by uranium atom fission on the cladding layer. The irradiation deformation is calculated using the following formula: (10) (11) In equations (10)-(11), and 1 / s represents the radial and tangential irradiation strain rates of dense pyrolytic carbon, respectively. This parameter describes the rate of volume change of pyrolytic carbon materials in the radial and tangential directions under neutron irradiation. Since pyrolytic carbon is anisotropic, its irradiation deformation behavior varies in different directions. This formula is obtained by fitting experimental data of neutron irradiation and can be expressed as a function of fast neutron flux, reflecting the evolution rate of the material's microstructure with cumulative irradiation. The irradiation creep of a material is calculated by the following formula: (12) (13) (14) In equations (12)-(14): The irradiation creep rate of dense pyrolytic carbon represents the rate at which pyrolytic carbon materials undergo irradiation-induced creep under external stress in an irradiated environment. In other words, it is the strain per unit time under constant stress and neutron irradiation. This parameter reflects the viscoplastic response characteristics of the material under irradiation. K is the creep constant, and σ... i Let υ be the i-th principal stress (i = 1, 2, 3) in the local coordinate system. c The Poisson's ratio for pyrolyzed carbon under irradiation conditions is 0.5, and T is the temperature. The irradiation swelling of silicon carbide is expressed by the following formula. (15) (16) In equations (11)-(12): σ represents the irradiation creep rate of silicon carbide, characterizing the creep deformation rate of silicon carbide materials under neutron irradiation. e This is the effective stress, i.e., the Mies stress; The bulk strain of silicon carbide caused by irradiation swelling is denoted as Φ0. This parameter characterizes the volume change of silicon carbide materials under neutron irradiation due to microstructural evolution such as point defect accumulation and lattice expansion. Φ0 is the neutron flux parameter, taken as 0.3396 x 10². 5 n / m². Due to nuclear fission and the special microstructure of uranium dioxide, the fission gas molecules and atoms are bound within the UO2 grains after being generated. Only after the bubbles aggregate, expand, and connect to a certain extent can they break free from the constraints and be released into free space, thereby generating gas pressure on the inner surface of the dense pyrolytic carbon layer. Before and after the concentration reaches saturation, the grain boundary gas concentration is calculated by the model proposed by Cui. This model is based on the diffusion equation of equivalent spherical grains, the incomplete sink boundary conditions of grain boundaries, and the gas mass conservation relationship. The analytical solutions of the grain boundary gas concentration before and after the concentration reaches saturation are obtained by Laplace transform and inverse transform, as shown in equations (13)-(17): (17) (18) (19) (20) (twenty one) In equations (13)-(17): N BS The concentration of the grain boundary gas before saturation is given, where 'a' is the average grain radius and 'Y' is the average number of fission gas atoms and molecules produced in a single fission event. λ is the fission rate, b' is the resolubility of gas at the grain boundary, and λ is the resolubility of gas at the grain boundary. b Let N be the full width at half maximum (FWHM) of the resoluble gas distribution function, D' be the effective diffusion coefficient, and N' be the effective diffusion coefficient. AS The concentration of the grain boundary gas after the concentration reaches saturation, t s This refers to the time it takes for the gas concentration to reach saturation; in this embodiment, the saturation concentration is calculated using the following formula: (twenty two) In the formula: Let be the radius of curvature of the lenticular grain boundary bubble when the grain boundary reaches saturation, taken as 0.5 μm. The angle between the two spherical surfaces of the lenticular bubble is 50°. , When grain boundaries reach saturation, the area fraction of the grain boundary surface occupied by intergranular bubbles is 0.5. The surface energy of the bubble is 0.626 J / m. 2 , , , These are Boltzmann's constant, temperature, and external pressure, respectively.

[0022] The saturation time is calculated using the following formula: (twenty three) To implement this in the program, this example solves the highly nonlinear equation that cannot be solved analytically using the Newton iteration method constructed by equations (24)-(26); (twenty four) (25) (26) After calculating the fission gas release, the gap width between the buffer layer and the inner dense pyrolytic carbon layer, and the average temperature of the air gap, the fission gas pressure acting on the inner surface of the inner dense pyrolytic carbon layer can be obtained by averaging the pore size of the buffer layer material and the temperature and volume of the gas in the gap between the buffer layer and the inner dense pyrolytic carbon layer using the ideal gas law, as shown in equation (18). (27) In the formula: and This refers to the amount of substance of the fission gas and CO. and This serves as a buffer layer and a space within the air gap to hold the gas. and This refers to the temperature of the gas in the buffer layer and air gap.

[0023] In this embodiment, the solid mechanics-phase-field model uses the solid mechanics governing equations and energy conservation equations. The solid mechanics governing equations are shown below. (28) In the formula: For Cauchy stress tensor, It is a volume force; The volumetric strain is calculated based on the small strain assumption, as shown in the following equation: (29) In the formula: Let be the initial volume of the domain. The volume of the domain after strain occurs. , , These are the strain tensors in the x, y, and z directions, respectively; Elastic strain is calculated using equations (21)-(23). (30) (31) (32) In equations (21)-(23): For total strain, and For elastic and inelastic strain, It is a displacement vector. The initial strain is determined by the initial conditions. External strain, which is related to boundary conditions. and For thermal strain and creep strain; To construct the energy conservation equation, we first consider the expression for the total potential energy. After considering the energy dissipation caused by fracture, when a crack Γ appears in a domain Ω made of a certain material, the total potential energy ψ... t As shown in the following formula: (33) In the formula: ψ ε Elastic strain energy density, ε strain tensor, G c Critical energy release rate, S crack surface area, b volume force tensor, u displacement tensor, f surface force tensor; To address the difficulty of integrating the fracture energy term, a phase field p is defined, and a crack surface density function γ(p, p) to approximate the fracture energy term in equation (24) by representing the crack region Г: (34) In the formula: l0 is the crack size parameter; The field variable p=1 indicates that the material at this point is in a completely fractured state, while p=0 indicates that the point is in an undamaged state. Since crack propagation is mainly driven by tensile load, in order to avoid the non-physical contribution of the compressive state to crack propagation, it is necessary to decompose the elastic strain energy density into tensile and compressive parts, and only the tensile part participates in the phase field evolution, while the compressive part does not contribute. Based on the positive and negative decomposition of the strain tensor, the elastic strain energy density decomposition is shown in equations (26)-(28): (35) (36) (37) In equations (26)-(28): ψ ε0 For elastic strain energy density without considering fracture, ψ ε0+ and ψ ε0- Let λ be the tensile and compressive components of the elastic strain energy density, λ and μ be the first and second Lamé constants, and tr[A] be the trace of matrix A. <x> ± This represents the operation (x ± |x|) / 2 on x, where ε + and ε - ψ represents the positive and negative components of the strain tensor. ε (ε,p) is the elastic strain energy density function considering fracture; The total Lagrange potential energy is the sum of the fracture energy term, the elastic strain energy term, the volumetric internal energy term, and the energy term caused by the external load: (38) In the formula: ρ is the material density, and t is time; By variational method, the strong formal equations of solid mechanics and phase-field models are obtained: (39) In the formula: σ ij Stress components, x-position components; If the above phase field model is used directly, when the load decreases or is removed, the phase field value at the location where the crack has formed will decrease directly, indicating that the crack has recovered or even disappeared, which contradicts the irreversibility of cracks in reality. To ensure the monotonically increasing phase field value, a historical strain energy field H(x,t) is defined to record the tensile component of the maximum elastic strain energy density in the historical solution time step: (40) Replacing the tensile component of the elastic strain energy in equation (30) with the historical strain energy H, the strong form equation of the phase field model is shown in equation (32): (41) In the solid mechanics equations, as shown in equation (33), the Cauchy stress σ ij It is also decomposed accordingly, and the Cauchy stress σ tensor is expressed as equation (34): (42) (43) In equation (34): I is the identity matrix; The weak form equations of solid mechanics and the phase field model are as shown in equations (35)-(36): (44) (45) In equations (35)-(36): δx represents a trial function of the variable x; To verify the reliability of the method, it should be compared with the results of classical fracture mechanics problems and other experimental and numerical methods studying crack propagation to analyze the model's reliability and differences. This embodiment uses the Kalthoff fracture experiment to model a notched plate, and applies inward-normal velocity boundary conditions at local locations on the side: (46) Transient calculations were performed from 0 to 90 μs. The crack propagation direction calculated by this model was approximately 65° from the notch. Figure 3 As shown, the crack morphology and propagation rate are consistent with the results of other numerical methods in fracture mechanics (such as the extended finite element method and the phantom node method), proving that the simulation of crack propagation by this model is accurate.

[0024] Step 6: Introduce the geometric defects to be studied and calculate First, defects commonly encountered in particle production and transportation are introduced into the coating layer. In this embodiment, initial cracks of different depths are introduced into the inner surface of the dense pyrolytic carbon layer to analyze its propagation behavior. Based on the outer coating layer model, the governing equations are discretized under the conditions of irradiation behavior and fission gas pressure. The weak form equations of solid mechanics and phase field model are solved using the finite element method to calculate the stress and phase field distribution of the coating layer at the current time t. As a continuous scalar between 0 and 1, the phase field can directly reflect the damage and crack propagation behavior in the coating layer.

[0025] Step 7: Calculation complete, output the calculation result. The calculation results are post-processed to obtain key mechanical parameters such as the proportion of interlayer separation region and the maximum shear stress of silicon carbide layer, and to draw stress distribution and phase field evolution images to evaluate crack propagation.

[0026] This embodiment calculates a single-particle power of 0.05 W and a service life of 1.5 × 10⁻⁶ W. 7 The results obtained under condition s are shown in Table 1. The calculations show that shallow initial cracks are insufficient to rapidly propagate to the silicon carbide layer interface and cause interface damage under service conditions. At this point, no significant interlayer separation has occurred, and the silicon carbide layer remains primarily under compressive stress, indicating that the overall cladding layer maintains good structural integrity, and the initial defects have a limited weakening effect on the safety barrier. When the initial crack depth exceeds 5.9 μm, the crack propagation rate significantly increases and begins to develop towards the interlayer interface, leading to interface separation. This indicates that the initial defects have reached the critical range for rapid evolution of cladding layer damage. Simultaneously, the large tensile stress in the silicon carbide layer indicates that cracks in the inner dense pyrolytic carbon layer have significantly affected the outer load-bearing structure, threatening the structural integrity of the cladding layer and significantly increasing the risk of failure.

[0027] Table 1. Calculated interlayer separation region ratio and maximum shear stress of silicon carbide layer under different initial crack conditions. In summary, the method established in this invention comprehensively utilizes the advantages of solid mechanics models in accurately analyzing stress-strain response and the fracture phase-field method in continuously characterizing crack initiation and propagation processes. For the crack propagation problem in TRISO particle coatings, it correlates material irradiation behavior, fission gas release, and geometric defects with the coating stress response and crack evolution process. Irradiation behavior and fission gas release models, solid mechanics governing equations, and fracture phase-field governing equations are established separately. The coupled model is solved using the finite element method. Factors such as irradiation deformation of pyrolytic carbon materials, irradiation creep, silicon carbide layer irradiation swelling, and the irreversibility of crack propagation are considered, which largely ensures the completeness of the crack propagation problem analysis, significantly improves model refinement, and ensures the accuracy of the calculation results. Furthermore, the calculation results of this embodiment also verify that this invention can accurately simulate the crack initiation and propagation process in TRISO particle coatings, ensuring the accuracy of the calculation results.< / x> < / x>

Claims

1. A method for simulating crack propagation in TRISO particle-coated layers, characterized in that: Includes the following steps: Step 1: Establishment of a two-dimensional rotationally axisymmetric geometric model of the three outer coating layers of TRISO particles; Step 2: Setting the physical properties of the three coating layer materials, including the density of pyrolytic carbon and silicon carbide. ρ Young's modulus E Poisson's ratio υ Crack size parameters l 0, Critical Energy Release Rate G c ; Step 3: Determine the service environment parameters of TRISO particles. The service environment parameters include core thermal power, service life, fuel temperature, fission rate, and fast neutron flux. Calculate the burnup depth and fast neutron flux of the fuel at a certain moment based on the service environment parameters to obtain the material properties under the corresponding environmental conditions and calculate the fission gas release. Step 4: Establish irradiation behavior and fission gas release models to simulate the possible irradiation deformation, irradiation creep, and irradiation swelling of nuclear fuel in the reactor due to neutron irradiation, as well as the release of fission gas generated by uranium atom fission and its effect on the cladding layer. Step 5: Establishment and verification of the fracture phase field method. Establish the governing equations and solve them using the finite element method. The governing equations include the solid mechanics governing equations and the energy conservation equations. Verify the reliability of the method using typical fracture mechanics problems. Step 6: Introduce the geometric imperfections to be studied into the two-dimensional rotational axisymmetric geometric model and solve the governing equations; Step 7: Calculation complete, output the calculation results, and evaluate the crack propagation.

2. The TRISO particle-coated layer crack propagation simulation method according to claim 1, characterized in that, The two-dimensional rotational axisymmetric geometric model in step 1 includes a buffer layer of a certain thickness, an inner dense pyrolytic carbon layer, a silicon carbide layer, and an outer dense pyrolytic carbon layer. To facilitate the introduction of anisotropic material properties and the calculation of radial and tangential stress and strain, while reducing unnecessary computational costs, a local spherical coordinate system is established. r, θ, φ ) and global two-dimensional axisymmetric coordinates ( ρ s , z, φ The transformation relationship is shown in equation (1): (1) In the formula, and These are the radial and tangential basis vectors in spherical coordinates. and These are the basis vectors for the polar radius and height directions in a two-dimensional axisymmetric coordinate system.

3. The TRISO particle-coated layer crack propagation simulation method according to claim 1, characterized in that, The service environment parameters mentioned in step 3 include core thermal power, service life, fuel temperature, fission rate, and fast neutron flux. Based on these parameters, the burnup depth and fast neutron flux at a given moment are calculated to obtain material properties under corresponding environmental conditions and to calculate fission gas release. The burnup depth is calculated using the following formula: (2) In the formula: Fuel consumption depth For fission rate, Initial uranium atom density; The fission rate and the initial uranium atom density are determined by equations (3)-(4): (3) (4) In equations (3)-(4): q For fuel enrichment, For effective fission cross section, For uranium atom density, Where is the fast neutron flux, D is the theoretical density ratio of the fuel (i.e., the ratio of actual density to theoretical density), and TD is the initial theoretical density of the fuel. Let Avogadro's constant be 1. is the molar mass of uranium dioxide; Fast neutron flux is calculated by the following formula: (5) In the formula: Let be the fast neutron fluence rate, and c be the conversion factor between line power and neutron fluence rate. The linear power density of the fuel.

4. The TRISO particle-coated layer crack propagation simulation method according to claim 1, characterized in that, The irradiation behavior in step 4 includes the irradiation deformation, irradiation creep, and irradiation swelling that may occur in the nuclear fuel due to neutron irradiation within the reactor, as well as the effect of fission gas generated by uranium atom fission on the cladding layer. The irradiation deformation is calculated using the formulas shown in (6)-(7): (6) (7) In equations (6)-(7), and The radial and tangential irradiation strain rates (1 / s) of dense pyrolytic carbon are respectively, and the irradiation creep is calculated using the formulas shown in (8)-(10): (8) (9) (10) In equations (8)-(10): The irradiation creep rate of dense pyrolytic carbon. K The creep constant is σ i The first in the local coordinate system i Principal stress ( i = 1,2,3), υ c The Poisson's ratio of pyrolyzed carbon under irradiation conditions. T For temperature; The irradiation swelling is calculated using the formulas shown in (11)-(12). (11) (12) In equations (11)-(12): The irradiation creep rate of silicon carbide. σ e For effective stress, This refers to the volumetric strain caused by irradiation swelling of silicon carbide. Ф 0 represents the neutron flux parameter.

5. The TRISO fuel particle coating crack propagation simulation method according to claim 1, characterized in that, The fission gas model in step 4 includes the release of fission gas and the gas pressure generated on the coating layer. Due to the occurrence of nuclear fission and the special microstructure of uranium dioxide material, the fission gas is generated into atoms and then bound within the UO2 grains. After the bubbles aggregate, expand, and connect to a certain extent, they can break through the constraints and be released into free space, thereby generating gas pressure on the inner surface of the dense pyrolytic carbon layer. Before and after the concentration reaches saturation, the grain boundary gas concentration is calculated by the model proposed by Cui. This model is based on the equivalent spherical grain diffusion equation, the incomplete sink boundary condition of the grain boundary, and the gas mass conservation relationship. The analytical solutions of the grain boundary gas concentration before and after the concentration reaches saturation are obtained by Laplace transform and inverse transform, as shown in equations (13)-(17): (13) (14) (15) (16) (17) In equations (13)-(17): N BS The concentration of the grain boundary gas before it reaches saturation. a The average grain radius is... Y The average number of fission gas atoms and molecules produced in a single fission event. For fission rate, b’ The resolubility of gas at grain boundaries. λ b Let be the full width at half maximum (FWHM) of the distribution function of the resoluble gas. D’ For the effective diffusion coefficient, N AS The concentration of the grain boundary gas after the concentration reaches saturation, t s The time it takes for the gas concentration to reach saturation; After calculating the fission gas release, the gap width between the buffer layer and the inner dense pyrolytic carbon layer, and the average temperature of the air gap, the fission gas pressure acting on the inner surface of the inner dense pyrolytic carbon layer can be obtained by averaging the pore size of the buffer layer material and the temperature and volume of the gas in the gap between the buffer layer and the inner dense pyrolytic carbon layer using the ideal gas law, as shown in equation (18). (18) In the formula: and This refers to the amount of substance of the fission gas and CO. and This serves as a buffer layer and a space within the air gap to hold the gas. and This refers to the temperature of the gas in the buffer layer and air gap.

6. The TRISO fuel particle coating crack propagation simulation method according to claim 1, characterized in that, The solid mechanics governing equation in step 5 is shown in equation (19). (19) In the formula: For Cauchy stress tensor, It is a volume force; The volumetric strain is calculated based on the small strain assumption, as shown in the following equation: (20) In the formula: Let be the initial volume of the domain. The volume of the domain after strain occurs. , , These are the strain tensors in the x, y, and z directions, respectively; Elastic strain is calculated using equations (21)-(23). (21) (22) (23) In equations (21)-(23): For total strain, and For elastic and inelastic strain, It is a displacement vector. The initial strain is determined by the initial conditions. External strain, which is related to boundary conditions. and For thermal strain and creep strain.

7. The TRISO fuel particle coating crack propagation simulation method according to claim 1, characterized in that, In step 5, to construct the energy conservation equation, an expression for the total potential energy is first established. Considering the energy dissipation caused by fracture, when a crack Г appears in a domain Ω made of a certain material, the total potential energy is... ψ t As shown in the following formula: (24) In the formula: ψ ε Elastic strain energy density, ε Strain tensor G c Critical energy release rate S Crack surface area b Volume force tensor u Displacement tensor f Surface external force tensor; Define phase field p and the introduction of crack surface density function γ(p, p) To approximate the crack region Г Thus, the fracture energy term in equation (24) can be approximated: (25) In the formula: l 0 represents the crack size parameter; field variables p=1 This indicates that the material at that point is in a state of complete fracture. p=0 This indicates that the point is in an undamaged state. Since crack propagation is mainly driven by tensile load, in order to avoid the non-physical contribution of the compressive state to crack propagation, the elastic strain energy density needs to be decomposed into tensile and compressive parts, and only the tensile part should participate in the phase field evolution, while the compressive part should not contribute. Based on the positive and negative decomposition of the strain tensor, the elastic strain energy density decomposition is shown in equations (26)-(28): (26) (27) (28) In equations (26)-(28): ψ ε0 For elastic strain energy density without considering fracture, ψ ε0+ and ψ ε0- These are the tensile and compressive components of the elastic strain energy density. λ and μ For the first and second Lamé constants, tr[A] For matrix A traces, <x>< / x> ± Indicates to x conduct (x±|x|) / 2 Operations, ε + and ε - These are the positive and negative components of the strain tensor. ψ ε ( ε,p () represents the elastic strain energy density function considering the fracture. The total Lagrange potential energy is the sum of the fracture energy term, the elastic strain energy term, the volumetric internal energy term, and the energy term caused by the external load: (29) In the formula: ρ For material density, t For time; By variational method, the strong formal equations of solid mechanics and phase-field models are obtained: (30) In the formula: σ ij Stress components, x Positional components; To ensure the monotonic increase of the phase field value, the history strain energy field is defined. H(x,t) The tensile component of the maximum elastic strain energy density recorded in the historical solution time step: (31) Replace the tensile component of the elastic strain energy in equation (30) with the historical strain energy. H The strong formal equation of the phase-field model is shown in equation (32): (32) In the solid mechanics equations, as shown in equation (33), Cauchy stress σ ij It is also decomposed accordingly, then the Cauchy stress σ The tensor is expressed as equation (34): (33) (34) In equation (34): I It is the identity matrix; The weak form equations of solid mechanics and the phase field model are as shown in equations (35)-(36): (35) (36) In equations (35)-(36): δx Representing variables x The trial function.

8. The TRISO fuel particle coating crack propagation simulation method according to claim 1, characterized in that, In step 6, the geometric defects to be studied are introduced into each coating layer and calculated. This includes first introducing defects that are common to particles during production and transportation into the coating layer, such as gaps in dense pyrolytic carbon layers and residual pores in silicon carbide layers.

9. The TRISO fuel particle coating crack propagation simulation method according to claim 1, characterized in that, After the calculation in step 7 is completed, the calculation results are output, stress distribution and phase field evolution images are plotted, and crack propagation is evaluated.

Images