A simplified calculation method of debris bed melt based on heat transfer rate

Abstract

This invention discloses a simplified calculation method for fragment bed melting based on heat transfer rate, comprising: solving the enthalpy conservation equation; calculating the liquid fraction and temperature gradient vector of the fragment bed region grid and storing them in the grid; calculating the geometric dimensions of the molten pool formed by the fragment bed melting; calculating the characteristic parameters of the molten pool in the fragment bed melting region; calculating the heat transfer characteristic scale within the molten pool; calculating the heat transfer rate in the liquid phase region; calculating the heat transfer rate components of the grid and the heat transfer rate in the mushy region; calculating the explicit convection term; uniformly removing the explicit convection term; calculating the source term after the energy conservation equation is corrected; iteratively calculating the energy conservation equation to obtain the latest temperature distribution and liquid fraction distribution, thus completing the calculation of the fragment bed melting process within a unit time step. This invention eliminates the need to solve the momentum equation in the fragment bed melting process, while still describing the heat dissipation and temperature changes of the fragment bed in different directions, greatly improving the computational efficiency of the fragment bed melting process and ensuring computational accuracy.

Description

Technical Field

[0001] This invention relates to the field of computational fluid dynamics, specifically to a simplified calculation method for the melting of a fragmented bed based on heat transfer velocity. Background Technology

[0002] Following a core meltdown in a severe accident, the molten material migrates to the lower head of the pressure vessel and is cooled by cooling water, forming a solid debris bed. If the debris bed is not effectively cooled, it will reheat and melt to form a molten pool as decay heat is released. The melting of the debris bed significantly impacts heat transfer within the molten pool and the heat flux density distribution on the walls, thus affecting the integrity of the pressure vessel and the state of molten material leakage; it is a crucial phenomenon in severe accidents. However, due to the long transient nature of the debris bed melting process, current research on the dynamic process of debris bed melting and molten pool formation is still limited. Therefore, finding an efficient and simple calculation method is of great significance for the study of severe reactor accidents.

[0003] In computational fluid dynamics, the enthalpy-porosity method is mainly used to solve solid-liquid phase transition processes. For molten pools formed by the melting of debris beds, the RANS equations often fail to describe the flow due to the strong turbulence within them; therefore, current research on heat transfer in molten pools primarily employs the LES method. The LES method directly solves for large-scale eddies, while the influence of small eddies on large eddies is considered through approximate models. This method typically requires high computer memory and speed. For transient processes like debris bed melting that can last for tens of hours, the computational efficiency of using the enthalpy-porosity method coupled with the LES method is low. Furthermore, the computational cost for three-dimensional simulations of full-scale pressure vessel head debris bed melting is extremely high. The direct solution of large eddies in the LES turbulence model is the main factor contributing to the high computational cost. However, for debris bed melting processes, we are more concerned with the heat dissipation and temperature changes of the debris bed in different directions during melting. If we could describe the temperature evolution and heat distribution during melting without solving for molten pool turbulence, it would be of great significance for improving the computational speed of numerical simulations. Summary of the Invention

[0004] To address the issue of high computational complexity in the fragment bed melting process caused by the use of the LES turbulence model, this invention aims to provide a simplified computational method for fragment bed melting based on heat transfer velocity. This method, based on the concept of heat transfer velocity, considers only the influence of natural convection within the molten pool on the melting process, thus avoiding the solution of the momentum equation in the fragment bed melting process. Using this method, only the energy conservation equation needs to be solved, effectively reducing the computational load while obtaining the transient temperature distribution during fragment bed melting. This method can be conveniently applied to two-dimensional and three-dimensional fragment bed melting calculations.

[0005] The objective of this invention is achieved through the following technical solution.

[0006] A simplified calculation method for melting of a fragmented bed based on heat transfer rate includes the following steps:

[0007] S1. Solve the enthalpy conservation equation, calculate the liquid fraction and temperature gradient vector of the fragment bed region grid, and store them in the corresponding grid.

[0008] S2. Calculate the geometric dimensions of the molten pool formed by the melting of the fragmented bed, and obtain the height and width of the molten pool;

