In the course of a postulated severe accident in an NPP, Direct Containment Heating (DCH) may occur after an eventual failure of the vessel. DCH is related to dynamical, thermal, and chemical phenomena involved by the eventual fine fragmentation and dispersal of the corium melt out of the vessel pit. It may threaten the integrity of the containment by pressurization of its atmosphere. Several simplified modellings have been proposed in the past but they require a very strong fitting which renders any extrapolation regarding geometry, material, and scales rather doubtful. With the development of multidimensional multiphase flow computer codes, it is now possible to investigate the phenomenon numerically with more details. We present an analysis of the potential of the MC3D code to support the analysis of this phenomenon, restricting our discussion to the dynamical processes. The analysis is applied to the case of French 1300 MWe PWR reactors for which we derive a correlation for the corium dispersal rate for application in a Probabilistic Safety Analysis (PSA) level 2 study.

1. Introduction

Direct Containment Heating (DCH) is a phenomenon that can potentially happen during a severe accident and may threaten the integrity of the containment by pressurization of its atmosphere. It can occur following a melting of the reactor fuel and its relocation in the bottom of the reactor vessel. Depending on the cooling conditions, the vessel might fail and the melt will be ejected out. If the primary circuit driving pressure is above some threshold, the melt will be, in a second stage, ejected out of the pit, going for a part into some small rooms (e.g., steam generator rooms, called for short hereafter subcompartments) and for the rest directly to the containment dome. If the melt can be oxidized (which is likely), then a very rapid oxidation of a large amount of the melt occurs. The hydrogen produced during this stage may subsequently burn altogether with the pre-existing one. Most of the past studies of this phenomenon were based on geometries of US NPP's (Surry, Zion, etc.). For a good synthesis of these works, the reader can consult a special issue of Nuclear Engineering and Design [1]. Regarding in particular the melt dispersion, several simplified models have been proposed in the past, mainly qualified on US reactor types, and implemented in lumped parameter codes as CONTAIN [2]. However, it is found that all models are highly dependent on the geometry with high variations of the main parameters [3]. Thus the extrapolations to other reactors, other materials, and other scales than those used for the fitting process of the parameters are doubtful. With the development of multidimensional multiphase flow computer codes, it is now possible to investigate the phenomenon numerically with more details, improving thus the comprehension, first, and prediction capabilities, second. At the IRSN, we are using for this purpose the 3-dimensional code MC3D, mainly developed for Fuel Coolant Interaction (FCI) analysis [4, 5].

The current analysis is realized in the frame of the development of a Probabilistic Safety Analysis (PSA) level 2 study on French 1300 MWe reactor. In order to minimize the uncertainties related to the use of unqualified models regarding the specific plant under consideration, it was decided to engage a program including specific experiments and the use of the MC3D code for the understanding and characterization of the dynamical aspects.

A priori, the DCH problem is mainly related to the evaluation of the following physical phenomena:

(i)melt ejection and dispersion to the containment;(ii)heat transfer from the melt to the containment atmosphere;(iii) melt oxidation, its energy, and consequences on material properties; (iv) hydrogen production coming from oxidation by steam; (v) combustion of the hydrogen (produced and pre-existing).

