High-burn-up foaming plate type nuclear fuel critical heat flux density numerical simulation method

Abstract

The invention discloses a numerical simulation method for critical heat flux density of high-burn-up foaming plate type nuclear fuel. The method comprises the following steps: 1, importing a geometric model and a grid model of the high-burn-up foaming plate type nuclear fuel into ANSYS fluent of commercial CFD (Computational Fluid Dynamics) software; 2, setting a calculation model and setting corresponding parameters; and 3, setting physical property parameters according to working conditions needing to be calculated. 4, setting boundary conditions of the model; and 5, starting calculation parameter monitoring. And 6, setting an initial volume heat source heating power value on the fuel pellet for heating, continuously increasing the volume heat source of the fuel pellet to increase the heat flux density of the wall surface of the fluid domain, and increasing the power density of the fuel pellet after the reading of the highest temperature and the highest void fraction of the heated wall surface converges. And 7, repeating the step 6 to continuously improve the power density of the heat source, and when it is monitored that the highest temperature of the heating wall surface is increased or the maximum void fraction of the heating wall surface is larger than 0.8, considering that critical occurs at the moment, and obtaining the obtained critical heat flux density. According to the numerical simulation method, the simulation result of existing test data is good, and the error ranges from-25% to + 10%.

Description

Technical Field

[0001] The invention belongs to the technical field of reactor safety analysis, and in particular relates to a numerical simulation method for critical heat flux density of high-burnup blister plate type nuclear fuel. Background Art

[0002] Plate-type fuel elements are a type of nuclear fuel that disperses fissile materials in the form of small particles in a non-fissile matrix of dispersed plate-type fuel. The shape of a single fuel is different from the rod-shaped fuel of a traditional pressurized water reactor, and is a long plate. Due to its low fuel element center temperature, high safety, large heat exchange area, large active area specific power, higher rigidity, high heat exchange efficiency, compact structure and other advantages, it is a relatively advanced compact core, and is currently widely used in research reactors and ship reactors. The coolant flow channel in the plate-type fuel assembly is a narrow rectangular flow channel formed between the plate fuel elements, and the narrow slit height of the narrow rectangular flow channel is generally 1 to 3 mm. From the inlet to the outlet, the coolant may experience several flow transitions, including single-phase water, steam-water two-phase flow including subcooled boiling and saturated boiling, and single-phase steam. The flow and heat transfer conditions are extremely complex.

[0003] The bubbling phenomenon is a unique failure form of dispersed plate fuel elements. Due to long-term in-core operation, cracks are generated inside the high-burnup fuel cores and fission gases such as krypton and xenon generated during nuclear fission burnup are released. The local accumulation of increasing gas pressure causes bubbling of the cladding, resulting in the distortion of the rectangular narrow slit flow channel. It is a key bottleneck technology in the thermal hydraulic design and safety analysis of dispersed fuel elements. The critical heat flux is the heat flux when the heat transfer on the surface of the nuclear fuel element deteriorates. It is the result of the transformation of the two-phase boiling heat transfer mechanism. At present, the research on boiling criticality and thermal hydraulic characteristics in rectangular narrow slit flow channels is mostly limited to normal rectangular narrow slit flow channels without bubbling. The research on gas-liquid two-phase boiling heat transfer and boiling criticality of rectangular narrow slit flow channels after bubbling is relatively lacking. In addition, the critical heat flux, an important safety criterion, should be focused on under normal operation and accident conditions of nuclear reactors.

[0004] After bubbling occurs in plate-type nuclear fuel elements, both boiling heat transfer characteristics and boiling criticality will be affected by changes in the coolant flow channel geometry. Therefore, numerical simulations are used to study important parameters such as the temperature distribution of the fuel element before the boiling criticality occurs, as well as the distribution of the near-wall cavitation fraction of the heated wall when the boiling criticality occurs; explore the effect of bubbling on the two-phase flow boiling mechanism inside the distorted rectangular narrow slit flow channel; analyze the boiling heat transfer characteristics after bubbling to provide a technical basis for the formulation of thermal safety criteria after bubbling. The above research objectives achieved by the present invention are of great significance for safety research under various operating conditions after high fuel consumption operation in the core of plate-type fuel elements. Summary of the invention

