F Guillard

Jacq

T Glantz
email: tony.glantz@irsn.fr

T Taurines

O De Luze

S Belon

G Guillard

F Jacq

DRACCAR: a multi-physics code for computational analysis of multi-rod ballooning and fuel relocation during LOCA transients Part one: General modeling description

Keywords:

HAL is

Introduction

Background and Context

The loss of coolant accident (LOCA) in a light water reactor subsequent to a break or a leakage on the primary circuit were and are still largely investigated in the framework of nuclear research due to their possible consequences regarding nuclear safety. For instance, core uncovery can lead to an increase of the fuel rod temperatures which can rise enough to lead to their deformation and possible damage. In the meantime, the progressive blockage of the core channel can seriously impact

The core management optimization of the French power reactors (increase of flexibility and availability, cost saving) leads the operators to increase progressively the burn-up of fuel and to introduce new types of fuel material. These evolutions involve the development of advanced cladding materials (M5 R , Zirlo TM ) and fuel micro-structures (gadolinium-doped fuel, MOX). In the meantime, efforts to improve simulation of LOCA lead to the development of more precise simulation tools (3D modeling, local models) and the use of calculation methods involving uncertainties to assess the safety demonstration.

IRSN conducted an extensive State-of-the-Art-Review relative to fuel behavior under Loss Of Coolant Accident (LOCA) conditions, covering the aspects of clad ballooning and flow blockage, coolability of partially blocked assemblies, clad oxidation and clad resistance to quench and post-quench loads [START_REF] Grandjean | A state-of-the-art review of past programs devoted to fuel behavior under LOCA conditions -part 1: Clad swelling and rupture -assembly flow blockage[END_REF][START_REF] Grandjean | A state-of-the-art review of past programs devoted to fuel behavior under LOCA conditions -part 2: Impact of clad swelling upon assembly cooling[END_REF]). Additionally, a review of existing computational tools devoted to the calculation of the fuel behavior under LOCA was performed [START_REF] Cunningham | FRAP-TRAN: A computer code for transient analysis of oxide fuel rods[END_REF]; [START_REF] Siefken | FRAP-T6: A computer code for the transient analysis of oxide fuel rod[END_REF]; [START_REF] Casadei | RODSWELL: A computer code for the thermo mechanical analysis of fuel rods under LOCA conditions[END_REF]; [START_REF] Haste | CANSWEL-2: A computer model for the creep deformation of Zircaloy cladding under Loss-of-Coolant Accident part 1 -model description[END_REF]; [START_REF] Uchida | Models of multi-rod code FRETA-B for transient fuel behavior analysis[END_REF]). Based on these reviews, IRSN has decided to launch a Research and Development program with the objectives of:

-Increasing the knowledge and reducing the uncertainties on the prediction of the thermo-mechanical behavior of the fuel rods and of the thermalhydraulics within an assembly in LOCA conditions, -Simulating the LOCA phenomenology taken into account all parameters: fuel management (burnup, linear power, ...), fuel assembly type (pellet, clad, grid, ...) and accident scenario (large break, intermediate break, ...).

The IRSN Research and Development program is composed of two parts, the first one dedicated to the collection and interpretation of complementary experimental data. The second part is dedicated to the development and validation (using experimental material obtained in the first part of the program) of a simulation tool, named DRACCAR and capable of simulating the reactor core during a LOCA. The experimental data interpretation and the use of the simulation tool will allow:

-Performing an analysis of the safety criteria relevance for the current and future fuel management, -Giving technical advice on revising the LOCA safety criteria, -Evaluating the safety margin compared with the current safety criteria.

The methodology widely used for the LOCA safety assessment is to evaluate the peak cladding temperature (PCT) using a system code whose description is limited to model the thermo-mechanical behavior of one average weighted fuel rod per core ring. In addition pessimistic conditions and assumptions are cumulated to ensure an over-estimation of the PCT.

However such approach cannot deal properly with the presence and the behavior of neighboring rods or with the complex flows induced by the progressive blockage of the fuel channels. The aim of the DRACCAR code is to realistically describe the 3D thermo-mechanical behavior and reflooding of a fuel assembly including its coolability as well as its embrittlement during a LOCA. Changes in the core geometry, as well as the ability to cool the assemblies are addressed and analyzed. This analysis is essential for a detailed understanding of circumstances and conditions for the assemblies quenching even in the case of a large flow blockage ratio due to clad ballooning. This implies complex computations involving two-phase flows coupled with heat transfers (decay heat, conduction, convection and radiation), mechanical behavior (creep behavior, thermal expansion, contact, oxidation and fuel relocation).

The flexibility of the DRACCAR code makes possible to model from a simple fuel rod up to a full core, and surrounding shrouds. In order to achieve this goal, DRACCAR allows the modeling of any kind of core components (pellet, cladding, grids, ...), fuel (UO 2 , MOX, ...), cladding (Zircaloy 4, M5 R , Zirlo TM , ...), core management (burn-up), and types of watercooled reactor including loops and safety systems.

Moreover, since the Fukushima accident, increased attention has been paid to the vulnerability of the Spent Fuel Pools (SFP). This vulnerability associated with the dramatic radiological consequences of an unmitigated accident is of a major concern for nuclear safety. Those consequences are related to the location of SFPs outside the reactor containment building and only the cladding acts as a confinment barrier. Accidents involving SFPs can be the result of a loss of heat removal system or a loss of coolant due to a breach in the pool structure or in a connected fluid circuit. To improve knowledge on the fuel behavior during these accident situations, IRSN has decided to complete its Research and Development program. Several experiments and specific developments in DRACCAR dedicated to SFP accident studies have been defined [START_REF] Mutelle | A new research program on accidents in spent fuel pools: the DENOPI project[END_REF]).

DRACCAR: a coupling platform

Modular approach

The main objective of DRACCAR is to describe the fuel assemblies' behavior during a LOCA on a reactor or on a spent fuel pool. Nevertheless, during such scenario, it is important in order to estimate well this behavior to:

-calculate the thermalhydraulics in all the reactor circuits in particular in the core or in the spent fuel pool, -simulate the auxiliary system which could be used during the accident, -calculate the residual fuel power evolution as a function of the fission products.

To achieve this main objective, DRACCAR couples two codes, the first one describes the thermo-mechanical behavior of the fuel assemblies and the second one the thermalhydraulics in these assemblies. The fluid code is also capable of simulating the thermalhydraulics in all the reactor circuits or the thermalhydraulics in a spent fuel pool (to some extent). Several additional modules are available in the DRACCAR platform. The -thermal, mechanical and chemical properties of all structural materials which could be found in a light water nuclear reactor, -all steam and liquid water properties, which are validated from 0 to 4000 bar for the pressure range and from 0 to 1000˚C for the temperature range. They can be extrapolated to 4000˚C.

-properties of fission products, -radioactive properties of all possible isotopes, -chemical reactions involving iodine in a containment, -properties of complex mixtures (corium/concrete interaction for instance) which could appear in the vessel (U-O-Zr-Fe) or in the cavity (U-O-Zr-Fe-Si-Ca-Al...). Moreover, MDB is sufficiently flexible to allow the modification of any properties law or to define a material not referenced in the data bank. In particular, it is possible for an existing material or a new one (for example, a specific cladding material) to modify or define the creep rate law (the anisotropy and the phase change can be taken into account), the burst criteria and the oxidation law. DRACCAR is then applicable to most of the existing and up-coming cladding materials.

In addition to the physical modelling, the management of engineered safety features can be simulated by SYSINT module, such as hydro-accumulators, valves or safety injection system as well as the management of pressurizer and steam generators.

It also allows the user to define in detail an accident scenario, with the succession of postulated events and operating safety procedures.

The decay of fission products and actinide isotopes in the reactor core is simulated by the ISODOP module. It starts the calculation using an initial isotope inventory (generated by another independent code) and estimates decay heat and activity in the core to be estimated. It is based on the JEFF database containing a description of all isotopes JEFF-3.1 (2006).