[0009] S3. Calculate the characteristic parameters of the molten pool in the fragmented bed melting region;

[0010] S4. Calculate the characteristic scale of heat transfer within the molten pool;

[0011] S5. Calculate the heat transfer rate in the liquid phase region where the fragmented bed grid is located;

[0012] S6. Calculate the heat transfer velocity components of the grid and the heat transfer velocity of the mushy region;

[0013] S7. Calculate the explicit convection term as the heat transferred to the boundary due to natural convection of the molten pool;

[0014] S8. Explicit convection terms are uniformly removed to achieve energy conservation throughout the entire fragmented bed molten pool;

[0015] S9. Calculate the source terms of the energy conservation equation after correction;

[0016] S10. Iteratively calculate the energy conservation equation to obtain the latest temperature distribution and liquid fraction distribution, and complete the calculation of the fragment bed melting process within a unit time step.

[0017] Furthermore, in step S2, the geometric parameters of the melting region of the debris bed are calculated based on the liquid fraction of each debris bed region grid; by traversing the grids with a liquid fraction greater than 0, the maximum and minimum values ​​of the grid's horizontal and vertical coordinates are obtained, and the height and width of the current melting region of the debris bed can be calculated.

[0018] Further, in step S3, based on the height of the molten pool in the fragmented bed obtained in step S2, the characteristic parameters of the molten pool are calculated. The characteristic parameters of the molten pool include the Rayleigh number Ra and the Nussel number Nu. The Rayleigh number Ra is calculated as shown in formula (1). When calculating the Nussel number Nu, the heat transfer correlation obtained from the molten pool experiment is used.

[0019]

[0020] In the formula, Ra represents the Rayleigh number of the molten pool after the fragmented bed begins to melt; g represents the acceleration due to gravity; β represents the coefficient of thermal expansion of the molten pool; q vν represents the decay heat power density of the fragmented bed; h represents the height of the fragmented bed molten pool calculated in step S2; α represents the thermal diffusivity of the melt; ν represents the kinematic viscosity of the melt; and λ represents the thermal conductivity of the melt.

[0021] Furthermore, in step S4, the heat transfer characteristic scale within the molten pool is obtained through the molten pool characteristic parameter Nusselt number, as follows:

[0022]

[0023]

[0024] In the formula, h upper Indicates the height of the mixing zone above the molten pool; h lower Indicates the height of the layered zone at the bottom of the molten pool; Nu up Nu represents the Nusselt number for upward heat transfer from the molten pool; dn This represents the Nusselt number, indicating the downward heat transfer from the molten pool.

[0025] Furthermore, in step S5, the heat transfer rate of the molten pool is derived based on the energy conservation equations for different heat transfer regions of the molten pool, and the upward, downward, and local heat transfer rates of the molten pool are calculated, as follows:

[0026] v up =αNu up / h-α / h upper (4)

[0027] v dn =αNu dn / h-α / h lower (5)