[0005] In order to solve the problems existing in the above-mentioned prior art, the purpose of the present invention is to provide a numerical simulation method for the critical heat flux density of high-burnup bubbled plate-type nuclear fuel. After the geometric model and grid of the plate-type fuel element after bubble formation are given, the fluid, solid properties and corresponding thermal-hydraulic parameters are set, and the critical heat flux density of high-burnup bubbled plate-type nuclear fuel can be calculated.

[0006] In order to achieve the above object, the present invention adopts the following technical solution:

[0007] A numerical simulation method for critical heat flux density of high-burnup blister plate nuclear fuel comprises the following steps:

[0008] Step 1: Establish a geometric model of high-burnup blister plate nuclear fuel and perform meshing. Import the geometric model and mesh file into the commercial CFD software ANSYS fluent through ANSYS workbench.

[0009] Step 2: Call the calculation model in ANSYS fluent. On the basis of the Euler-Euler two-fluid model and the RPI wall boiling model, add the interphase action model, critical heat flux density model, turbulent Reynolds pressure model RSM and auxiliary models. The auxiliary models include the wall nucleation density model, bubble detachment diameter model, bubble detachment frequency model and bubble diameter model, and set the corresponding parameters according to the pressure, temperature and flow rate of the calculation conditions.

[0010] Step 3: Set the physical property parameters. Set the material physical property parameters according to the numerical simulation conditions. In the numerical simulation process of critical heat flux density, the coolant involves a phase change process from liquid to gas. The coolant physical property parameters are set using a linear function relationship to characterize the physical property changes. The physical property parameters include density, pressure heat capacity, thermal conductivity, viscosity and enthalpy value of the fluid, and density, specific heat and thermal conductivity of the solid;

[0011] Step 4: Set boundary conditions, set the fluid domain inlet to mass flow rate inlet, the outlet to pressure outlet, the fuel pellet to volume heat source, and use standard wall function for fluid domain-solid domain wall treatment;

[0012] Step 5: Use the monitors tab provided by Fluent to add thermal parameter monitoring during the numerical simulation process, monitor the maximum temperature of the heated wall of the fluid domain and the maximum cavitation fraction of the gas phase in the grid near the heated wall, and present them as a two-dimensional curve, where the horizontal axis is the number of iterative calculation steps of Fluent;

[0013] Step 6: Set the initial volume heat source heating power value on the fuel core block for heating. The initial volume heat source heating power value is set to a heat flux density of the volume heat source through the plate-type nuclear fuel cladding to the fluid heating surface of 30% of the critical heat flux density CHF experimental value. Observe the monitoring chart. After the maximum temperature and the maximum cavitation fraction readings of the heated wall converge, increase the volume heat source power density. The power density increment is 5% of the experimental value. When the power density is close to the experimental value, reduce the power density increment to 1% of the experimental value.

[0014] Step 7: Repeat step 6 to continuously increase the volume heat source power density. When the maximum temperature of the heated wall is monitored to rise or the maximum cavitation fraction of the heated wall is greater than 0.8, it is considered that the critical state has occurred. The wall heat flux density of the heated wall in the fluid domain at this time is the critical heat flux density under this working condition.

[0015] The high-burnup bubble plate nuclear fuel geometric model described in step 1 is modeled using SolidWorks modeling software. The fuel core block, fuel cladding and coolant fluid domain are established, and fission gas air gaps are established in the bubbled cladding bubbles. The internal air gap thermal resistance of the bubble is equivalently replaced, and the fluidity of the fission gas air gap is ignored. The air gap in the cladding bubble is equivalent to a small solid thermal resistance, which simplifies the geometric model and reduces the calculation cost.