ODESSA is a set of tools dedicated to the data managements in scientific softwares. It offers a specific data structure (ODESSA database) designed for efficient and fast exchanges between parts of codes. It is composed of several tools useful for reading, completing and checking users input data decks also. ODESSA fulfills too tasks dedicated to the processing of simulation results, in particular it is possible to visualize on-line the results during the calculation performance.

The several tools coupled in DRACCAR can also run in stand-alone mode. Each of them has its own syntax which can be very different as well as the input data decks. In order to make DRACCAR more user-friendly, the coupling platform has a special interface, which role is to generate the appropriate input data decks of each code from a unique input data deck, written in the ODESSA language. This method presents several advantages:

-the use of a single language in the input data deck to perform an ICARE3D/CESAR or an ICARE3D/CATHARE-3 calculation reduces the knowledge necessary for the user, -the use of the ODESSA language offers compatibility to the peripheric tools, in particular for pre-and postprocessing, -the structure and the fluid meshings are automatically generated together in the same single mesh creating phase even for an ICARE3D/CATHARE-3 simulation. This ensures the consistency and relation between the two meshings needed to perform the coupling whatever the thermalhydraulic code used.

Modeling concept

Concerning the modeling, only the core part is concerned by the coupling of structure code ICARE3D to the thermalhydraulic code. In case of reactor simulation, the circuits (fluid and structures) are treated only by the thermalhydraulic code.

The circuits discretization is an assembly of zero-dimensional (0-D) elements (homogeneous or with a swollen level for twophase flow situations) and one-dimensional (1-D) elements (axial mesh) connected with junction elements.

In the core, the flexibility of DRACCAR allows modeling precisely from one rod to a fuel assembly, the rest of the core or the pool being simulated by an average weighted modeling. According to the existing symmetry in fuel assemblies, a reasonable target for fuel assembly calculations is to model an eighth of a fuel assembly including fuel rods (stack of fuel pellets, claddings), control rods, guide tubes, instrumentation tubes, spacer grid. Each structure is in mechanical and thermal interaction with other ones, including contacts between fuel rods and eventually with guide tubes or spacer grids. The structure geometry evolution during the transient has directly an impact on the fluid (for example, mechanical deformation of a rod leads to a reduction of the fluid flow section). Two meshings are generated by the DRACCAR code from the ODESSA input data deck: the structure meshing (figure 2) and the sub-channel meshing (figure 3):

-The structure has a 3D description and a r-θ-z meshing type, -For the sub-channel, a 3D unstructured meshing is used (x-y-z meshing type).

One constraint imposed by the coupling is only to have the same axial meshing in the core domain for the structure and the fluid-subchannels. A proper use of symmetries reduces the size of the numerical system to solve and thus saves the CPU time. All these interactions must be defined in the input data deck.

The thermo-mechanical code: ICARE3D

Mechanical module

During LOCA, the drop of pressure inside the primary circuit associated with the heat up of fuel rods lead to the creep and burst of the fuel cladding. In the meantime, the deformation of fuel cladding blocks progressively the flow channel which impacts the coolability of fuel bundles. The contact between neighbouring rods can be observed and results in complex configurations. The mechanical module deals with strains due to thermal expansion and creep. PCT as well as the burst time and location is predicted by the code.

Basic assumptions

The structure deformation is modelled in quasi-3D at each axial level with a simplified 3D model (model assumes that the axial slices do not interact with each other and the shear or the bending stresses are not taken into account). The stress tensor is assumed diagonal:

σ =                                     σ r = - 1 2 • (P rod + P coolant ) σ θ = (P rod -P coolant ) • r e + σ r σ z = S int • P rod -S ext • P coolant S ext -S int                    (3) with r = 1 2 • (r ext + r int ) e = r ext -r int S ... = π • r 2 ...
With P coolant denoting the coolant pressure, P rod the internal pressure of the rod and r the cladding mid-radius. To deal with non-circular cases, for instance due to the azimuthal gradients or the contact with neighbor rods, the curvature radius is used instead of the simple radius. The internal pressure is calculated with the ideal gas law at each elevation. The free volumes used to calculate the pressure are the gap between the fuel and the cladding, the fuel porosity and plenums. Currently, the fission gas release in fuel pellets and gas flow within the gap are not taken into account in the internal pressure calculation.

To obtain the creep rate, the effective stress σ e can be calculated with the Hill criteria to carry on the possible zirconium anisotropy into the creep rate evaluation.

σ e = (F(σ θ -σ z ) 2 + G(σ z -σ r ) 2 + H(σ r -σ θ ) 2 ) (4)
Only the secondary creep is taken into account and modelled by default by an empirical Norton-like law to relate the creep rate to the effective stress:

d e dt = Aσ e n exp(-Q/RT ) (5)
Classically, the A, n and Q coefficients were obtained from steady and ramp tests (pressure and/or temperature). Those parameters depend on the crystallographic phase of the zirconium alloy (α, α + β, β). Near phase transitions, such as α to α + β and α + β to β, a specific approach is used. It consists in considering two components that behave as parallel systems with the same creep rate (α//α + β or α + β//β). The creep rates in each direction are then calculated from the effective creep rate (from (5)) as follows:

                 ˙ cr-r ˙ cr -θ ˙ cr-z                  = ˙ e ∂σ e ∂σ (6)
with

∂σ e ∂σ = 1 σ e •                  H • (σ r -σ θ ) + G • (σ r -σ z ) F • (σ θ -σ z ) + H • (σ θ -σ r ) G • (σ z -σ r ) + F • (σ z -σ θ )                  (7)
A special algorithm has been developped to adapt the mechanical time step to the creep rate. Indeed, it is important to reduce the time step to predict burst when creep rate is high or when the contact with surrounding structures is close.

Burst criteria

Several burst criteria can be defined with respect to:

-maximum strain, -maximum stress, -cumulative damage on strain.

The strain criteria are derived from two experimental databases described hereafter. The first one is the well-known NUREG-0630 database [START_REF] Owers | Cladding swelling and rupture models for LOCA analysis[END_REF] This criterion gives the maximum strain at the burst elevation depending of the heating rate (0 K/s, above 25 K/s with a linear interpolation in between). These two criteria are also available in their integral form (damage concept [START_REF] Raff | Entwicklung eines Deformations und Versagenmodells für Zircaloy in Hochtemperaturbereich zur Anwendung bei Kühlmittelverluststorfalluntersuchungen an Leichtwasserreaktoren[END_REF])) which is particularly well suited for transient scenarios since the strain criteria depend on the temperature. The burst is assumed to occur when the relative cumulated strain reaches unity. A burst criterion on stress is also derived from the EDGAR program, it gives the burst stress versus temperature as follows:

σ θburst ( f α , T ) = (k exp(-qT )) f α (k exp(-q T )) 1-f α (9)
The stress at burst is evaluated above and below the tear and not at the maximum circumferential strain elevation. In DRAC-CAR, the criteria can be activated at each node, or at a given elevation or simultaneously at two elevations separated by a length (chosen by the user).

Contact between structures

When the model is activated, a search for contact is performed inside all meshes and between all the structures.

The contact may happen:

-initially when the user defines structures in contact, -due to the thermal expansion, -due to the oxidation, -due to the creep.

To detect the contact between structures, the external contour of structures is first discretized. When the external surface is defined with a specific shape (structure 3 with center O3 in figure 4), this shape is used for the discretization. When the external surface is a cylinder and not managed by the mechanical creep model, a regular discretization of the surface is performed (structure 1 with center O1 in figure 4). When the external structure is managed by the creep model, the discretization of the external surface is the one computed by the creep model (structure 2 with center O2 in figure 4). Inside each mesh, a contact fraction is computed as depicted on figure 4 by the f1, f2 and f3 fractions.

For structures in contact, some physical models are then automatically adapted:

-the wall to fluid heat transfers as well as the oxidation of walls by steam and/or air are limited to free surfaces, -the radiative heat transfer between structures is reduced to free surfaces and conductive heat transfer is applied in contact areas,

-the mechanical creep model stops the deformation of wall in contact zones.

Fuel relocation management

The axial fuel relocation in the cladding balloon can be represented using a simple parametric model. Thus, the user has to provide the main characteristics of the fuel fragments able to move: balloon filling ratio, fraction of lost fuel through the break, residual gap thickness, onset of fuel relocation (strain threshold or burst). The model is represented on figure 5.

The physical phenomenon is modeled as depicted on figure 5: -phase 1: radial expansion of fuel fragments, -phase 2: axial relocation of fuel fragments.

The first phase (optional) allows the user to impose an higher porosity in the fuel to simulate a possible radial relocation of fuel fragments. The second phase consists in axially relocating fuel fragments, the model works as follows:

-At cladding burst or on strain threshold, the axial relocation is triggered, -The fuel is then moved from the stacks above to fill the free space consistent with the filling ratio chosen by the user.

The filling ratio is only applied to relocated fuel, therefore the filling ratio varies axially.

Relocated fuel can be regarded as a macroscopic porous media composed by the fuel fragments and the gas. Its equivalent thermal conductivity is evaluated accounting for the conduction in the solid, in the gas and the radiation between fuel fragments [START_REF] Repetto | Application of the porous medium heat transfer model of ICARE/CATHARE code against debris bed and bundle experiments[END_REF]; [START_REF] Yagi | Studies on effective thermal conductivities in packed beds[END_REF]; [START_REF] Imura | Effect of gas pressure on the effective thermal conductivity of packed beds[END_REF]). That equivalent conductivity can be expressed as :

λ = ψ(λ g + α r .D) + 1 -ψ ω λ g + α r .D.ω + 1 -ω λ s (10)
with:

ω = 0.3[(1 -FR)] 1.6 λ s λ g -0.044 ψ = 1 -ω -FR 1 -ω and α r = 4 σT 3
With α r the radiative heat transfer coefficient ( emissivity, σ Stefan-Boltzmann constant), FR the filling ratio, λ g (W/mK) gas conductivity, λ s (W/mK) solid conductivity, D(m) particle mean diameter (considered spherical). The conductivity of fuel fragments is calculated using the approach proposed by [START_REF] Lucuta | A pragmatic approach to modelling thermal conductivity of irradiated UO2 fuel -review and recommendations[END_REF].

Thermal module

The thermo-mechanical code ICARE3D calculates the heat exchange by conduction inside a structure and between structures in contact, through the gaseous gap, between structures by radiation and between structures and fluid by convection. Several options in the code are available to define a power source in a structure including decay heat evolution associated to the FP inventory (for reactor analysis) or electrical heating (for experiment simulation).

Conduction inside a structure mesh

Assuming a steady-state regime for heat diffusion inside a structure mesh, after some simplifications, the solution of the heat balance equation leads to the following equation giving the conductive flux inside a structure mesh:

T s = T av + R s • Φ s (11) with:
-T s surface or mesh face temperature [K],

-T av average temperature of the structure mesh [K],

-R s thermal resistance of the structure depending on con-

ductivity [K/W],
-Φ s total power exchanged by the structure mesh through the face [W].

Conduction between structure meshes

Heat transfer by conduction between a structure mesh and another is done through faces. Currently a structure mesh has a maximum of six faces with the three-dimensional modeling of a solid component (axial, radial and azimuthal). The conductive power exchanged between two structure meshes (figure 6) is equal to:

Φ 1→2 = T 2 av -T 1 av R 1 s + R cont + R 2 s ( 12
)
with:

-R cont contact thermal resistance (possibly defined by the user) [K/W], -R k s thermal resistance associated with the structure mesh k [K/W].

Heat exchange between structures meshes through a gap

In a fuel rod modeling with DRACCAR, a gap can be defined between the fuel pellet and the cladding. The thermal heat transfer model used between those structures (figure 7) is either a simple conductive heat transfer (when the gap is closed) or both radiative and conductive heat transfers in the gap. In this last case, the convective exchange between the gap gas and the structures is neglected (in fact, this assumption is justified because it is possible to show, from the kinetic theory of gases and gaseous flow, that the magnitude of the Grashof number in the fuel-clad gap is less than 2000 and so the heat exchanges due to free convection are small enough to be neglected). So, the heat exchange between the two structures separated by a gap is equal to:

Φ 1→2 = H gap • (T 2 av -T 1 av ) (13)
where H gap is the global exchange coefficient.

When the gap is collapsed (case of contact between the fuel pellets and the cladding for example), the heat exchange is determined as pure conduction between two meshes of structures (equation 12), the contact thermal resistance being set in the input data deck.

If it is not the case, there are two parallel heat transfers in the gap in series with heat transfer conduction in each struc- So, assuming a steady thermal diffusion, the global exchange coefficient is given by:

1

H gap = R 1 s + R gap + R 2 s ( 14
)
with:

-R k s thermal resistance associated with the structure mesh k

[K/W], -R gap thermal resistance of the gap [K/W],.
The thermal resistance of the gap R gap includes the thermal resistance of the gap (R gas ) defining the heat transfer by conduction and the heat transfer coefficient by radiation (H rad ):

R gap = 1 1 R gas + H rad (15)
The thermal resistance of the gas is given by:

R gas = r 2 • ln p 2 p 1 λ gas • S 2 (16)
with:

- The gas conductivity depends on the gas composition and the average temperature in the gap which is assumed equal to:

r
T gas = T 1 s • S 1 + T 2 s • S 2 S 1 + S 2
The radiative heat transfer coefficient is calculated by using Stefan's law assuming that the gas is transparent and surfaces are grey:

H rad = σ • S 1 S 2 • T 3 s12 (17)
where σ is Stefan's constant, S 1 S 2 is the radiation total exchange area calculated by using the classical formula for enclosed structures and T 3 s12 is equal to T 1

s 2 + T 2 s 2 • T 1 s + T 2 s .

Radiation

A radiation heat transfer model is available in ICARE3D to describe this important mode of energy exchange dominant at high temperatures. An important hypothesis is made: the radiation fraction leaving a fluid mesh is assumed to be negligible.

In fact, a part of the photons emitted by the objects within the considered fluid mesh could leave that mesh either radially (entering a neighbor mesh located at the same elevation) or axially (entering a mesh above or below by direct emission or by reflection). Anyway, this hypothesis is rather well verified in the case of DRACCAR simulations because radiation plays a dominant role only in zones at rather high temperatures (above 1000˚C).

Yet in this zones, rod claddings are largely deformed and close enough together to consider that the areas available for a photon leak are small enough and so photon leak are negligible.

As a conclusion, it's reasonable to consider that radiative heat transfer occurs only between structures located in a same fluid mesh.

Nevertheless it appears that the axial radiation has to be taken into account in case of simulations exhibiting a sharp axial temperature profile, which is generally the case in spent fuel pool severe accidents. So a complementary model has been introduced in ICARE3D to calculate the axial radiative heat transfer.

Radiative heat transfer inside a fluid mesh

The implemented model has been adjusted to deal with the geometry depicted in figure 9. It is based on the net radiation enclosure model which considers surface to surface radiation and the interaction between radiation and a gas filling the enclosure. The radiative heat transfer is determined by equation 18

Φ i→ j = σS i S j (T 4 i -T 4 j ) ( 18
)
where σ is Stefan's constant, S i S j is the radiation total exchange area between structure mesh i and structure mesh j, T i and T j the average temperatures of the structure meshes i and j. The two-dimensionnal Hottel's crossed string method is applied to determine the view factors. This method assumes that the object height is much larger than the distance between the objects. The view factor from a first object A to a second object B (figure 10) is computed as follows: For the transmittance calculation, it is supposed that:

(A1B1 + A2B2) -(A1B2 + A2B1) 2 = A1A2 F A→B (19)
-Radiative properties of the gas in the enclosure are taken into account only for a mixture of water vapor and a noncondensable gas, -The gas is considered as an optically thin grey medium, that is to say that the transmittance of the gas is supposed to be wavelength independent. When water is present, no radiation between structures is assumed to occur.

-The geometric mean beam length between each pair of surfaces in the enclosure is very complicated to compute in details. In ICARE3D, it is simply equal to the distance between the chord middles (figure 10).

So the transmittance of the gas is:

τ AB = e -K p (T g ) P s L AB (20)
with K p (T g ) the Planck mean absorption coefficient, P s the steam pressure and L AB the geometric mean beam length.

Axial radiative heat transfer

The axial radiative flux is often neglected for two main reasons:

-the use of a rather coarse axial meshing pointing out a small ratio of the fluid cross section over the vertical struc- Considering all these assumptions, the heat flux from mesh i to mesh j, with j = i -1 or j = i + 1 is:

φ i→ j = H i j • (T 4 i -T 4 j ) (21)
with H i j given by the following equation:

H i j = σ 1 i S i f i→ j + 1 j S j f j→i - 2 i S i + j S j (22)
where:

-S i and S j are the surfaces of the two structure meshes, i and j are the emissivities of the structures, -σ is Stefan-Boltzmann's constant.

The horizontal view factor of the structure mesh i ( f i→ j is computed assuming that the geometric mean beam length is equal to p/2, p being the pitch between rods). If the height of the structure mesh i is ∆z i then: The flow regime is automatically determined by the code itself according to criteria on structure temperatures and fluid void fraction. Two independent flow regions can be chosen. In the first one, the liquid water is always in contact with the structure. On the contrary, in the second region, there is no more liquid in contact with the structure. The choice between the two regions is done through the comparison between the structure temperature and the minimum film stable temperature.

f i→ j = 1 - 2 π atan(α) - log(α 2 + 1) 2α 2 with α = 2∆z i p (
The radiative heat transfer is not mentioned hereafter. In fact, the radiation heat transfer model between the structure and the fluid (calculated through the equation 18) is activated only when the void fraction is above 0.9999.

In reflooding conditions (typical phase of a LOCA in a reactor), a special treatment is done for the heat fluxes in order to calculate better the heat transfer in the mesh where the quench front is located.

Minimum film boiling temperature

The minimum film boiling temperature allows a choice between a pre-critical heat flux regime (liquid convection or nucleate boiling) and a post-critical heat flux regime (gas convection or film boiling) as it can be seen in figure 12. The minimum film boiling temperature is based on the Groeneveld-Stewart corre-lation [START_REF] Groeneveld | The minimum film boiling temperature for water during film boiling collapse[END_REF]) and is equal to:

           f (P) if P < 9.MPa T sat (P) + ( f (P 0 ) -T sat (P 0 )) • P cr -P P cr -P 0 if P ≥ 9.MPa (24)
where:

f (P) = 557.85 + 4.411 • 10 5 • P -3.72 • 10 -12 • P 2 -10 4 • X eq 2.82 + 1.22 • 10 -6 • P -P 0 = 90 • 10 5 Pa,

-P cr = 221.2 • 10 5 Pa,
-X eq is the subcooled quality.

Structure to fluid heat transfer in liquid convection regime

The heat transfer between structure and liquid phase by convection can only occur if the structure temperature (T w ) is lower than the saturation temperature at a given pressure (T sat (P)).

Otherwise, it is equal to 0. The heat transfer by liquid convection is:

φ CFL = λ l • Nu l D h • (T w -T l ) (25)
In equation 25, the Nusselt number Nu l can be given by the user in the DRACCAR input data deck. Otherwise this number describes a laminar flow regime or a turbulent one. So it is equal to:

Nu l = max(3.66, Nu t f c ) (26)
Nu t f c is given by the Dittus-Boelter correlation:

Nu t f c = 0.023 • Re 0.8 l • Pr 0.4 l ( 27
)
The liquid properties used in these previous equations are taken in the fluid mesh in contact with the structure mesh:

-D h is the hydraulic diameter of the fluid mesh [m],

-T l is the liquid temperature

[K], -λ l is the liquid conductivity [W.m -1 .K -1 ],
-Re l is the liquid Reynolds number, -Pr l is the liquid Prandtl number.

Structure to fluid heat transfer in steam convection regime

Heat transfer between structure and gas phase occurs when the structure temperature (T w ) is greater than the saturation temperature. The structure to steam heat flux depends on the steam flow regime (natural or forced convection) and on the geometry simulated (pipe or bundle). It is given by:

φ CG = λ g D h • Nu g • (T w -T g ) (28)
where Nu g is the Nusselt number which is determined as follows:

Nu g = max(Nu geom , Nu f c1 , Nu f c2 , Nu nc ) (29)
where:

-Nu f c1 is given by the Groeneveld-Delorme correlation [START_REF] Groeneveld | Prediction of thermal non-equilibrium in the post dry-out regime[END_REF]):

Nu f c1 = 0.0083 • Re 0.877 g • Pr 0.611 g -Nu f c2
is given by the Kays-Crawford correlation [START_REF] Rohsenow | Hand book of heat transfer fundamentals -second edition -forced convection[END_REF]):

Nu f c2 = 0.0210 • Re 0.8 g • Pr 0.5 g -Nu nc is given by Rohsenow-Harnett correlation [START_REF] Rohsenow | Hand book of heat transfer[END_REF]):

Nu nc =                  0.135(Gr g • Pr g ) 0.
33 if Gr g • Pr g > 2.10 7 0.54(Gr g • Pr g ) 0.25 if Gr g • Pr g > 5.10 2 1.18(Gr g • Pr g ) 0.125 else -Nu geom depends on the geometry. In case of heat exchange in a pipe, it is equal to 3.66. In case of heat exchange in a bundle, it is given by Sparrow's correlation [START_REF] Sparrow | Heat transfer to laminar flow between cylinders[END_REF]): The gas properties used in these previous equations are taken in the fluid mesh in contact with the structure mesh:

Nu geom =                          30.0 • psd -28.5 if psd < 1.2 2.
-D h is the hydraulic diameter of the fluid mesh [m],

-T g is the gas temperature [K],

λ g is the gas conductivity [W.m -1 .K -1 ],

-Re g is the gas Reynolds number, -Pr g is the gas Prandtl number, -Gr g is the gas Grashof number. Film boiling heat flux used both for the inverted annular/slug and dispersed flow regime is given by the following correlation:

φ bo = H ber • (T w -T sat ) (31)
with H ber the heat transfer coefficient from the original Berenson correlation [START_REF] Berenson | Film boiling heat transfer from an horizontal surface[END_REF]).

Reflooding model

Reflooding is a key point of reactor accident management.

Physically speaking, it is characterized by the injection of large amounts of water into the very hot structures of the core.

The resulting flow pattern is usually split into four zones (fig-

ure 13):

-below the quench front, all the temperatures are under the burn-out temperature and the structures are wet, -just above the quench front, a non-stable two-phase flow establishes at the structure side; there, very large axial temperature gradients can drive the cooling (top-bottom conductive heat transfers),

-further above, a stable gaseous film establishes along the structures, limiting drastically the coolability of the structures and decreasing the liquid-gas interfacial area (inverted annular flow),

-far above the quench front, the gas phase is predominant, possibly entraining droplets and/or water slugs.

CHF = F 0.25 CHF • ρ g • (H v -H l ) • Csub • 0.149 (32)
where:

-

F CHF = g • σ l • ρ l -ρ g ρ 2 g -Csub =              1 if T sat ≤ T l 1 + 0.1 • C p l • (T sat -T l ) H v -H l • ρ l ρ g 0.75 if T sat > T l
ρ k is the density of the fluid phase k (liquid or gas), -H k is the enthalpy of the fluid phase k (liquid or gas), -σ l is the liquid surface tension, -C p l is the heat capacity for liquid.

