Uniform debris bed reflooding numerical simulation method based on porous medium thermal unbalance model

Abstract

The invention discloses a uniform debris bed reflooding numerical simulation method based on a porous medium thermal unbalance model. The method comprises the following steps: 1, determining components and forms of a debris bed, coolant flow and an injection mode according to actual conditions; 2, establishing a simplified geometric model of the debris bed according to the debris bed form determined in the step 1, and performing mesh generation on the simplified geometric model to obtain a mesh model; 3, compiling a UDF based on the porous medium thermal unbalance model; 4, importing the grid model obtained in the step 2, setting boundary conditions according to the coolant injection mode and flow determined in the step 1, and setting specific parameters of a solver; 5, performing transient calculation until the whole debris bed re-submerging process is completed; and 6, performing post-processing on a numerical simulation result by utilizing a fluent post-processing tool. According to the method, the physical phenomenon, the re-submerging time and shock cooling leading edge migration of the debris bed re-submerging process are accurately simulated, the calculation speed is high, and the method has guiding significance on debris bed coolability prediction and a melt pressure vessel inner retention strategy.

Description

Technical Field

[0001] The present invention relates to the technical field of research on the cooling heat transfer characteristics of the debris bed at the lower head of a pressure vessel after a severe accident in a nuclear reactor, and particularly to a numerical simulation method for debris bed reflooding that couples internal flow and heat transfer in the debris bed based on a thermal non-equilibrium model of porous media. Background Art

[0002] When a nuclear reactor experiences an accident beyond the design basis, the core materials may melt, and the high-temperature melt will form a debris bed under the action of the core channels and coolant, or migrate to the lower head of the pressure vessel, forming a loose and porous debris bed under the action of the coolant inside the lower head. The cooling process of the debris bed involves complex single-phase (water), two-phase (water and steam) mixed flow states and coupled heat transfer, flow heat transfer, and phase change heat transfer processes among multiple phases (heating particles, cooling water, and escaping steam), with numerous influencing factors. The debris bed is the most easily cooled form of radioactive melt. If the debris bed in the core channels or the lower head can be effectively cooled, it can prevent the further development of nuclear power plant accidents, reduce the release of radioactive fission products, and lower the accident risk. The process of cooling the high-temperature debris bed by continuously injecting coolant is called debris bed reflooding. Therefore, studying the reflooding process of the nuclear reactor debris bed is of great significance.

[0003] With the improvement of computer technology and the booming development of Computational Fluid Dynamics (CFD) and numerical heat transfer, numerical simulation methods have become an important means for studying the coolability of debris beds. Through CFD technology, the operating state of the calculation model under various working conditions can be clearly understood. Many scholars have used numerical simulation methods for the reflooding process of debris beds.

[0004] The thermal non-equilibrium model of porous media in existing commercial CFD software is incompatible with the multiphase flow model. Abandoning either one will seriously affect the calculation accuracy. Therefore, it is of great significance to study numerical simulation methods that are compatible with the thermal non-equilibrium model of porous media and the multiphase flow model. Summary of the Invention

[0005] In order to overcome the problems existing in the above-mentioned prior art, the purpose of the present invention is to provide a numerical simulation method for uniform debris bed reflooding based on the thermal non-equilibrium model of porous media, which can accurately calculate the reflooding process of debris beds with different sizes and porosities under different coolant flow rates, subcooling degrees, and injection modes.

[0006] To achieve the above purpose, the present invention adopts the following technical solutions:

[0007] A numerical simulation method for uniform debris bed reflooding based on the thermal non-equilibrium model of porous media includes the following steps:

[0008] Step 1: Determine the composition and morphology of the debris bed, coolant flow rate, and coolant injection mode according to the actual situation; among them, the composition and morphology of the debris bed include the composition material of the debris bed particles, the diameter of the debris bed particles, the porosity of the debris bed, the height and diameter of the debris bed; the coolant flow rate is given in terms of velocity flow rate or mass flow rate, and the coolant injection mode includes injection from the top of the debris bed and injection from the bottom of the debris bed;

[0009] Step 2: Geometric modeling and mesh generation: First, simplify the debris bed into a cylinder for geometric modeling according to the height and diameter of the debris bed determined in Step 1, and generate a mesh to obtain a mesh model for simulating the reflooding process of the debris bed;