However, it is important to stress at this stage that this should be considered only as what can be called the “academic” DCH problem. In this viewpoint, the boundary conditions (i.e., the global geometry), are not affected by the phenomenon itself. The break will be supposed to be central and with a constant diameter, although some models exist for the break enlargement due to the melting. The break enlargement is not considered here due to the very high uncertainty regarding the break formation itself (other arguments that support this choice will appear further from the analysis). All in all, this has to be considered as an uncertain parameter. The vessel is supposed to have the shape of an intact one, that is, with no enlargement due to the internal pressurization (which induces the break). Also, the presence of the vessel insulator is not taken into account. The behaviour of this equipment is very uncertain. It should at least loose strongly its mechanical properties but should not disappear. It might then reduce the flow path area to the containment. At first, considering the vessel as undeformed and the vessel insulator as absent could be considered as a conservative point of view. This conclusion is however not totally clear if one considers that a flow path restriction tends to increase the mechanical effort onto the vessel itself and then the mechanical stress onto the attachment points of the vessel. This might involve an important displacement of the vessel that could in turn open some new paths. A second limiting feature of the academic DCH problem is that it considers that no water is present, neither in the pit nor in the vessel (laying over the melt or trapped somewhere in the Reactor Coolant System (RCS)). Both cases are in fact quite unlikely. The presence of water adds a very important additional difficulty. This is not taken into account in this study (neither investigated in the experiments). The impact of remaining water in the RCS was investigated in few CE/CES tests (Calvert Cliffs geometry) [6]. This was possible because in these experiments the melt charge was initially placed in the pit. These tests revealed an important increase of the duration of the vessel pressurization, due to water flashing. However, the overall pressurization was slightly smaller, due to the coejection of water with the melt. This is not sufficient to have firm conclusions since the results should strongly depend on the characteristics of the remaining water (mass, temperature). This effect could in principle be investigated, at least qualitatively, with tools like MC3D.

2. Experimental Database

The analysis is based on an experimental program investigating the phenomenon with different geometries representative of European reactors such as the EPR (DISCO series, performed by FZK, Germany [7]), the French 1300 MWe P4 (DISCO-F), VVER, and Konvoï German reactors. A general presentation of the program and the main outcomes can be found in [8]. We will focus our present study on the 1300 MWe P4 geometry although similar results are available for the EPR geometry.

For the P4 geometry the scale is 1 : 16th. Figure 1 gives a global sketch of the geometry as represented for a 3D calculation with MC3D. The geometry is complex, 3-dimensional, with an important path section towards the bottom access (path 3), and allows a distribution of the fuel ejected upwards to different spaces: containment dome (path 6) and subcompartments (path 7). A closer view of exit pathways towards the dome and subcompartments is given in Figure 2.

Some of the tests, denoted “2D,” were carried out with a simplified geometry with a cylindrical cavity without the access path. The purpose of this simplification is to allow an easier code qualification and analysis. The program comprised:

(i)10 water-tests called DISCO-F(X), including 2 tests in “2D” geometry;(ii) 4 cold tests with a heavy metal, made from a gallium alloy, called DISCO-FM(X), of which one is “2D”;(iii) 6 cold tests with a second heavy metal, the Wood's Metal, called DISCO-FW(X), of which 3 are “2D”;(iv) 5 hot integral tests with an alumina-iron thermite, called DISCO-FH(X) with(a) 1 test in neutral environment, avoiding chemical effects (FH02),(b) 1 test in “2D” geometry (FH05).

The pressure differences between the vessel and the pit were varied from 10 to 22 bars approximately. Three diameters of vessel breaks (i.e., vessel outlet nozzle) were used: 30, 45, and 60 mm, corresponding to 0.5, 0.75, and 1 m break diameters at full scale. Note that only tests with a central break were performed for P4 geometry but previous experiments with EPR geometry showed that this situation was bounding.

In the rest of the paper, we will firstly recall the major characteristics of the MC3D code and show the status of the analysis that has been performed thanks to its use. As MC3D does not treat combustion, we will focus the analysis on the dynamical aspects. Indeed, the dynamical features drives the thermal and chemical phenomena and it is then of first importance to describe this stage accurately.

3. MC3D: Analysis of the Potential andLimitations of the Code to Describe the Main Features of the Phenomenon

MC3D is a computer code devoted to study multiphase and multiconstituent flows in the field of nuclear safety. It is one of the reference tools for the evaluation of Fuel Coolant Interaction (FCI) but it is more generally a Computational Multifluid Dynamics (CMFD) tool that can treat various problems [5]. One of the most important particularities of MC3D is that the fuel is modelled with two fields which allow describing the two possible states of the fuel: either continuous or discontinuous (Figure 3). These two fields are connected by mass transfer constitutive laws (i.e., fragmentation or coalescence). By this way, the drop field is supposed to describe only discontinuous fuel state. The continuous fuel field is modelled with an original Volume of Fluid, Piece Wise Linear Interface Construction (VOF-PLIC) method, avoiding then numerical dispersion (Figure 4).