Downstream, the flux decreases exponentially with respect to the distance from the quench front.

Φ T B = CHF • e ZQF-Z L 1 (33)
The length L 1 is obtained in the Ishii experiments [START_REF] Ishii | Flow visualization study of inverted annular flow of post dry-out heat transfer region[END_REF]). It is proportional to the square root of the capillary number:

L 1 = 5 8 • Ca 0,5 for Re l < 2000 and L 1 = Ca 0,5 for Re l ≥ 2000 where Ca, the capillary number, is expressed as follows:

Ca = µ l • V L σ l

Power generation

Because of the origin of the nuclear power, there are two models in DRACCAR allowing the simulation of experimental conditions as well as accident transient scenarios of nuclear power plants. Moreover, those models can be activated at the same time during a calculation. The total nuclear power is the result of:

-the power delivered by neutrons (and delayed neutrons after reactor scram),

-the decay heat of fission products (alpha, beta and gamma radiations).

In the model for neutron power, the following assumption is made: no relocation of heating material occurs till the end of the delayed neutron power generation.

The nuclear power generated in a structure mesh i is expressed as follows:

P nucl,i = P neut nucl,i + P FP nucl,i (34)
with P neut nucl,i the power delivered by neutrons directly given by the user in the input deck or calculated by a classical neutron point kinetics model. P FP nucl,i is the decay heat delivered by fission products evaluated by the DRACCAR module ISODOP or imposed by the user in the input deck.

Chemical model and DIFFOX module

The fuel rod cladding oxidation is one key phenomenon involved during LOCA especially as the oxygen diffusion within cladding material favours its embrittlement. That risk is associated to a safety criteria on the maximum equivalent cladding reacted (ECR) which imposes a limit of oxygen gained by the Zircaloy cladding. As a consequence to assess that risk, codes

Zircaloy oxidation by steam

The oxidation of Zircaloy (Zry) in steam is a complex process, modeled in a simple way using experimental correlations.

The following schematic equation for the reaction shows that both heat and hydrogen gas are liberated by the reaction:

Zr + 2H 2 O → ZrO 2 + 2H 2 + ∆heat (35)
with: ∆heat=526 kJ/mol-Zr at 25˚C

The reaction of Zry with steam at elevated temperatures involves the growth of discrete layers of oxides and oxygen-rich phases (illustrated on figure 15). Numerous experiments show that under unlimited steam conditions, the oxidation process is an oxygen diffusion-controlled process: the chemical reaction kinetics for the layers growth and the total oxygen mass gain can be therefore modeled using parabolic laws in the following differential form:

dX 2 dt = K p (T ) (36)
with X being either the total oxygen mass gain in both α-Zr(O)

and ZrO 2 layers [kg/m 2 ] or the ZrO 2 layer thickness [m]

Experimentally measured parabolic correlations K are strongly temperature dependent, and are presented as Arrhenius type functions with constant coefficients :

K p = A exp(-B/(R.T )) (37)
where:

-T is the wall temperature [K],

-A is the constant of the Arrhenius formulation in [kg 2 /m 4 /s] or [m 2 /s] and B is the activation energy in

[J/mol],
-R is the ideal gas constant.

In the case of limited steam supply conditions (starvation) Reactions between Zircaloy and steam depends on the kind of zirconium alloy. In ICARE3D, four oxidation laws are available for Zircaloy-4 obtained for different temperature ranges.

and
The recommended correlation (default option) is the BEST-FIT correlation introduced by A. Volchek (NSI NRC KI) in 2004 [START_REF] Volchek | Advanced treatment of zircaloy cladding high-temperature oxidation in severe accident code calculations: Part ii. best-fitted parabolic correlations[END_REF]). This correlation has been evaluated

analytically on the basis of oxygen diffusion theory and of an experimental database.

Zircaloy oxidation by air

The presence of air -in a SFP accident for instance -can lead to vigorous oxidation of the cladding resulting in an additional heat generation. In comparison with steam, Zr oxidation by air yields about twice more heat, driving the degradation further. These considerations have encouraged IRSN to launch the MOZART research program [START_REF] Duriez | Zircaloy-4 and M5 R high temperature oxidation and nitriding in air[END_REF]) consisting of experimental investigations and development of an accurate modeling of air oxidation [START_REF] Coindreau | Air oxidation of Zircaloy-4 in the 600-1000˚C temperature range: Modeling for ASTEC code application[END_REF]). The schematic equation for the reaction is the following: 38) with: ∆heat=1100 kJ/mol-Zr at 25˚C

Zr + O 2 → ZrO 2 + ∆heat (
The main difference with the previous steam oxidation modeling is that when Zircaloy-4 is scaled in air, two different kinetic regimes unfold, separated by a kinetic transition named breakaway. Only the total oxygen weight gain is quantified in these available data. Accordingly, it is assumed that 90% of the weight gain correspond to the oxide layer growth (consistent with some recommendations made by OPSA partners [START_REF] Shepherd | Oxidation phenomena in severe accidents (opsa)final report[END_REF])).

Pre A = 45.9(kg/m 2 ) 0.5 s -1 and B = 1.32510 5 J/mol

This correlation is used as long as the temperature is lower than 1200˚C. For higher temperatures, a parabolic correlation is reapplied. -Nitriding :

ZrO x + y 2 N 2 → ZrO x N y -Nitriding oxidation by O 2 : ZrO x N y + (1 -x 2 )O 2 → ZrO 2 + y 2 N 2
-Nitriding oxidation by steam :

ZrO x N y + (2 -x)H 2 O → ZrO 2 + y 2 N 2 + (2 -x)H 2
The model is rather simple since it's assumed that :

-the kinetics of these reactions are modeled by linear laws (using Arrhenius constant coefficients),

-the three reactions are not simultaneous and happen in so far as the chemical reagents are present. In that case, gaseous pressure thresholds defined as input data (parametric model) are used to trigger the reactions, -the fluid involved in the chemical reactions is evaluated by the thermal-hydraulic code (CESAR or CATHARE-3).

Oxygen diffusion: DIFFOX module