[0010] Step 3: Write a user-defined function UDF based on the thermal non-equilibrium model of porous media: The thermal non-equilibrium model of porous media includes a porous media flow resistance model and a porous media flow heat transfer model; the two models have different flow patterns in different flow pattern regions, and the flow patterns are divided according to the void fraction in the mesh. The void fraction is divided from high to low as follows: the void fraction in the single-phase gas phase region is 1, the void fraction in the film boiling region is between 0.9 and 1, the void fraction in the nucleate boiling region is between 0 and 0.9, and the void fraction in the single-phase water region is 0;

[0011] Step 4: Import the mesh model obtained in Step 2 and the user-defined function UDF obtained in Step 3 into fluent, set the boundary conditions according to the coolant injection mode and coolant flow rate determined in Step 1, set the specific parameters of the solver, and start the calculation;

[0012] Step 5: Based on the thermal non-equilibrium model of porous media, perform a transient simulation on the reflooding process of the debris bed to be simulated, iterate until the reflooding process ends and obtain the corresponding numerical simulation results;

[0013] Step 6: Use the fluent post-processing tool to post-process the numerical simulation results, obtain the phenomena and specific data of the debris bed reflooding, and predict the coolability of the debris bed.

[0014] The porous media flow resistance model and the porous media flow heat transfer model included in the thermal non-equilibrium model of porous media are as follows:

[0015] The porous media flow resistance model: The pressure drop is the sum of the pressure loss caused by viscous resistance and the inertial energy loss. The relationship is given in the following form, and this relationship is applicable to laminar flow, turbulent flow, and transitional flow patterns;

[0016]

[0017] In the formula:

[0018] ΔP——pressure drop / Pa;

[0019] L—the length of the pressure measurement section / m;

[0020] ε—the porosity of the porous medium;

[0021] μ—the dynamic viscosity of the fluid / kg·m -1 ·s -1 ;

[0022] d p —the particle diameter of the porous medium / m;

[0023] ρ—the density of the fluid / m;

[0024] U—the velocity of the fluid / m;

[0025] In the flow and heat transfer model of the porous medium, the temperatures of the solid phase and the fluid phase in the same grid are different, and there is heat transfer. In each time step, each grid is recalculated, and its heat transfer is given in the following form:

[0026] S = A fs h(T s -T f ) (2)

[0027] In the formula:

[0028] S—the heat transfer between the solid phase and the fluid phase of the porous medium / W·m -3 ;

[0029] A fs —the heat transfer area between the solid phase and the liquid phase per unit volume / m -1 ;

[0030] h—the convective heat transfer coefficient / W·m -2 ·K -1 ;

[0031] T s —the temperature of the solid phase of the porous medium / K;

[0032] T f —the temperature of the fluid phase of the porous medium / W·m -2 ·K -1 ;

[0033] The form of the convective heat transfer coefficient in formula (2) is different in different heat transfer regions. In the single-phase gas region, it is given in the form of formula (3), which is used to calculate the heat transfer between the debris bed particles and the vapor:

[0034] Nu1 = 0.27Re 0.8 Pr 0.4 (3)

[0035] In the formula:

[0036] Nu1——Nusselt number for convective heat transfer in the single-phase gas phase region;

[0037] Re——Reynolds number of the vapor;

[0038] Pr——Prandtl number of the vapor;

[0039] The Nusselt number for convective heat transfer in the single-phase region obtained by calculating through formula (3), and the relationship between the Nusselt number of convective heat transfer and the convective heat transfer coefficient given by formula (4) are used to obtain the convective heat transfer coefficient required by formula (2):

[0040] Nu1 = (hd p ) / k g (4)

[0041] Wherein:

[0042] k g ——Thermal conductivity of the vapor phase / WmK;

[0043] The Nusselt number for convective heat transfer in the single-phase liquid phase region is given by the following formula:

[0044]

[0045] Wherein:

[0046] Nu2——Nusselt number for convective heat transfer in the single-phase liquid phase region;

[0047] Re——Reynolds number of the coolant liquid phase;

[0048] The convective heat transfer coefficient in the nucleate boiling region is given by the following formula:

[0049] h smuc =4.63×10 6 f(prop)(T s -T sat ) m (6)

[0050] Wherein:

[0051] h smuc ——Convective heat transfer coefficient in nucleate boiling / W·m -2 ·K -1 ;

[0052] T s ——Solid phase temperature of the porous medium / K;

[0053] T sat ——Saturation temperature of the coolant / K;

[0054] f(prop) is a function of the fluid physical properties and is given by the following formula:

[0055]

[0056] In the formula:

[0057] μ f —— Dynamic viscosity of the liquid-phase coolant / kg·m -1 ·s -1 ;

[0058] c pf —— Specific heat at constant pressure of the liquid-phase coolant / J·kg·K;

[0059] h fg —— Latent heat of vaporization / J·kg -1 ;

[0060] σ —— Surface tension / kg·s -2 ;

[0061] g —— Acceleration due to gravity / m·s -2 ;

[0062] ρ f —— Density of the saturated liquid-phase coolant / kg·m -3 ;

[0063] ρ g —— Density of the saturated gas-phase coolant / kg·m -3 ;

[0064] k f —— Thermal conductivity of the liquid-phase coolant / W·m -1 ·K -1 ;

[0065] m in formula (6) is given by the following formula

[0066] m = 3.3 - 9.0e -d (8)

[0067] Where:

[0068]

[0069] The convective heat transfer coefficient in the film boiling region is calculated by the following formula:

[0070]

[0071] In the formula:

[0072] Nu3 —— Convective heat transfer coefficient in the film boiling region;

[0073] h sfb —— Convective heat transfer coefficient of film boiling / W·m -2 ·K -1 ;

[0074] D p —— Effective diameter of the particles / m;

[0075] k g —— Thermal conductivity of the gas coolant / W·m -1 ·K -1 ;

[0076] k f —— Thermal conductivity of the liquid coolant / W·m -1 ·K -1 ;

[0077] σ b —— Boltzmann constant;

[0078] —— Average Nusselt number of the particles for the saturated film boiling heat transfer coefficient;

[0079] —— Average Nusselt number of the particles based on the natural convection heat transfer coefficient;

[0080] Pr g —— Prandtl number of the gas coolant;

[0081] Pr f —— Prandtl number of the liquid coolant;

[0082] Sc —— Subcooling parameter of the liquid coolant; Sh —— Superheat parameter of the gas coolant;

[0083] h fg —— Latent heat of vaporization / J·kg -1 ;

[0084] μ g —— Dynamic viscosity of the gas coolant / kg·m -1 ·s -1 ;

[0085] μ f —— Dynamic viscosity of the liquid coolant / kg·m -1 ·s -1 ;

[0086] ρ f —— Density of the saturated liquid coolant / kg·m -3 ;

[0087] ρ g —— Density of the saturated gas coolant / kg·m -3 ;

[0088] T sat —— Saturation temperature / K;

[0089] T w —— The surface temperature of the particle / K;

[0090] ΔT sub —— The difference between the saturation temperature and the mainstream liquid temperature;

[0091] ΔT w —— The difference between the surface temperature of the particle and the saturation temperature;

[0092] β —— The thermal expansion coefficient of the liquid.

[0093] In the second step, the geometric modeling and meshing of the debris bed are three-dimensional. The parameter characteristics concerned in this step include the diameter of the debris bed particles, the height of the debris bed, the inlet diameter, and the coolant injection mode.

[0094] The fourth step is specifically as follows:

[0095] a) Import the mesh model in the second step and the user-defined function UDF in the third step into fluent;

[0096] b) Initialize the multiphase flow model: In fluent, set the multiphase flow model to the mixture model and turn on the energy equation;

[0097] c) Select the debris bed material, set the entire debris bed area as a porous medium area, and set its porosity and resistance coefficient;

[0098] d) Select the coolant material, use the physical property query tool to obtain the density, enthalpy, specific heat, viscosity, and thermal conductivity of the coolant at different temperatures and 0.1 MPa for sampling. Use data processing software to fit the above thermal physical parameters into a function of temperature and write them in the form of UDF and import them into fluent;

[0099] e) Set the boundary conditions, including the coolant inlet mass flow rate and subcooling degree, the outlet pressure, and the initial temperature of the debris bed. Set the specific parameters of the solver, including the relaxation factor, the velocity-pressure coupling method, the spatial discretization method, the time step, and the number of iterations.