For the DCH evaluations, the fragmentation process is obtained through a model extending the Kelvin-Helmholtz (KH) instability model to multiphase fragmentation, taking into account local velocities. In this model, we write that the volumetric rate of fragmentation of the jet is given by

Γ𝑓=𝑁𝑓𝑐𝑖,(1) where

𝑐𝑖=1𝜌𝑓+𝜌𝑎𝜌𝑓𝜌𝑎𝑉𝑓𝑉𝑎2𝜌𝜎𝑘𝑓+𝜌𝑎||||𝑘=𝑘max(2) is the standard instability growth rate of KH model. Index 𝑓 stands for the continuous fuel, and 𝑎 stands for the ambient fluid. In a fully local model, the difficulty is to provide

(i) an adequate spatial averaging of the properties in order to limit the effect of mesh size,(ii) an averaging of the fluid properties and particularly a definition of the ambient fluid.

These numerical aspects cannot be detailed here. Let us simply say that the ambient fluid is related to a mixture formed by liquid water and the gases. Fuel droplets themselves are not taken into account. Then, for DCH simulation, the ambient fluid definition is straightforward since, in general, we do not consider the presence of liquid water.

In (1), 𝑁𝑓 is called the fragmentation parameter. It is related to the relative height (regarding the wavelength) of the instability when the fragmentation occurs. It also includes eventual nonlinear effects. The expected value stands between 1 and 6. In the calculations, 𝑁𝑓=2 was used, as for FCI evaluations. The wave number corresponding to the maximum growth is given by


In this model, the diameter of the created drops is related to the wavelength 𝜆 of the instability so that


𝑁𝑑 is a parameter with expected value between 0.1 and 0.5. In calculations, we still use the recommended parameter as for FCI calculations, that is, 𝑁𝑑=0.2.

Coalescence occurs from drop collision with the continuous fuel. Coalescence happens mainly at the impact along a wall.

Moreover, an oxidation model is available. Due to the specific conditions that MC3D has generally to face (high temperature melts, high pressure), the model is however quite simple and parametric. Melt oxidation in such conditions is a particularly complex problem. It is currently under investigation at IRSN in order to improve the code with more predictive constitutive laws in future versions. Currently, we simply set an amount of potential oxidation per surface area of the melt (called capacity), calculated with a transport equation. The oxidation itself occurs at a rate depending on the capacity, the time scale, and the vapour content. Up to now, the approach for this dealing with this problem for DCH evaluations has been mainly heuristic and empirical, based on the experimental evidences of strong and very rapid oxidation. We estimate that some questions regarding the actual characteristics of prototypical corium oxidation are still pending, for example, oxidation level of uranium, as well as the question of scaling. Improvements regarding this issue in codes like MC3D should then give important benefits regarding the understanding and evaluations.

The major key limitation in MC3D is then the evaluation of combustion which is currently not possible. Experiments have shown that this phenomenon drives the essential of the “heating” problem. This means that the treatment of combustion needs a particular attention. However, an important number of situations involve already high containment pressure at the vessel break. The high content of vapour in these cases might inert partially the combustion and then simple thermal heat transfers might play an important role. In the reactor case, oxidation of zirconium will also add an important amount of energy. For the issue of combustion, the treatment up to now has also been mainly empirical. To the knowledge of the authors, only few simulations with dedicated combustion codes were done and failed to describe with high accuracy the characteristics of combustion during a DCH event. Large uncertainties exist due to important energy losses in experiments, particularly if one considers the problem of extrapolating and scaling. However, taking into account the various strong uncertainties and the fact that we are seeking for a practical evaluation for nuclear safety, margins can be taken and conservative approaches might be sufficient. This topic will not then be addressed further in the following of the paper.

Note also that although this could be done quite easily, the code cannot compute friction and heat transfer with structures. In the present application, this should have only a minor impact, particularly at reactor scale.