[0028] v dn,local =v dn (Nu dn,local / Nu dn (6)

[0029] In the formula, v up This represents the upward heat transfer rate of the molten pool; v dn This indicates the downward heat transfer rate of the molten pool; v dn,local Indicates the downward local heat transfer rate of the molten pool; Nu dn,local This represents the Nusselt number, indicating the localized downward heat transfer within the molten pool.

[0030] Furthermore, based on the liquid fraction stored in the fragmented bed region grid in step S1, it is determined whether the region where the grid is located is a solid region, a phase change paste region, or a liquid phase region. The heat transfer rate in the phase change paste region is then corrected to describe the influence of the phase change process on heat transfer. Simultaneously, the direction of the heat transfer rate component is determined based on the sign of the temperature gradient component stored in the fragmented bed region grid in step S1, thereby obtaining the heat transfer rate in the x, y, and z directions for each grid, as detailed below:

[0031]

[0032]

[0033]

[0034] In the formula, u i This represents the heat transfer velocity in the x-direction of grid i; v i w represents the heat transfer velocity in the y-direction of grid i; i f represents the heat transfer velocity in the z-direction of grid i; l,i Represents the liquid fraction at grid i; gradT i,x This represents the component of the temperature gradient in the x-direction of grid i; gradT i,y This represents the component of the temperature gradient in the y-direction of grid i; gradT i,z This represents the component of the temperature gradient in the z-direction of grid i.

[0035] Furthermore, in step S7, the heat transfer rate is replaced with the velocity in the convection term of the energy conservation equation, and the temperature gradient components stored in the fragmented bed region grid in step S1 are used to calculate the explicit convection term as the heat transferred to the boundary due to natural convection of the molten pool, as follows:

[0036] q i,con =-ρ i c p,i (u i gradT i,x +v i gradT i,y +w i gradT i,z (10)

[0037] In the formula, q i,con ρ represents the convective heat transfer obtained by grid i based on the heat transfer rate; i c represents the density of the melt in grid i; p,i This represents the specific heat of the molten material in grid i.

[0038] Further, in step S8, the volume-weighted average of the explicit convection term is calculated so that this portion of heat is subtracted from the entire molten pool in each grid, thereby achieving the energy conservation relationship within the entire molten pool, as detailed below:

[0039]

[0040] In the formula, q ave V represents the volume-weighted average of the explicit convection terms; i,l This represents the volume of grid i located in the liquid phase region or the phase transition paste region.

[0041] Further, in step S9, the explicit convection term and the heat that needs to be uniformly removed from each grid are treated as source terms, and combined with the decay heat power density of the debris bed, the final corrected source terms of the energy conservation equation are calculated, as follows:

[0042] S i =q v +q i,con -q ave (12)

[0043] In the formula, S i This represents the corrected source term for grid i.

[0044] Compared with existing technologies, the calculation method proposed in this invention has the following advantages:

[0045] 1. This invention uses heat transfer rate instead of actual fluid velocity for calculation, omitting the solution of momentum equations during the melting of the fragmented bed, thus greatly reducing the computational load. Simultaneously, the heat transfer rate obtained based on the heat transfer correlation can better calculate the heat dissipation and temperature changes of the molten pool in different directions.

[0046] 2. This invention calculates the heat transfer rate based on existing experimental correlations. Depending on the size of the molten pool, different Ra ranges of heat transfer correlations can be selected, thereby making the calculation results more accurate.

[0047] 3. The simplified calculation method for fragment bed melting based on heat transfer rate proposed in this invention can be conveniently applied to two-dimensional and three-dimensional calculations. Attached Figure Description

[0048] Figure 1 This is a flowchart of a simplified calculation method for fragment bed melting based on heat transfer rate in an embodiment of the present invention.

[0049] Figure 2 This is a schematic diagram of the heat transfer characteristic scale and heat transfer rate within the fragment bed melting area in an embodiment of the present invention.

[0050] Figure 3This invention compares the results of the simplified calculation method based on heat transfer rate in this embodiment with the results of the maximum temperature change over time during the melting process of the fragment bed obtained by the LES method.

[0051] Figure 4 This invention compares the results of the simplified calculation method based on heat transfer rate in this embodiment with the results of the change of heat dissipation in different directions over time during the melting process of the fragment bed obtained by the LES method. Detailed Implementation

[0052] The present invention will be further described below with reference to the accompanying drawings and embodiments.

[0053] Example:

[0054] A simplified calculation method for fragment bed melting based on heat transfer rate, such as Figure 1 As shown, the main steps include:

[0055] S1. Solve the enthalpy conservation equation, calculate the liquid fraction and temperature gradient vector of the fragment bed region grid, and store them in the corresponding grid.

[0056] S2. Calculate the geometric parameters of the melting region of the fragmented bed based on the liquid fraction of each fragmented bed grid. By traversing the grids with a liquid fraction greater than 0, the maximum and minimum values ​​of the grid's horizontal and vertical coordinates can be obtained, thus calculating the height and width of the current melting region.

[0057] S3. Based on the geometric dimensions of the molten pool obtained in step S2, calculate the characteristic parameters of the molten pool, mainly including the Rayleigh number Ra and the Nusselt number Nu. The Rayleigh number Ra is calculated as shown in formula (1). When using the experimental correlation of the molten pool to obtain the Nusselt number Nu, in this embodiment, the heat transfer correlation obtained by the COPRA molten pool experiment is used (L. Zhang, Y. Zhang, Y. Zhou, et al. COPRA experiments on natural convection heat transfer in a volumetrically heated slice pool with high Rayleigh numbers. Annals of nuclear energy, 2016, 87(JAN.PT.2):81-88.), specifically as shown in formulas (2)(3)(4):

[0058]

[0059] Nu up =0.384Ra 0.233 (2)

[0060] Nudn =0.00453Ra 0.32393 (3)

[0061]

[0062] In the formula, Ra represents the Rayleigh number of the molten pool after the fragmented bed begins to melt; Nu up Nu represents the Nusselt number for upward heat transfer from the molten pool; dn Nu represents the Nusselt number, indicating the number of heat transfers downwards from the molten pool. dn,local The Nusselt number represents the downward local heat transfer from the molten pool; g represents the acceleration due to gravity; β represents the coefficient of thermal expansion of the molten pool; q v ν represents the decay heat power density of the fragmented bed; h represents the height of the fragmented bed molten pool calculated in step S2; α represents the thermal diffusivity of the melt; ν represents the kinematic viscosity of the melt; λ represents the thermal conductivity of the melt; and θ represents the radial angle.

[0063] S4, such as Figure 2 As shown, the characteristic scale of heat transfer within the molten pool is calculated using the Nusselt number (Nu), a characteristic parameter of heat transfer within the molten pool, as detailed below:

[0064]

[0065]

[0066] In the formula, h upper Indicates the height of the well-mixed zone at the top of the molten pool; h lower This indicates the height of the layered zone at the bottom of the molten pool.

[0067] S5. The heat transfer rate of the molten pool is derived based on the energy conservation equations for different heat transfer regions of the molten pool, such as... Figure 2 As shown. The upward, downward, and local heat transfer velocities of the molten pool are calculated as follows:

[0068] v up =αNu up / h-α / h upper (7)

[0069] v dn =αNu dn / h-α / h lower (8)

[0070] v dn,local =v dn,local (Nu dn,lcoal / Nu dn (9)

[0071] In the formula, v up This represents the upward heat transfer rate of the molten pool; v dn This indicates the downward heat transfer rate of the molten pool; vdn,local This indicates the rate of localized downward heat transfer within the molten pool.

[0072] S6. Based on the liquid fraction stored in the fragmented bed region grid in step S1, determine whether the region where the grid is located is a solid region, a phase change paste region, or a liquid phase region, and correct the heat transfer rate in the phase change paste region to describe the influence of the phase change process on heat transfer; at the same time, determine the direction of the heat transfer rate component based on the sign of the temperature gradient component stored in the fragmented bed region grid in step S1, thereby obtaining the heat transfer rate in the x, y, and z directions for each grid, as detailed below:

[0073]

[0074]

[0075]

[0076] In the formula, u i This represents the heat transfer velocity in the x-direction of grid i; v i w represents the heat transfer velocity in the y-direction of grid i; i f represents the heat transfer velocity in the z-direction of grid i; l,i θ represents the liquid fraction of grid i; i Represents the radial angle of grid i; gradT i,x This represents the component of the temperature gradient in the x-direction of grid i; gradT i,y This represents the component of the temperature gradient in the y-direction of grid i; gradT i,z This represents the component of the temperature gradient in the z-direction of grid i.

[0077] S7. Replace the velocity in the convection term of the energy conservation equation with the heat transfer rate, and use the temperature gradient components stored in the fragmented bed region grid in step S1 to calculate the explicit convection term as the heat transferred to the boundary due to natural convection of the molten pool, as follows:

[0078] q i,con =-ρ i c p,i (u i gradT i,x +v i gradT i,y +w i gradT i,z (13)

[0079] In the formula, q i,con ρ represents the convective heat transfer obtained by grid i based on the heat transfer rate; i c represents the density of the melt in grid i; p,i This represents the specific heat of the molten material in grid i.

[0080] S8. Calculate the volume-weighted average of the explicit convection terms so that this portion of heat in each grid can be subtracted from the entire molten pool, thereby achieving the energy conservation relationship within the entire molten pool, as follows:

[0081]

[0082] In the formula, q ave V represents the volume-weighted average of the explicit convection terms; i,l This represents the volume of grid i located in the liquid phase region or the phase transition paste region.

[0083] S9. Treating the explicit convection term and the heat that needs to be uniformly removed from each grid as source terms, and combining this with the decay heat power density of the debris bed, the final corrected source terms of the energy conservation equation are calculated, as follows:

[0084] S i =q v +q i,con -q ave (15)

[0085] In the formula, S i This represents the corrected source term for grid i.

[0086] S10. Iteratively calculate the energy conservation equation to obtain the latest temperature distribution and liquid fraction distribution, and complete the calculation of the fragment bed melting process within a unit time step.

[0087] To verify the effectiveness of the simplified calculation method for fragment bed melting based on heat transfer rate proposed in this invention, the traditional LES method and the simplified calculation method for fragment bed melting based on heat transfer rate proposed in this invention were used to calculate the fragment bed melting process in the lower head of the pressure vessel.

[0088] Figure 3 The results of calculating the maximum temperature change over time during the melting process of a fragment bed are compared between the traditional LES method and the method proposed in this invention. It can be seen that the transient temperature changes calculated by both methods agree well.

[0089] Figure 4 This paper compares the calculation results of heat dissipation over time in different directions during the melting of a fragmented bed using the traditional LES method and the method proposed in this invention. It can be seen that the heat dissipation values ​​simulated using the method proposed in this invention are in good overall agreement with the values ​​calculated by the traditional LES method, with a maximum deviation of less than 10%, which is within an acceptable range.

[0090] The traditional LES method requires over 30 days to calculate the transient process of a bed of debris melting over 6 hours, while the simplified method proposed in this invention requires only 1 day. This demonstrates that the simplified calculation method significantly improves the computational efficiency of the bed of debris melting process while maintaining good accuracy.

Claims

1. A simplified calculation method for fragment bed melting based on heat transfer rate, characterized in that: Includes the following steps: S1. Solve the enthalpy conservation equation, calculate the liquid fraction and temperature gradient vector of the fragment bed region grid, and store them in the corresponding grid. S2. Calculate the geometric dimensions of the molten pool formed by the melting of the fragmented bed, and obtain the height and width of the molten pool; S3. Calculate the characteristic parameters of the molten pool in the fragmented bed melting region; S4. Calculate the characteristic scale of heat transfer within the molten pool; S5. Calculate the heat transfer rate in the liquid phase region where the fragmented bed grid is located; S6. Calculate the heat transfer velocity components of the grid and the heat transfer velocity of the mushy region; S7. Calculate the explicit convection term as the heat transferred to the boundary due to natural convection of the molten pool; S8. Explicit convection terms are uniformly removed to achieve energy conservation throughout the entire fragmented bed molten pool; S9. Calculate the source terms of the energy conservation equation after correction; S10. Iteratively calculate the energy conservation equation to obtain the latest temperature distribution and liquid fraction distribution, and complete the calculation of the fragment bed melting process within a unit time step. In step S5, the heat transfer rate of the molten pool is derived based on the energy conservation equations for different heat transfer regions of the molten pool, and the upward, downward, and local heat transfer rates of the molten pool are calculated, as follows: （4） （5） （6） In the formula, This indicates the upward heat transfer rate of the molten pool; This indicates the rate of downward heat transfer from the molten pool; This indicates the rate of localized downward heat transfer within the molten pool; The Nusselt number represents the number of localized downward heat transfer within the molten pool. Indicates the height of the mixing zone above the molten pool; Indicates the height of the layered zone at the bottom of the molten pool; This indicates the Nusselt number representing the upward heat transfer from the molten pool; The Nusselt number represents the number of heat transferred downwards from the molten pool. In step S6, the liquid fraction stored in the fragmented bed region grid in step S1 is used to determine whether the region where the grid is located is a solid region, a phase change paste region, or a liquid phase region. The heat transfer rate in the phase change paste region is then corrected to describe the influence of the phase change process on heat transfer. Simultaneously, the sign of the temperature gradient component stored in the fragmented bed region grid in step S1 is used to determine the direction of the heat transfer rate component, thereby obtaining the heat transfer rate in the x, y, and z directions for each grid, as detailed below: （7） （8） （9） In the formula, Represents a grid i Heat transfer rate in the x-direction; Represents a grid i Heat transfer rate in the y-direction; Represents a grid i The heat transfer rate in the z-direction; Represents a grid i The liquid fraction; Represents a grid i The component of the temperature gradient in the x-direction; Represents a grid i The component of the temperature gradient in the y-direction; Represents a grid i The component of the temperature gradient in the z-direction; Represents a grid i The radial angle; In step S7, the heat transfer rate is replaced with the velocity in the convection term of the energy conservation equation, and the temperature gradient components stored in the fragmented bed region grid in step S1 are used to calculate the explicit convection term as the heat transferred to the boundary due to natural convection of the molten pool, as follows: （10） In the formula, Represents a grid i Convective heat transfer based on heat transfer rate; Represents a grid i Medium melt density; Represents a grid i Specific heat of the molten material.

2. The simplified calculation method for fragment bed melting based on heat transfer rate according to claim 1, characterized in that: In step S2, the geometric parameters of the melting region of the fragment bed are calculated based on the liquid fraction of each fragment bed region grid. By traversing the grids where the liquid fraction is greater than 0, the maximum and minimum values ​​of the grid's horizontal and vertical coordinates can be obtained, thus enabling the calculation of the height and width of the current melting region of the fragment bed.

3. The simplified calculation method for fragment bed melting based on heat transfer rate according to claim 1, characterized in that: In step S3, based on the height of the molten pool in the fragmented bed obtained in step S2, the characteristic parameters of the molten pool are calculated. The characteristic parameters of the molten pool include the Rayleigh number Ra and the Nussel number Nu. The Rayleigh number Ra is calculated as shown in formula (1). When calculating the Nussel number Nu, the heat transfer correlation obtained from the molten pool experiment is used. （1） In the formula, This indicates the Rayleigh number of the molten pool after the fragment bed begins to melt; Represents gravitational acceleration; Indicates the coefficient of thermal expansion of the molten pool; This represents the decay heat power density of the fragment bed; This represents the height of the molten pool in the fragment bed calculated in step S2; Indicates the thermal diffusivity of the molten material; Indicates the kinematic viscosity of the melt; This represents the thermal conductivity of the molten material.

4. The simplified calculation method for fragment bed melting based on heat transfer rate according to claim 1, characterized in that: In step S4, the heat transfer characteristic scale within the molten pool is obtained through the molten pool characteristic parameter Nusselt number, as follows: （2） （3） In the formula, Indicates the height of the mixing zone above the molten pool; Indicates the height of the layered zone at the bottom of the molten pool; This indicates the Nusselt number representing the upward heat transfer from the molten pool; This represents the Nusselt number, indicating the downward heat transfer from the molten pool.

5. The simplified calculation method for fragment bed melting based on heat transfer rate according to claim 1, characterized in that: In step S8, the volume-weighted average of the explicit convection term is calculated so that this portion of heat is subtracted from the entire molten pool in each grid, thereby achieving the energy conservation relationship within the entire molten pool, as detailed below: （11） In the formula, This represents the volume-weighted average of the explicit convection terms; Mesh representing the liquid phase or phase transition paste region i The volume.

6. The simplified calculation method for fragment bed melting based on heat transfer rate according to claim 1, characterized in that: In step S9, the explicit convection term and the heat that needs to be uniformly removed from each grid are treated as source terms. Combined with the decay heat power density of the debris bed, the final corrected source terms of the energy conservation equation are calculated, as follows: （12） In the formula, Represents a grid i The corrected source term.

Images