[0100] The specific sign that the reflooding process in the fifth step ends is that the temperature of the entire debris bed is lower than the saturation temperature of the coolant.

[0101] In the sixth step, the phenomena of debris bed reflooding include the shape and migration process of the quench migration during the debris bed reflooding process; the specific data includes the flow field characteristic distribution, the gas volume fraction distribution, the phase change rate distribution, the debris bed fluid and particle temperature distribution, the debris bed pressure distribution, the quench front migration characteristics, and the debris bed reflooding time. Compared with the prior art, the present invention has the following advantages:

[0102] The numerical simulation method for the reflooding of a uniform debris bed based on the thermal non-equilibrium model of porous media of the present invention comprehensively considers various influencing parameters in the debris bed reflooding process, including the composition and morphology of the debris bed, the coolant flow rate, and the coolant injection mode, and has wide applicability; the method of the present invention can accurately simulate the physical phenomena, reflooding time, and quench front migration during the debris bed reflooding process, and realizes the visualization of the results by drawing contour maps, and has a fast calculation speed. In summary, the present invention has guiding significance for the prediction of the coolability of the debris bed and the retention strategy of the melt in the pressure vessel. Description of the Drawings

[0103] Figure 1 It is a schematic flow chart of the numerical simulation method for the reflooding of a uniform debris bed based on the thermal non-equilibrium model of porous media of the present invention.

[0104] Figure 2 It is the contour map of the solid phase temperature, void fraction, and phase change rate in the first stage of the top injection reflooding obtained by this numerical simulation method.

[0105] Figure 3 It is the contour map of the solid phase temperature, void fraction, and phase change rate in the second stage of the top injection reflooding obtained by this numerical simulation method.

[0106] Figure 4 It is the comparison between the quench front of the top injection reflooding obtained by this numerical simulation method and the experimental value.

[0107] Figure 5 It is the comparison between the quench front of the top injection reflooding obtained by this numerical simulation method and the experimental value. Detailed Embodiments

[0108] In order to introduce the present invention more clearly, the present invention will be further described in detail below with reference to the accompanying drawings. The introduction of the implementation scheme of the present invention is only for explaining the advantages of the present invention and does not constitute a limitation on the present invention.

[0109] As Figure 1 shown, the numerical simulation method for the reflooding of a uniform debris bed based on the thermal non-equilibrium model of porous media of the present invention takes the working condition of top injection with a diameter of 200 cm, a height of 600 mm, a porosity of 0.37, a water injection flow rate of 70 kg / h, an injection temperature of 30 °C, and an initial temperature of the debris bed of 230 °C as an example; it includes the following steps:

[0110] Step 1, determine according to the actual situation that the debris bed is a cylinder with a diameter of 200 cm and a height of 600 mm, its composition is iron balls with a diameter of 3 mm, and the porosity of its porous medium is 0.37. The coolant is liquid water, the flow rate is 70 kg / h, and the top injection mode.

[0111] Step 2: Geometric modeling and mesh generation. First, simplify the debris bed into a cylinder for geometric modeling according to the height and diameter of the debris bed determined in Step 1, and generate a three-dimensional mesh model for simulating the debris bed reflooding process by meshing in Fluent Meshing.

[0112] Step 3: Write a UDF based on the thermal non-equilibrium multiphase flow model in porous media. The thermal non-equilibrium model in porous media includes a flow resistance model and a flow heat transfer model in porous media. The above models are different in different flow pattern regions, and the flow patterns are divided according to the void fraction in the mesh. The void fraction is divided from high to low into: single-phase gas region (void fraction is 1), film boiling region (void fraction is between 0.9 and 1), nucleate boiling region (void fraction is between 0 and 0.9), and single-phase water region (void fraction is 0).

[0113] The flow resistance model and the flow heat transfer model in porous media are as follows:

[0114] For the flow resistance model in porous media, the pressure drop is the sum of the pressure loss caused by viscous resistance and the inertial energy loss. The relationship is given in the following form, and this relationship is applicable to laminar, turbulent, and transitional flow patterns;

[0115]

[0116] In the formula:

[0117] ΔP - Pressure drop / Pa;

[0118] L - Length of the pressure measurement section / m;

[0119] ε - Porosity of the porous media;

[0120] μ - Fluid dynamic viscosity / kg·m -1 ·s -1 ;