One of the advantages of the MC3D code is its recognized ability to treat the interaction between the fuel and the water, either for slow mixing, either for steam explosion. It is then possible to combine with the same tool the effect of DCH and FCI. This will be exemplified in Section 8. The qualification of the tool for the specific DCH problem allows us to go far beyond what can be obtained with analytical models with a more comfortable confidence regarding extensions to, for example, different scales and material or different initial conditions as including the effect of the presence of water. Varying the mesh precision allows us to make several kinds of studies. With an intermediate mesh size, we will discuss the flow patterns calculated with the code. With a somewhat rough mesh, one can make a lot of calculations and thus having a statistical treatment to build correlations for application in simplified models for, for example, system codes or PSA studies. In contrast, very fine meshes can allow analysis of particularly complex features of the geometry as the upper part of the annular space where occurs the separation of the melt between the subcompartments and the containment. This last point is however still under investigation and will not be discussed further in this paper. We will also restrict the discussion to 2D geometries although some 3D situations have been studied with the code. However, the most important features are already present with 2D situations and it is sufficient to focus on this limiting case.

4. Flow Patterns Analyzed with MC3D

4.1. Pure Gas Flow

Investigating the pure gas flow has a limited practical interested. However, this is a necessary step particularly if one intends to build analytical simplified models as for system codes. Indeed, the gas is the carrier fluid and it is important to characterize its flow. It is interesting to notice the existence of a single intense convection roll occupying the whole pit (Figure 5). Then, the definition of a characteristic velocity in the pit (where occurs the entrainment), as required by simple lumped parameter models, is not an easy task. Along the wall, the velocity shows a very strong gradient. On the mean, we find that the upcoming velocity is in the same order as the one at the annular space. The velocity at the break is obviously sonic, but it might be noted that this is not the location of the maximum speed. We obtain a situation which is similar to a Venturi tube with supersonic flow after the contraction.

The vessel depressurization is computed on the basis of an adiabatic flow since no heat transfer is computed between the gas and the vessel. Experiments show that the situation is closer to an isothermal flow. However, the impact on the overall process is very limited. At reactor scale, the impact should be even lower.

4.2. Melt Ejection and Upward Dispersion

The qualification of MC3D regarding the melt and gas ejection out of the vessel have already been discussed in previous analyses [8, 10]. We propose now a visualization for the purpose of illustration and completeness of analysis in Figure 6. At first, the melt flows as an undisturbed jet. Despite the high ejection velocity, the jet fragmentation is limited: due to the small distance between the vessel and the floor in the experiment, the travelling time scale is smaller than the fragmentation time scale. This in fact is not true at full scale: the ejection velocity and then the fragmentation time scale are still the same but the travelling time is increased accordingly to the scale. However, we find that this modification of the flow pattern does not have an important influence on the overall process of dispersion out of the pit. Regarding the impact of the jet with the floor, the code does not (cannot) represent an eventual splashing effect as it was observed in some experiments. However, this splashing does not seem to have a significant effect. More, the effect should also diminish at full scale. Rapidly, a hole is formed inside the melt pool (third picture) and the gas begins to flow out of the vessel with a high velocity. As the gas flows inside the melt, forming first some kind of expanding bubble in the case of large breaks, this one is totally fragmented and dispersed in the pit. However, the drops rapidly hit the floor or the wall and undergo a recoalescence. Before this two-phase flow stage starts, the melt can be accumulated along the wall simply by the effect of inertia. Depending on the vessel pressure, it can either flow back and return on the floor, either be entrained and flow along the wall.

The gas velocity becomes really important only once most of the melt is ejected. During the stage of two-phase flow, the gas velocity is strongly reduced by friction with the dispersed melt. We find that in MC3D calculations, the melt along the wall is really entrained and fragmented only once this essentially single-phase gas flow starts. In MC3D, it is found that the fragmentation mostly occurs as the continuous melt approaches the annular space. It is also seen from Figure 6 that gas flow patterns seen with pure gas flow in Figure 5 are strongly disturbed. We often see the development of two circular rolls (best seen on Figures 8 and 9).