[0016] In step 2, the mass, energy and momentum equations of each phase of the gas and liquid phases were solved respectively in the two-phase boiling heat transfer process of water as a nuclear fuel coolant heated by the wall by using the Euler-Euler two-fluid model, the RPI wall boiling model and the turbulent Reynolds pressure model RSM; the numerical solution of the phase interface characteristics and the mass, energy and momentum transfer process between the phases were realized by adding the interphase action model and the auxiliary model, and the solution parameters included the phase interface density, drag, lift, wall lubrication force and turbulent dissipation force; the critical heat flux density model ensured the convergence of the numerical simulation calculation under the condition of boiling critical approximate burning.

[0017] In step 3, the physical property parameters are set. The density, thermal conductivity and viscosity of the fluid are characterized by the linear interpolation method with respect to temperature. The pressure heat capacity and enthalpy values ​​are implemented by the UDF-based interpolation method. The pressure heat capacity selects the pressure heat capacity of the inlet supercooled water and the saturated temperature water at the corresponding pressure for temperature interpolation. The enthalpy value is obtained by integrating the pressure heat capacity. The density, specific heat and thermal conductivity of the fuel core and fuel cladding are interpolated at the corresponding temperature. The thermal resistance of the fission gas in the cladding bubble is characterized by temperature interpolation of density, specific heat and thermal conductivity based on the mass fraction of krypton and xenon of 90% and 10% respectively.

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

[0019] The present invention realizes the numerical simulation of the critical heat flux density of high-burnup blister plate type nuclear fuel through the Euler-Euler two-fluid model CFD fluid-solid coupling simulation method, which has guiding significance for the steady-state and accident conditions of this type of nuclear fuel under the condition of high operating life, and provides support for the design of the distortion heating plate test section of similar geometric models. The numerical simulation method of the present invention has good simulation results for existing test data, with an error of -25% to +10%. BRIEF DESCRIPTION OF THE DRAWINGS

[0020] Figure 1 It is a flow chart of the numerical simulation method of the critical heat flux density of high-burnup blister plate type nuclear fuel of the present invention.

[0021] Figure 2 It is an error analysis diagram of the numerical simulation method of the critical heat flux density of the high-burnup blister plate type nuclear fuel compared with experimental data of the present invention. DETAILED DESCRIPTION

[0022] The present invention will be further described in detail below in conjunction with the accompanying drawings and specific embodiments. Figure 1 As shown, the present invention is a numerical simulation method for the critical heat flux density of high-burnup blister plate nuclear fuel, which mainly comprises the following steps:

[0023] Step 1: Establish the geometric model of high-burnup bubbler plate nuclear fuel. The geometric model of high-burnup bubbler plate nuclear fuel is modeled using SolidWorks 3D modeling software, and the solid and liquid phases are distinguished at the same time. The geometric structure parameters include the length (mm), width (mm), and height (mm) of the rectangular narrow slit flow channel in the fluid domain, the thickness (mm) of the fuel cladding in the solid domain, and the thickness (mm) of the fuel core block. The bubble parameters include the bubble size and specific distribution position. In ANSYS workbench, the geometric model is imported into Meshing for meshing, and the hexahedral mesh structure is used. The mesh quality is controlled to be above 0.8, and the maximum aspect ratio of the mesh is controlled to be below 30 to avoid the problem of non-convergence in the numerical simulation. The same boundary conditions are used to perform mesh independence analysis and verification, calculate and determine the optimal number of meshes, and generate mesh files. In ANSYS workbench, the geometric model and mesh files are imported into ANSYS fluent.

[0024] Step 2: Call the calculation model in the models tab of ANSYS fluent, including the Euler-Euler two-fluid model, the RPI wall boiling model, and the turbulent Reynolds pressure model RSM and auxiliary models.

[0025] The Euler-Euler two-fluid model is a numerical simulation model that solves the mass conservation, momentum conservation, and energy conservation control equations of each phase separately, and solves the phase interface characteristics and the heat and mass transfer process between the phases by introducing appropriate interphase interaction models and auxiliary models. The interphase energy exchange in the present invention is realized by the Ranz-Marshall formula. The heat exchange between the liquid phase and the gas phase Q gl for:

[0026] Q gl =-Q lg =h l (T g -T l )A i

[0027]

[0028] In the formula, Q lg ——heat exchange between liquid and gas phase; T g ——gas phase temperature; R1——liquid phase temperature; A i ——The total area of ​​the vapor-liquid interface per unit volume (interface density); d g ——bubble diameter; h l ——Heat transfer coefficient / W·m -2 ·K -1 ; l ——Liquid phase thermal conductivity / W·m -1 ·K -1 ;Re p ——Reynolds number with gas phase diameter as reference length and relative velocity between phases as reference velocity; Pr——Prandtl number of liquid phase.

[0029] The phase exchange added by the present invention is the evaporation of the liquid phase Condensation of the gas phase In two-phase flow heat transfer, it is expressed as follows:

[0030]

[0031] Among them, q e ——Evaporation heat flux / W·m -2 , extracted from the RPI wall boiling model, F / V——the ratio of the wall area to the volume of the near-wall grid / m 2 ·m -3 , which is used to convert the surface heat flux density into the volume heat source of the first layer of grid; h g ——Gas phase enthalpy / J·kg -1 ;h l ——Liquid phase enthalpy / J·kg -1 .

[0032] The interphase interaction model added by the present invention includes drag force, lift force, turbulent dissipation force and wall lubrication force.

[0033] The drag force refers to the force acting in the opposite direction of the relative motion on the bubble due to the viscosity when the bubble is dispersed in the continuous fluid and there is relative motion between the two. The drag force acting on a single bubble can be expressed as:

[0034]

[0035] In the formula: C d - Drag coefficient, calculated using the relationship recommended by Schiller-Naumann:

[0036]

[0037] In the formula, Re p ——Reynolds number with the gas phase diameter as the reference length and the relative velocity between phases as the reference velocity; ρ l ——Liquid phase density; ——Liquid phase velocity; ——Gas phase velocity.

[0038] The lift is caused by the unevenness of the flow field. The velocity gradient of the liquid phase will inevitably exist in the direction perpendicular to the flow direction. When the dispersed bubbles move in the liquid shear flow field, they will be affected by the lift. The direction of the lift is perpendicular to the relative velocity of the vapor and liquid phases, which can be expressed as the expression of the relative velocity of the two phases and the liquid phase velocity curl, that is:

[0039]

[0040] Where: C L ——lift coefficient; ρ l ——Liquid phase density; α g - cavitation fraction; ——Liquid phase velocity; ——Gas phase velocity.

[0041] The present invention selects the Moraga model to calculate the phase lift, and the lift coefficient calculation model is:

[0042]

[0043] Where: φ is the product of the bubble Reynolds number and the vortex Reynolds number, which is calculated by the bubble diameter d g Calculate using the following formula:

[0044]

[0045]

[0046] Where: l ——Liquid phase density; ——Liquid phase velocity; ——Gas phase velocity; μ l - viscosity;

[0047] In the Euler-Euler two-fluid model, the eddy viscosity turbulence model eliminates the additional force of the pulsation on the discrete phase by the ensemble average of the turbulent pulsation, so additional turbulent dissipation force is required. The turbulent dissipation force itself is a force generated by modeling, and the turbulent dissipation force mainly depends on the gradient of the vapor phase cavitation fraction in the flow field. The present invention uses the calculation model based on phase drag proposed by Burns et al. after mass averaging, and the expression is as follows:

[0048]

[0049] Where: ——turbulent dissipation force; C TD ——turbulence dissipation coefficient, the value is 1.0; σ lg ——Constant, 0.9; C d ——Drag coefficient; μ t,l ——Liquid turbulent viscosity / Pa·s; μ l ——Liquid phase molecular viscosity / Pa·s; the meanings of other symbols are the same as above.

[0050] The wall lubrication force refers to the large difference in velocity distribution around the bubble near the wall, which leads to a higher near-wall liquid phase velocity gradient. Due to the no-slip condition, the continuous phase flow between the bubble and the wall decreases while the flow on the other side of the bubble increases, thus causing an additional force on the bubble that pushes the bubble away from the wall. The wall lubrication force is calculated by the following formula:

[0051]

[0052] Where: ——unit vector normal to the wall; α g - cavitation fraction; ——Liquid phase velocity; ——gas phase velocity; C W ——Wall lubrication coefficient, calculated by Antal model:

[0053]

[0054] Where: C1, C2——dimensionless constants, respectively -0.01 and 0.05; y w ——Distance from the wall / m. The wall lubrication force only acts on the area with liquid phase velocity gradient, that is, the fluid domain near the wall. When the distance from the wall meets this condition: When , the wall lubrication force is 0, indicating that the position is in the fully developed turbulence area, the radial gradient of the liquid phase velocity is very small, and the wall lubrication force can be ignored.

[0055] The RPI wall boiling model is the most widely used wall boiling model. It divides the heat transfer process into three heat transfer modes, namely single convection heat flow, evaporation heat flow, and quenching heat flow. The sum of the three is the total wall heat flow. The turbulence model uses the turbulent Reynolds pressure model RSM. This model is the best for simulating flows in complex geometries in the Reynolds time-averaged method RANS turbulence model category. Although the solution cost increases, the simulation results of complex flows are better than other simple RANS turbulence models, such as the standard k-ε model.

[0056] Finally, in order to close the equation and realize the numerical solution of flow boiling, some auxiliary models need to be introduced, including the heating surface nucleation density model, the bubble detachment diameter model, the bubble detachment frequency model and the bubble diameter model.

[0057] The bubble detachment diameter model is the diameter of the bubble when it detaches from the heated wall. Its size depends on the force balance during the bubble growth process. The bubble detachment model can be obtained through bubble force analysis or through data summary of experimental observations. The present invention uses the empirical model proposed by Tolubinsky et al. based on experimental observations to calculate the bubble detachment diameter:

[0058]

[0059] Where: ΔT sub ——Subcooling degree of fluid near the wall / K; d ref =0.6mm; ΔT ref =45K; d max =1.4mm.

[0060] The wall nucleation density model is the number of vaporization cores per unit heated wall area, usually expressed as N w The present invention adopts the empirical formula proposed by Lemmert-Chawla to calculate the nucleation density on the heating surface. In this model, the nucleation density is only related to the wall superheat, that is:

[0061] N w =C n ΔT sup n

[0062] Where: C = 210; n = 1.805; ΔT sup ——Wall superheat / K, ΔT sup =T w -T sat .

[0063] The bubble detachment frequency is a measure of the bubble detachment speed on the heated wall and is an important parameter in bubble dynamics. For a certain nucleation point, assuming that the bubble waiting time is t wt , the bubble growth time is t ct , then the bubble detachment frequency can be expressed as:

[0064]

[0065] Cole et al. recommended using the following relationship for calculation:

[0066]

[0067] Where: C f ——Bubble detachment wall resistance factor, the default value is 1; g——gravitational acceleration; ρ l ——Liquid phase density; ρ g ——Gas phase density, d w ——Bubble detachment diameter The present invention performs calculations based on this model.

[0068] The bubble diameter model adopts an improved model of the bubble diameter model proposed by Anglart et al. This model sets upper and lower limits for the bubble diameter. It is believed that when the fluid supercooling is higher than ΔT0, the bubble diameter in the mainstream is d0, and when the fluid supercooling is less than ΔT1, the bubble diameter is d1, as shown in the following formula:

[0069]

[0070] Where: d g ——77 bubble diameter; d0 = 0.0001 m; d1 = 0.0015 m; ΔT0 = 13.5 K; ΔT1 = 0 K; ΔT sup ——Wall superheat / K.

[0071] The above models can all be selected in the fluent model library.