[0121] d p - Diameter of the porous media particles / m;

[0122] ρ - Fluid density / m;

[0123] U - Fluid velocity / m;

[0124] In the flow heat transfer model in porous media, the temperatures of the solid phase and the fluid phase are different in the same mesh, and there is heat transfer. In each time step, each mesh is recalculated, and its heat transfer is given in the following form:

[0125] S = A fs h(T s - T f ) (2)

[0126] In the formula:

[0127] S —— Heat transfer between the solid phase and the fluid phase of the porous medium / W·m -3 ;

[0128] A fs —— Heat transfer area between the solid phase and the liquid phase per unit volume / m -1 ;

[0129] h —— Convective heat transfer coefficient / W·m -2 ·K -1 ;

[0130] T s —— Solid phase temperature of the porous medium / K;

[0131] T f —— Fluid phase temperature of the porous medium / W·m -2 ·K -1 ;

[0132] The convective heat transfer coefficient in formula (2) has different forms in different heat transfer regions. In the single-phase gas phase region, it is given in the form of formula (3), which is used to calculate the heat transfer between the debris bed particles and the vapor:

[0133] Nu1 = 0.27Re 0.8 Pr 0.4 (3)

[0134] In the formula:

[0135] Nu1 —— Nusselt number of convective heat transfer in the single-phase gas phase region;

[0136] Re —— Reynolds number of the vapor;

[0137] Pr —— Prandtl number of the vapor;

[0138] By calculating the Nusselt number of convective heat transfer in the single-phase gas phase region obtained from formula (3), and the relationship between the Nusselt number of convective heat transfer and the convective heat transfer coefficient given by formula (4), the convective heat transfer coefficient required by formula (2) is obtained:

[0139] Nu1 = (hd p ) / k g (4)

[0140] In the formula:

[0141] k g —— Thermal conductivity of the vapor phase / WmK;

[0142] The Nusselt number of convective heat transfer in the single-phase liquid phase region is given by the following formula:

[0143]

[0144] In the formula:

[0145] Nu2——Nusselt number of convective heat transfer in the single-phase liquid phase region;

[0146] Re——Reynolds number of the coolant liquid phase;

[0147] The convective heat transfer coefficient in the nucleate boiling region is given by the following formula:

[0148] h smuc =4.63×10 6 f(prop)(T s -T sat ) m (6)

[0149] In the formula:

[0150] h smuc ——Convective heat transfer coefficient in nucleate boiling / W·m -2 ·K -1 ;

[0151] T s ——Solid phase temperature of the porous medium / K;

[0152] T sat ——Saturation temperature of the coolant / K;

[0153] f(prop) is a function of fluid physical properties and is given by the following formula:

[0154]

[0155] In the formula:

[0156] μ f ——Dynamic viscosity of the liquid-phase coolant / kg·m -1 ·s -1 ;

[0157] c pf ——Specific heat at constant pressure of the liquid-phase coolant / J·kg·K;

[0158] h fg ——Latent heat of vaporization / J·kg -1 ;

[0159] σ——Surface tension / kg·s -2 ;

[0160] g——Gravitational acceleration / m·s -2 ;

[0161] ρ f ——Density of the saturated liquid-phase coolant / kg·m -3 ;

[0162] ρ g —— Density of saturated vapor coolant / kg·m -3 ;

[0163] k f —— Thermal conductivity of liquid coolant / W·m -1 ·K -1 ;

[0164] m in formula (6) is given by the following formula

[0165] m = 3.3 - 9.0e -d (8)

[0166] where:

[0167]

[0168] The convective heat transfer coefficient in the film boiling region is calculated by the following formula:

[0169]

[0170] In the formula:

[0171] Nu3 —— Convective heat transfer coefficient in the film boiling region;

[0172] h sfb —— Convective heat transfer coefficient of film boiling / W·m -2 ·K -1 ;

[0173] D p —— Effective diameter of particles / m;

[0174] k g —— Thermal conductivity of vapor coolant / W·m -1 ·K -1 ;

[0175] k f —— Thermal conductivity of liquid coolant / W·m -1 ·K -1 ;

[0176] σ b —— Boltzmann constant;

[0177] —— Average Nusselt number of particles for saturated film boiling heat transfer coefficient;

