Junying Xu ¹ Xiaomeng Dong ^2* Huiyong Zhang ¹ Xiangyu Yun ¹ Lei Zhang ¹ Yuzhao Ye ² 刘 湘 ² Jinhong Mo ² Ming Yang ²

• ¹ China Nuclear Power Technology Research Institute Co., Ltd, Shenzhen, China
• ² Shenzhen University, Shenzhen, China

After the severe accident of nuclear reactor occurs, the fuel element may melt if the decay heat is not removed in time. The molten corium is expected to be ejected into the lower head of the Reactor Pressure Vessel (RPV) and reacts violently with the residual coolant. In some extremely conditions, the reactor corium is expected to discharge from the RPV into a deep-water pool in the cavity below the RPV, which is called the reactor pit. As a result, a debris bed is formed after the initial quenching. The debris bed is a porous structure composed of particles in different sizes. Due to the continuous fission reaction in the debris bed, it continues to release decay heat in the water pool. If no useful cooling method is carried out, the decay heat is large enough to result in the re-melting of the debris bed and a further damage to the RPV structure. In general, the pool boiling of residual coolant in the lower head is the sufficient way to remove the generated decay heat. However, the boiling phenomenon is limited by the Critical Heat Flux (CHF) which also brings potential threat to the debris bed and associated structure. So, it is of great significance to investigate the pool boiling phenomenon and the CHF in the porous debris bed after severe accident.In general, the boiling crisis in debris bed is of 'dry-out' type because of the not very high decay heat and the porous structure. So, the dryout heat flux (DHF) or limited power density is regarded as the most important parameter to describe the heat loading capacity of the porous debris bed. Furthermore, the flow resistance is another key parameter which has a great influence on the formation of natural convection in porous media. The value and position of DHF are also changed according to the strength of natural convection.After a wide literature review, numerous researches have been found around the two key parameters in porous media. Without heating, the mechanistic studies are carried out around flow resistance by many scientists in simple geometry. The first mechanistic model is developed by Ergun (Ergun, 1952) for the single-phase flow through packed columns. The correlation (equation (1-1)) among the pressure drop, superficial velocity and structure of porous media is established through theoretical analysis and data fitting. (1-1)Where parameter K and η are called permeability and passability, respectively.They are the key parameters to describe the structure of porous media. They are calculated by equation (1-2). (1-2)Where the □□ □□ and □□ □□ refer to the effective particle diameter and porosity respectively. When it comes to two-phase boiling flow, the Ergun equation is extended through including the effect of relative permeability and relative passability. The model is provided by Lipinski (Lipinski, R.J., 1982) and many researchers (Reed, 1982;Schulenberg and Muller, 1987;Hu and Theofanous, 1991) provides lots of transformation by using different correlations. All the modifications are calculated according to the void fraction of fluid and vapor. As a validation, experiments of Li et al. (Liangxing Li, 2015) proves that Reed's model has a high fidelity for the two-phase flow resistance. In the validation (Huang, Z., 2018) of 2D thermal hydraulic code MEWA, Reed's model (Reed, 1982) shows best performance in the validation of frictional pressure drop while Schulenberg&Müller's model (Schulenberg and Muller, 1987) is more suitable for the prediction of DHF. However, it is still necessary to do the validation and verification in the research.Apart from the flow resistance, numbers of experiments and numerical simulations (Trenberth and Stevens, 1980;Barleon and Werle, 1981;Lipinski, 1982;Hofmann, 1984) have been carried out around the debris cooling ability in the past thirty years. Among all the factors, properties of structure, shape of outer contour, velocity of fluid and operating pressure are the hot points of the research on porous debris bed. In general, the properties of structure refer to the porosity and effective particle diameter which are the main parameters to describe the porous debris bed. Many experiments (Kulkarni, P.P., 2010;Leininger, S., 2014;Takasuo, E., 2016;Gourbil, A., 2019) and simulations (Chakravarty, A., 2020) have shown that the value of DHF increases along with porosity and effective particle diameter. However, most of the researches used uniform values which is difficult to gain a deeper comprehension on the relationship between DHF and structural properties. In a heat-generating porous media, boiling crisis always occurs on the top of the geometry. So, non-uniform distributions are studied to investigate the determining factor of DHF in this paper.During the process of reactor corium ejection, kinds of debris bed's outer contour may form due to the violent effect with water. In the classical analyses, the conical, cylindric and truncated cone shaped debris bed are widely investigated by many researchers (Takasuo, E., 2016;Chakravarty, A., 2020) through experiments and simulations. In these investigation, volume of the debris bed should be kept consistent to guarantee the same decay heat density level while the other parameters, such as height, basal area and surface area, could be changed to investigate the effect on DHF. However, the shape of realistic debris bed is seldom investigated in existing researches. According to the realistic physical process, considering the influence of molten corium injection, the volcanic type bed with a downcomer is possible in severe accident. W.M. Ma et al. (Huang, Z., 2018) have studied the dryout phenomenon in the debris bed with downcomer by using the 2D thermal-hydraulic code MEWA. The relative error of dryout power density is under 20%. In addition, a smaller scale of downcomer delivers a larger error. Research on the shape of debris bed also shows that the dryout power density is inversely proportional to the bed's height regardless of the bed's shape. The MEWA code is a system code whose accuracy is still not enough for the analysis. For the debris bed with complex shape, Computational Fluid Dynamics (CFD) method is more suitable (Wang, M.J., 2021).Except for the value of DHF, the position where boiling crisis occurs is also an important parameter in the safety analysis of porous debris bed. For the realistic severe accident process, it is difficult to determine the true position of boiling crisis. After a wide investigation, the experimental data (Takasuo, E., 2016) shows that the boiling crisis always occurs close to the center of the geometry while results of the numerical simulation (Chakravarty, A., 2020) indicates that the position locates at the edge of the top surface. No explanation is given in the existing literature. This phenomenon is not obvious for the conical porous debris bed but it influenced a lot for the cylinder type. One of the differences between experiment and simulation is the heating-process. For the experiment, the DHF is obtained through the 'step-wise heating', i.e., the heating power is increased step by step until the temperature excursion. In contrast, the heating power is usually set as a constant value directly to check the dryout phenomenon in the simulation. However, the influence of different heating process is not clear which should be investigated to give more useful information for the analysis of severe accident in PWR.In this paper, the two-phase Eulerian model with internal heat source is established for porous medium. The additional correlations for the flow resistance, heat transfer coefficient and mass exchange are also provided for the calculation of pressure drop and limited power density. All the calculation are carried out based on the CFD method which could give a detailed distribution for each parameter. Results are discussed around the effect of heating type, non-uniform distribution of structural parameters, and shape of geometry. In this paper, the debris bed is simplified as a porous media which consists of many solid particles with internal heat source. The whole porous media is immersed in the water pool to simulate the cooling process after severe accident which is shown in Fig. 1. During all the boiling process, the solid particles should be kept in the original positions although they are impacted by the water flow. The porous media is fully-flooded by the water and pool boiling occurs under the effect of decay heat. So, three phases, solid, liquid and vapor, exists in the geometry zone. Porous region (debris bed)Fig. 1 the geometry schematic of porous media Eulerian two-fluid model is used to describe the fluid phase which consists of mass, momentum and energy governing equations. The equations are nearly the same for the fluid in the continuous fluid region and porous region. The difference is from the effect of the porosity in porous region. For the solid phase, additional energy equation is applied to solve the temperature distribution of solid particles. The governing equations are shown below.For the fluid in continuous fluid region, the governing equations are shown as below.(1) mass conservation equation (2-10)(6) energy conservation equation □□ □□□□ (□□ □□ □□ □□ ℎ □□ ) + □□ ⋅ (□□ □□ □□ □□ □□ □□ ℎ □□ ) = □□ □□ □□ □□ □□ □□ ∇ 2 □□ □□ + □□ □□□□ + □□ □□□□ + ∆□□ □□□□ ℎ □□□□ (2-11) □□ □□□□ (□□ □□ □□ □□ ℎ □□ ) + □□ ⋅ (□□ □□ □□ □□ □□ □□ ℎ □□ ) = □□ □□ □□ □□ □□ □□ ∇ 2 □□ □□ + □□ □□□□ + □□ □□□□ + ∆□□ □□□□ ℎ □□□□ (2-12)In addition, the energy equation for the solid phase is described below. It is assumed that the solid phase would not melt or move during the process. So, the mass and momentum conservation equations are neglected.□□((1-□□ □□ )□□ □□ □□ □□,□□ □□ □□ ) □□□□ = (1 -□□ □□ )□□ □□ ∇ 2 □□ □□ + □□ -□□ □□□□ -□□ □□□□ (2-13)2.2 Closure equations Based on the basic conservation equations, additional correlations should be given for the solution of the interfacial mass, momentum and energy exchange associated with the phenomena. The correlations include the model of flow resistance for the porous media, the heat exchange among solid phase, liquid phase and vapor phase, and the mass exchange between liquid and vapor. They are discussed in detail in the following sections.In the multiphase flow, the flow regime is an important parameter in the correlation. In the present study, fluid flow within the porous bed as well as the clear fluid region is assumed to be divided into three different flow regimes. Depending on volume fraction of the fluid phases, they are described as liquid continuous (or bubbly) regime, transition regime and vapor continuous (or droplet) regime. In this paper, the divided rule of Koushik Ghosh (Chakravarty, A., 2020) is used which is shown in Table 1. In the following correlations, the subscript j inside refers to the continuous phase while k represents the dispersed phase.Table 1 The main source of flow resistance considered in this paper is the interfacial drag force among the three phases. It is different in different regions. For the continuous fluid region, it only exists the interfacial drag force between water and vapor which comes from the velocity difference of two phases. For the porous region, the drag forces between solid particle and fluids occupy the major fraction of flow resistance other than the interfacial drag force.(1) Continuous fluid region In the continuous fluid region, the interfacial momentum exchange □□ □□□□ is determined by the following equations.□□ □□□□ = □□ □□ □□ 6□□ □□ □□ □□ □□ □□□□ ∆□□ (2-14)In the equation (2-14), the subscript k refers to the dispersed phase. □□ □□ is the effective diameter of the dispersed phase. ∆□□ refers to the velocity difference of two phases. □□ □□ represents the relaxation time and is defined as□□ □□ = □□ □□ □□ □□ 2 18□□ □□ (2-15)The interfacial drag coefficient f is determined by the model of Schiller and Naumann (Schiller, L., 1935) which is the widely-used in the analysis of boiling phenomenon. It is a usually defaults model and adequate for all fluid-fluid phases in the CFD codes (Khan, I., 2020). The expression is shown as equation (2-16). □□ □□ □□□□ □□ (2-16) Where the coefficient □□ □□ is expressed as□□ □□ = { 24(1 + 0.15□□□□ □□ 0.687 ) □□□□ □□ 0.44 □□□□ □□ ≤ 1000 □□□□ □□ ≥ 1000(2) Porous region For the porous region, the existing three phases would generate three kinds of force which include the interfacial drag between the solid phase and the saturating fluid phases as well as the interfacial drag between the liquid phase and vapor phase. All the drag forces are described by the terms □□ □□□□ and □□ □□□□ in the equation (2-9) and (2-10).After a wide literature review, the interfacial drag between the saturating fluid phases is ignored by many researchers except Schulenberg and Muller (Schulenberg and Muller, 1987). They considered the interfacial drag force □□ □□□□ to be invariant with flow regime change. The expression is shown in equation (2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16)(17).□□ □□□□ = 350□□ □□ 7 □□ □□ □□ □□ (□□ □□ -□□ □□ ) □□□□ □□□□ ( □□ □□ □□ □□ - □□ □□ □□ □□ ) 2 (2-17)Compared to the interfacial drag force, the drag forces □□ □□□□ between the solid phase and fluid phase play a more significant role in the simulation. The value is determined by the permeability K and passability η. The expression of □□ □□□□ is described as□□ □□□□ = -□□ □□ □□ □□ ( □□ □□ □□□□ □□.□□ □□ □□ + □□ □□ □□□□ □□,□□ |□□ □□ |□□ □□ ) (2-18)Where the permeability K and passability η are calculated according to equation (1-2).The coefficient of correction □□ □□.□□ and □□ □□,□□ is determined by different models which are shown in Table 2. According to the conclusion of Weimin Ma (Huang, Z., 2018), Reed's model shows the best performance in the validation of frictional pressure drop while Schulenberg&Müller's model is more suitable for the prediction of DHF. So, it is necessary to do another validation for the model of flow resistance in porous media. During the cooling process after severe accident, the major heat source is the decay heat from the debris bed. The heat exchange exists among the three phases and different mechanisms are taken into account. In the continuous fluid region, the interfacial liquidvapor heat transfer is the main type of heat exchanging. In the porous region, the convective heat transfer and boiling heat transfer between solid phase and fluid phases should be added as the major types of heat exchanging. All of these mechanisms are discussed in the following sections with the correlations for the calculations.(1) Continuous fluid region In the continuous fluid region, the interfacial heat transfer is especially important for the two-phase flow with 0 < □□ □□ < 1. It is described as□□ □□□□ = □□ □□□□ ℎ □□□□ (□□ □□ -□□ □□ ) (2-19)Where the subscript j refers to the continuous phase. The heat transfer coefficient ℎ □□□□ is determined based on the existing flow regime.For the liquid continuous regime ℎ □□□□ = (2 + 0.6□□□□ □□ 0.5 □□□□ □□ 0.33 ) □□ □□ □□ □□ For the vapor continuous regime ℎ □□□□ = (2 + 0.6□□□□ □□ 0.5 □□□□ □□ 0.33 ) □□ □□ □□ □□ For the transition regime, a linear interpolation between continuous and vapor continuous regimes is used which is ℎ □□□□ = (1 -□□)ℎ □□□□ + □□ℎ □□□□ Where W is a factor of weight which is calculated according to the void fraction and its limitation in Table 1 (Chakravarty, A., 2020).(2) Porous region For the interfacial liquid-vapor heat transfer, the model in porous region is the same as that in the continuous fluid region. However, it occupies small fraction of the heat exchange in porous region. The heat transfer between solid phase and fluid phases are the major source which are generated from the decay heat. This heat exchange can be divided to two types which are convective heat transfer and boiling heat transfer. Both of these two types occur on the surface of the porous debris bed. But difference exists when the flow regimes vary. The heat transfer is estimated as below □□ □□□□,□□ = □□ □□□□ ℎ □□□□,□□ (□□ □□ -□□ □□ )(2-20) The heat transfer coefficient is determined by using the Ranz and Marshall (Ranz, W.E., 1952) correlation which can be expressed as ℎ □□□□,□□ = (2 + 0.6□□□□ □□ 0.5 □□□□ □□ 0.33 )□□ □□ □□ □□ (2-21)The coefficient of convective heat transfer above is only used for the continuous fluid phase while the □□□□ □□ and □□ □□ are calculated by using the parameters of dispersed phase. The convective coefficient between solid phase and dispersed phase is assumed to be zero in the fluid continuous regime.In the transition regime, the heat transfer coefficient is evaluated by a linear interpolation between the nucleate and film boiling regimes using a weighing parameter w. The coefficient correlation proposed by Rohsenow (Rohsenow, W.M., 1951) and Bromley (Bromley, L., 1950) For the coefficient ℎ □□□□,□□ of the vapor phase, the thermal properties of liquid are replaced by that of equation .Except the convective heat transfer, the boiling heat transfer is more important in which the influence of latent heat is far more than convective heat transfer in the boiling process. In this paper, the method of RPI model (Kurul N., 1991) is used which is described as□□ □□□□,□□ = □□ □□ □□ □□ □□ □□ ℎ □□□□ □□ (2-23)Where □□ □□ is the volume of the bubble based on the bubble departure diameter, □□ □□ is the active nucleate site density, and ℎ □□□□ is the latent heat of evaporation, and □□ is the bubble departure frequency.Because of the narrow space in the porous region, the quenching heat flux is neglected. For the boiling phenomenon under low pressure, the bubble departure diameter is always 1mm which is larger than the scale of flow space in the porous region. So, the action of bubble lift-off is suppressed, as well as the quenching heat flux. The interfacial mass exchange only occurs between the liquid and vapor phases. It is assessed by the boiling heat transfer and interfacial heat transfer. The former value is used for the vaporization and the latter one is used for the condensation. The expression is described as∆□□ □□□□ = ∆□□ □□□□ = □□ □□□□,□□ +□□ □□□□ ℎ □□□□(2-24) The calculations are carried out based on the commercial computational fluid dynamics platform ANSYS FLUENT. The simplified geometry and boundary conditions are shown in Fig. 2. The investigated geometry is simplified to 2-D axissymmetry body and the axis is set to X-direction because of the rule of platform. The bottom and outer wall surface are considered as adiabatic wall and the upper water surface is set as a constant pressure boundary. The shape of porous region varies according to the setting which could be cylindrical, conical and truncated conical bodies. The Eulerian multiphase model is utilized for the solution of governing conservation equations. Furthermore, the packed bed method is used to integrate the above settings. In detail, the porous debris bed is set as packed bed whose velocity is limited to zero strictly. The volume fraction of solid phase is set to 1 -□□ □□ through the PATCH setting.The RPI boiling model is used as the basic method for the calculation of porous media. The mentioned heat coefficients among three phases are utilized and coupled to the corresponding transport equations by using the User-Defined Functions (UDFs).For the boiling heat transfer in equation (2-23), the models provided by platform are used directly which include Tolubinski-Kostanchuk model (Tolubinski, V.i., 1970) for the bubble departure diameter, Cole model (Cole, R., 1960) for the bubble departure frequency, Lemmert-chawla model (Lemmert, M., 1977) for the nucleation site density and Delvalle-kenning model (Del, V.V.H., 1985) for the area influence coefficient. The symmetric model is used for the interfacial area □□ □□□□ among the three phases.In the correlation related to flow resistance, SST k-ω model is used for the calculation of turbulence in the present study. The integrated drag model of Schiller-Naumann is utilized directly for the continuous fluid region. In the porous region, the additional force □□ □□□□ is calculated through the UDF of modification. The drag effect between solid phase and fluid phases are calculated by the UDFs in the porous zone. The parameter □□ and □□ □□.□□ are used to calculate the viscous resistance while the parameter □□ and □□ □□,□□ are used for the inertial resistance. In addition, another modification should be applied for the flow resistance □□ □□□□ in porous region because of the porosity and the related superficial velocity of fluid phases. These modifications are loaded through the UDFs of □□ □□.□□ and □□ □□,□□ in equation (2-18). In general, the pressure drop and limited power density (or dryout heat flux) are the key parameters to be used to validate the accuracy of the developed model (Liangxing Li, 2015). In the porous region, bubble generates at the whole-body surface and aggregates into the continuous vapor regime. When the void fraction exceeds the limited value, dryout phenomenon occurs which bring a potential damage to the debris bed. However, the process of bubble aggregation is influenced by the flow resistance. The equilibrium relationship between bubble generation, movement and aggregation is the main factor to determine the limited power density.3.1 grid-independent validation Before the other validations, grid-independent calculation has been carried out to select the best grid scale. The bed 3 of DEBECO facility (Liangxing Li, 2015) is used as the test case which considers isothermal air-water flow. It is simplified as a 2-D geometry with 60mm × 600mm with the porosity of 0.4 and the effective particle diameter of 1.44mm. The material properties assumed for the solid phase are taken from Takasuo et al. (Takasuo, E., 2016) ( □□ □□ = 4200 kg/m 3 , □□ □□ = 2 W/mK , □□ □□,□□ = 775 J/kg. K). The test matrix of structured grid includes 0.3mm, 0.5mm, 1mm, 2mm and 3mm. Result of Fig. 3 shows that the pressure drop is nearly consistent when the grid scale is lower than 1mm. So, 1mm is chosen as the grid scale which could meet the requirement of accuracy. Pressure drop is the main parameter to reflect the flow resistance. The same bed of experimental facility (DEBECO, (Liangxing Li, 2015)) is used for the validation. In addition, single-phase and air-water two-phase flow are validated among the experimental data, theoretical results and CFD calculations. For the single-phase flow, Ergun model (Ergun, 1952) which is described in equation (1-1) is used for the simulation. The compared results are shown in Fig. 4 under different superficial velocity of liquid. The analytic result is gained through solve the equation (1-1). Results have shown good agreement with the relative error under 0.1% for the single-phase flow resistance. Fig. 5 The comparison of results for the two-phase flow resistance Compared to the experimental data, the analytic result of Reed's model shows the best performance while the CFD results of S&M and H&T's models are proximate. All the relative errors are less than 13% which meet the requirement of accuracy. For the validation of heat and mass transfer correlation, a fully-flooded porous media from the COOLOCE facility (Takasuo, E., 2016) is used as the integrated case. According to the experimental facility, the scale of whole zone could be simplified to 306.5mm × 500mm. Inside, the porous region is a 2D geometry with 152.5mm × 269.6mm which means a cylinder-type porous media. The porosity is set to 0.392 and the effective particle diameter is 0.97mm. This case operates under 1.3bar pressure which is a bit more than the atmospheric pressure.During the experimental process, a step-wise heating method is utilized to determine the limited power density. The heating power is set to a relative low value at the initial condition. The internal thermocouples would feedback the temperature at the measured points whose growth determines the occurrence of the dryout phenomenon. If no temperature rise is monitored, the heating power would increase step by step every five minutes. Through this method, the measured power level is 34.1kW in the experiment while the limited power density is 1728.6kW/m 3 . Based on the ANSYS FLUENT platform, the same method is used in the calculation to determine the limited power density. The initial power density to set to be 25% lower than experimental value. The maximum temperature of porous media is also monitored at the same time. The limited power density is determined when the obvious increase of temperature occurs. Under the four correlations of flow resistance, the calculated limited power density is shown in Table 3.Table 3 1.24% According to the validation above, results show that Schulenberg&Muller drag model have both good performance in the calculation of flow resistance and limited power density. So, it is chosen as the man drag model in the following discussion. The validation also proves a certain accuracy of the proposed method of heat and mass transfer which is suitable for the analysis of pool boiling in porous media. In the following discussion, the influence of natural convection, structural parameters and shape are analyzed based on the proposed model above. The same fullyflooded porous media from the Takasuo's facility (Takasuo, E., 2016) is used as the main research object in this section. In addition, the determined factor of limited power density is also discussed under different boundary conditions and settings.4.1 The effect of heating type The literature review shows that the position where boiling crisis occurs is difficult to determine and non-general conclusion is provided. the experimental data (Takasuo, E., 2016) showed that the boiling crisis always occurs close to the center of the geometry while results of the numerical simulation (Chakravarty, A., 2020) indicated that the position located at the edge of the top surface. So in this section, the position of maximum temperature is investigated based on the proposed model.After comparing the experiment and existing simulation, the heating type is found to be the main difference during the research process. In general, the step-wise heating method is used in the experiment. However, a given value of heating source is used in the simulation directly. If it is not the limited power density, all the setting of calculation would be initialized before the next step. Fig. 6 shows the variation of temperature distribution under two different heating types. Along with the time, the position of maximum temperature suffers a movement from the edge of top surface to the center under the setting of the direct heating source. However, the maximum temperature appears at the same position under the setting of step-wise heating source. Fig. 6 The time variation of temperature under two heating types In order to investigate the reason of the location migration, Fig. 7 and Fig. 8 shows the void fraction and velocity vectors of liquid and vapor separately which are at the same moment with Fig. 6. It is concluded that the velocity vectors are of great significance in the distribution of vapor phase and temperature. In the setting of direct heating source, the bubble moves along the vertical direction under the effect of buoyancy force at the initial time. However, it is difficult to form a strong natural convection because of the low temperature difference and high flow resistance in porous region. So, the generated bubbles would gather at the edge of the top surface which may lead to a potential boiling crisis. Along with the develop of natural convection, the flow towards to center of the porous region is enhanced which push the bubbles move to the center. When the bubble aggregation occurs at the center, the high temperature appears as a consequence.The phenomenon is different in the setting of step-wise heating type. The initial power level is not high enough to bring about dryout phenomenon. But the waiting time (not less than five minutes) is enough to form natural convection at a certain level. When the additional power is added step by step, the generated bubbles would have a high radius velocity and move to the center of porous region. So, the bubble aggregation would occur at the center at the first time which deliver the maximum temperature. All in all, the strengthen of natural convection is the determined factor on the position of boiling crisis. It is also the reason for the different performance of two heating types during the process of realistic severe accident, Fig. 7 The time variation of vapor phase fraction under two heating types Fig. 8 The time variation of velocity vector under two heating types 4.2 The effect of operating conditions Under different boundary conditions, the limited power density may have a large difference. The effect of pressure and structural parameters are investigated in this section. Furthermore, the non-uniform distributions of structural parameters are also studied which is more similar to the real porous debris bed.(1) Operating pressure The operating pressure is an important parameter during the severe accident process. The realistic pressure of debris bed may in the range of 0.1MPa to several atmospheres which depends on the realistic process. So, several conditions should be investigated to study the relationship between the limited power density and pressure. Fig. 9 shows the trend of limited power density under different pressure. It is concluded that a linear correlation exists between the limited power density and pressure. This trend is also consistent with the experimental value. It proves a confidence of the proposed model again with a relative error lower than 5%. Fig. 9 The trend of limited power density under different pressure(2) Structural parameters The structural parameters discussed in this section refer to the porosity and effective particle diameter after the formation of debris bed. In the existing literatures, uniform distributions are considered for these structural parameters under different conditions. From the view of Solid Geometry, the porosity has a minimum and maximum value if the uniform spherical particles are assumed to contact with each other in all three dimensions. The two contact modes are shown in Fig. 10. According to the mathematic calculation, the minimum and maximum porosity are approximate 0.108 and 1 -π/6, respectively. So, the range from 0.2 to 0.5 is discussed the following content. For the effective particle diameter, a test is carried out for the potential range. It is found that too large diameter is not suitable for the discussion because non limited power density is gained. So, the range from 0.8mm to 1.1mm is used to investigate the influence of diameter.(a) minimum porosity (b) maximum porosity Fig. 10 The solid geometry of minimum and maximum porosity Fig. 11 shows the trend of limited power density under different structural parameters. Higher limited power density is affordable when the porosity increases. The analysis of results also shows the great effect on the reduction of flow resistance when the porosity increases. The flow passability of fluid increases under the lower flow resistance which improves the heat transfer coefficient upon the heat surface of debris bed. It is also difficult for the bubble to accumulate at the top of the porous region which have a positive influence on the heat removal ability. In addition, the rise is more obvious when the porosity increases from 0.4 to 0.5 which is very approach to the limited value of porosity. This growth phenomenon is consistent with the result of Koushik Ghosh et al. (Chakravarty, A., 2020). However, a higher porosity than 0.5 is difficult to be used in this geometry and no limited power density is gained based on the proposed model. Fig. 11b shows the effect of particle diameter. During the severe accident process, the whole mass of porous debris remains unchanged. The volume of the debris bed is also the same if the density of reactor corium is assumed constant. So lower limited power density or worse cool ability is found under the condition of smaller particles, which means too many fragments and pellets are formed after the vigorous react between corium and coolant.(3) Non-uniform distribution of effective particle diameter During the process of realistic severe accident, it may exist multi-scale fragments at different regions which leads to a variation of structural parameters. So, the nonuniform distributions of structural parameters are also discussed in this section. Simple cases with difference from top to bottom are used to investigate the key factor of limited power density. As shown in Fig. 12a, the structural parameter in porous region 2 remain 1.0mm in all the cases while that in porous region 1 varies from 0.9mm to 1.1mm. Compared with the value in Fig. 11b, the results show that the limited power density is determined by the structural parameter of top half geometry rather than the whole region. (a) the geometry (b) effective particle diameter Fig. 12 Limited power density under non-uniform distribution of effective particle diameter Porosity is also an important variable parameter in the analysis of porous debris bed. However, the influence of inhomogeneous distribution is not easy to investigate through the simulation. In the ANSYS FLUENT platform, the porosity is calculated based on the equation (4-1). □□ inside refers to the number of particles.□□ □□ = 1 - □□□□ □□□□□□□□□□□□□□□□ □□ □□□□□□□□□□□□ □□□□□□ (4-1)Based on the assumption that the whole mass of reactor corium and its composition are constant, the volume of all the solid phase would not change during the process of severe accident. So, the porosity is determined by the volume of space occupied by debris bed. So, the shape and volume of debris bed are investigated in the section 4.3 which could also reflect the effect of porosity.4.3 The effect of shape During the process of severe accident, kinds of shape may form based on the reactor corium of same mass. In this paper, three kinds of heap-like geometry are used for the investigation which include conical, truncated conical and cylindrical geometries. During the comparison, the volume of debris bed is kept constant which may influence the total decay heat. Fig. 13 shows the limited power density under different boundary setting. For the truncated cone, the geometries with same height, same bottom radius and same top radius are investigated respectively. For the conical geometry, the debris bed with same height is discussed while the height of case with same bottom radius is too higher than the total scale. In addition, the cylindrical geometry with downcomer is also discussed in this section. The comparison of limited power density is shown in Fig. 13. It is concluded that the conical geometry has the minimum limited power density which is far lower than the others, which can be explained by the aggregation of bubbles in the tip area. Fig. 14 shows the distribution of temperature and void fraction. The maximum temperature appears in the top of the conical geometry which is also the remelting start position. The cylinder with downcomer has a stronger cooling ability than that of the other geometries. As a result, it is concluded that the volcanic type with a downcomer would be safer during the realistic process. Compared the results of different geometries with same height, we can also know that the height is not the unique factor to determine the limited power density, even though the boundary conditions and structural parameters are kept the same. In order to investigate the cooling ability of porous debris bed after the severe accident, the two-phase conservation equations with closure correlations are proposed for the pool boiling phenomenon inside. All the calculation are carried out based on the CFD method. The analysis and discussion are around the effect of heating type, nonuniform distribution of structural parameters, and shape of geometry. Result shows the key effect of natural convection between the different boundary settings of heating type. The time series of strong natural convection formation and decay power heat are the determined factor for the position of boiling crisis. In addition, the limited power density is determined by the top half of debris bed. The increase of structural parameters and operating pressure would lead to a better cooling ability. For the shape of debris bed, regular cylinder is a better structure for the heat removal while the conical shape would reduce the limited power density a lot which is dangerous for the long-term cooling of the debris bed. The cooling ability would be better if a downcomer exists in the porous debris bed.