The qualification of the code regarding the experimental results can be shortly discussed with the help of Figure 7 which illustrates the results of a calculation of the test F07 with water as simulant in a 2D geometry (burst pressure =1.6 MPa, break diameter =60 mm). There is a mist which is formed in the experimentation due to fragmentation at the edge of the break. This has no identified effect on the dispersion, even for cases with hot simulant. We observe again the formation of a hole inside the jet. This leads to the formation of a kind of pressurized bubble inside the jet. In the present simulation, the bubble breaks just before the jet reaches the floor. This explains why no single-phase liquid flow pattern is detected from the pressure traces in F07. The pictures are taken during the test show that the upward dispersion starts quite lately after about 50–60 milliseconds. For both the calculation and the experiment, we see that the upward ejection begins at approximately the same time as the transition to single-phase flow. This behaviour is confirmed by all MC3D calculations with both the P4 and EPR geometries. It also confirms the representativity of experiments as CE/CES tests [6] where the melt is initially placed on the pit floor.

4.3. Gas Flow Characteristics during Melt Ejection and Dispersion Steps

Some further characteristics of the gas flow can be qualitatively discussed with the support of Figures 8 and 9. For a characteristic calculation with hot reacting material, both the gas temperature (Figure 8) and the vapour partial pressure (Figure 9) are strongly affected by the fragmentation stage at the break. In this calculation, we set the time scale for melt oxidation at the relatively small value of 10 ms. The real characteristics of thermite oxidation kinetic has not been yet investigated but it is likely that the phenomenon is very rapid due to the high temperature. We however are only seeking here for a qualitative description. As soon as the two-phase flow appears at the break, the gas is strongly heated.

At the experimental scale, the pit is entirely heated at a temperature close to the melt one in a few milliseconds. Meanwhile (but independently in fact) the vapour pressure has mostly disappeared from the pit (Figure 9). The minimum residual value of the vapour partial pressure is due to the limitation of oxidation set in the calculation close to 1000 pa. When the flow becomes essentially a single-phase gas flow, then the pit gas starts cooling down and fresh vapour can re-enter in the pit. As the dispersion process occurs mostly during this single-phase gas flow, it is seen that dispersion occurs with very transient conditions regarding the gas properties.

5. Dispersion Correlation

5.1. Method and Limitations

Finding correlations that would work for any case and any geometry has been the subject of work of many researchers in the past (see the important number of correlations in the CONTAIN code). Our purpose is to concentrate on a particular geometry and correlate parameters regarding the initial physical conditions (gas, melt) and the break section. The limited number of experiments makes this task difficult. One of the obvious interests of Computational Fluid Dynamics (CFD) tools (when qualified) is to allow a substitution of the experiments with the calculation. The experiments at our disposal are mainly used to qualify the codes in some particular simulating situations and to verify the correlation. Then we propose a mixed method where we will seek for a correlation with parameters that need to be adjusted thanks to CMFD calculations and validated with experimental results.

At first, we have to admit that, unless using fine 3D meshes, it is necessary to treat the dispersion globally, that is, without taking into account the particular paths to the containment and subcompartments. The particular problem of separation of melt can be investigated only with a different method. Also, the bottom exit path cannot be treated easily unless 3D calculations are performed. Although this is feasible, previous estimations tend to show that the 3D complete geometry behaves similarly to the 2D case [10]: the inceptions of entrainment were found very close and only the maximum dispersed fraction were changed. Then, it is likely that we can use the 2D geometry to build a general correlation that will be calibrated for the 3D case with experimental results and some additional particular 3D calculations. It is also rather clear that we cannot have a very fine precision regarding the maximum possible dispersal. In the code, this depends mostly on the meshing and numerical methods. In the experiment and real reactor case, it depends on local complex geometries particularities that cannot be reproduced.

5.2. Theoretical Analysis

The first idea that comes in mind to correlate the upward dispersion is to relate it with the ability of the gas to entrain melt droplets. The relative terminal velocity of a monodispersed droplet population in an ambient fluid can approximately be given by