[0178] —— Average Nusselt number of particles based on natural convection heat transfer coefficient;

[0179] Pr g —— Prandtl number of vapor coolant;

[0180] Pr f —— Prandtl number of the liquid coolant;

[0181] Sc—— Subcooling parameter of the liquid coolant; Sh—— Superheat parameter of the gas coolant;

[0182] h fg —— Latent heat of vaporization / J·kg -1 ;

[0183] μ g —— Dynamic viscosity of the gas coolant / kg·m -1 ·s -1 ;

[0184] μ f —— Dynamic viscosity of the liquid coolant / kg·m -1 ·s -1 ;

[0185] ρ f —— Density of the saturated liquid coolant / kg·m -3 ;

[0186] ρ g —— Density of the saturated gas coolant / kg·m -3 ;

[0187] T sat —— Saturation temperature / K;

[0188] T w —— Temperature of the particle surface / K;

[0189] ΔT sub —— Difference between the saturation temperature and the mainstream liquid temperature;

[0190] ΔT w —— Difference between the particle surface temperature and the saturation temperature;

[0191] β—— Thermal expansion coefficient of the liquid.

[0192] Step 4: Import the mesh model obtained in Step 2 and the UDF obtained in Step 3 into the solver, set the boundary conditions according to the coolant injection mode and flow rate determined in Step 1, set the specific parameters of the solver, and start the calculation. Step 4 is specifically as follows:

[0193] a) Import the mesh model in Step 2 and the user-defined function UDF in Step 3 into fluent;

[0194] b) Initialize the multiphase flow model: In fluent, set the multiphase flow model to the mixture model and turn on the energy equation;

[0195] c) Select the debris bed material as iron, set the entire debris bed area as a porous medium area, and set its porosity to 0.37 and the resistance coefficient;

[0196] d) Select the coolant material, use the physical property query tool to sample the density, enthalpy, specific heat, viscosity, and thermal conductivity of the coolant at different temperatures and 0.1 MPa. Use data processing software to fit the above thermal physical parameters into a function of temperature and write them in the form of a UDF and import them into fluent;

[0197] e) Set the boundary conditions, including the coolant inlet mass flow rate of 70 kg / h and temperature of 303 K, the outlet pressure of 0, and the initial temperature of the debris bed. Set the specific parameters of the solver, including the relaxation factor (pressure 0.3, velocity 0.7), the velocity-pressure coupling method as SIMPLE, the spatial discretization method as first-order upwind, the time step of 0.1 s, and the number of iterations of 5000 times.

[0198] Step 5: Based on the porous medium thermal non-equilibrium model, perform a transient simulation on the reflooding process of the debris bed to be simulated, iterate until the reflooding process ends and obtain the corresponding numerical simulation results. The specific sign of the end of reflooding is that the temperature of the entire debris bed is lower than the saturation temperature of the coolant, and this moment is 450 s.

[0199] Step 6: Use the fluent post-processing tool to post-process the numerical simulation results to obtain the phenomena and specific data of the debris bed reflooding. Draw a contour map according to the simulation results to predict the coolability of the debris bed. The phenomena include the shape and migration process of the quench front migration during the debris bed reflooding process, Figure 2 and Figure 3 The contour map can well reflect this process; the specific data includes the flow field characteristic distribution, gas volume fraction distribution, phase change rate distribution, debris bed fluid and particle temperature distribution, debris bed pressure distribution, quench front migration characteristics, and debris bed reflooding time. As Figure 4 and Figure 5 shown, the numerical simulation results of the reflooding time are in good agreement with the experiments.

Claims

