- Email: [email protected]
The simulation of melt spreading with THEMA code
Nuclear Engineering and Design 236 (2006) 415–424
The simulation of melt spreading with THEMA code Part 1: Model, assessment strategy and assessment against analytical and numerical solutions B. Spindler ∗ , J.M. Veteau Commissariat a` l’Energie Atomique, CEA-Grenoble, 17 rue des Martyrs, 38054 Grenoble, France Received 28 September 2004; accepted 16 September 2005
Abstract In the framework of R&D on severe accidents for PWRs , definition of mitigation means for cooling of the molten core after a vessel melt-through is of utmost importance. In the case of spreading-based concepts, the THEMA code, which considers integrated balance equations over the fluid depth, aims at predicting the spreading extent of non-eutectic mixtures under defined inlet conditions (pouring rate, temperature and composition). In this first part, the models and the assessment strategy of the code are presented. The satisfactory running of the main moduli of the code is then verified by comparing the numerical results against known solutions of basic problems: dam break problem, conduction with or without melting in the basemat, crust formation kinetics on a constant or “free” bottom plate temperature. © 2005 Elsevier B.V. All rights reserved.
1. Introduction The safety authorities of several countries have requested the industries and utilities to take into account severe accidents in the design of the next generation of nuclear power plants. Severe accidents refer to situations where high temperatures are reached in the core, leading to the formation of high temperature mixture of metals and oxides (corium). Among the mitigation means, which can be envisaged, special devices (core catchers) are considered, in which the corium is collected and cooled in order to insure the integrity of the vessel containment. A peculiar concept of core catcher, retained for instance for the European Pressurized Water Reactor (Weissh¨aupl, 1999), is based on the cooling of corium spread on a large area. Within this context, large research efforts are conducted throughout the world in order to get better understandings of the phenomena occurring during spreading of high temperature oxide and metal mixtures. CEA has undertaken a large program on severe accidents (Cognet et al., 1997), which aims at providing the tools dedicated to understand and simulate these phenomena.
∗
DOI of original article:10.1016/j.nucengdes.2005.09.016. Corresponding author. Tel.: +33 438 78 46 87; fax: +33 438 78 52 51. E-mail address: [email protected] (B. Spindler).
0029-5493/$ – see front matter © 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.nucengdes.2005.09.009
The development of THEMA code is part of this work. The THEMA code is devoted to the simulation of spreading, taking into account the melt solidification, as well as the possible melting of the basemat. The main characteristic of THEMA code is the use of balance equations integrated over the melt thickness. Integrated equations for the melt were chosen because it was stated that a very accurate formulation is of weak interest, since large uncertainties remain concerning the physical properties of real material on one hand, and for purpose of numerical simplicity and robustness, on the other hand. This model is also used in the MELTSPREAD code (Farmer et al., 1993), but limited to 1D geometry, whereas THEMA is able to simulate 1D or 2D geometry. The LAVA code also uses integrated equations, but inertia terms are not taken into account in the momentum equations (Allelein et al., 1999). Beside theses codes, some other codes consider a meshing within the melt thickness: CORFLOW developed by Framatome-ANP (Wittmaack, 1999), and CROCO developed by IRSN (Piar et al., 1999). The run time of such codes is very large. Nevertheless, computing time of THEMA is also increased for a 2D geometry, due to the 3D meshing used for computing heat conduction and melting in the basemat. In this first part, the models and the assessment strategy of THEMA code are first presented, followed by the results obtained in the frame work of comparison with known analytical
416
B. Spindler, J.M. Veteau / Nuclear Engineering and Design 236 (2006) 415–424
Nomenclature A coefficient C coefficient C0 coefficient Cp heat capacity (J kg−1 K−1 ) du/dy velocity gradient (s−1 ) e thickness (m) f friction factor or volumetric solid fraction F imposed heat flux (W m−2 ) H enthalpy (J kg−1 ) f0 volumetric solid fraction g gravity (m s−2 ) h heat transfer coefficient (W m−2 K) H enthalpy (J kg−1 ) J superficial velocity (m s−1 ) K heat conductivity L latent heat (J kg−1 ) S surface (m2 ) t time (s) T, T0 , T1 temperature (K) liquidus temperature (K) Tliq V velocity (m s−1 ) Vol volume (m3 ) drift flux velocity (m s−1 ) Vgj x abscissa (m) Greek letters α void fraction or diffusivity γ melting constant δ, δ0 crust thickness (m) Γ mass transfer term (kg s−1 ) µ viscosity (Pa s) µ0 viscosity for f = 0 σ surface tension (Pa m) τ0, τ shear stress (Pa)
or numerical solutions of basic problems, then ensuring a correct running of the main moduli of the code. The second part is devoted to the assessment results based on various spreading experiments: CORINE, KATS, ECOKATS, COMAS, FARO and VULCANO. Details of the assessment synthesis of THEMA can be found in a CEA report by Spindler and Veteau (2004). 2. THEMA models 2.1. Spatial distribution of the phases and temperature profile The spatial distribution of the solid and liquid phases considered in THEMA is presented in Fig. 1, together with the corresponding temperature profiles. The use of integrated equations corresponds to a single temperature for the melt. Transient temperature profiles are calculated in the bottom and the top crust and in the basemat. 2.2. Equations 2.2.1. Flow of a fluid with a floating crust The integrated equations were obtained from the well-posed complete system of equations related to the flow of a fluid composed with several phases (metal and oxide), with a solid crust above it. The basic integrated equations are those of the classical two-phase flows (Ishii, 1975) and the detailed process of the simplification of the system of equations is presented in the thesis of Eberl´e (1997). The objective of this work was to describe, together with the liquid flow, the flow of the solid crust at the upper surface, entrained by the liquid. However, for practical applications concerning corium, the mechanical properties of the crust at high temperature, which are necessary to close the model, are not known. Hence, it was decided to simplify the system equations, discarding the mechanical effects of the upper crust. Another reason for this simplification is that the effects of the upper crust in the spreading experiment are not clear. One
Non-dimensional numbers Fr Froude number Nu Nusselt number Pr Prandtl numer Re Reynolds number Subscripts b basemat bcr bottom crust c collapsed fus fusion g gas i interface L liquid S solid u upper
Fig. 1. Spatial distribution of the solid and the liquid phases and corresponding temperature profiles. Ts is the solidus temperature of the melt.
B. Spindler, J.M. Veteau / Nuclear Engineering and Design 236 (2006) 415–424
conclusion of the European CSC project (Cognet et al., 2001) was that, for the melts with a large solidification interval, there are no crust effects, while for the melts with a low solidification interval, the melt seems to be stopped by the formation of the crust. The investigations made in the KATS program (Engel et al., 1999) conclude that the upper crust effects occur only at the end of the spreading transient. 2.2.2. Main hypotheses The thickness of the melt is considered as low compared to the extent of the melt in the main direction of the flow. This hypothesis, known as the shallow water hypothesis, means that the pressure term in the momentum equation reduces to the gravity term, and hence the pressure term no more appears in the equation as a main variable. An other hypothesis which is just mentioned here, because practical applications were only performed with a single oxide or metal melt, is that the same velocity is used for the different liquid phases. 2.2.3. Set of equations A 1D or 2D, rectangular or cylindrical geometry can be used. After simplifications corresponding to the hypotheses stated above, the final set of equations is presented here, in case of a 1D rectangular geometry. Mass equation: ∂(e · ρ) ∂(e.ρ.V ) + = ∂t ∂X
Γ S
where e, ρ and V are, respectively, the thickness, the density and the velocity of the melt. S is the exchange area between the melt and the basemat and Γ is the mass transfer terms (from the melt to the crusts and from the ablated basemat to the melt). Momentum equation: ρ·
∂V ∂V ∂[ρ · (e + eb + ebcr )] +ρ·V · +g· ∂t ∂X ∂X
1 ρ · f · V2 + · =0 2 e where es and ebcr are the basemat and the bottom crust thickness, while f denotes the friction factor. Energy equation: ∂(e · ρ · H) ∂(e · ρ · V · H) + ∂t ∂x = hbcr · (T − Ti ) + hu · (T − Tu ) +
(Γ · H) S
where H is enthalpy of the melt, hb and hu the heat transfer coefficients, respectively, from the melt to the bottom crust and from the melt to the upper crust (or to the surface itself if no upper crust exist). Ti is the interfacial temperature between melt and crust, and Tu is either equal to Ti if there is an upper crust, or to the surface temperature of the melt. The last term corresponds to the enthalpy associated with the mass transfer
417
from the melt to the crusts, and from the ablated basemat to the melt. In case of gas issuing from the concrete basemat decomposition, a void fraction α is calculated. It is taken into account in the momentum equation, where the density becomes the two-phase density [=α·ρg + (1 − α)·ρL ], and e becomes the height of the swollen level. 2.2.4. Discretization and resolution The numerical integration of the finite difference set of equations uses a semi-implicit method with a Newton–Raphson iteration scheme with the melt height as main variable. A secondorder discretization in space is used, in order to get a better accuracy of the melt front progression. The time step is first limited by the Courant condition but is also controlled by the time variations of the wall heat flux and the crust thickness. Hence, the sharp variations of theses variables encountered during the thermal shock, occurring when the melt comes into contact with the cold basemat at the front, can be predicted with a controlled accuracy. 2.3. Limiting height at the leading edge The surface tension term is taken into account through a limiting height at the spreading front, obtained through the hypotheses of equilibrium between gravity and surface tension forces, and of a radius of curvature at the front equal to the height of the melt. Hence, the melt progression at the front is allowed only if the height of the melt is higher than the minimum height emin : 2·σ emin = ρ·g 2.4. Solidification and crust models Part of the melt is encrusted on the basemat on one hand, and at the upper surface on the other hand. The interface temperature between melt and crusts is the solidus temperature of the melt, since no thermodynamic equilibrium may occur within the short time of spreading. The crust has the same composition as the melt. Moreover, for non-eutectic mixtures, a specific viscosity model applies when the melt temperature decreases below the liquidus temperature, in order to depict the resulting semi-solid flow. 2.4.1. Crust models The transient conduction equation in the crust together with the energy balance at the interfaces (melt–basemat interface for the bottom crust; melt–surroundings for the upper crust) is solved using, for the temperature profile in the crust, an analytical iteration procedure from the steady state solution, as already proposed by Savino and Siegel (1968). For the bottom crust, the crust thickness and the heat flux from the crust to the basemat is obtained. For the surface crust, a set of two equations is
418
B. Spindler, J.M. Veteau / Nuclear Engineering and Design 236 (2006) 415–424
solved with two unknowns: the crust thickness and the surface temperature at the crust surface. As already pointed out, the present validation of THEMA is made without mechanical effects of the crusts. It appears that it is sufficient with regard to the accuracy that is required in the frame of the severe accident studies. 2.4.2. Viscosity models in the solidification interval A complete methodology for the calculation of the viscosity in the solidification interval is given by Ramacciotti (1999) and by Ramacciotti et al. (2001). The method of Ramacciotti gives a viscosity depending on the volumetric solid fraction f and a coefficient C laying between 4 and 8: µ = µ0 · exp(2.5 · C · f ) where µ0 is the viscosity of the liquid phase (f = 0). In the full set of THEMA calculations, the value used for µ0 is the viscosity of the liquid phase given by the experimentalists. It may also be predicted, for instance, with the Andrade model (Ramacciotti et al., 2001). Values of 6.1 and 6.3 for coefficient C were obtained, respectively, from the VULCANO tests VEU1 and VEU7 (Journeau et al., 2003). This formula is sometimes fitted (coefficient A in the formula hereunder) by the Shaw formula (Shaw, 1969) involving the melt temperature: µ = µ0 · exp[−A · (T − Tliq )] where Tliq is the liquidus temperature of the melt. The formulation of Ramacciotti using the solid fraction seems to be better, since for some mixtures the solid fraction variation versus temperature is not linear, and the equivalent viscosity is rather directly connected to the solid fraction than to the temperature. The first validations of THEMA were made with the model of Chong et al. (1971) better known as the Stedman model (Stedman et al., 1990): 0.75f µ = µ0 · 1 + f0 − f
with f0 = 0.5. That means that the value f = 0.5 corresponds to an infinite viscosity, i.e. to an immobilization of the melt. The Pinkerton–Stevenson viscosity model (Pinkerton and Stevenson, 1992) is a mix of the Stedman model (with f0 = 0.5), and the Shaw model: [−0.04 · (T − Tliq )] µ = µ0 · exp (1 − f/f0 )2.5 The comparison of the viscosity factor µ/µ0 given by the Stedman correlation, the Ramacciotti correlation (with C = 6) and the Pinkerton–Stevenson correlation (with f0 = 0.5) with a typical solidification interval of 500 K and with a linear relationship between temperature and solid fraction for simplification) is shown in Fig. 2. The Stedman correlation undergoes a slow increase of the viscosity for a low solid fraction, and a large increase when the value 0.5 is approached, while the Ramacciotti correlation provides also a large increase for small values of f. The Pinkerton–Stevenson is a combination of the two others. Concerning the solid fraction at which immobilization occurs, it is imposed with the parameter f0 for the Stedman and Pinkerton–Stevenson correlations. For the Ramacciotti correlation, immobilization occurs at a solid fraction not precisely determined, depending on the pouring rate. The first assessment of THEMA was made in the past using the Stedman correlation (Spindler et al., 1999). In the synthesis presented here, the calculations are made with the Stedman correlation with f0 = 0.5, and with the Ramacciotti correlation with C = 6 (corresponding to the mean value between 4 and 8, and to the value obtained by Journeau et al. (2003). The Stedman or Ramacciotti models are used for the mixtures of oxides only, while a specific model was developed for the non-eutectic mixtures of metals used in the CORINE experiments. For the spreading of pure metals (KATS iron tests) no semi-solid viscosity model is used, since there is no solidification range. 2.4.3. Simulation of Bingham fluids The validation of THEMA is made with the viscosity model described above, corresponding to a Newtonian fluid. Nevertheless, it is also possible to make simulations with THEMA using a Bingham fluid model, as in the CORFLOW and LAVA codes.
Fig. 2. Viscosity factor µ/µ0 vs. solid fraction, calculated with the Stedman correlation (with f0 = 0.5), the Ramacciotti correlation (with C = 6) and the Pinkerton–Stevenson correlation [with f0 = 0.5 and f = (Tliq − T)/500].
B. Spindler, J.M. Veteau / Nuclear Engineering and Design 236 (2006) 415–424
For Bingham fluids, the relationship between the shear stress τ and the velocity gradient du/dy is the following: τ = τ0 + du = 0, dy
µ du , dy
if τ > τ0 if τ ≤ τ0
where µ is the viscosity, as defined and used for a Newtonian fluid and τ 0 is the shear stress threshold. Practically, for a Bingham fluid, the wall shear stress coefficient is modified taken into account the supplementary parameter τ 0 . The complete formulation of the modified wall shear stress coefficient, in terms of the depth integrated fluid velocity, was introduced in THEMA. The comparison of simulations of the same spreading experiment, performed with the Newtonian model and with the Bingham model leads to the conclusion that the two results are very close. Since the parameter τ 0 is not well known, this kind of model cannot give better results than the Newtonian model. Furthermore, the selection of the viscosity model is often a way to match satisfactorily calculations with experiments. For this purpose, the Bingham model provides two fitting parameters, which is not desirable, even though corium is probably not a pure Newtonian fluid. Finally, the use of the Bingham model is not recommended. 2.5. Wall friction coefficient A classical correlation is used for the wall friction coefficient, based on a Reynolds number Re calculated with a hydraulic diameter equal to four times the melt thickness (geometry of a free surface flow, with a small thickness compared to the width). For low Reynolds numbers, the laminar friction factor is 24/Re, whereas the coefficient in the turbulent regime is 0.083 Re−0.25 . A more complex model, based on the boundary layer theory, was also implemented in THEMA (Eberl´e, 1997). But generally, the boundary layer thickness rapidly reaches the melt thickness, and the results are then quite similar (Spindler and Veteau, 1998). As already pointed out, the gas effect is taken into account through an increase of the melt thickness. Furthermore, a modification of the wall shear stress coefficient is approached by the concept of two-phase friction multipliers, already worked out for two-phase flows in pipes. Details of the model used for spreading with gas percolation are given by Veteau et al. (2003). The two-phase friction multipliers depend on the ratio of the fluid to the gas viscosity. The following correlations are used: • For
µ µgas
≤ 103 : Friedel correlation (Hewitt, 1982).
• For 103 < µµgas ≤ 105 : Lockhart–Martinelli correlation (Hewitt, 1982). • For µµgas > 105 : multiplier is (1 − α)2 deduced from the CORINE experiments. Details of the adaptation for free surface flows are given by Veteau et al. (2003).
419
2.6. Heat transfer coefficients Classical correlations are also used for the heat transfer coefficients. For laminar flows the Nusselt number is equal to 7.6, and for turbulent flow the well-known Colburn correlation Nu = 0.023 Re0.8 Pr0.4 is used for oxide melts. For metallic melts, specific correlations are used (Duchatelle and Vautrey, 1964): Nu = 6.0 for laminar flows, and Nu = 5.14 + 0.0127 Re0.8 Pr0.8 for turbulent flows. No specific effect of the presence of gas is considered in the heat transfer coefficients, despite it is known that gas bubbling increases the heat transfer. 2.7. Void fraction model A void fraction model is required in the two-phase friction multiplier approach. This quantity is deduced from special experiments conducted with the simulating materials used in the isothermal CORINE spreading tests with gas fed from the bottom. The data were analysed in the frame of the drift-flux void fraction model: α=
1 C0 +
Vgj Jg
where α is the volumetric void fraction and Jg is the superficial gas velocity. The parameters C0 and Vgj have been correlated against the collapsed melt thickness ec (Veteau et al., 2003). For water as working fluid, the following correlation was found: C0 = 26.8ec + 0.8224
with ec in m, and Vgj = 0.043 m/s.
These values are used not only for the simulation of the CORINE tests with water, but also for high temperature tests with freezing (oxide and iron KATS tests on a concrete basemat, VULCANO tests) and are recommended for reactor applications. Similar but different correlations have been developed to simulate the CORINE isothermal tests conducted with viscous fluids. 2.8. Heat transfer in the basemat A specific 3D heat transfer module is used, solving the heat transfer equation in the solid basemat. A fine meshing of the basemat is used near the melt–basemat or crust–basemat interface, a region where the temperature variations in the basemat are large. The heat of fusion is taken into account in the specific energy of the solid material, which enables the calculation of the basemat ablation. The gas release in a concrete basemat is taken into account: bounded water and carbon dioxide release when the temperature of the meshes reaches 373 and 746 K, respectively. The corresponding latent heat and decomposition heat are considered with a modification of the concrete specific heat at these temperatures. A thermal resistance may also be included between the basemat and the melt or crust, but in absence of a sufficiently accurate value of this parameter, it is recommended to set this value
420
B. Spindler, J.M. Veteau / Nuclear Engineering and Design 236 (2006) 415–424
at zero. Furthermore, the sensitivity studies that were made show that the effect of the thermal resistance on the spreading length is low for the spreading of oxide mixtures, but large for iron or eutectic metallic melts. Some sensitivity studies are presented in the second part, which show the effect of the thermal resistance, and how it can be (improperly) used as a fitting parameter. The determination of the thermal resistance from matching the computed and experimental temperature profiles in the basemat from standard spreading experiments is quite hazardous, especially with low conductivity materials like ceramic or concrete, owing to the very large temperature gradient in the neighbourhood of the interface between the crust and the support. For instance, given a typical heat flux density of 1 × 105 W/m2 , the resulting temperature difference for a conductivity of 1 W/m/K is 100 K along 1 mm, which is, at least, the size of a thermocouple. Furthermore, the accuracy on the position of the thermocouple in the support is generally 1 mm or more.
4. Assessment against analytical and numerical solutions 4.1. Dam break for inertial flows This kind of calculation validates the basic mass and momentum equations of THEMA without viscous effects, and are representative of spreading at high mass flow rates. 4.1.1. 1D rectangular geometry The 1D solution of this problem is given by Stoker (1957). At time zero fluid with an initial thickness of 0.1 m starts to slump. The analytical solutions of the inertial flows provide successive thickness profiles versus time, with a common point at 4/9 = 0.444 of the initial thickness, and a constant velocity at this point. Figs. 3 and 4 show the consistency of THEMA results with theory. The differences with theory should be attributed to the space discretization in THEMA, and to the applicability of the shallow water hypotheses.
3. Assessment strategy The THEMA code was developed for safety purpose, in the frame of severe accidents of nuclear power plants. The assessment of the code must prove that it can be used with a sufficient level of confidence as a predictive tool for reactor cases. This goal implies the use of a frozen version of the code, with a fixed set of models, and no tuning parameters. It is believed that, even though all the phenomena that occur during the spreading of a mixture melt with solidification are not totally known, the major part of them is sufficiently understood and modelled in the code, so that predictive calculations for the reactor can be performed. With this respect, no tuning parameters optimised from experiments with simulating materials can be tolerated in the code, unless sensitivity studies are done within a wide range of these parameters. The assessment work includes comparisons with analytical or numerical solutions of basic problems which are presented in this first part. In a second part are given the results of the comparison of computations against experiments performed with simulating materials: the CORINE tests with water, water and viscous additives as well as low fusion temperature eutectic or non-eutectic mixtures of metals or oxide salts (up to 473 K). Comparisons with simulating material at high temperature (KATS and ECOKATS experiments at 2000 K) are also presented. Finally, results of assessment against experiments conducted with prototypic material (COMAS, FARO and VULCANO) are reported. All the presented results are obtained with the same code version (THEMA V2.5), and the same models. These models are consequently recommended for reactor calculations. Nevertheless, a specific study concerning the viscosity model in the solidification interval is generally made in the assessment. The results presented in previous papers (Spindler and Veteau, 1998; Spindler et al., 1999; Veteau et al., 2003) were obtained with previous versions of THEMA. This is the reason why they may slightly differ from those obtained with THEMA V2.5.
4.1.2. 2D cylindrical geometry At time 0, a cylindrical volume of fluid slumps. For an inertial flow, the self-similar solution by Fannelop and Waldman (1972) can be derived in terms of the Froude number at the spreading front. The spreading length is then given by the following relationship: Fr · g · Vol 0.25 √ r(t) = 2 · · t π · (4 − Fr) where Vol denotes the initial volume of the fluid, Fr = √Vg·e the Froude number and t is the time. V, g and e denote, respectively, the front velocity, the gravity and the depth of the fluid at the leading edge. Calculation with THEMA gives a mean value of Fr = 1.32 at the spreading front. The comparison on spreading length vs time with the analytical result (with Fr = 1.32) is presented in Fig. 5 and exhibits a very satisfactory agreement.
Fig. 3. 1D dam break problem. Free surface at different times computed with THEMA.
B. Spindler, J.M. Veteau / Nuclear Engineering and Design 236 (2006) 415–424
Fig. 4. 1D dam break problem. Velocity at x = 0 vs. time. Theory and THEMA calculations.
4.2. Conduction in the basemat This kind of conduction calculation validates the heat transfer module, which calculates the temperature field in the basemat. Refined meshes of 0.5 mm are used for all the THEMA simulations. 4.2.1. Conduction with imposed temperature The first considered analytical solution, given by Carslaw and Jaeger (1993), is the temperature field, versus abscissa x in a solid of thickness l and thermal diffusivity α at initial temperature T0 , with an imposed temperature T1 on one side: x 2 + · (T1 − T0 ) · l π ∞ n · π · x cos(n · π) × · sin · n l n=1 α · n2 · π 2 · t × exp − l2
T (x) = T0 + (T1 − T0 ) ·
421
Fig. 6. Conduction in a refractory basemat (d = 0.5 m) with imposed temperature on one side. Temperature vs. time at depths 0., 0.005, 0.010, 0.015, 0.020 and 0.050 m. Symbols refer to THEMA calculations, lines to the reference solution.
THEMA calculations presented in Fig. 6 correspond to a refractory basemat, for characteristic times of spreading. Similar computations with a metallic plate lead to similar agreement between the code and the reference solution. 4.2.2. Conduction and ablation with imposed temperature The second considered analytical solution, given by Carslaw and Jaeger (1993), is the temperature field in a semi-infinite solid at initial temperature T0 , with an imposed temperature T1 on one side, greater than the fusion temperature Tfus of the solid. An ablation front is then observed given by the following equation: X(t) = 2 · γ ·
√ αL · t
where γ is the solution of the following equation: √ kS · Tfus · αL · exp −γ 2 · ααLS exp(−γ 2 ) − √ erfc(γ) kL · (T1 − Tfus ) · αS · erfc γ · αL √
−
αS
γ · Lfus · π =0 CpL · (T1 − Tfus )
Lfus is the latent heat, k the conductivity, Cp the heat capacity, α the diffusivity and subscripts L and S refer, respectively, to the liquid and solid phases. Furthermore, the temperature versus abscissa x in the solid is given by: T x 1 · erfc T (x, t) = √ 2 · αL · t erfc γ · ααLS
Fig. 5. 2D dam break problem. Spreading length vs. time. Theory and THEMA calculations.
THEMA calculations presented here (Figs. 7 and 8) correspond to a refractory basemat with a diffusivity of 1 × 10−6 m2 /s, for times representative of spreading. As the analytical solutions are given for a semi-infinite solid, the simulations with THEMA are performed using a basemat with very large thickness (3 m). The same level of agreement was found for a metallic basemat.
422
B. Spindler, J.M. Veteau / Nuclear Engineering and Design 236 (2006) 415–424
Fig. 7. Conduction and ablation in a refractory basemat. Ablation front position vs. time. THEMA calculation and reference solution.
4.2.3. Conduction with imposed heat flux The third analytical solution regards the temperature field in a semi-infinite solid of diffusivity α and conductivity k, at initial temperature T0 , with an imposed heat flux Φ on one side obeys the following relationship (Carslaw and Jaeger, 1993): 2·Φ α·t x2 T (x) = T0 + · · exp − k π 4·α·t
x x √ − · erfc 2 2· α·t THEMA calculations presented at (Fig. 9) correspond to a ceramic basemat, for times and heat flux representative of spreading. Similar agreement with the reference solution was also found for a metallic basemat.
Fig. 9. Conduction with imposed heat flux in a refractory basemat. Temperature vs. time at depths 0, 0.01, 0.02, 0.03 and 0.05 m. Symbols refer to THEMA calculations, lines to the reference solution.
4.3.1. Crust growth on a basemat at constant temperature Savino and Siegel (1968) give the solution of the heat equation in the crust of conductivity k, heat capacity Cp and density ρ, in case of a basemat with constant temperature Tb . The crust thickness δ(t) corresponding to the solution at order 2 is obtained from the following equation: Cp · (TS − Tb ) h · (TL − TS ) · t = 1+ · ρ · L · δ0 3·L δ(t) δ(t) × − − ln 1 − δ0 δ0 where δ0 is the crust thickness corresponding to the steady state: δ0 =
k · (TS − Tb ) h · (TL − TS )
This kind of calculation validates the crust model used for the crust formed at the interface between melt and basemat.
where h is the heat transfer coefficient between the liquid at temperature TL and the crust, with an interfacial temperature TS and L is the latent heat. The comparison with THEMA (Fig. 10) gives very good agreement, since the same resolution approach is followed in the code and in the reference solution (Section 2.4.1).
Fig. 8. Conduction and ablation in a refractory basemat. Temperature of the basemat vs. time at 8 mm under the surface. THEMA calculated and reference solution.
Fig. 10. Crust growth on a basemat at constant temperature. Crust thickness vs. time for a reactor material. Symbols refer to THEMA calculations, lines to the reference solution.
4.3. Crust growth
B. Spindler, J.M. Veteau / Nuclear Engineering and Design 236 (2006) 415–424
423
temperature, necessary to get stable results in THEMA, as well as to reasonably small time steps (to ensure reasonable computing times). The basemat cell size is 0.5 mm. The time step ranges from 10−5 to 1 s depending of the variation of the heat flux and crust thickness. More accurate results could be obtained with a reduced time step and a reduced basemat cell size, but only the results with the standard numerical parameters are presented.
5. Conclusion
Fig. 11. Crust growth with conduction in the basemat. Crust thickness vs. time for a reactor melt on a refractory basemat. High and medium heat flux. Comparison between THEMA (T) calculations and the reference solution (RS).
4.3.2. Crust growth on a basemat with conduction Epstein (1976) gives the system of equations of the problem of the crust behaviour with conduction in the basemat with temperature varying with time. The system is solved numerically, assuming second-order polynomials for the temperature distribution in the crust and in the basemat. Then important quantities such as the maximum crust thickness and time of occurrence as well as time life of the crust are provided on charts as a function of non-dimensional numbers. For the purpose of the present comparison, the system was, again, numerically solved to provide the full data regarding the evolution of crust thickness versus time. The results presented here (Figs. 11 and 12) correspond to a reactor case (oxide melt on a refractory basemat). Several calculations were made with different values of the heat flux density between the fluid and the crust. THEMA generally overestimates the crust thickness by a factor of about 15%. This discrepancy is attributed to smoothing versus time of the crust-to-basemat heat flux and of the basemat
The THEMA code is devoted to the simulation of spreading of melts accompanied by solidification. Basically, the balance equations are integrated over the melt thickness. Crust formation on the basemat as well as at the melt surface is taken into account. The interface temperature between the melt and the crust is the solidus temperature of the melt mixture. A viscosity model depending on the solid fraction in the solidification interval is used. It is possible to simulate Bingham fluids, but the simulation of the melts as Newtonian fluids is preferred, since no data are available for a more complex model. The effect of gas issued from a concrete basemat is taken into account through a void fraction model depending on the superficial gas velocity, and a two-phase multiplier of the wall shear stress. The strategy of the assessment of THEMA code is to select a defined set of models (THEMA version 2.5), and to perform the simulations without any fitting of some parameters (for instance coefficients of the viscosity model or the contact resistance), because, by nature, these parameters cannot be fitted for a reactor case simulation. The analytical and numerical solutions are useful to assess separately the basic models of the code. In the case of the dam break problem, comparisons of the results obtained with analytical solutions for an inertial flow are made for a 1D and a 2D geometry. The 3D heat transfer module in the basemat was assessed using analytical solutions of heat conduction with or without melting of the solid. The crust growth on a basemat at constant temperature is compared with the solution of Savino and Siegel, and the problem of crust growth on a conducting basemat, is compared with the numerical solution of Epstein. Generally speaking, agreement between THEMA and the reference solutions can be perfect, but compromises between run time and accuracy of the results have to be done. This agreement remains fully acceptable when the numerical parameters (mesh size, time step) typical of a spreading calculation are considered.
Acknowledgements
Fig. 12. Crust growth with conduction in the basemat. Crust thickness vs. time for a reactor melt on a refractory basemat. Low value of the heat flux (0.25 MW/m2 ). Comparison between THEMA calculations and the reference solution.
THEMA code was developed within an agreement first between CEA, EDF and FRAMATOME, and then between CEA and EDF. A large part of the assessment work was performed in the frame of the Corium Spreading and Coolability (CSC) project of the 4th Framework Programme on Nuclear Fission of the European Community.
424
B. Spindler, J.M. Veteau / Nuclear Engineering and Design 236 (2006) 415–424
References Allelein, H.J., Breest, A., Spengler, C., 1999. Simulation of core melt spreading with LAVA: theoretical background and status of validation. In: OECD Workshop on Ex-Vessel Debris Coolability, Karlsruhe, Germany, 15–18 November. Carslaw, H.S., Jaeger, J.C., 1993. Conduction of Heat in Solid, second ed. Oxford University Press. Chong, J.S., Christiansen, E.B., Baer, A.D., 1971. Rheology of concentrated suspensions. J. Appl. Polym. Sci. 15, 2007–2021. Cognet, G., Seiler, J.M., Szabo, I., Latch´e, J.C., Spindler, B., Humbert, J.M., 1997. La r´ecup´eration du corium hors cuve. Revue G´en´erale Nucl´eaire 1, 38–43. Cognet, G., Alsmeyer, H., Tromm, W., Magallon, D., Wittmaack, R., Sehgal, B.R., Widmann, W., De Cecco, L., Ocelli, R., Azarian, G., Pineau, D., Spindler, B., Fieg, G., Werle, H., Journeau, C., Cranga, M., Laffont, G., 2001. Corium spreading and coolability CSC project. Nucl. Eng. Des. 209, 127–138. Duchatelle, L., Vautrey, J., 1964. D´etermination des coefficients de convection d’un alliage en e´ coulement turbulent entre plaques parall`eles. Int. J. Heat Mass Transfer 7, 1017–1031. Eberl´e, P., 1997. Mod´elisation physique et num´erique de l’´etalement d’un fluide avec solidification dans le cadre des e´ tudes de sˆuret´e pour les r´eacteurs a` eau sous pression. Th`ese de l’Universit´e Joseph Fourier Grenoble I. Engel, G., Fieg, G., Massier, H., Stegmaier, U., Sch¨utz, W., 1999. KATS experiments to simulate corium spreading in the EPR core catcher concept. In: OECD Workshop on Ex-Vessel Debris Coolability, Karlsruhe, Germany, 15–18 November. Epstein, M., 1976. The growth and decay of a frozen layer in forced flow. Int. J. Heat Mass Transfer 19, 1281–1288. Fannelop, T.K., Waldman, G.D., 1972. Dynamics of Oil Slicks. AIAA J. 10 (4), 506–510. Farmer, M.T., Sienicki, J.J., Chu, C.C., Spencer, B.W., 1993. The MELTSPREAD-1 Computer Code for the Analysis of Transient Spreading and Cooling of High Temperature Melts. EPRI Report TR-103413. Hewitt, G.F., 1982. In: Hetsroni, G. (Ed.), Handbook of Multiphase Systems. Chapter 2. Pressure Drop. Hemisphere Publishing Corporation, Mc GrawHill Book Company. Ishii, M., 1975. Thermo-fluid dynamic theory of two-phase flow, Eyrolles. Journeau, C., Boccaccio, E., Brayer, C., Cognet, G., Haquet, J.F., J´egou, C., Piluso, P., Monerris, J., 2003. Ex-vessel corium spreading: results
from the VULCANO spreading tests. Nucl. Eng. Des. 223, 75– 102. Piar, B., Michel, B.D., Babik, F., Latch´e, J.C., Guillard, G., Ruggi´eri, J.M., 1999. CROCO: a computer code for corium spreading. In: 9th International Topical Meeting on Nuclear Reactor Thermal Hydraulics (NURETH-9), San Francisco, California, 3–8 October. Pinkerton, H., Stevenson, R.J., 1992. Methods of determining the rheological properties of magmas at sub-liquidus temperatures. J. Volcanol. Geotherm. Res. 53, 47. Ramacciotti, M., 1999. Etude du comportement rh´eologique de m´elanges issus de l’interaction corium/b´eton. Th`ese de doctorat. Universit´e de Provence (Aix-Marseille I). Ramacciotti, M., Journeau, C., Sudreau, F., Cognet, G., 2001. Viscosity models for corium melts. Nucl. Eng. Des. 204, 377–389. Savino, J.M., Siegel, R., 1968. An analytical solution for solidification of a moving warm liquid onto an isothermal cold wall. Int. J. Heat Mass Transfer 12, 803–809. Shaw, H.R., 1969. Rheology of basalt in the melting range. J. Petrol. 10 (Part 3), 510–535. Spindler, B., Veteau, J.M., 1998. Status of the assessment of the spreading code THEMA against the CORINE experiments. In: SARJ Meeting, 4–6 November, Tokyo, Japan. Spindler, B., Veteau, J.M., 2004. Simulation of Spreading with Solidification: Assessment Synthesis of THEMA Code. CEA Report R-6053. Spindler, B., Veteau, J.M., Brayer, C., Cranga, M., De Cecco, L., Montanelli, P., Pineau, D., 1999. Assessment of THEMA code against spreading experiments. In: OECD Workshop on Ex-Vessel Debris Coolability, Karlsruhe, Germany, 15–18 November. Stedman, S.J., Evans, J.R.G., Woodthorpe, J., 1990. Rheology of composite ceramic injection moulding suspensions. J. Mater. Sci. 25, 1833–1841. Stoker, J.J., 1957. Water Waves. The Mathematical Theory with Applications. Interscience Publishers, Inc., New York, p. 313. Veteau, J.M., Spindler, B., Daum, G., 2003. Modelling of two-phase friction from isothermal spreading experiments with gas fed from the bottom and application to spreading accompanied by solidification. In: 10th International Topical Meeting on Nuclear Reactor Thermal Hydraulics (NURETH-10), Seoul, Korea, 5–9 October. Weissh¨aupl, H.A., 1999. Severe accident mitigation concept of EPR. Nucl. Eng. Des. 187, 35–45. Wittmaack, R., 1999. Numerical simulation of corium spreading in the EPR with CORFLOW. In: OECD Workshop on Ex-Vessel Debris Coolability, Karlsruhe, Germany, 15–18 November 1999.
The simulation of melt spreading with THEMA code Part 1: Model, assessment strategy and assessment against analytical and numerical solutions B. Spindler ∗ , J.M. Veteau Commissariat a` l’Energie Atomique, CEA-Grenoble, 17 rue des Martyrs, 38054 Grenoble, France Received 28 September 2004; accepted 16 September 2005
Abstract In the framework of R&D on severe accidents for PWRs , definition of mitigation means for cooling of the molten core after a vessel melt-through is of utmost importance. In the case of spreading-based concepts, the THEMA code, which considers integrated balance equations over the fluid depth, aims at predicting the spreading extent of non-eutectic mixtures under defined inlet conditions (pouring rate, temperature and composition). In this first part, the models and the assessment strategy of the code are presented. The satisfactory running of the main moduli of the code is then verified by comparing the numerical results against known solutions of basic problems: dam break problem, conduction with or without melting in the basemat, crust formation kinetics on a constant or “free” bottom plate temperature. © 2005 Elsevier B.V. All rights reserved.
1. Introduction The safety authorities of several countries have requested the industries and utilities to take into account severe accidents in the design of the next generation of nuclear power plants. Severe accidents refer to situations where high temperatures are reached in the core, leading to the formation of high temperature mixture of metals and oxides (corium). Among the mitigation means, which can be envisaged, special devices (core catchers) are considered, in which the corium is collected and cooled in order to insure the integrity of the vessel containment. A peculiar concept of core catcher, retained for instance for the European Pressurized Water Reactor (Weissh¨aupl, 1999), is based on the cooling of corium spread on a large area. Within this context, large research efforts are conducted throughout the world in order to get better understandings of the phenomena occurring during spreading of high temperature oxide and metal mixtures. CEA has undertaken a large program on severe accidents (Cognet et al., 1997), which aims at providing the tools dedicated to understand and simulate these phenomena.
∗
DOI of original article:10.1016/j.nucengdes.2005.09.016. Corresponding author. Tel.: +33 438 78 46 87; fax: +33 438 78 52 51. E-mail address: [email protected] (B. Spindler).
0029-5493/$ – see front matter © 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.nucengdes.2005.09.009
The development of THEMA code is part of this work. The THEMA code is devoted to the simulation of spreading, taking into account the melt solidification, as well as the possible melting of the basemat. The main characteristic of THEMA code is the use of balance equations integrated over the melt thickness. Integrated equations for the melt were chosen because it was stated that a very accurate formulation is of weak interest, since large uncertainties remain concerning the physical properties of real material on one hand, and for purpose of numerical simplicity and robustness, on the other hand. This model is also used in the MELTSPREAD code (Farmer et al., 1993), but limited to 1D geometry, whereas THEMA is able to simulate 1D or 2D geometry. The LAVA code also uses integrated equations, but inertia terms are not taken into account in the momentum equations (Allelein et al., 1999). Beside theses codes, some other codes consider a meshing within the melt thickness: CORFLOW developed by Framatome-ANP (Wittmaack, 1999), and CROCO developed by IRSN (Piar et al., 1999). The run time of such codes is very large. Nevertheless, computing time of THEMA is also increased for a 2D geometry, due to the 3D meshing used for computing heat conduction and melting in the basemat. In this first part, the models and the assessment strategy of THEMA code are first presented, followed by the results obtained in the frame work of comparison with known analytical
416
B. Spindler, J.M. Veteau / Nuclear Engineering and Design 236 (2006) 415–424
Nomenclature A coefficient C coefficient C0 coefficient Cp heat capacity (J kg−1 K−1 ) du/dy velocity gradient (s−1 ) e thickness (m) f friction factor or volumetric solid fraction F imposed heat flux (W m−2 ) H enthalpy (J kg−1 ) f0 volumetric solid fraction g gravity (m s−2 ) h heat transfer coefficient (W m−2 K) H enthalpy (J kg−1 ) J superficial velocity (m s−1 ) K heat conductivity L latent heat (J kg−1 ) S surface (m2 ) t time (s) T, T0 , T1 temperature (K) liquidus temperature (K) Tliq V velocity (m s−1 ) Vol volume (m3 ) drift flux velocity (m s−1 ) Vgj x abscissa (m) Greek letters α void fraction or diffusivity γ melting constant δ, δ0 crust thickness (m) Γ mass transfer term (kg s−1 ) µ viscosity (Pa s) µ0 viscosity for f = 0 σ surface tension (Pa m) τ0, τ shear stress (Pa)
or numerical solutions of basic problems, then ensuring a correct running of the main moduli of the code. The second part is devoted to the assessment results based on various spreading experiments: CORINE, KATS, ECOKATS, COMAS, FARO and VULCANO. Details of the assessment synthesis of THEMA can be found in a CEA report by Spindler and Veteau (2004). 2. THEMA models 2.1. Spatial distribution of the phases and temperature profile The spatial distribution of the solid and liquid phases considered in THEMA is presented in Fig. 1, together with the corresponding temperature profiles. The use of integrated equations corresponds to a single temperature for the melt. Transient temperature profiles are calculated in the bottom and the top crust and in the basemat. 2.2. Equations 2.2.1. Flow of a fluid with a floating crust The integrated equations were obtained from the well-posed complete system of equations related to the flow of a fluid composed with several phases (metal and oxide), with a solid crust above it. The basic integrated equations are those of the classical two-phase flows (Ishii, 1975) and the detailed process of the simplification of the system of equations is presented in the thesis of Eberl´e (1997). The objective of this work was to describe, together with the liquid flow, the flow of the solid crust at the upper surface, entrained by the liquid. However, for practical applications concerning corium, the mechanical properties of the crust at high temperature, which are necessary to close the model, are not known. Hence, it was decided to simplify the system equations, discarding the mechanical effects of the upper crust. Another reason for this simplification is that the effects of the upper crust in the spreading experiment are not clear. One
Non-dimensional numbers Fr Froude number Nu Nusselt number Pr Prandtl numer Re Reynolds number Subscripts b basemat bcr bottom crust c collapsed fus fusion g gas i interface L liquid S solid u upper
Fig. 1. Spatial distribution of the solid and the liquid phases and corresponding temperature profiles. Ts is the solidus temperature of the melt.
B. Spindler, J.M. Veteau / Nuclear Engineering and Design 236 (2006) 415–424
conclusion of the European CSC project (Cognet et al., 2001) was that, for the melts with a large solidification interval, there are no crust effects, while for the melts with a low solidification interval, the melt seems to be stopped by the formation of the crust. The investigations made in the KATS program (Engel et al., 1999) conclude that the upper crust effects occur only at the end of the spreading transient. 2.2.2. Main hypotheses The thickness of the melt is considered as low compared to the extent of the melt in the main direction of the flow. This hypothesis, known as the shallow water hypothesis, means that the pressure term in the momentum equation reduces to the gravity term, and hence the pressure term no more appears in the equation as a main variable. An other hypothesis which is just mentioned here, because practical applications were only performed with a single oxide or metal melt, is that the same velocity is used for the different liquid phases. 2.2.3. Set of equations A 1D or 2D, rectangular or cylindrical geometry can be used. After simplifications corresponding to the hypotheses stated above, the final set of equations is presented here, in case of a 1D rectangular geometry. Mass equation: ∂(e · ρ) ∂(e.ρ.V ) + = ∂t ∂X
Γ S
where e, ρ and V are, respectively, the thickness, the density and the velocity of the melt. S is the exchange area between the melt and the basemat and Γ is the mass transfer terms (from the melt to the crusts and from the ablated basemat to the melt). Momentum equation: ρ·
∂V ∂V ∂[ρ · (e + eb + ebcr )] +ρ·V · +g· ∂t ∂X ∂X
1 ρ · f · V2 + · =0 2 e where es and ebcr are the basemat and the bottom crust thickness, while f denotes the friction factor. Energy equation: ∂(e · ρ · H) ∂(e · ρ · V · H) + ∂t ∂x = hbcr · (T − Ti ) + hu · (T − Tu ) +
(Γ · H) S
where H is enthalpy of the melt, hb and hu the heat transfer coefficients, respectively, from the melt to the bottom crust and from the melt to the upper crust (or to the surface itself if no upper crust exist). Ti is the interfacial temperature between melt and crust, and Tu is either equal to Ti if there is an upper crust, or to the surface temperature of the melt. The last term corresponds to the enthalpy associated with the mass transfer
417
from the melt to the crusts, and from the ablated basemat to the melt. In case of gas issuing from the concrete basemat decomposition, a void fraction α is calculated. It is taken into account in the momentum equation, where the density becomes the two-phase density [=α·ρg + (1 − α)·ρL ], and e becomes the height of the swollen level. 2.2.4. Discretization and resolution The numerical integration of the finite difference set of equations uses a semi-implicit method with a Newton–Raphson iteration scheme with the melt height as main variable. A secondorder discretization in space is used, in order to get a better accuracy of the melt front progression. The time step is first limited by the Courant condition but is also controlled by the time variations of the wall heat flux and the crust thickness. Hence, the sharp variations of theses variables encountered during the thermal shock, occurring when the melt comes into contact with the cold basemat at the front, can be predicted with a controlled accuracy. 2.3. Limiting height at the leading edge The surface tension term is taken into account through a limiting height at the spreading front, obtained through the hypotheses of equilibrium between gravity and surface tension forces, and of a radius of curvature at the front equal to the height of the melt. Hence, the melt progression at the front is allowed only if the height of the melt is higher than the minimum height emin : 2·σ emin = ρ·g 2.4. Solidification and crust models Part of the melt is encrusted on the basemat on one hand, and at the upper surface on the other hand. The interface temperature between melt and crusts is the solidus temperature of the melt, since no thermodynamic equilibrium may occur within the short time of spreading. The crust has the same composition as the melt. Moreover, for non-eutectic mixtures, a specific viscosity model applies when the melt temperature decreases below the liquidus temperature, in order to depict the resulting semi-solid flow. 2.4.1. Crust models The transient conduction equation in the crust together with the energy balance at the interfaces (melt–basemat interface for the bottom crust; melt–surroundings for the upper crust) is solved using, for the temperature profile in the crust, an analytical iteration procedure from the steady state solution, as already proposed by Savino and Siegel (1968). For the bottom crust, the crust thickness and the heat flux from the crust to the basemat is obtained. For the surface crust, a set of two equations is
418
B. Spindler, J.M. Veteau / Nuclear Engineering and Design 236 (2006) 415–424
solved with two unknowns: the crust thickness and the surface temperature at the crust surface. As already pointed out, the present validation of THEMA is made without mechanical effects of the crusts. It appears that it is sufficient with regard to the accuracy that is required in the frame of the severe accident studies. 2.4.2. Viscosity models in the solidification interval A complete methodology for the calculation of the viscosity in the solidification interval is given by Ramacciotti (1999) and by Ramacciotti et al. (2001). The method of Ramacciotti gives a viscosity depending on the volumetric solid fraction f and a coefficient C laying between 4 and 8: µ = µ0 · exp(2.5 · C · f ) where µ0 is the viscosity of the liquid phase (f = 0). In the full set of THEMA calculations, the value used for µ0 is the viscosity of the liquid phase given by the experimentalists. It may also be predicted, for instance, with the Andrade model (Ramacciotti et al., 2001). Values of 6.1 and 6.3 for coefficient C were obtained, respectively, from the VULCANO tests VEU1 and VEU7 (Journeau et al., 2003). This formula is sometimes fitted (coefficient A in the formula hereunder) by the Shaw formula (Shaw, 1969) involving the melt temperature: µ = µ0 · exp[−A · (T − Tliq )] where Tliq is the liquidus temperature of the melt. The formulation of Ramacciotti using the solid fraction seems to be better, since for some mixtures the solid fraction variation versus temperature is not linear, and the equivalent viscosity is rather directly connected to the solid fraction than to the temperature. The first validations of THEMA were made with the model of Chong et al. (1971) better known as the Stedman model (Stedman et al., 1990): 0.75f µ = µ0 · 1 + f0 − f
with f0 = 0.5. That means that the value f = 0.5 corresponds to an infinite viscosity, i.e. to an immobilization of the melt. The Pinkerton–Stevenson viscosity model (Pinkerton and Stevenson, 1992) is a mix of the Stedman model (with f0 = 0.5), and the Shaw model: [−0.04 · (T − Tliq )] µ = µ0 · exp (1 − f/f0 )2.5 The comparison of the viscosity factor µ/µ0 given by the Stedman correlation, the Ramacciotti correlation (with C = 6) and the Pinkerton–Stevenson correlation (with f0 = 0.5) with a typical solidification interval of 500 K and with a linear relationship between temperature and solid fraction for simplification) is shown in Fig. 2. The Stedman correlation undergoes a slow increase of the viscosity for a low solid fraction, and a large increase when the value 0.5 is approached, while the Ramacciotti correlation provides also a large increase for small values of f. The Pinkerton–Stevenson is a combination of the two others. Concerning the solid fraction at which immobilization occurs, it is imposed with the parameter f0 for the Stedman and Pinkerton–Stevenson correlations. For the Ramacciotti correlation, immobilization occurs at a solid fraction not precisely determined, depending on the pouring rate. The first assessment of THEMA was made in the past using the Stedman correlation (Spindler et al., 1999). In the synthesis presented here, the calculations are made with the Stedman correlation with f0 = 0.5, and with the Ramacciotti correlation with C = 6 (corresponding to the mean value between 4 and 8, and to the value obtained by Journeau et al. (2003). The Stedman or Ramacciotti models are used for the mixtures of oxides only, while a specific model was developed for the non-eutectic mixtures of metals used in the CORINE experiments. For the spreading of pure metals (KATS iron tests) no semi-solid viscosity model is used, since there is no solidification range. 2.4.3. Simulation of Bingham fluids The validation of THEMA is made with the viscosity model described above, corresponding to a Newtonian fluid. Nevertheless, it is also possible to make simulations with THEMA using a Bingham fluid model, as in the CORFLOW and LAVA codes.
Fig. 2. Viscosity factor µ/µ0 vs. solid fraction, calculated with the Stedman correlation (with f0 = 0.5), the Ramacciotti correlation (with C = 6) and the Pinkerton–Stevenson correlation [with f0 = 0.5 and f = (Tliq − T)/500].
B. Spindler, J.M. Veteau / Nuclear Engineering and Design 236 (2006) 415–424
For Bingham fluids, the relationship between the shear stress τ and the velocity gradient du/dy is the following: τ = τ0 + du = 0, dy
µ du , dy
if τ > τ0 if τ ≤ τ0
where µ is the viscosity, as defined and used for a Newtonian fluid and τ 0 is the shear stress threshold. Practically, for a Bingham fluid, the wall shear stress coefficient is modified taken into account the supplementary parameter τ 0 . The complete formulation of the modified wall shear stress coefficient, in terms of the depth integrated fluid velocity, was introduced in THEMA. The comparison of simulations of the same spreading experiment, performed with the Newtonian model and with the Bingham model leads to the conclusion that the two results are very close. Since the parameter τ 0 is not well known, this kind of model cannot give better results than the Newtonian model. Furthermore, the selection of the viscosity model is often a way to match satisfactorily calculations with experiments. For this purpose, the Bingham model provides two fitting parameters, which is not desirable, even though corium is probably not a pure Newtonian fluid. Finally, the use of the Bingham model is not recommended. 2.5. Wall friction coefficient A classical correlation is used for the wall friction coefficient, based on a Reynolds number Re calculated with a hydraulic diameter equal to four times the melt thickness (geometry of a free surface flow, with a small thickness compared to the width). For low Reynolds numbers, the laminar friction factor is 24/Re, whereas the coefficient in the turbulent regime is 0.083 Re−0.25 . A more complex model, based on the boundary layer theory, was also implemented in THEMA (Eberl´e, 1997). But generally, the boundary layer thickness rapidly reaches the melt thickness, and the results are then quite similar (Spindler and Veteau, 1998). As already pointed out, the gas effect is taken into account through an increase of the melt thickness. Furthermore, a modification of the wall shear stress coefficient is approached by the concept of two-phase friction multipliers, already worked out for two-phase flows in pipes. Details of the model used for spreading with gas percolation are given by Veteau et al. (2003). The two-phase friction multipliers depend on the ratio of the fluid to the gas viscosity. The following correlations are used: • For
µ µgas
≤ 103 : Friedel correlation (Hewitt, 1982).
• For 103 < µµgas ≤ 105 : Lockhart–Martinelli correlation (Hewitt, 1982). • For µµgas > 105 : multiplier is (1 − α)2 deduced from the CORINE experiments. Details of the adaptation for free surface flows are given by Veteau et al. (2003).
419
2.6. Heat transfer coefficients Classical correlations are also used for the heat transfer coefficients. For laminar flows the Nusselt number is equal to 7.6, and for turbulent flow the well-known Colburn correlation Nu = 0.023 Re0.8 Pr0.4 is used for oxide melts. For metallic melts, specific correlations are used (Duchatelle and Vautrey, 1964): Nu = 6.0 for laminar flows, and Nu = 5.14 + 0.0127 Re0.8 Pr0.8 for turbulent flows. No specific effect of the presence of gas is considered in the heat transfer coefficients, despite it is known that gas bubbling increases the heat transfer. 2.7. Void fraction model A void fraction model is required in the two-phase friction multiplier approach. This quantity is deduced from special experiments conducted with the simulating materials used in the isothermal CORINE spreading tests with gas fed from the bottom. The data were analysed in the frame of the drift-flux void fraction model: α=
1 C0 +
Vgj Jg
where α is the volumetric void fraction and Jg is the superficial gas velocity. The parameters C0 and Vgj have been correlated against the collapsed melt thickness ec (Veteau et al., 2003). For water as working fluid, the following correlation was found: C0 = 26.8ec + 0.8224
with ec in m, and Vgj = 0.043 m/s.
These values are used not only for the simulation of the CORINE tests with water, but also for high temperature tests with freezing (oxide and iron KATS tests on a concrete basemat, VULCANO tests) and are recommended for reactor applications. Similar but different correlations have been developed to simulate the CORINE isothermal tests conducted with viscous fluids. 2.8. Heat transfer in the basemat A specific 3D heat transfer module is used, solving the heat transfer equation in the solid basemat. A fine meshing of the basemat is used near the melt–basemat or crust–basemat interface, a region where the temperature variations in the basemat are large. The heat of fusion is taken into account in the specific energy of the solid material, which enables the calculation of the basemat ablation. The gas release in a concrete basemat is taken into account: bounded water and carbon dioxide release when the temperature of the meshes reaches 373 and 746 K, respectively. The corresponding latent heat and decomposition heat are considered with a modification of the concrete specific heat at these temperatures. A thermal resistance may also be included between the basemat and the melt or crust, but in absence of a sufficiently accurate value of this parameter, it is recommended to set this value
420
B. Spindler, J.M. Veteau / Nuclear Engineering and Design 236 (2006) 415–424
at zero. Furthermore, the sensitivity studies that were made show that the effect of the thermal resistance on the spreading length is low for the spreading of oxide mixtures, but large for iron or eutectic metallic melts. Some sensitivity studies are presented in the second part, which show the effect of the thermal resistance, and how it can be (improperly) used as a fitting parameter. The determination of the thermal resistance from matching the computed and experimental temperature profiles in the basemat from standard spreading experiments is quite hazardous, especially with low conductivity materials like ceramic or concrete, owing to the very large temperature gradient in the neighbourhood of the interface between the crust and the support. For instance, given a typical heat flux density of 1 × 105 W/m2 , the resulting temperature difference for a conductivity of 1 W/m/K is 100 K along 1 mm, which is, at least, the size of a thermocouple. Furthermore, the accuracy on the position of the thermocouple in the support is generally 1 mm or more.
4. Assessment against analytical and numerical solutions 4.1. Dam break for inertial flows This kind of calculation validates the basic mass and momentum equations of THEMA without viscous effects, and are representative of spreading at high mass flow rates. 4.1.1. 1D rectangular geometry The 1D solution of this problem is given by Stoker (1957). At time zero fluid with an initial thickness of 0.1 m starts to slump. The analytical solutions of the inertial flows provide successive thickness profiles versus time, with a common point at 4/9 = 0.444 of the initial thickness, and a constant velocity at this point. Figs. 3 and 4 show the consistency of THEMA results with theory. The differences with theory should be attributed to the space discretization in THEMA, and to the applicability of the shallow water hypotheses.
3. Assessment strategy The THEMA code was developed for safety purpose, in the frame of severe accidents of nuclear power plants. The assessment of the code must prove that it can be used with a sufficient level of confidence as a predictive tool for reactor cases. This goal implies the use of a frozen version of the code, with a fixed set of models, and no tuning parameters. It is believed that, even though all the phenomena that occur during the spreading of a mixture melt with solidification are not totally known, the major part of them is sufficiently understood and modelled in the code, so that predictive calculations for the reactor can be performed. With this respect, no tuning parameters optimised from experiments with simulating materials can be tolerated in the code, unless sensitivity studies are done within a wide range of these parameters. The assessment work includes comparisons with analytical or numerical solutions of basic problems which are presented in this first part. In a second part are given the results of the comparison of computations against experiments performed with simulating materials: the CORINE tests with water, water and viscous additives as well as low fusion temperature eutectic or non-eutectic mixtures of metals or oxide salts (up to 473 K). Comparisons with simulating material at high temperature (KATS and ECOKATS experiments at 2000 K) are also presented. Finally, results of assessment against experiments conducted with prototypic material (COMAS, FARO and VULCANO) are reported. All the presented results are obtained with the same code version (THEMA V2.5), and the same models. These models are consequently recommended for reactor calculations. Nevertheless, a specific study concerning the viscosity model in the solidification interval is generally made in the assessment. The results presented in previous papers (Spindler and Veteau, 1998; Spindler et al., 1999; Veteau et al., 2003) were obtained with previous versions of THEMA. This is the reason why they may slightly differ from those obtained with THEMA V2.5.
4.1.2. 2D cylindrical geometry At time 0, a cylindrical volume of fluid slumps. For an inertial flow, the self-similar solution by Fannelop and Waldman (1972) can be derived in terms of the Froude number at the spreading front. The spreading length is then given by the following relationship: Fr · g · Vol 0.25 √ r(t) = 2 · · t π · (4 − Fr) where Vol denotes the initial volume of the fluid, Fr = √Vg·e the Froude number and t is the time. V, g and e denote, respectively, the front velocity, the gravity and the depth of the fluid at the leading edge. Calculation with THEMA gives a mean value of Fr = 1.32 at the spreading front. The comparison on spreading length vs time with the analytical result (with Fr = 1.32) is presented in Fig. 5 and exhibits a very satisfactory agreement.
Fig. 3. 1D dam break problem. Free surface at different times computed with THEMA.
B. Spindler, J.M. Veteau / Nuclear Engineering and Design 236 (2006) 415–424
Fig. 4. 1D dam break problem. Velocity at x = 0 vs. time. Theory and THEMA calculations.
4.2. Conduction in the basemat This kind of conduction calculation validates the heat transfer module, which calculates the temperature field in the basemat. Refined meshes of 0.5 mm are used for all the THEMA simulations. 4.2.1. Conduction with imposed temperature The first considered analytical solution, given by Carslaw and Jaeger (1993), is the temperature field, versus abscissa x in a solid of thickness l and thermal diffusivity α at initial temperature T0 , with an imposed temperature T1 on one side: x 2 + · (T1 − T0 ) · l π ∞ n · π · x cos(n · π) × · sin · n l n=1 α · n2 · π 2 · t × exp − l2
T (x) = T0 + (T1 − T0 ) ·
421
Fig. 6. Conduction in a refractory basemat (d = 0.5 m) with imposed temperature on one side. Temperature vs. time at depths 0., 0.005, 0.010, 0.015, 0.020 and 0.050 m. Symbols refer to THEMA calculations, lines to the reference solution.
THEMA calculations presented in Fig. 6 correspond to a refractory basemat, for characteristic times of spreading. Similar computations with a metallic plate lead to similar agreement between the code and the reference solution. 4.2.2. Conduction and ablation with imposed temperature The second considered analytical solution, given by Carslaw and Jaeger (1993), is the temperature field in a semi-infinite solid at initial temperature T0 , with an imposed temperature T1 on one side, greater than the fusion temperature Tfus of the solid. An ablation front is then observed given by the following equation: X(t) = 2 · γ ·
√ αL · t
where γ is the solution of the following equation: √ kS · Tfus · αL · exp −γ 2 · ααLS exp(−γ 2 ) − √ erfc(γ) kL · (T1 − Tfus ) · αS · erfc γ · αL √
−
αS
γ · Lfus · π =0 CpL · (T1 − Tfus )
Lfus is the latent heat, k the conductivity, Cp the heat capacity, α the diffusivity and subscripts L and S refer, respectively, to the liquid and solid phases. Furthermore, the temperature versus abscissa x in the solid is given by: T x 1 · erfc T (x, t) = √ 2 · αL · t erfc γ · ααLS
Fig. 5. 2D dam break problem. Spreading length vs. time. Theory and THEMA calculations.
THEMA calculations presented here (Figs. 7 and 8) correspond to a refractory basemat with a diffusivity of 1 × 10−6 m2 /s, for times representative of spreading. As the analytical solutions are given for a semi-infinite solid, the simulations with THEMA are performed using a basemat with very large thickness (3 m). The same level of agreement was found for a metallic basemat.
422
B. Spindler, J.M. Veteau / Nuclear Engineering and Design 236 (2006) 415–424
Fig. 7. Conduction and ablation in a refractory basemat. Ablation front position vs. time. THEMA calculation and reference solution.
4.2.3. Conduction with imposed heat flux The third analytical solution regards the temperature field in a semi-infinite solid of diffusivity α and conductivity k, at initial temperature T0 , with an imposed heat flux Φ on one side obeys the following relationship (Carslaw and Jaeger, 1993): 2·Φ α·t x2 T (x) = T0 + · · exp − k π 4·α·t
x x √ − · erfc 2 2· α·t THEMA calculations presented at (Fig. 9) correspond to a ceramic basemat, for times and heat flux representative of spreading. Similar agreement with the reference solution was also found for a metallic basemat.
Fig. 9. Conduction with imposed heat flux in a refractory basemat. Temperature vs. time at depths 0, 0.01, 0.02, 0.03 and 0.05 m. Symbols refer to THEMA calculations, lines to the reference solution.
4.3.1. Crust growth on a basemat at constant temperature Savino and Siegel (1968) give the solution of the heat equation in the crust of conductivity k, heat capacity Cp and density ρ, in case of a basemat with constant temperature Tb . The crust thickness δ(t) corresponding to the solution at order 2 is obtained from the following equation: Cp · (TS − Tb ) h · (TL − TS ) · t = 1+ · ρ · L · δ0 3·L δ(t) δ(t) × − − ln 1 − δ0 δ0 where δ0 is the crust thickness corresponding to the steady state: δ0 =
k · (TS − Tb ) h · (TL − TS )
This kind of calculation validates the crust model used for the crust formed at the interface between melt and basemat.
where h is the heat transfer coefficient between the liquid at temperature TL and the crust, with an interfacial temperature TS and L is the latent heat. The comparison with THEMA (Fig. 10) gives very good agreement, since the same resolution approach is followed in the code and in the reference solution (Section 2.4.1).
Fig. 8. Conduction and ablation in a refractory basemat. Temperature of the basemat vs. time at 8 mm under the surface. THEMA calculated and reference solution.
Fig. 10. Crust growth on a basemat at constant temperature. Crust thickness vs. time for a reactor material. Symbols refer to THEMA calculations, lines to the reference solution.
4.3. Crust growth
B. Spindler, J.M. Veteau / Nuclear Engineering and Design 236 (2006) 415–424
423
temperature, necessary to get stable results in THEMA, as well as to reasonably small time steps (to ensure reasonable computing times). The basemat cell size is 0.5 mm. The time step ranges from 10−5 to 1 s depending of the variation of the heat flux and crust thickness. More accurate results could be obtained with a reduced time step and a reduced basemat cell size, but only the results with the standard numerical parameters are presented.
5. Conclusion
Fig. 11. Crust growth with conduction in the basemat. Crust thickness vs. time for a reactor melt on a refractory basemat. High and medium heat flux. Comparison between THEMA (T) calculations and the reference solution (RS).
4.3.2. Crust growth on a basemat with conduction Epstein (1976) gives the system of equations of the problem of the crust behaviour with conduction in the basemat with temperature varying with time. The system is solved numerically, assuming second-order polynomials for the temperature distribution in the crust and in the basemat. Then important quantities such as the maximum crust thickness and time of occurrence as well as time life of the crust are provided on charts as a function of non-dimensional numbers. For the purpose of the present comparison, the system was, again, numerically solved to provide the full data regarding the evolution of crust thickness versus time. The results presented here (Figs. 11 and 12) correspond to a reactor case (oxide melt on a refractory basemat). Several calculations were made with different values of the heat flux density between the fluid and the crust. THEMA generally overestimates the crust thickness by a factor of about 15%. This discrepancy is attributed to smoothing versus time of the crust-to-basemat heat flux and of the basemat
The THEMA code is devoted to the simulation of spreading of melts accompanied by solidification. Basically, the balance equations are integrated over the melt thickness. Crust formation on the basemat as well as at the melt surface is taken into account. The interface temperature between the melt and the crust is the solidus temperature of the melt mixture. A viscosity model depending on the solid fraction in the solidification interval is used. It is possible to simulate Bingham fluids, but the simulation of the melts as Newtonian fluids is preferred, since no data are available for a more complex model. The effect of gas issued from a concrete basemat is taken into account through a void fraction model depending on the superficial gas velocity, and a two-phase multiplier of the wall shear stress. The strategy of the assessment of THEMA code is to select a defined set of models (THEMA version 2.5), and to perform the simulations without any fitting of some parameters (for instance coefficients of the viscosity model or the contact resistance), because, by nature, these parameters cannot be fitted for a reactor case simulation. The analytical and numerical solutions are useful to assess separately the basic models of the code. In the case of the dam break problem, comparisons of the results obtained with analytical solutions for an inertial flow are made for a 1D and a 2D geometry. The 3D heat transfer module in the basemat was assessed using analytical solutions of heat conduction with or without melting of the solid. The crust growth on a basemat at constant temperature is compared with the solution of Savino and Siegel, and the problem of crust growth on a conducting basemat, is compared with the numerical solution of Epstein. Generally speaking, agreement between THEMA and the reference solutions can be perfect, but compromises between run time and accuracy of the results have to be done. This agreement remains fully acceptable when the numerical parameters (mesh size, time step) typical of a spreading calculation are considered.
Acknowledgements
Fig. 12. Crust growth with conduction in the basemat. Crust thickness vs. time for a reactor melt on a refractory basemat. Low value of the heat flux (0.25 MW/m2 ). Comparison between THEMA calculations and the reference solution.
THEMA code was developed within an agreement first between CEA, EDF and FRAMATOME, and then between CEA and EDF. A large part of the assessment work was performed in the frame of the Corium Spreading and Coolability (CSC) project of the 4th Framework Programme on Nuclear Fission of the European Community.
424
B. Spindler, J.M. Veteau / Nuclear Engineering and Design 236 (2006) 415–424
References Allelein, H.J., Breest, A., Spengler, C., 1999. Simulation of core melt spreading with LAVA: theoretical background and status of validation. In: OECD Workshop on Ex-Vessel Debris Coolability, Karlsruhe, Germany, 15–18 November. Carslaw, H.S., Jaeger, J.C., 1993. Conduction of Heat in Solid, second ed. Oxford University Press. Chong, J.S., Christiansen, E.B., Baer, A.D., 1971. Rheology of concentrated suspensions. J. Appl. Polym. Sci. 15, 2007–2021. Cognet, G., Seiler, J.M., Szabo, I., Latch´e, J.C., Spindler, B., Humbert, J.M., 1997. La r´ecup´eration du corium hors cuve. Revue G´en´erale Nucl´eaire 1, 38–43. Cognet, G., Alsmeyer, H., Tromm, W., Magallon, D., Wittmaack, R., Sehgal, B.R., Widmann, W., De Cecco, L., Ocelli, R., Azarian, G., Pineau, D., Spindler, B., Fieg, G., Werle, H., Journeau, C., Cranga, M., Laffont, G., 2001. Corium spreading and coolability CSC project. Nucl. Eng. Des. 209, 127–138. Duchatelle, L., Vautrey, J., 1964. D´etermination des coefficients de convection d’un alliage en e´ coulement turbulent entre plaques parall`eles. Int. J. Heat Mass Transfer 7, 1017–1031. Eberl´e, P., 1997. Mod´elisation physique et num´erique de l’´etalement d’un fluide avec solidification dans le cadre des e´ tudes de sˆuret´e pour les r´eacteurs a` eau sous pression. Th`ese de l’Universit´e Joseph Fourier Grenoble I. Engel, G., Fieg, G., Massier, H., Stegmaier, U., Sch¨utz, W., 1999. KATS experiments to simulate corium spreading in the EPR core catcher concept. In: OECD Workshop on Ex-Vessel Debris Coolability, Karlsruhe, Germany, 15–18 November. Epstein, M., 1976. The growth and decay of a frozen layer in forced flow. Int. J. Heat Mass Transfer 19, 1281–1288. Fannelop, T.K., Waldman, G.D., 1972. Dynamics of Oil Slicks. AIAA J. 10 (4), 506–510. Farmer, M.T., Sienicki, J.J., Chu, C.C., Spencer, B.W., 1993. The MELTSPREAD-1 Computer Code for the Analysis of Transient Spreading and Cooling of High Temperature Melts. EPRI Report TR-103413. Hewitt, G.F., 1982. In: Hetsroni, G. (Ed.), Handbook of Multiphase Systems. Chapter 2. Pressure Drop. Hemisphere Publishing Corporation, Mc GrawHill Book Company. Ishii, M., 1975. Thermo-fluid dynamic theory of two-phase flow, Eyrolles. Journeau, C., Boccaccio, E., Brayer, C., Cognet, G., Haquet, J.F., J´egou, C., Piluso, P., Monerris, J., 2003. Ex-vessel corium spreading: results
from the VULCANO spreading tests. Nucl. Eng. Des. 223, 75– 102. Piar, B., Michel, B.D., Babik, F., Latch´e, J.C., Guillard, G., Ruggi´eri, J.M., 1999. CROCO: a computer code for corium spreading. In: 9th International Topical Meeting on Nuclear Reactor Thermal Hydraulics (NURETH-9), San Francisco, California, 3–8 October. Pinkerton, H., Stevenson, R.J., 1992. Methods of determining the rheological properties of magmas at sub-liquidus temperatures. J. Volcanol. Geotherm. Res. 53, 47. Ramacciotti, M., 1999. Etude du comportement rh´eologique de m´elanges issus de l’interaction corium/b´eton. Th`ese de doctorat. Universit´e de Provence (Aix-Marseille I). Ramacciotti, M., Journeau, C., Sudreau, F., Cognet, G., 2001. Viscosity models for corium melts. Nucl. Eng. Des. 204, 377–389. Savino, J.M., Siegel, R., 1968. An analytical solution for solidification of a moving warm liquid onto an isothermal cold wall. Int. J. Heat Mass Transfer 12, 803–809. Shaw, H.R., 1969. Rheology of basalt in the melting range. J. Petrol. 10 (Part 3), 510–535. Spindler, B., Veteau, J.M., 1998. Status of the assessment of the spreading code THEMA against the CORINE experiments. In: SARJ Meeting, 4–6 November, Tokyo, Japan. Spindler, B., Veteau, J.M., 2004. Simulation of Spreading with Solidification: Assessment Synthesis of THEMA Code. CEA Report R-6053. Spindler, B., Veteau, J.M., Brayer, C., Cranga, M., De Cecco, L., Montanelli, P., Pineau, D., 1999. Assessment of THEMA code against spreading experiments. In: OECD Workshop on Ex-Vessel Debris Coolability, Karlsruhe, Germany, 15–18 November. Stedman, S.J., Evans, J.R.G., Woodthorpe, J., 1990. Rheology of composite ceramic injection moulding suspensions. J. Mater. Sci. 25, 1833–1841. Stoker, J.J., 1957. Water Waves. The Mathematical Theory with Applications. Interscience Publishers, Inc., New York, p. 313. Veteau, J.M., Spindler, B., Daum, G., 2003. Modelling of two-phase friction from isothermal spreading experiments with gas fed from the bottom and application to spreading accompanied by solidification. In: 10th International Topical Meeting on Nuclear Reactor Thermal Hydraulics (NURETH-10), Seoul, Korea, 5–9 October. Weissh¨aupl, H.A., 1999. Severe accident mitigation concept of EPR. Nucl. Eng. Des. 187, 35–45. Wittmaack, R., 1999. Numerical simulation of corium spreading in the EPR with CORFLOW. In: OECD Workshop on Ex-Vessel Debris Coolability, Karlsruhe, Germany, 15–18 November 1999.