𝑉𝑉01𝛼𝑑=1𝛼𝑑83𝑅𝑑𝐶𝑑0||𝜌𝑑𝜌𝑎||𝜌𝑐𝑔,(5) where the subscript 0 stands for one single droplet.

On the other hand, the relative flow tends to fragment the drops up to equilibrium. If one single drop is stable for a relative Weber number of, say, 12, it is found experimentally [11] that the mean Weber number of the created drops (relatively to the initial velocity) is lower. More precisely, the drops are not created such that the Weber is equal to the critical one but can be significantly smaller. We will call this number the equilibrium Weber number and take it as an uncertain parameter that should range from 1 to 6. The equilibrium Weber number gives us a second relationship between the velocity and the drop mean diameter:


In MC3D, the equilibrium Weber number gives its actual limit for secondary fragmentation. Considering only dispersed flow with low drop volume fractions and drop densities far higher than the ambient one, one has an expression of the relative velocity as


This gives a lower limit for the gas velocity that can lead to melt entrainment. This in fact defines a dimensionless number called the Kutadeladze number:

𝜌𝐾=𝑎𝑣2𝜎𝜌𝑑𝑔,(8) and a criteria for entrainement can be expressed as


If the drop diameter is in the inertial range, 𝐶𝑑0 is about 0.5 and we obtain a critical Kutadeladze of about 3. The Kutadeladze number gives a criterion for inception of entrainment but our purpose is to find a correlation giving the upward dispersion fraction. It is however likely that we can use this number for this purpose. The dispersion is a function of the mass flow rate and thus we will seek for a relation using the square root of the Kutadeladze number, that is,


The problem now is to have a characteristic gas velocity in our system. The fragmentation and entrainment occurs in the reactor pit where no obvious characteristic velocity can be defined. However, it is likely that the velocity along the wall can be related to the velocity in the annular path. Then, for the purpose of a correlation, one can use the annular path velocity as a characteristic velocity to compute the Kutadeladze number. For the case of pure gas flow, all evaluations show that the pressure inside the pit comes very rapidly to equilibrium and the flow is rapidly in a quasisteady state. In that case, the mass flow at the entrance and outlet of the pit are balanced. At the vessel break, the flow is, in all cases of interest, sonic. The classical theory of compressible gas flow gives us the mass flow rate as


Then, the velocity in the annular space is


The correlation will then take the form


In this relation, we identify two difficulties. The first difficulty comes from the influence of the break diameter. A small break will increase the duration of the gas flow, having then an influence on dispersion opposite to the one in (13) (i.e., inversely proportional). Also, for a small break, the two-phase flow stage is longer and there might be an important pressure loss in the vessel during this stage. Then the influence is not so simple. On the overall, the influence should not be linear but rather proportional to 𝐴𝑛𝑏 with a small exponent 𝑛. The second difficulty comes from the influence of gas density in the pit. This one is not known a priori and we have seen previously that it should change drastically during the transient.

5.3. Evaluation with MC3D and Comparison to Experimental Data

In order to be able to have a statistical treatment, we use the rough mesh shown in Figure 10, the pit being represented by a 8×15 grid. Such rough mesh has proven to be sufficient to predict the upward dispersion for all tested geometries [10]. Unless indicated, all 2D calculations presented hereafter are done with an equilibrium Weber number of 1. The parameters that have been checked to build the correlation are listed in Table 1 with the range of variation.

Figure 11 shows the influence of break area (scaled with the annular path area). We find that the best fit is obtained with an exponent equal to 1/4 instead of 1, probably due to the longer ejection time for small breaks, as expected. The scatter of the data is, as already explained, mostly due to the very rough mesh used. This is a characteristic due to the continuous fuel field numerical treatment (VOF-PLIC), which behaves only approximately with rough meshes. However, in most applications, it is not necessary to have a precise result for a particular situation, but rather a global tendency. On the overall, it is then found that the break diameter is not a very critical parameter, thus limiting the impact of the important uncertainty regarding the break characteristic during a severe accident. We also see on this figure the good correlation with the function 𝑃𝑣/𝑇𝑣. Then, the temperature inside the vessel is also an important parameter.