1. A numerical simulation method for the reflooding of a uniform fragmented bed based on a thermal non-equilibrium model of porous media, characterized in that, It includes the following steps: Step 1: Determine the composition and morphology of the debris bed, coolant flow rate, and coolant injection mode according to the actual situation. Among them, the composition and morphology of the debris bed include the composition material of the debris bed particles, the diameter of the debris bed particles, the porosity of the debris bed, the height and diameter of the debris bed; the coolant flow rate is given in terms of velocity flow rate or mass flow rate, and the coolant injection mode includes injecting from the top of the debris bed and injecting from the bottom of the debris bed; Step 2: Geometric modeling and mesh generation: First, simplify the debris bed into a cylinder for geometric modeling according to the height and diameter of the debris bed determined in Step 1, and generate a mesh to obtain a mesh model for simulating the reflooding process of the debris bed; Step 3: Write a user-defined function UDF based on the thermal non-equilibrium model of porous media. The thermal non-equilibrium model of porous media includes a porous media flow resistance model and a porous media flow heat transfer model; the two models have different flow patterns in different flow pattern regions, and the flow patterns are divided according to the void fraction in the grid. The void fraction is divided from high to low as follows: the void fraction in the single-phase gas phase region is 1, the void fraction in the film boiling region is between 0.9 and 1, the void fraction in the nucleate boiling region is between 0 and 0.9, and the void fraction in the single-phase water region is 0; Step 4: Import the mesh model obtained in Step 2 and the user-defined function UDF obtained in Step 3 into fluent, set the boundary conditions according to the coolant injection mode and coolant flow rate determined in Step 1, set the specific parameters of the solver, and start the calculation; Step 5: Based on the thermal non-equilibrium model of porous media, perform a transient simulation on the reflooding process of the debris bed to be simulated, iterate until the reflooding process ends and obtain the corresponding numerical simulation results; Step 6: Use the fluent post-processing tool to post-process the numerical simulation results, obtain the phenomena and specific data of the debris bed reflooding, and predict the coolability of the debris bed.

2. The numerical simulation method for the reflooding of a uniform debris bed based on the thermal non-equilibrium model of porous media according to claim 1, wherein: The porous media flow resistance model and the porous media flow heat transfer model included in the thermal non-equilibrium model of porous media are specifically as follows: For the porous media flow resistance model, the pressure drop is the sum of the pressure loss caused by viscous resistance and the inertial energy loss, and the relationship is given in the following form, and this relationship is applicable to laminar flow, turbulent flow, and transitional flow patterns; In the formula: ΔP——pressure drop / Pa; L——length of the pressure measurement section / m; ε——porosity of the porous media; μ —— Kinematic viscosity / kg·m -1 ·s -1 ; d p ——Diameter of porous medium particle / m; ρ——fluid density / m; U——fluid velocity / m; In the porous media flow heat transfer model, the temperatures of the solid phase and the fluid phase in the same grid are different and there is heat transfer. In each time step, each grid is recalculated, and its heat transfer is given in the following form: S = A fs h(T s -T f ) (2) In the formula: S——Heat transfer between the solid phase and the fluid phase of the porous medium / W·m -3 ; A fs —— Heat transfer area between solid phase and liquid phase per unit volume / m -1 ; h——Convective heat transfer coefficient / W·m -2 ·K -1 ; T s ——Solid temperature of porous medium / K; T f —— Temperature of fluid phase in porous medium / W·m -2 ·K -1 ; The convective heat transfer coefficient in formula (2) has different forms in different heat transfer regions. In the single-phase gas phase region, it is given in the form of formula (3), which is used to calculate the heat transfer between the debris bed particles and the vapor: Nu1 = 0.27Re 0.8 Pr 0.4 (3) In the formula: Nu1——Nusselt number of convective heat transfer in the single-phase gas phase region; Re——Reynolds number of vapor; Pr——Prandtl number of vapor; The Nusselt number of gas-phase convective heat transfer in the single-phase region calculated by formula (3), and the relationship between the Nusselt number of convective heat transfer and the convective heat transfer coefficient given by formula (4) are used to obtain the convective heat transfer coefficient required by formula (2): Nu1 = (hd p ) / k g (4) Where: k g ——Vapor thermal conductivity / WmK; The Nusselt number of convective heat transfer in the single-phase liquid region is given by the following formula: Where: Nu2——The Nusselt number of convective heat transfer in the single-phase liquid region; Re——The Reynolds number of the coolant liquid phase; The convective heat transfer coefficient in the nucleate boiling region is given by the following formula: h smuc = 4.63×10 6 f(prop)(T s -T sat ) m (6) Where: h smuc —— Nucleate boiling convective heat transfer coefficient / W·m -2 ·K -1 ; T s ——Solid temperature of porous medium / K; T sat ——Coolant saturation temperature / K; f(prop) is a function of the physical properties of the fluid and is given by the following formula: Where: μ f —— Dynamic viscosity of the liquid coolant / kg·m -1 ·s -1 ; c pf —— Specific heat at constant pressure of the liquid coolant / J·kg·K; h fg —— Latent heat of vaporization / J·kg -1 ; σ —— surface tension / kg·s -2 ; g——acceleration of gravity / m·s -2 ; ρ f —— Saturated liquid coolant density / kg·m -3 ; ρ g —— Saturated vapor coolant density / kg·m -3 ; k f ——Thermal conductivity of the liquid coolant / W·m -1 ·K -1 ; m in formula (6) is given by the following formula m = 3.3 - 9.0e -d (8) Where: The convective heat transfer coefficient in the film boiling region is calculated using the following formula: Where: Nu3——The convective heat transfer coefficient in the film boiling region; h sfb —— Convective heat transfer coefficient of film boiling / W·m -2 ·K -1 ; D p —— effective diameter of particle / m; k g ——Thermal conductivity of the gas coolant / W·m -1 ·K -1 ; k f —— Thermal conductivity of the liquid coolant / W·m -1 ·K -1 ; σ b —— Boltzmann constant; ——Average Nusselt number of particles for saturated film boiling heat transfer coefficient; —— Average Nusselt number of particles based on the natural convection heat transfer coefficient; Pr g —— Prandtl number of the gas coolant; Pr f —— Prandtl number of the liquid coolant; Sc——The subcooling parameter of the liquid-phase coolant; Sh——The superheat parameter of the gas-phase coolant; h fg —— Latent heat of vaporization / J·kg -1 ; μ g —— Dynamic viscosity of the gas coolant / kg·m -1 ·s -1 ; μ f —— Dynamic viscosity of the liquid coolant / kg·m -1 ·s -1 ; ρ f —— Saturated liquid coolant density / kg·m -3 ; ρ g —— Saturated vapor coolant density / kg·m -3 ; T sat —— Saturation temperature / K; T w —— The surface temperature of the particle / K; ΔT sub —— The difference between the saturation temperature and the mainstream liquid temperature; ΔT w —— The difference between the particle surface temperature and the saturation temperature; β——The thermal expansion coefficient of the liquid.