[0072] Step 3: Set the material properties. Set the physical property parameters according to the working conditions to be calculated. Set the fluid domain system pressure in the cell zone condition tab. For example, the saturation temperature of the coolant at this pressure is 13.8MPa, which is 608.68K at 13.8MPa. In the materials tab, the changes in physical property parameters caused by temperature changes are expressed by linear interpolation. The pressure heat capacity and enthalpy values ​​of gas and liquid phase fluids are expressed by inserting UDF compilation. Density, thermal conductivity, and viscosity are interpolated between the inlet subcooled fluid temperature and the saturation temperature through the piece linear method. The materials in the solid domain include fuel cores and cladding, and the density, specific heat, and thermal conductivity are interpolated with respect to temperature. Specific materials include water, water vapor, aluminum alloy, U10Mo alloy, etc.

[0073] Step 4: Set the model boundary conditions, use the mass flow rate inlet and pressure outlet, add the source term to the entire solid domain of the fuel pellet as a volume heat source through fluent, realize the fluid-solid coupling heat transfer calculation of the fuel pellet-fuel cladding-fluid domain, and use the standard wall function for wall treatment. Enter the initial conditions calculated by the corresponding working conditions in the initialization interface, such as inlet temperature, fluid velocity, turbulent kinetic energy, turbulent dissipation rate, and saturation temperature.

[0074] Step 5: Use the report plots function in the monitors tab of Fluent to add thermal parameter monitoring during the numerical simulation calculation process. Select the heating boundary surface of the solid to the fluid as the monitoring surface. Use the facet-maximum option to add the maximum temperature wall temperature and the maximum gas phase cavitation volume fraction to the monitoring parameters. This enables monitoring of the maximum temperature of the heated wall of the fluid domain and the maximum gas phase cavitation volume fraction of the near-wall grid of the heated wall, and presents them as a two-dimensional curve. The horizontal coordinate of the curve image is the number of Fluent iterative calculation steps.

[0075] Step 6: Set the initial volume heat source heating power value in the fuel pellets for heating. Set the initial volume heat source heating power value in the source term of the solid domain of the fuel pellets to 30% of the experimental value of the critical heat flux density CHF through the plate-type nuclear fuel cladding on the fluid heating surface, and start the simulation calculation until convergence. During the calculation, observe the monitoring chart, and increase the volume heat source power density after the maximum temperature and the maximum cavitation fraction readings of the heated wall converge. The power density increment is 5% of the experimental value. When the power density is close to the experimental value, reduce the power density increment to 1% of the experimental value.

[0076] Step 7: Repeat step 6 to continuously increase the volume heat source power density of the fuel core. When the maximum temperature of the heated wall is monitored to rise or the maximum cavitation fraction of the heated wall is greater than 0.8, it is considered that the critical state has occurred. The wall heat flux density of the heated wall in the fluid domain at this time is the critical heat flux density under this working condition. Figure 2 The numerical simulation method of the critical heat flux density of high-burnup bubbled plate type nuclear fuel of the present invention simulates the critical heat flux density in a rectangular narrow slit channel under normal non-bubbling conditions, and compares the error of the critical heat flux density experimental data of a rectangular narrow slit channel with double heating surfaces in the Bettis laboratory of Westinghouse Electric Company. The simulation error of the critical heat flux density in a normal rectangular narrow slit channel is -24.92% to 9.39%, and the overall simulation result is good.

Claims