The impact of the other parameters not explicitly appearing in (13) (but eventually through the gas density) is found quite mild. We do not see a significant impact of the initial containment pressure (if limited to 3 bars). In contrast, the melt temperature seems to increase noticeably the dispersion. However, in the hot calculations, the combustion effects are not included and it is likely that the increase in containment pressure would lead to a decrease in dispersion. Due to the various uncertainties, it did not seem relevant to include additional parameters in (13) and the gas density was simply changed to a reference density 𝜌0 (1 kg/m3).

The overall result of all the 2-D calculations is plotted in Figure 12 with the available experimental data and the best fit correlation of the form


We find, for the “2D” DISCO 𝑃4 experimental facility: 𝐾50=170, 𝐴=5, 𝐹𝑑,max=0.9.

It will be emphasized that the experiments confirms that integral tests behave approximately in the same way as cold tests.

The correlation can now be compared to the experimental 3D data, just changing the maximum dispersion in (14). This is done in Figure 13. In such case, we only modify the maximum dispersion which from experimental results but also MC3D calculations (see [10]), is about 0.6. There is a somewhat important scatter of the data around the correlation. However, considering the very important variations of the physical conditions between a nonreacting water flow and an integral test with chemical reactions, it can be estimated that the scatter is largely in the uncertainty range of actual conditions during a severe accident and probably of some of the experimental conditions. For the use in a particular problem, it is envisaged to fit more particularly on both the calculations and experimental results (if any) with the particular conditions of consideration. The fit might then be calibrated using corium physical properties, vapour as a driving gas, and high melting temperature.

Nonetheless, with the help of MC3D, it was possible to find a best fit by modifying slightly the influence of melt density (exponent 1/2 instead of 1/4) and suppressing the surface tension. Since no rational can be found for these modifications, and due to the high uncertainty attached to this evaluation, we estimate that correlation (14) is sufficeint for pratical evaluations.

Using the derived correlation, we can now already fix ideas of what can be expected for a representative DCH in a severe accident with however the scale of DISCO tests. This is presented in Figure 14 (see legend for conditions). The inception of entrainment is expected to occur approximately at 20 bars. A full dispersion is obtained at approximately 40 bars with a large break (1 m at full scale) and 60 bars with an intermediate one.

6. Scaling

The scaling is of course a key issue. A large uncertainty can be attached to this problem if one relies only on experimental results and simplified models. The use of CMFD codes as MC3D allows investigating this issue with an increased confidence, although scaling can have an important effect on some specific constitutive laws, particularly when it is obtained mostly with a simple mesh size modification. In the particular case of DCH, the primary fragmentation (i.e., continuous fuel to drops) seems the most critical model regarding mesh sensitivity. However, our evaluations showed a strong sensitivity to the secondary fragmentation model which led us to consider an equilibrium Weber number of the order of the unity. Then, the primary fragmentation model is expected to have, on the overall, a limited impact. Thus, the impact on scaling on the models is expected also to be limited.

Regarding dispersion, the calculations did not show a noticeable effect regarding the global results of melt dispersion. Some of the flow patterns are different particularly regarding melt ejection. As an example, at full scale, the nearly intact jet issuing the vessel is not visible at full scale. In most cases, the aerodynamic processes of fragmentation have sufficient time to break the jet before hitting the floor. For large break, the two-phase flow also occurs before the melt-bottom contact. Then on the overall, the initial melt fragmentation is more intense. However, due to the recoalescence along wall and floor, there is no really clear effect in the calculations. Figure 15 gives an example of comparison of scaling effect. The dispersion is plotted as a function of the modified Kutadeladze number (see previous section). What has been found from the various trials of calculation at full scale is that the inception of dispersion is not modified. For conditions with strong dispersal, the effect is not clear and, on the mean, no important effect has been detected.