3. The numerical simulation method for the reflooding of a uniform debris bed based on the thermal non-equilibrium model of porous media according to claim 1, characterized in that: In the second step described above, the geometric modeling and meshing of the debris bed are three-dimensional. The parameter characteristics concerned in this step include the diameter of the debris bed particles, the height of the debris bed, the inlet diameter, and the coolant injection mode.

4. The numerical simulation method for the reflooding of a uniform debris bed based on the thermal non-equilibrium model of porous media according to claim 1, characterized in that: The fourth step is specifically as follows: a) Import the mesh model in the second step and the user-defined function UDF in the third step into fluent; b) Initialize the multiphase flow model: In fluent, set the multiphase flow model to the mixture model and turn on the energy equation; c) Select the debris bed material, set the entire debris bed area as a porous medium area, and set its porosity and resistance coefficient; d) Select the coolant material, use the physical property query tool to obtain the density, enthalpy, specific heat, viscosity, and thermal conductivity of the coolant at different temperatures and 0.1 MPa for sampling. Use data processing software to fit the above thermal physical parameters as a function of temperature and write them in the form of UDF and import them into fluent; e) Set the boundary conditions, including the coolant inlet mass flow rate and subcooling degree, the outlet pressure, and the initial temperature of the debris bed. Set the specific parameters of the solver, including the relaxation factor, the velocity-pressure coupling method, the spatial discretization method, the time step, and the number of iterations.

5. The numerical simulation method for the reflooding of a uniform debris bed based on the thermal non-equilibrium model of porous media according to claim 1, characterized in that: The specific sign of the end of the reflooding process in the fifth step described above is that the temperature of the entire debris bed is lower than the saturation temperature of the coolant.

6. The numerical simulation method for the reflooding of a uniform debris bed based on the thermal non-equilibrium model of porous media according to claim 1, characterized in that: In the sixth step described above, the phenomena of debris bed reflooding include the shape and migration process of quenching migration during the debris bed reflooding process; The specific data includes the distribution of flow field characteristics, the distribution of gas volume fraction, the distribution of phase change rate, the distribution of debris bed fluid and particle temperature, the distribution of debris bed pressure, the migration characteristics of the quench front, and the debris bed reflooding time.

Images