In order to calculate oxygen diffusion through the cladding The oxygen concentration profile is calculated in each layer by solving the 1-D diffusion equation (Fick's second law) using the implicit finite difference method (FDM) in cylindrical coordinates:

∂C(r, t) ∂t = - 1 r • ∂ ∂r (r • F r ) (40)
where the left side is the oxygen accumulation and the right side is the space change of the diffusion flux F r . The oxygen flux is assumed to be proportional to the oxygen gradient (first Fick's law):

F r = -D(T ) • ∂C(r, t) ∂r [kg/m 2 .s] (41)
where C(r,t) is the oxygen concentration [kg/m 3 ], D(T) is the oxygen diffusion coefficient as function of temperature [m 2 /s],

t is the current time [s], r is the radius [m].
The layer interface movement at the end of the time step is determined from the condition to preserve both the oxygen and the zirconium mass in the control volume at the end of the time step before and after layer interface movement.

The simplifying assumptions are listed below:

-The cladding tube is a one-dimensional system with finite geometry.

-At the beginning of the oxidation transient the initial temperature must be in the α-temperature range, i.e. below 800˚C. The system consists of the α -Zr phase possibly covered by the corrosion layer on the outer (or on both the outer and the inner) cladding surface. The subsequent creation or disappearance of reaction layers due to temperature change (or due to layer consumption) is performed automatically.

-The oxidation is controlled by oxygen diffusion processes (Fick's first law) combined with α to β phase transformation during heating and cooling.

-The mass balance at the layer interfaces (the Stefan equation) and in the volume of the reaction layers is performed for oxygen only.

-All the hydrogen initially present is supposed to move into the β -Zr layer when it appears (hydrogen pick-up is not considered).

-Oxygen concentrations at the layer interfaces are taken from Zr-O and Zircaloy-O equilibrium phase diagrams.

The influence of hydrogen on the equilibrium concentrations can be taken into account. To avoid an unphysical state under temperature transients the change of the equilibrium values with temperature is limited in some cases.

-Oxygen diffusion coefficients in ZrO 2 , α -Zr and β -Zr depend only on temperature and not on concentration, following an Arrhenius law :

D = D 0 exp( -Q RT ).
-The equivalent oxygen diffusion coefficient in the twophase layer (α + β) -Zr is built from the individual diffusion coefficients for each phase and the corresponding volume fractions. This assumption was adopted with the aim to keep the model as simple as possible.

-The volume expansion related to the metal-to-oxide transformation (the Pilling-Bedworth ratio) is 1.5. The volume change related to the Zr to ZrO 2 phase transformation is neglected.

-The movement of the layer interfaces is calculated before solving the diffusion equation (the second Fick's law) in each reaction layer.

The oxide layer in models for Zircaloy oxidation in LOCA conditions is usually approximated by a homogeneous mate-rial. This is acceptable for as-received Zircaloy and for temperatures well above 1000˚C, where the high-temperature oxide is tetragonal and the oxidation kinetics is parabolic. The presence of both the tetragonal sub-layer at the oxide/metal interface below 1000-1100˚C (where the monoclinic phase is stable) and the monoclinic pre-oxides at the cladding surface(s) complicate the situation, because the material properties (diffusion coefficient) of the monoclinic and tetragonal phases differ from each other. Thus, the oxide model in DIFFOX distinguishes up to three zones (sub-layers) in the oxide layer: tetragonal sub-layer, bulk oxide and pre-oxide. Closure laws of such a model is under development and some validation has been performed.

Balance equations of a 6-equation model

The equations of the thermalhydraulic codes are the timeaveraged equations of the two-fluid model [START_REF] Ishii | Thermo-fluid dynamic theory of two-phase flow[END_REF]). CE-SAR and CATHARE-3 use the same form for the mass and momentum balance equations. The mass balance equations are given by:

∂α k ρ k ∂t + ∇ • (α k ρ k -→ V k ) = Γ k (42)
where Γ k is the volumetric production rate of phase k due to phase change depending in particular on convective heat flux calculated by ICARE3D.

The momentum balance equations are:

α k ρ k         ∂ -→ V k ∂t + -→ V k • ∇ -→ V k         = -α k ∇P + (P i -P)∇α k + (-1) k-→ τ i + --→ τ kw + α k ρ k - → g (43)
where P, P i , -→ τ i and --→ τ kw are respectively the pressure common to the two phases, the interfacial pressure, the interfacial friction per unit volume and the wall friction per unit volume. Moreover, a simplification has been done in the CESAR momentum balance equations: all the crossed contributions ( ∂V x ∂y , ∂V x ∂z , ...) in the momentum convective terms are not taken into account. It is not the case in CATHARE-3.

For the energy balance equations, CESAR uses the enthalpy form (equation 44). CATHARE-3 uses the internal energy form (equation 45) for the vessel modeling and the enthalpy form for the circuits modeling.

∂ ∂t

      α k ρ k       H k + V 2 k 2             + ∇(α k ρ k -→ V k (H k + V 2 k /2)) = (-1) k Γ k       H k + W 2 i 2       + q ki + q kw -α k ∂P ∂t + α k ρ k -→ V k • - → g (44)
where q ki and q kw denote respectively the volumetric energy exchange between phase k and interface and the one between phase k and walls (which is calculated by ICARE3D)

∂α k ρ k e k ∂t + ∇(α k -→ V k (ρ k e k + P)) = Γ k h ki + q ki + q kw -P ∂α k ∂t (45)
where h ki is the mean enthalpy of phase k near the interface.

It can be noticed that the molecular and turbulent diffusion terms are not taken into account in the momentum and enthalpy/energy equations. Otherwise the presence of structures in the hydraulic meshes is taken into account by introducing volume and surface porosities when integrating the system of equations in space according to the finite volume approach.

Numerics and coupling scheme

An implicit coupling between the structure and the fluid codes is needed to prevent artificial gain or loss of energy related to the convective exchanges, as well as gain or loss of oxygen mass in case of starvation condition variations. Independently from any coupling, each code or module may use, internally, its own implicit technique.

In a standard coupling, the thermo-mechanical code and the thermal-hydraulic one have the same time step in order to calculate the evolution of structure and fluid unknown variables.

But considering the consequent number of meshes (3D mesh-

dH dt = F(T p , t) -MR • E (46) with:
-H: the enthalpy of the structure, -F: all the energy exchange inside and between structures.

dT dt = 1 ∆t • [δT 0 + A • E] (47)
with:

-A the inverse matrix of the system evaluated for the converged temperature field T * , -δT 0 the structure temperature variation after a time ∆t when the external flux are exactly equal to E 0 .

Considering this fact, a predictive/corrective method is applied to perform the coupling:

-the behavior of the structure is predicted by the thermo-

Validation and Application

The description of the validation status of DRACCAR is out of the scope from this present article. Precise information on this topic can be found in [START_REF] Bascou | Development and validation of the multi-physics DRACCAR code[END_REF]) and will be extensively presented in a dedicated up-coming article.

At this moment, the code DRACCAR is intensively used to define precisely the appropriate test conditions in several experimental programs. For example, in the framework of the PERFROI project (Repetto et al. (2015a)), DRACCAR has been used to define some parameters of the thermalhydraulic COAL experiment (dedicated to the study of ballooned bundles reflooding study) such as the blocked bundle geometries, the heat-up scenario or the test matrix. It is also used to perform pre-calculation of all the tests foreseen in this experiment.

The code is also used to perform LOCA calculations on existing or future nuclear power plant in support of safety demontration analysis.

In the framework of the European project, NUclear REactor SAFEty simulation platform (NURESAFE SP3) project, apart from IRSN, AREVA [START_REF] Fargette | Simulation of REBEKA-6 with DRACCAR v2.1[END_REF]) and ENEA [START_REF] Bascou | Development and validation of the multi-physics DRACCAR code[END_REF]) used the code to perform some validation work

Perspectives

In a near future, the modeling efforts will be enhanced to better simulate the mechanical contact between rods and neighboring structures. One of the main ideas is to address the issue of the coolability, as a safety criterion, of a partially blocked bundle resulting from the dynamic ballooning of rod claddings and sub-channel blockage during a loss of coolant accident. In this context, the experimental program called PERFROI (Repetto et al. (2015b)) (in particular, the COCAGNE and ELFE experiments), supported by Agence National pour la Recherche (ANR) and Electricité de France (EDF), is underway in order to help to validate and improve the thermo-mechanics models of ICARE 3D dealing with cladding contact, creep and burst.

In the same program, the COAL experiment will be used to improve structure to fluid heat transfer models (ICARE3D) and

  codes and modules of the DRACCAR platform are presented on figure 1. Two different thermalhydraulic codes (CESAR and CATHARE-3) can currently be coupled to the thermomechanical code ICARE3D. CESAR and ICARE3D are two codes developed by IRSN. CATHARE-3 is the advanced

Figure 1 :

 1 Figure 1: Codes and modules available in DRACCAR platform

Figure 2 :

 2 Figure 2: Example of a structure meshing

Figure

  Figure 4: Contact detection between structures

Figure 5 :

 5 Figure 5: Principle of the radial and axial relocation model

Figure 6 :

 6 Figure 6: Temperature profile through 2 structure meshes

Figure 7 :

 7 Figure 7: Temperature profile through meshes of 2 structures separated by a gap

Figure 8 :

 8 Figure 8: Thermal resistances linking two structures meshes separated by a gas filled gap

  2 the external radius of the gap [m]; p 1 and p 2 the internal and external perimeter of the gap [m]; -S 2 the external heat exchange area of the gap [m 2 ]; -λ gas the conductivity of the gas [W/m/K].

Figure 9 :

 9 Figure 9: Geometry on which the radiative heat transfer model is applied

Figure 10 :

 10 Figure 10: View factor determination

  ture surface: the view factor of a vertical structure in the horizontal slice if much larger than its view factors upward and downward, especially because structures are close together, -a smooth axial temperature gradient. But if the axial discretization is quite fine (typically close to the hydraulic diameter), and if the temperature gradient is steep, then the axial radiative heat transfer can become significant and even dominent. In the scope of DRACCAR, this situation is only met during the late phase of spent fuel pool uncovery when the assemblies are oxidized by air. In order to well reproduce these situations, a rather simple radiative axial model has been chosen. It is based on an axial conduction evaluating coarsely the axial radiative heat flux. The main assumptions are: -in a fluid channel, the difference of temperatures between two wet structures in two neighboring axial fluid meshes is greater than the difference of two wet structures belonging to the same fluid mesh, -the heat transfer is performed locally, to the meshes just above and just below, without trying to propagate the radiative flux to the next meshes.The first assumption is used to simplify strongly the heat transfer: instead of computing a heat transfer for all the possible couple of wet structures (one in each considered neighboring meshes), a wet structure mesh will exchange energy only with the wet structure meshes just above and just below which belong to the same structure (figure11). That assumption is valid if the equivalent weighted rod represents a large enough number of rods.

Figure 11 :

 11 Figure 11: Structure meshes involved in the axial radiative heat exchange model

  23)3.2.5. Convective heat transferThis mode of heat transfer is the major one during a LOCA that drives the coolability of the core. Particular attention was paid to its modelling in DRACCAR including the reflooding configuration. Structure to fluid heat transfer models are necessary both in the thermo-mechanical code and thermalhydraulic code coupled in DRACCAR. In the thermo-mechanical code ICARE3D, these models are used in the energy balance equation (46). In the thermalhydraulic code, these models are closure laws of the mass and energy fluid balances. In order to well simulate the core behavior during a reactor or spent fuel pool LOCA, it is needed to describe all structure heat transfers of each flow regime met in such accidents. As a consequence, all the models implemented in ICARE3D allow the description of the boiling curve (figure12). The flow regimes which can be simulated are:

Figure 12 :

 12 Figure 12: Structure to fluid heat flux

  3.2.5.4. Structure to fluid heat transfer in nucleate boiling regime Nucleate boiling heat transfer occurs only if the structure temperature is greater than the saturation temperature (otherwise it is set to 0). It is based on Thom's correlation (Thom et al. (1965)): φ NB = 1.97 • 10 3 • exp(0.23 • 10 -6 • P) • (T w -T sat ) 2 (30) 3.2.5.5. Structure to fluid heat transfer in film boiling regime

Figure 13 :

 13 Figure 13: Flow regimes during reflooding and associated flux

Figure 14 :

 14 Figure 14: Heat flux imposed in the meshes closed to the quench front during reflooding

simulating

  LOCAs require dedicated models to evaluate the development of the oxidation. Moreover, at high temperature, the cladding oxidation is one of the key phenomena influencing the core behavior. It is well established that this chemical reaction results in a very large exothermal heat release, determines the hydrogen source term and can cause fuel rod degradation. The chemical model of ICARE3D deals with Zircaloy oxidation by steam and by air at low and high temperatures. Even if high temperature scenarios are not strictly in the scope of DRACCAR code, oxidation and accelerated kinetics are available and useful for SFP accident simulation or in the demonstration of coolability of the core in LOCA conditions. In addition to study the cladding embrittlement, a model of oxygen diffusion in zirconium has been developed to determine the oxygen profile within the oxidized cladding (DIFFOX module).

Figure 15 :

 15 Figure 15: Schematic view of the radial oxygen concentration profile in the three layer at the beginning of oxidation

  when the 3 layers are present, an additional assumption is needed. It is assumed that the ZrO 2 layer can grow if the oxygen available allows at least the increase of the α-Zr(O) layer as in non-starvation conditions (assumption partially verified using diffusion equation codes for low starvation rate). Double-sided oxidation after cladding burst may be considered in the model. The interaction surface is updated in case of ballooning or contact between structures.

=

  -breakaway kinetics: The model calculates the interaction between Zr and O 2 with parabolic kinetics. Many parabolic kinetic relations are available in the code (correlation from Hayes and Roberson recommended by OPSA partners and various NUREG correlations[START_REF] Basu | Air oxidation kinetics for Zr-based alloys[END_REF])). However, the default correlation in ICARE3D was built from the MOZART program, with coefficients in eq. 37 for the oxygen mass gain: A = 696.75(kg/m 2 ) 2 .s -1 and B = 1.94910 5 J/mol Remark: Those coefficients are equivalent to the one published[START_REF] Coindreau | Air oxidation of Zircaloy-4 in the 600-1000˚C temperature range: Modeling for ASTEC code application[END_REF]) and provide the same simulation results.Breakaway transition: breakaway transitions were observed in the MOZART experiments for samples oxidized in isothermal conditions below 1200˚C. These transitions occurred when a critical oxygen weight gain, depending on the temperature, was reached. It was assumed that the breakaway transition is linked to the "tetragonal to monoclinic" zirconia transformation, the latter phase being thermodynamically stable below temperature Tb=1174˚C. Based on that assumption, the weight gain for Zircaloy-4 is correlated with the temperature by an hyperbolic law(Coindreau et al. (2010)) for T<Tb. Post-Breakaway kinetics: After the kinetics transition, the oxidation rate increases due to breakaway propagation and it is modeled by an accelerated law as it was observed for sample oxidized in air at moderate temperatures (<1200˚C) in the framework of the MOZART program. The transition to accelerated kinetics can be simply modeled by the following equation: Ka(T ) = A exp(-B/(RT )) (39) with X being the total oxygen weight gain in both α-Zr(O) and ZrO 2 layers [kg/m 2 ]. The weight gain correlation for the accelerated constant was established through the MOZART program. For Zircaloy-4, it was shown that:

3. 3

 3 .3. Zircaloy nitriding and nitride re-oxidation From 2009 to 2013, the US-NRC and 12 partners from the OECD/NEA undertook an experimental program, called the Sandia Fuel Project, to address thermalhydraulic conditions and zirconium fire propagation during a complete loss of coolant in a spent fuel pool from a pressurized-water reactor (Durbin et al. (2016a,b)). From that project, it was shown that nitriding reactions occurred above the thermal front, for cladding above 1200˚C and in a zone free of oxygen (starvation). From the experimental studies performed at KIT on high-temperature reaction of α-Zr(O) with nitrogen (Steinbrück (2014)), an empirical model was built in ICARE3D to reproduce these experiments. The possible re-oxidation of the nitrides by oxygen, or by steam, is also evaluated by ICARE3D. The nitriding model aims to represent the following chemical equations (the values x=0.25 and y =1. are taken in the code):

  the DIFFOX module can be activated locally in a region of interest. This module numerically solves a typical non-linear multilayer Stefan problem with moving boundaries within the Zircaloy cladding tube. Up to seven reaction layers can be simulated in the cladding wall: [InnerZrO 2 ; α -Zr ; (α + β) -Zr ; β -Zr ; (α + β) -Zr ; α-Zr ; outerZrO 2 ]. During oxidation or temperature transients, some layers may be consumed by other layers or dissolve and some layers may form again.

  4. Thermal-hydraulic codes: CESAR and CATHARE 3During LOCA, complex flow dynamics occurs with phenomena involving fast kinetics such as the blowdown phase and complex phenomenology such as reflooding. Safety requirements need demonstration of the coolability of the core even in deformed geometries with blockage of the core channels. To achieve that objective, DRACCAR requires modelling of the two-phase flow within the core channel by taking into account the complex geometrical configuration generated by rod deformation. When dealing with the two-phase fluid dynamics simulation, the choice of the appropriate physical model is crucial.The main types of two-phase flow models are the following:-Models treating the two-phase mixture as a whole (one fluid model):• homogeneous-equilibrium mixture (HEM) model is a 3-equation model where one assumes that the velocity, temperature and pressure between the phases are equal. The two-phase mixture is treated as a pseudo single-phase fluid supposing that the phases are at thermodynamic equilibrium and that they move with the same velocity. Thus, one mean mass balance equation, one mean energy balance equation and one mean momentum balance equation are used.• drift-flux model (DFM) is a 4-equation model. The two-phase mixture is considered as a whole, whereas each phase has its own velocity. It relates the gasliquid velocity difference to the drift flux (or socalled drift velocity).-Two-fluid model: • 5-equation model involves two mass balance equations, two energy balance equations and one mean momentum balance equation. Thermal nonequilibrium can thus be simulated. A phase slip model, which is closely related to the interfacial friction, is added in order to model the mechanical nonequilibrium between the liquid and the gas phases. • 6-equation model is a more general model since it includes two mass balance equations, two energy balance equations and two momentum balance equations for the liquid and the gas phases. The adapted selection of the physical model depends basically on the balance between the required accuracy, the accepted simulation time, and the character of the physical problem at hand. In nuclear applications, the HEM model is suitable for instance for the computation of pressure losses in the core channels in nominal state, whereas for the calculation of void fraction, the DFM is necessary to use since the relative velocity between the liquid and the gas has to be taken into account. In the case of fast pressure transient conditions, the two-fluid models are adequate. The advantage of a 5-equation model over a 6-equation model consists of simplifying the momentum equation since there is no interfacial friction involved and this leads to a higher computing efficiency because of the suppression of one momentum balance equation. The drawback is the relative loss of representativity of the velocity field particularly when countercurrent flows of vapor and liquid occur. Multiple liquid field models distinguishing the dispersed liquid from the continuum can also be mentioned as an evolution of the 6-equation model. Regarding the LOCA, such a model could describe the transport and behavior of liquid droplets during the reflooding phase. CESAR and CATHARE 3, the two codes included in the DRACCAR platform, are both based on a two-fluid model. In CESAR, it is possible to use a 5-equation model or a 6-equation model depending on the type of simulation and accuracy desired. In CATHARE 3, a classical 6-equation model is available and a multi-field model is under development but this last model cannot yet be activated in the DRACCAR platform. 4.1. Presentation of the CESAR and the CATHARE-3 codes CESAR is the fluid module of the ASTEC code developped by IRSN (Chatelard et al. (2014). It is designed to simulate the two-phase thermalhydraulics of reactor coolant system during the different phases of reactor accident. CESAR was initially based on a 5-equation model completed by a phase slip correlation. It provides a good compromise between the need of the representativity and resolution time in particular in case of simulation of severe accidents, framework of the ASTEC code. With the development of the DRACCAR code which is dedicated, in particular to simulate precisely the core behavior during a loss of coolant accident, it has been decided to use CESAR to calculate the thermalhydraulic. However, it appears that the 5-equation model cannot well represent the fluid behavior in the core. Then a 6-equation model has been implemented in CE-SAR. So, it is possible to use a 5-equation model or a 6-equation model in the CESAR code through selection in the input data deck. The Code for Analysis of THermalhydraulics during an Accident of Reactor and safety Evaluation (CATHARE) is a system code for PWR safety analysis, accident management, definition of plant operating procedures and for research and development. It is also used to quantify conservative analysis margins and for licensing. CATHARE is the result of a joint effort of AREVA, CEA, EDF and IRSN. CATHARE-3 (Emonot et al. (2011)) is in continuity with the CATHARE 2 code which is the current industrial version of CATHARE and internationally used for nuclear power plant safety analysis. CATHARE-3 uses the classical two-fluid 6-equation model which has been extensively validated. One limitation of the 6-equation model is the use of a flow regime map based on simple local algebraic criteria which do not describe the dynamics of transitions. Moreover when aphase is split into separate fields (e.g. continuous liquid plus droplets) the dynamics of the two fields cannot be well represented by one field. That is why, an approach to overcome these difficulties is to use a multi-fields model. For CATHARE-3, a three-field model consisting of a continuous liquid plus a droplet field is considered as the most promising way of improvement. Such a model is already available in CATHARE-3.

  ing) used in DRACCAR simulation (e.g. for modeling only 1/8 of an assembly, nearly 100000 structure meshes and 10000 fluid meshes are frequently used), the time duration of such calculation would be too long principally because the thermalhydraulic time step have to respect the Courant/Friedrich/Levy criteria generally imposed by the numerical scheme of the thermal-hydraulic code. As a consequence, an advanced coupling technique between ICARE3D and any thermal-hydraulic code has been developed in order to overcome this difficulty. In fact, this technique allows each code to use its own normal time step: a relatively large one for ICARE 3D (called macro time step) and a small one for thermal-hydraulic codes (called time step). Moreover, it has the advantage to reduce the frequency of the data exchange between codes.The coupling is based on the fact that the behavior of the structures, in the absence of fluid, is usually almost linear unless phase changes occur (in fact in case of melt of material, which is out of the scope of DRACCAR simulations). Considering the general energy balance equation, which described all the physical phenomena treated by the thermo-mechanical code, only one leads to introduce a potentially non-linear term in the structure energy balance equation, the convective flux itself :

  This term includes the energy exchange by conduction and radiation in one hand, and the energy transfer due to the structure changes (chemical reactions, relocation) in another hand, -T p : the averaged structure temperature, -E: the heat exchange between the structure and the fluid, -MR: a matrix that distinguish wet structure elements among the structure meshes -only these structure elements can exchange with the fluid by convection. So the key idea of the coupling is to eliminate the explicit use of the convective flux and its partial derivatives in the structure energy balance equations. It is done by considering that this term is completely unknown. Moreover, if a reasonable estimation of E is known (named E 0 ), it is possible to solve the previous equation 46, the solution being denoted T * . As a consequence, it is possible to express the structure temperature as a function of the convective flux term during a structure macro time step:

  mechanical during a macro time step using the convective flux of the previous time step (E 0 ). During the prediction step, the expression of the structure temperature is determined by DRACCAR as a function of the convective flux, its value at the previous time step and the estimated sturcture temperatures (T * ). -then, DRACCAR calculates at each thermalhydraulic time step the structure temperature using the expression 47 and deduces the energy exchanged by convection during a time step. So the evolutions of the fluid variables can be calculated by the thermalhydraulic code. The mean convective flux during the macro-time step is determined too.-in the correction step phase, DRACCAR calculates the structure temperature evolution using the mean convective flux, -Then, the good agreement between the prediction step and the correction one is checked. This checking consists in the comparison between the structure temperature evaluated during the prediction step by the linear expression and the structure temperature calculated at the end of the correction step. If the verification test is not conclusive, the predictive step is started again using the mean convective flux obtained in the previous aborded step.

Figure 16 :

 16 Figure 16: Scheme of the fluid/structure algorithm

  3D thermal-hydraulics (CATHARE-3 and CESAR) in reflooding conditions. Otherwise, several developments are planned for the DRAC-CAR code including: -a better definition of the fuel rod initial state (oxidation and hydrogen concentration in cladding, fuel burn up, gas inventory and fuel micro-structure) by linking DRACCAR to a performance fuel code like FRAPCON (Geelhood et al. (2015)), -an adaptation of the DRACCAR models in order to deal with multi-scale problems allowing to simulate an entire reactor core in LOCA conditions in an acceptable computation time, -a generalization of parametric and sensitivitiy studies already performed with DRACCAR in order to evaluate the model uncertainties, -an increase of the code performance (CPU time, portability) through numerical developments.

Acknowledgements The work described in the present paper has been performed in the continuity of many years of development of the DRACCAR code at IRSN with the support of EDF.

The authors are thankful to the former or current colleagues from IRSN for their fruitful contribution to the physical modeling. In particular, thanks are due to Sophie Bascou, Tim Haste, Franc ¸ois Lamare and Jean-Marc Ricaud. Some valuable contributions to the code improvement have been provided by Stefano Ederli from ENEA as well as by members from CS-SI company involved both in the code maintenance and development.