It should be recalled that MC3D does not compute effects related to the presence of the walls such as friction or heat transfer. At large scale, the impact of wall should be smaller due to a smaller ratio surface/volume. We did not notice important effects at DISCO scale so it is likely that these effects are also minor at reactor scale. On the overall, we conclude that the scaling has no noticeable effect regarding the dispersal processes of the melt.

7. Further Potential Capabilities of MC3D

The “academic” problem of DCH does not take into account the presence of water. However, the probability for presence of water either in the vessel or in the pit is rather high for various reasons. MC3D seems a particularly well-suited code for taking into account the presence of water at least regarding the problem of dispersion. The case of water in the pit is quite easy, as MC3D is already used for the evaluation of Fuel Coolant Interaction in the so-called exvessel situation. This is nothing else, for the code, than a DCH with water in the pit. The case of water in the vessel either just above the melt, either somewhere trapped in the lines of the Reactor Cooling System is more interesting since the flashing might induce a longer pressurization in the pit and then an enhanced dispersion. The case of water lying above the melt has been investigated recently for demonstration purpose and is illustrated in Figure 16.

It is very likely that the water in the vessel will be close to saturation. Then, at the exit from the vessel, a flashing phenomenon will occur, where some part of the water will evaporate quasi-instantaneously. In the example given, the overpressure is only 5 bars, so far below the inception of entrainment without water. Another test with a smaller mass of water and 8 bars as overpressure gives a full dispersion. We find that the flashing increases strongly the pressure in the pit and then the dispersion. In case of large masses of water, there is a large remaining part of water which is expelled together with the melt. This expelled water might inert the combustion and reduce the pressurization. In other situations, the water might flash completely thus avoiding these limiting effects. The impact of water is still under investigation at IRSN.

8. Conclusion and Perspectives for Further Developments

We highlighted the interest of the use of qualified multiphase flow codes as MC3D for the evaluation of the DCH phenomenon during the course of a hypothetical severe accident. Due to the complex geometries of NPPs, this problem has for a long time been faced to the technical difficulties in representing adequately the phenomenon with simplified analytical models. Correlations were largely inspired by experimental results with only restrained applications and weak confidence in extrapolations to real material and real scale.

With the use of the computer code MC3D with relatively rough meshes, we could present a method for developing a dispersion correlation. For each geometry (EPR, large pit reactors), the calculations can be used to calibrate precisely the correlation, as it has been presently done for the case of French P4 reactors.

For a full evaluation of the phenomenon, one needs a model of combustion. However, even codes dedicated to combustion can have some difficulties to evaluate precisely the phenomenon. Then, despite the progress made and presented here, we are still faced to a complex problem. In the case of MC3D, it is then envisaged to seek only for a simple model, dedicated to the problem.

Despite this difficulty we have elaborated some patterns of explanation for the dispersion of the melt, we find that 2D and 3D situations might not be so different. We are also able to specify with a good confidence the range of physical conditions for which the dispersion can occur and then propose a range of conditions for which the phenomenon cannot threaten the containment. We can also integrate some supplementary conditions that are usually not taken into account as the presence of water.


𝐴:Section, area (m2), coefficient in correlation(14)
𝑐𝑖:Characteristic velocity of instability (ms1, (2))
𝐶𝑑:Friction coefficient
𝐷:Diameter (m)
𝐹:Dispersed fraction
𝑔:Gravitational constant (ms2)
𝑘:Wave number (m1)
𝐾:Kutateladze number
𝑀:Molar mass (kgmol1)
𝑁𝑓:Fragmentation parameter(1)
𝑃:Pressure (pa)
𝑄:Mass flow rate (kgs1m2)
𝑅:Radius (m), gas constant
𝑣,𝑉:Velocity (ms1)
We:Weber number
Γ𝑓:Volumetric fragmentation rate (ms1, (1))
𝛾:Adiabatic index
𝜌:Density (kgm3)
𝜎:Surface tension (Nm1 or Jm2)
𝑎:Ambient or annular space
𝑑:Drop or dispersed


The authors gratefully acknowledge Dr. Leo Meyer from the Forschungszentrum Karlsruhe (FzK) for successfully conducting the DISCO program.