1. A numerical simulation method for the critical heat flux density of a high-burnup blister plate type nuclear fuel, characterized in that: The following steps are involved: Step 1: Establish a geometric model of high-burnup blister plate nuclear fuel and perform meshing. Import the geometric model and mesh file into the commercial CFD software ANSYS fluent through ANSYS workbench. Step 2: Call the calculation model in ANSYS fluent. On the basis of the Euler-Euler two-fluid model and the RPI wall boiling model, add the interphase action model, critical heat flux density model, turbulent Reynolds pressure model RSM and auxiliary models. The auxiliary models include the wall nucleation density model, bubble detachment diameter model, bubble detachment frequency model and bubble diameter model, and set the corresponding parameters according to the pressure, temperature and flow rate of the calculation conditions. Step 3: Set the physical property parameters. Set the material physical property parameters according to the numerical simulation conditions. In the numerical simulation of the critical heat flux density, the coolant involves a phase change process from liquid to gas. The coolant physical property parameters are set using a linear function relationship to characterize the physical property changes. The physical property parameters include density, pressure heat capacity, thermal conductivity, viscosity and enthalpy of the fluid, and density, specific heat and thermal conductivity of the solid; Step 4: Set boundary conditions, set the fluid domain inlet to mass flow rate inlet, the outlet to pressure outlet, the fuel pellet to volume heat source, and use standard wall function for fluid domain-solid domain wall treatment; Step 5: Use the monitors tab provided by Fluent to add relevant thermal parameters monitoring during the numerical simulation process, monitor the maximum temperature of the heated wall of the fluid domain and the maximum cavitation fraction of the gas phase in the grid near the heated wall, and present them as a two-dimensional curve, where the horizontal axis is the number of iterative calculation steps of Fluent; Step 6: Set the initial volume heat source heating power value for heating the fuel core block. The initial volume heat source heating power value is set to a heat flux density of the volume heat source through the plate-type nuclear fuel cladding to the fluid heating surface of 30% of the critical heat flux density CHF experimental value. Start the simulation calculation and observe the monitoring chart. After the maximum temperature and the maximum cavitation fraction readings of the heated wall converge, increase the volume heat source power density. The power density increment is 5% of the experimental value. When the power density is close to the experimental value, reduce the power density increment to 1% of the experimental value. Step 7: Repeat step 6 to continuously increase the volume heat source power density. When the maximum temperature of the heated wall is monitored to rise or the maximum cavitation fraction of the heated wall is greater than 0.8, it is considered that the critical state has occurred. The wall heat flux density of the heated wall in the fluid domain at this time is the critical heat flux density under this working condition.

2. The numerical simulation method for the critical heat flux density of a high-burnup blister plate type nuclear fuel according to claim 1 is characterized in that: The high-burnup bubble plate nuclear fuel geometric model described in step 1 is modeled using SolidWorks modeling software. The fuel core block, fuel cladding and coolant fluid domain are established, and fission gas air gaps are established in the bubbled cladding bubbles. The internal air gap thermal resistance of the bubble is equivalently replaced, and the fluidity of the fission gas air gap is ignored. The air gap in the cladding bubble is equivalent to a small solid thermal resistance, which simplifies the geometric model and reduces the calculation cost.

3. The numerical simulation method for the critical heat flux density of a high burnup blister plate type nuclear fuel according to claim 1 is characterized in that: In step 2, the mass, energy and momentum equations of each phase of the gas phase and the liquid phase are solved respectively in the two-phase boiling heat transfer process of water as a nuclear fuel coolant heated by the wall by using the Euler-Euler two-fluid model, the RPI wall boiling model and the turbulent Reynolds pressure model RSM; the numerical solution of the phase interface characteristics and the mass, energy and momentum transfer process between the phases are realized by adding the interphase action model and the auxiliary model, and the solution parameters include the phase interface density, drag, lift, wall lubrication force and turbulent dissipation force; The critical heat flux model ensures the convergence of numerical simulation calculations under the condition of boiling critical approximation burnout.

4. The numerical simulation method for the critical heat flux density of a high burnup blister plate type nuclear fuel according to claim 1 is characterized in that: In step 3, the physical property parameters are set. The density, thermal conductivity and viscosity of the fluid are characterized by the linear interpolation method with respect to temperature. The pressure heat capacity and enthalpy values ​​are implemented by the UDF-based interpolation method. The pressure heat capacity selects the pressure heat capacity of the inlet supercooled water and the saturated temperature water at the corresponding pressure for temperature interpolation. The enthalpy value is obtained by integrating the pressure heat capacity. The density, specific heat and thermal conductivity of the fuel core and fuel cladding are interpolated at the corresponding temperature. The thermal resistance of the fission gas in the cladding bubble is characterized by temperature interpolation of density, specific heat and thermal conductivity based on the mass fraction of krypton and xenon of 90% and 10% respectively.

Images