Thierry Beck
email: thierry.beck@cea.fr

Victor Blanc

Jean-Michel Escleine

David Haubensack

Michel Pelletier

Mayeul Phelip

Benoît Perrin

Christophe Venard

Conceptual design of ASTRID fuel sub-assemblies

Keywords: ASTRID, SFR, fuel, sub-assembly

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The French 600 MWe sodium-cooled fast reactor demonstrator called ASTRID has just reached the end of its Conceptual Design phase (AVP2). The aim of the past three years (2013)(2014)(2015) of studies was to propose a consistent, innovative preliminary design [START_REF] Rouault | ASTRID, the SFR GEN IV Technology Demonstrator Project: Where Are We, Where Do We Stand For?[END_REF]. The core design studies are being conducted by the CEA with support from AREVA and EDF. Innovative design choices for the core have been made to comply with the GEN IV reactor objectives, thereby marking a break with the former Phénix and SuperPhénix French Sodium Fast Reactors.

The main objective to improve safety levels compared with current GEN II or III reactors has led to a core design that shows intrinsically safe behaviour in transients of accident situations. A negative or null sodium void worth can be obtained by using a heterogeneous core (CFV). Other challenging objectives such as minimising the secondary sodium activity level, preventing reactivity insertions due to core compaction and applying costkilling measures have led to a rather new fuel sub-assembly design compared with former SFRs.

This paper describes the design of the fuel sub-assembly (S/A) for the ASTRID CFV v4 core (version 4) at the end of the Conceptual Design phase with focus on innovations and specificities.

CFV Core

The layout of the CFV v4 core remains the same as the v3 core previously described in [START_REF] Venard | The ASTRID Core at the Midterm of the Conceptual Design Phase (AVP2)[END_REF] The major changes integrated into v4 concern some specific details of the sub-assembly design, with no significant impact on the core performance levels.

The CFV v4 core architecture consists of 180 inner fuel sub-assemblies and 108 outer fuel sub-assemblies, surrounded by 3 rows of reflector sub-assemblies with MgO, 3 rows of shielding sub-assemblies with B 4 C, 3 rows of shielding sub-assemblies with MgO, and another 2 rows of shielding sub-assemblies with B 4 C. The lateral reflector and shielding arrangement was optimised to provide enhanced radial neutron absorption [START_REF] Chapoutier | ASTRID Core Shielding -Design Studies and Benchmark Analysis[END_REF]. The core's absorber device relies on an innovative architecture composed of 9 control/shutdown rods and 9 diverse control/shutdown rods (Guénot-Delahaie et al., 2013), which are used both to manage core reactivity during the cycle, and to shut down the reactor. The other main characteristics of the ASTRID core include the integration of Complementary Safety Devices (CSD) for severe accident prevention [START_REF] Guénot-Delahaie | The Innovative RBH Complementary Safety Device for ASTRID to Address Unprotected Loss of Flow Transients: From Design to Qualification[END_REF] and mitigation, and an internal storage and debugging positions in the lateral shielding area.

The safety and performance goals [START_REF] Varaine | Pre-Conceptual Design Study of ASTRID Core[END_REF] assigned to the CFV core can be summarised by:

• Intrinsically safe behaviour, possibly supplemented by CSD during Unprotected Loss-Of-Flow (ULOF) transients and unprotected Control Rod Withdrawal (CRW) events. • Negative sodium void worth (CFV core).

• Search for the best performance in terms of fuel burn-up and S/A residence times.

The core is managed on the basis of 4 batches of fuel sub-assemblies. The objective fuel residence time is 1440 EFPD and the equilibrium cycle length is 360 EFPD. The fuel S/A positions are divided into 5 different flow zones to get a flat temperature distribution. The sodium flow rate required through each zone, varying between 22.5 and 26.7 kg/s, is calculated according to the fission heat generation and design criteria of a nominal cladding temperature lower than 620°C (best estimate) and an outlet temperature dispersion between two neighbouring sub-assemblies lower than 50°C. The whole core pressure drop is close to 3 bars.

Overview of the fuel sub-assembly

Although the overall design of the ASTRID fuel S/A mainly relies on feedback from French SFRs Phénix and SuperPhénix, the CFV core concept involves some appreciable differences [START_REF] Varaine | Pre-Conceptual Design Study of ASTRID Core[END_REF]:

• Axially heterogeneous fuel pins with a thick internal fertile plate for the inner core.

• Longer fissile zone for the outer core.

• Large cladding/small spacer wire fuel bundle.

• Sodium plenum (30 to 40 cm) above the fissile area.

• Enriched B 4 C zone at the bottom of upper neutron shielding. Due to washing considerations, this component is made removable through the S/A head.

A schematic diagram of the CFV v4 fuel S/A is given in Fig. 1. The main parts are detailed in the following paragraphs. Cost-killing processes performed at the end of the project's pre-conceptual design phase (AVP1) [START_REF] Chenaud | State of the Art of the Conceptual Designs for ASTRID Control and Shutdown Rods[END_REF] confirmed the objective of minimising the S/A height to decrease the overall costs of the reactor and the fuel cycle. A core height reduction of about 0.5 m was possible for ASTRID thanks to removal of the upper axial blanket used in former SFRs, coupled with the choice of a compact pressure drop device making it possible to shorten the spike by another 0.4 m. This resulted in a total S/A height of 4.50 m which has been reduced by 0.9 m compared with SuperPhénix (5.40 m), and which is not much higher than Phénix S/A (4.30 m).

Detailed description

Fuel Pins

One major innovation that has helped reducing the core's sodium void effect by about 5$ [START_REF] Varaine | Pre-Conceptual Design Study of ASTRID Core[END_REF][START_REF] Chenaud | State of the Art of the Conceptual Designs for ASTRID Control and Shutdown Rods[END_REF] compared with previous reactors is the wide fuel pin/small wire concept. The ASTRID design consists of 217 fuel pins in each S/A with a 9.70 mm outer diameter separated by a 1 mmdiameter spacer wire helically wound around the pins. The pin's lower plugs are mounted on the rails of a stainless-steel single-part grid to form a bundle (Fig. 2) that is vertically held within the hexagonal wrapper tube. The inner and outer core fuels pins differ slightly by their composition and their length (Fig. 3). The CFV core is characterised by axially heterogeneous fuel pins, with a UO 2 fertile zone inside the (U,Pu)O 2 fissile area for the inner fuel sub-assemblies. A 10 cm-longer fissile zone for the outer fuel sub-assemblies provides a global asymmetrical, crucible-shaped core. Fuel pins in the core are provided with a lower fertile axial blanket of 30 cm. Fissile pellets have a central hole to improve margins with respect to the fuel melting temperature, while fertile pellets are solid to increase the breeding gain, to reduce the loss of reactivity and to make it possible to be differentiated from fissile pellets. Fuel pins are hermetically sealed once filled with a helium-based mixture at 1 bar. Thermomechanical analysis based on finite-element computations of the fuel pins during irradiation have started with the CEA code called GERMINAL V2 [START_REF] Lainet | Recent Modelling Improvements in Fuel Performance Code GERMINAL for SFR Oxide Fuel Pins[END_REF] within the PLEIADES fuel simulation platform (Helfer et al. 2015). The aim of this analysis is to check that the design criteria (maximal temperatures, mechanical constraints, cladding swelling, etc.) are met in nominal conditions and incident transients (CRW, etc.), and to estimate a number of pin failures to be expected during severe accident sequences using a statistical analysis technique. By avoiding excessive pressurisation in the pins, the lower plenum length will undoubtedly be the key point of these in-depth studies. In nominal conditions, the end-of-life pressure in the pins reaches about 40 bars. Life duration of the fuel element is limited by the cladding embrittlement due to the high radiation damage rate. The reference cladding material for ASTRID is composed of special 15%Cr-15%Ni austenitic steel [START_REF] Flem | Status of the French R&D on ASTRID Core Materials[END_REF] which was tested in Phénix: AIM1 (Austenitic Improved Material #1). The level of cold-work (~20%), stabilisation by titanium, and control over the amounts of minor elements (Si and P) have made it possible to delay the beginning of swelling under irradiation, which makes AIM1 suitable for use up to high dose rates. AIM1 has also been chosen for the spacer wire to prevent differential behaviour under irradiation. A dedicated qualification programme is required on the upper plug/cladding welding which is subjected to a higher dose rate in the CFV core.

Wrapper Tube and Spacer Pads

The closed leaktight hexagonal wrapper tube (hexcan) is made of EM10 steel which has been extensively tested in Phénix. This martensitic 9%Cr-1%Mo steel exhibits excellent dimensional stability under irradiation (Le [START_REF] Flem | Status of the French R&D on ASTRID Core Materials[END_REF]. A record dose of 155 dpa was reached on the EM10 hexcan during the BOITIX-9 experiment in Phénix, with a resulting irradiation swelling limited to 0.5%. At the beginning of the ASTRID project, the slight geometrical variation expected for the EM10 wrapper tubes at the target dose has made it possible to optimise the core by reducing the duct wall thickness and the inter-assembly gap compared with austenitic 316Ti SuperPhénix hexcans. The external and internal widths of the ASTRID hexcan are respectively 168.7 mm and 161.5 mm.

To ensure suitable clearance between the sub-assemblies in the core, a rectangular spacer pad is stamped through each side of the hexcan just above the fuel pins (Fig. 4). This pad design was originally used in Phénix and SuperPhénix, and renewed for ASTRID because of its reliability and its simplicity. Not only do the pads play a central role in the mechanical equilibrium of the core, they also enhance safety by preventing reactivity increases due to a core compaction that follows a flowering, which may occur when core is subjected to dynamic stresses such as earthquake or internal pulse load release. Following the safety requirements for ASTRID wherein core compaction must be minimised, studies have been engaged to analyse the effect of a pad stiffness increase on the core mechanical behaviour. For this purpose, extensive studies have been performed to increase the stiffness of the stamped pads [START_REF] Blanc | Characterization, Simulation and Improvement of Spacer Pads Mechanical Behaviour for Sodium Fast Reactor Fuel Subassemblies[END_REF].

Finite element calculations have been done using the LICOS code (Helfer et al. 2015b) based on the CAST3M solver within the PLEIADES fuel simulation platform (Helfer et al. 2015a). A model was first developed to reproduce the cold stamping process and obtain representative pad geometry, and then was validated on the SuperPhénix pads geometry thanks to comparison with available experimental results of duct crushing tests at 550°C. A parametric study considering variations in the pad dimensions (height, width, radius in the corners) was performed for the 6-loaded-pad configuration using an elastoplastic law of EM10 at 550°C. More than 200 plots were generated to cover the range of parameters (Fig. 5).

Thickness

Wrapper tube

Stamped spacer pads Stiffness is defined by the slope of the curves in the elastic zone (small displacements). The highest stiffness of 192 kN/mm was obtained for the extremal parameters: maximal width, maximal height and minimal radius of the pads. These computations showed that pad stiffness could be increased by a maximal ratio of ~7 compared with the previous reference design.

The final pad design for ASTRID will be defined following the core static and dynamic mechanical behaviour studies that are currently in progress.

Upper Neutron Shielding

Up to now, the upper neutron shielding is by far the component that has required the most brainstorming and iterative studies. The result is an innovative design covered by a patent application [START_REF] Lorenzo | Assemblage Combustible pour Réacteur Nucléaire de Type RNR-Na, à Boîtier Logeant un Dispositif de Protection Neutronique Solidarisé de Manière Amovible[END_REF].

ASTRID Requirements

The CFV core's requirement for ASTRID makes it impossible to reuse the upper neutron shielding designs of Phénix and SuperPhénix. Indeed, the new requirement specifies a negative sodium void worth. For this reason, the bottom of upper neutron shielding, located just above the sodium plenum, should reflect neutrons as little as possible. This required choosing the best neutron absorber -boron carbide (B 4 C) -which benefits from extensive feedback since it has been used in SFR absorber rods. Core studies [START_REF] Venard | The ASTRID Core at the Midterm of the Conceptual Design Phase (AVP2)[END_REF] concluded that the optimum design for reaching the CFV effect would be a lower part made of 7.5 cm-thick B 4 C enriched to 90% in 10 B for enhanced absorption, together with the smallest steel structure thickness to minimise neutrons reflection.

The second main function assigned to the upper neutron shielding is obviously to provide proper axial neutron shielding to minimise damage to reactor internals and reduce the activation of surrounding Na circuits. The specification for ASTRID is a target 24 Na activity level less than 10 Bq/cm 3 (best-estimate calculation) for the secondary sodium in the intermediate heat exchangers (IHX) to ensure that the ASTRID steam generator building remains classified as a non-controlled area. This objective is much more severe than for Phénix, and also more difficult to meet than for SuperPhénix because of the distance between the IHXs and the core which is reduced. Moreover, considering the inherent tendency of the CFV core to amplify upper neutron leakages (due to the suppression of the upper fertile axial blanket), and the limited available space (about 90 cm) in the top the S/A (to meet a total height of 4.50 m), it is clear that the upper neutron shielding is a key issue.

Extensive 3D Monte-Carlo calculations were conducted to optimise the ASTRID core shielding, considering the best "performance vs. cost" ratio [START_REF] Chapoutier | ASTRID Core Shielding -Design Studies and Benchmark Analysis[END_REF]. This consisted in identifying the best arrangement for the reflector, moderator and absorber materials in the lateral shielding sub-assemblies, as well as the material in the top part of the upper neutron shielding. Concerning the lateral shielding, successive alternations of MgO and B 4 C sub-assemblies were adopted (see Section 0). Concerning the upper neutron shielding, the top part was filled with natural B 4 C (19.78% of 10 B), the best neutron absorber, which gave almost as good results as 10 Benriched B 4 C, while being much less expensive. Additional neutron shields made of borated steel sleeves were added in parallel to the IHXs to decrease the secondary sodium activation to an acceptable level of about 8 Bq/cm 3 .

Upper Neutron Shielding Design

Knowing that the upper shield needed to be made of both enriched and natural B 4 C (see previous section), the other complex part of the task consisted in proposing a design compatible with the design criteria and operational constraints.

The axially heterogeneous fuel pins provide a maximal flux shifted towards the top of the sub-assemblies compared with former SFRs. This means that the B 4 C-enriched lower part of the neutron shielding is subjected to a higher flux (3.0E+14 n/cm 2 /s) which produces a large volume of helium (~330 dm 3 ) and dissipates a great deal of power (up to 50 W/cm 3 ) due to the exothermal 10 B(n,α) 7 Li neutron capture reactions. The latter aspect led to the choice of B 4 C in the form of a 19-pin bundle to meet the maximal temperature criterion, together with a Na-bond to take advantage of the high thermal conductivity of sodium. Non-leaktight pins design also provided a solution to the problem of significant helium production. By allowing He to be released from the pins, any unacceptable pressurisation is avoided inside the pins and the height of the pins -i.e. the fuel S/A -is by the way minimised, because axial plenum is no longer necessary.

But Na-bonded pins maybe raised the most important issue that needed to be resolved. Absorber pins comprise a lower and an upper porous vent through which helium is released and sodium flows inside the pins. Due to the fragmentation of B 4 C under irradiation, small fragments may plug the micrometric pores of the vents and hinder pin draining. This is made even more difficult when considering the sub-millimetre cladding-pellet gap inside the pins. It is well-known that the sodium inside the Na-bonded absorber pins cannot be totally removed after draining. Yet this non-negligible amount of Na trapped inside the pins engenders unacceptable safety risks during the fuel S/A washing or underwater storage phases due to possible uncontrolled exothermal sodium-water reactions The solution to this problem was found following a brainstorming and value engineering process: the whole upper neutron shielding would be made removable on-line through the assembly head just before the washing operations.

The main difficulty was to design a device which enables the upper shielding extraction and that fulfilled the following requirements:

• The whole upper shielding must be removable on-line from the S/A head while in the washing pit.

• The smallest possible axial dimension must be applied to maximise the B 4 C pin height within the overall available space. • The upper neutron shielding must be securely locked in position in the S/A, by acting against the hydraulic lift force.

The extracting device is located at the top of the upper neutron shielding. It was designed on the basis of the experimental capsules head from SuperPhénix which were removable from the S/A in the outer Na storage area. A strict design process was carried out to reduce its length by a factor of ~4. The upward axial locking of neutron shielding is ensured by three lateral fingers spread out inside a groove in the S/A head and maintained thanks to the weight of a movable solid stainless steel part (Fig. 6). The downward axial locking is ensured by stopping fins located in the hexcan corners and mounted on a cylindrical skirt fixed to the S/A head structure. The three locking fingers are withdrawn by the upward movement of the movable part, allowing the upper neutron shielding to be extracted by the top through the S/A head. The translation motion of the movable part is ensured by a specific mechanical grapnel actuated in the washing pit and whose fingers are inserted in an internal groove. This grapnel is very similar to the one used for S/A lifting but with smaller dimensions.

The removable upper neutron shielding consists of a 14.5 cm-diameter and 94 cm-long cylinder. These dimensions comply with a S/A height of 4.50 m while fulfilling all ASTRID requirements. Fig. 6. Removable upper neutron shielding inside the S/A (section view).

The lower part of the neutron shield, in its first design, comprises the 19-pin bundle equipped with flow equalisers, while the upper part comprises the extracting head as described above. The pins are attached to the fixed part of the extracting head thanks to a grid support composed of parallel rails, like the fuel pins (see Section 4.1). The pin design is equivalent to the one used for the absorber rods. It comprises a cladding made of 15-15Ti AIM1 stainless steel, capped by two plugs equipped with a porous vent. From bottom to top inside the pins, there are 7.5 cm of B 4 C enriched to 90% in 10 B, 69.5 cm of natural B 4 C, and an upper plenum comprising a spring. The pins have a diameter of 25.8 mm and are separated by a 2 mm-diameter helical spacer wire. The surface fraction of B 4 C is 26.5%. B 4 C pellets are confined by a stainless steel shroud and are fabricated according to a "carbothermal" process instead of the "magnesiothermal" process (Phénix feedback). Like the control and shutdown rods (Guénot-Delahaie et al., 2013), the shroud and the "carbothermal" B 4 C are both the subject of in-depth R&D actions that aim to double the pin lifetime from start-up value which is based on feedback (two cycles of 360 EFPD) to the required full lifetime at equilibrium (four cycles of 360 EFPD).

Thermal-Hydraulic CFD Computations

CFD computations were performed first on a preliminary design, then on an optimised design. The movable part in the extracting head does not only ensure the locking/unlocking and extraction of the upper neutron shielding, but also the sodium flow path to the S/A outlet. A small central hole in the preliminary design was necessary to house the actuating mechanism of the three fingers. This sodium section reduction was expected to increase the pressure drop in the S/A. Three peripheral bean-shaped holes were consequently added between the fingers to increase the hydraulic diameter and reduce the total pressure drop (Fig. 7). Preliminary thermal-hydraulic calculations on this first design were performed to evaluate the pressure drop and to characterise the sodium flow fields (velocity and temperature) above the S/A outlet at the core monitoring system level (thermocouples and flowmeter). It is required to have the more homogeneous flow at the monitoring level to ensure a proper protection of the core, with measured values as close as possible to the mean values at the S/A output. Simplified but representative modelling of the S/A upper zone was used with the CFD STAR-CCM+ code, starting from fuel pins output up to the monitoring system in the hot pool. A constant sodium flow rate was considered at the inlet of the modelled area, corresponding to the maximum value in the core (26.7 kg/s), but with the heterogeneous temperature map that was separately computed at the fuel bundle output. The mean temperature is 555°C at the S/A outlet.

Wrapper tube

Removable upper neutron shielding

Movable part

The sodium velocity and temperature fields at the sub-assembly outlet are shown in Fig. 8. This figure shows a heterogeneous flow in temperature with a gradient of ~13°C between the center and the periphery in the monitoring section plane, with a measured value 9°C higher than the mean temperature. The velocity field shows that the flow is much more disturbed as it goes through the extracting head of the upper shielding. The central and the three peripheral holes in the movable part create four separated flows at the S/A outlet. The peripheral flows have a velocity higher than the central one, which may alter the measure by the flowmeter. The total pressure drop in the upper shielding is ~0.6 bar, 90% of which is located in the absorber pins bundle because of the wires and flow equalisers that reduce the sodium section. These preliminary computations revealed several ways of optimising the pins bundle and extracting head design in order to improve the temperature and velocity fields at the sub-assembly outlet and to decrease the total pressure drop.

Movable part

Sodium flow

First, the wires and flow equalizer were removed, in order to enlarge the sodium section in the pins bundle. The parallel grid support was replaced by a symmetrical star-shaped one, to avoid heterogeneous flow separation, and a similar star-shaped holder was included at the bottom of the pins bundle, to ensure its stability now that wires do not endorse this role anymore. At last, the movable part in the extracting head was simplified: the three peripheral holes were suppressed while the central hole was enlarged with a smoother entrance to reduce the pressure drop (Fig. 9). This figure shows a more symmetrical and flat response in temperature with a gradient of ~5°C between the center and the periphery in the monitoring section plane, with a measured value just 4°C higher than the mean temperature at S/A outlet. The field shows that the flow is much more homogeneous and acceptable too with a single flow at the monitoring system level. The velocity is higher in the center of the flow which is supposed to increase the accuracy of the measurement. The total pressure drop in the upper shielding is decreased to ~0.3 bar, which represents a reduction by ~0.3 bar comparatively to the previous design.

Lifting Head

The S/A lifting head was designed in coherence with the upper neutron shielding (see previous section). The objective of maximising the B 4 C volume in the removable upper shielding led to increasing the head inner diameter (15 cm) through which it is extracted. The sub-assembly head is a ~1 meter-long single part made of AISI 316L stainless steel. The top part is cylindrical while the bottom part is hexagonal (Fig. 11).

The cylindrical top part of the head is 27 cm-long. It comprises an internal groove at the top for subassembly lifting using the handling grapnel. Several radial holes will be included to ensure proper S/A cooling in the case of exceptional situations where heads are uncovered for visual inspection. The top outer surface is used for surface detection purposes by the under-sodium visualisation system (VISUS), which is designed to ensure that the S/A head does not stick out further than the others. To prevent handling errors, a bar code is engraved on the edge of the head for sub-assembly identification by under-sodium ultrasonic readings. The bottom part of the head has an outer hexagonal section while retaining the inner cylindrical hole. It is extended downwards to the interface between the enriched B 4 C and the natural B 4 C of the upper neutron shielding (Fig. 6). This design was necessary to provide the negative sodium void effect by limiting the amount of steel (neutron reflector) in the bottom of the upper shielding. Moreover, the steel in the upper part of the head, which is next to the natural B 4 C, improves the global upper neutron shielding and helps to achieve a low secondary sodium activity level. The hexagonal part of the head is inserted inside the wrapper tube with a gap. The binding with the wrapper tube is ensured by cold stamping after mounting of the fuel pins bundle.

Na velocity field

Spike

All core sub-assemblies are supported vertically with the spike inserted in the shroud tubes provided in the grid plate. The bearing surface is spherical on the spike side and conical on the grid plate side (Fig. 12). The shroud tubes comprise multiple holes to allow cold sodium from the grid plate to flow into the S/A spike through oblong slots.

Two helicoid labyrinths are provided on the spike with a reduced gap with the shroud tube. The function of the top labyrinth is to control the leakage flow of cold sodium -~200 kg/s for the whole core -into the hot pool in the interspace between sub-assemblies. The spherical bearing is provided with six slots to prevent cavitation and erosion at this level due to the leakage.

The function of the bottom labyrinth is to ensure the minimum flow required for main vessel cooling. It also prevents the S/A from hydraulic lifting or floating, which may occur due to high coolant drag pressure, by generating a hold-down force owing to the pressure drop developed here. The choice of a compact pressure drop device for ASTRID made it possible to reduce the spike length by about 0.4 m compared with the SuperPhénix design equipped with diaphragm-plates device. Sub-assemblies of each flow zone in the core are equipped with a compact pressure drop in the spike to achieve the required flow rate through the S/A (Fig. 1). This compact device consists of multiple aligned or crossed honeycomb orifice plates (Fig. 13), within a total length of less than 15 cm. This device has been successfully experimented with some sub-assemblies in Phénix and SuperPhénix. Additional protection against handling errors during refueling is provided thanks to discriminatory features located at the bottom of the spike. Those discriminators have a specific design for each flow zone in the core. They prevent a fuel S/A from being completely lowered into the wrong position in the core, i.e. in a position where the S/A would be under-cooled, or in a position of an absorber rod that would lead to a reactivity increase. This protection involves a combination of inner and outer diameters at the spike's end that are reproduced on a lock located inside the shroud tubes (Fig. 14).

Though this discriminator concept comes from SuperPhénix, the novelty for ASTRID is that the spike end is provided with self-orientation devices. These devices are arranged in parallel to the discriminator at the bottom of the spike, instead of being located at the top just above the grid plate as was the case for former SFRs. The discriminator consists of two cams at the spike end coupled with six vertical grooves inside the shroud tubes. The top of each groove is provided with conical ramps to rectify the angular mis-orientation. While refueling, self-orientation guarantees the angular position of each S/A with a good accuracy with respect to the grid plate. Such an orientation with respect to the fixed structure of the grid plate prevents the S/A network from angular drift that may cause S/A blockages while refueling, and guarantees the proper S/A orientation even in case of refueling after a WCD, i.e. without adjacent sub-assemblies. The spike is made of AISI 316L austenitic steel. Its thickness will be defined according to the results of the core dynamic mechanical studies that are currently underway.

All core sub-assemblies have the same spike design except for the discriminator which depends on the flow zone, and for the mitigation CSD which have a specific spike.

Conclusion

The fuel sub-assembly design for the ASTRID CFV v4 core at the end of the Conceptual Design phase (AVP2) is described.

Innovative design choices have been made to meet the ASTRID project requirements, marking a break with the former Phénix and SuperPhénix French SFRs. A negative sodium void worth (CFV core) will be achieved thanks to axially heterogeneous fuel pins, a wide cladding/small spacer wire bundle, a sodium plenum above the fuel pins, and upper neutron shielding with both enriched and natural B 4 C sodium-bonded pins.

The upper neutron shielding will help to reach a low secondary sodium activity level and will be removable on-line through the assembly head so as to meet washing constraints. Studies have been performed to increase the stiffness of the stamped spacer pads in order to analyse its effect on the core mechanical behaviour during hypothetical radial core flowering and compaction events. To prevent handling errors, an under-sodium ultrasonic bar code identification system will equip the lifting head, and the spike end will be fitted with discriminatory features arranged in parallel to self-orientation devices.

The beginning of the next Basic Design phase now focuses on consolidation of main design choices, performance levels confirmation, design optimization, as well as the start of qualification programme.

Fig. 1 .

 1 Fig. 1. ASTRID CFV v4 fuel subassembly.

Fig. 2 .

 2 Fig. 2. Fuel pin bundle with the grid support.

Fig. 3 .

 3 Fig. 3. Fuel pins composition (dimensions in cm).

Fig. 4 .

 4 Fig. 4. Section of a hexagonal wrapper tube with stamped spacer pads.

Fig. 5 .

 5 Fig. 5. Pad design optimisation -Crushing computations.

Finger

  Fig. 7. Simplified cross-section of the upper shielding extracting head (fingers not represented). First design.

  Fig. 8. Primary CFD computations at the S/A outlet (first design). Left: along Z-axis; Right: in monitoring section plane.

Fig. 9 .Fig. 10 .

 910 Fig. 9. Simplified cross-section of the upper shielding extracting head (fingers not represented). Second design. The sodium velocity and temperature fields at the sub-assembly outlet are shown in Fig. 10.

  Fig. 11. Lifting head.

Fig

  Fig. 12. Spike.

Fig. 13 .

 13 Fig. 13. Grid plate of the pressure drop device.

Fig. 14 .

 14 Fig. 14. Discriminatory and self-orientation features.

Acknowledgments

The authors wish to thank all the teams involved in the fuel S/A design studies at CEA and AREVA-NP for their support.


R Le Tellier
email: romain.le-tellier@cea.fr

E Skrzypek

L Saas

On the treatment of plane fusion front in lumped parameter thermal models with convection

Keywords: plane fusion front, two-phase Stefan problem, lumped parameter model, focusing effect Research highlights

Within the framework of lumped parameter models for integral codes, this paper focus on the modelling of a two-phase Stefan fusion problem with natural convection in the liquid phase. In particular, this specific Stefan problem is of interest when studying corium pool behavior in the framework of light water reactor severe accident analysis. The objective of this research is to analyze the applicability of different approximations related to the modelling of the solid phase in terms of boundary heat flux closure relations. Three different approximations are considered: a quadratic profile based model, a model where a parameter controls the power partitioning at the interface and the steady state conduction assumption. These models are compared with an accurate fronttracking solution of this plane fusion front problem. This "reference" is obtained by combining the same integral conservation equations as the approximate models with a mesh-based solution of the 1D heat equation. Numerical results are

discussed for a typical configuration of interest for corium pool analysis. Different fusion transients (constructed from nondimensionalization considerations in terms of Biot and Stefan numbers) are used in order to highlight the potential and limitations of the different approximations.

Introduction

Mathematical and numerical description of solidification and melting of materials is a substantial issue arising in numerous engineering disciplines, including power and manufacturing engineering. For instance, in the context of nuclear power engineering, melting processes are of prime interest during a severe accident progression as a consequence of the insufficiency of the nuclear reactor core cooling. During this kind of event, the residual heat produced in the nuclear fuel coming from the fission products decay, is transferred depending on the system state and boundary conditions by convection, conduction or radiation across internal materials of the reactor pressure vessel.

For the special case of single pure material phase change, the moving boundary problem reduces to the so-called Stefan problem. In the Stefan problem, the moving phases interface and the associated conditions are non-linear with unknown location of the interface changing in time. The source of non-linearity lies in three aspects: the phase change itself at the interface, the properties of the material that are usually temperature dependent and lastly the boundary conditions of the domain that can be non-linear. This leads to a problem that is space and time dependent, which numerical solution aims at determining the interface position in the analysed domain (see, for instance, the review in [START_REF] Hu | Mathematical modelling of solidification and melting: a review[END_REF]).

Modelling of the solidification and melting becomes a subject of many studies in which the behaviour of the material under various conditions is predicted.

The problem was investigated in different ways that involved both analytical, numerical and experimental research.

In this paper, we are interested in the numerical simulation of such problems.

From this computational point of view, the different approaches found in the litterature can be classified into two broad categories: mesh-based discretizations of the local conservation equations (including fixed-grid and moving-grid methods) and lumped parameter models.

On the one hand, within the framework of mesh-based approaches, fixedgrid methods have received a lot of attention for their computational efficiency when 2D or 3D geometries are involved. The various methods differ in terms of the formulation for the energy conservation equation (see, for instance, the review in [START_REF] Voller | An overview of numerical methods for solving phase change problems[END_REF]). When convection in the liquid is to be taken into account, such methods becomes computationally demanding as the momentum conservation equation has to be solved; Computational Fluid Dynamics (CFD) tools are then used. For instance, in [START_REF] Assis | Numerical and experimental study of melting in a spherical shell[END_REF], CFD calculations were performed and compared with experimental results in order to validate the numerical approach based on the "Volume Of Fluid" method and an enthalpy formulation. In [START_REF] Wang | A comprehensive numerical model for melting with natural convection[END_REF], a temperature transforming model is used in order to model melting with natural convection.

On the other hand, lumped parameter models for phase change problems are developed and used in so-called "0D" or integral codes that require fastrunning models, in particular when Monte-Carlo based sensitivity analyses are performed. In this context, one important case of interest is when the solid can be treated as one-dimensional. Many different approximations and simplifications have been proposed in this case (see [START_REF] An | Lumped parameter model for one-dimensional melting in a slab with volumetric heat generation[END_REF] and citations therein). For instance, in [START_REF] Crepeau | Approximate solutions to the Stefan problem with internal heat generation[END_REF], under the hypothesis that the Stefan number is lower than one, a quasi-static approximation was discussed and compared in [START_REF] Crepeau | On the Stefan problem with volumetric energy generation[END_REF] with CFD calculations with and without convection; a good agreement was reported. In [START_REF] An | Lumped parameter model for one-dimensional melting in a slab with volumetric heat generation[END_REF], closures for the complete set of integral energy equations (liquid, solid and interface) are obtained based on Hermite approximations for integrals that define the average temperatures and boundary heat fluxes. Good agreement was shown with respect to both analytical formula for a semi-infinite medium without internal heat generation and a fixed-grid enthalpy method for internally heated slab configurations.

In this paper, the specific Stefan problem we are interested in is a twophase problem where convection in the liquid is taken into account through correlation-based closures and conduction in the solid can be approximated by the heat equation in a 1D slab. This specific Stefan problem is of interest when studying corium pool behavior in the framework of light water reactor severe accident analyses. In this context, integral models are used and the liquid phase modelling heavily relies on correlations obtained for natural convection in different configurations and for different regimes [START_REF] Bonnet | Thermohydraulic phenomena in corium pool: the bali experiment[END_REF]. To the authors knowledge, a comprehensive analysis of the possible approximations of this Stefan problem from the point of view of solid conduction treatment is missing. This is the object of this paper where different approximations for the solid phase in such a lumped parameter modelling are compared in detail. We have considered different approximations that are in use in severe accident analysis comprising a quadratic profile based model (that is shown to be equivalent to the approach in [START_REF] An | Lumped parameter model for one-dimensional melting in a slab with volumetric heat generation[END_REF]), a parametric model reported in [START_REF] Tellier | Phenomenological analyses of corium propagation in LWRs: the PROCOR software platform[END_REF] and the steady state conduction assumption (similar to [START_REF] Crepeau | Approximate solutions to the Stefan problem with internal heat generation[END_REF]). The important feature of this comparison is that it is carried out with respect to a "reference" model that combines the same integral formulation of the conservation equations as the approximate models with a mesh-based solution of the 1D heat equation in the solid. In this way, an accurate front-tracking solution of this Stefan problem is obtained. It has the exact same modelling features as the approximate models under consideration except only for the solid boundary heat flux closure relations that are the object of our analysis.

The paper is structured as follows. Section 2 gives a brief presentation of the mathematical settings for the rest of the paper, both in terms of local and macroscopic conservation equations. From there, the different models to be compared are discussed in Section 3. The "reference" model is given in Section 3.1 while Section 3.2 details the different approximate models under consideration.

Numerical results and associated discussion are reported in Section 4 where a typical configuration of interest for light water reactor severe accident analysis is used along with nondimensionalization considerations (see Section 4.1) for constructing fusion transient cases. Finally, concluding remarks are given in Section 5.

Mathematical settings

This section briefly presents the mathematical settings of the physical model that represents a liquid domain Ω l and a solid domain Ω s separated by a sharp interface Γ ls associated with a plane fusion or solidification front. For the sake of clarity, the macroscopic balance equations are obtained from the local conservation equations and, then, they are simplified under the approximation of of a 1D slab solid.

Local conservation equations

The set of equations describing the system is composed of local conservation equations with the Stefan condition related to the moving interface.

For the liquid domain Ω l (t), the mass, momentum and energy conservation equations are written under the Boussinesq approximation in terms of the liquid velocity v l and the liquid temperature T l as:

∇ • v l ( r, t) = 0 for r ∈ Ω l (t) (1) ∂ v l ∂t ( r, t) + v l • ∇ v l ( r, t) = -1 ρ ref l ∇p( r, t) + ν l ∆ v l ( r, t) -β l T l ( r, t) -T ref l g for r ∈ Ω l (t) (2)
ρ ref l C p,l ∂T l ∂t ( r, t) + v l • ∇T l ( r, t) = λ l ∆T l ( r, t) for r ∈ Ω l (t) (3) v l ( r, t) • n = 0 for r ∈ ∂Ω l (4)
conditions in T l ( r, t) and/orλ l ∇T l ( r, t)

• n for r ∈ ∂Ω l (t) \ Γ ls (t) (5)
T l ( r, t) = T fus. for r ∈ Γ ls (t) [START_REF] Crepeau | Approximate solutions to the Stefan problem with internal heat generation[END_REF] with appropriate initial conditions (in particular, in terms of temperature T l ( r, 0) = T 0 l ( r)). The liquid properties are: ρ ref l the mass density at a reference tem-perature T ref l , ν l the kinematic viscosity, β l the volumetric thermal expansion coefficient, C p,l the specific heat capacity and λ l the thermal conductivity.

In the solid domain Ω s (t), the temperature T s is governed by the heat conduction equation:

ρ s C p,s ∂T s ∂t ( r, t) = λ s ∆T s ( r, t) for r ∈ Ω s (t) (7)
conditions in T s ( r, t) and/or -

λ s ∇T s ( r, t) for r ∈ ∂Ω s \ Γ ls (t) (8)
T s ( r, t) = T fus. for r ∈ Γ ls (t) [START_REF] Tellier | Phenomenological analyses of corium propagation in LWRs: the PROCOR software platform[END_REF] with the initial condition T s ( r, 0) = T 0 s ( r). The solid properties are: ρ s the mass density, C p,s the specific heat capacity and λ s the thermal conductivity.

Finally, the liquid-solid interface Γ ls velocity is described by a plane fusion front equation:

v ls ( r, t) = 1 ρ s ∆H fus. -λ l ∇T l ( r, t) • n ls + λ s ∇T s ( r, t) • n ls n ls for r ∈ Γ ls (t) (10)
with n ls is oriented from the liquid to the solid. The fusion properties are: T fus.

the fusion temperature, ∆H fus. the specific latent heat of fusion.

Note, that for the sake of conciseness and because of the selected applications in the numerical results section, no volumetric heat source has been introduced in these equations but the models presented hereafter are not restrictive with respect to this simplification.

Lumped parameter formulation

From there, letting aside the momentum equation in the liquid, the problem in both phases can be cast into an integral form in terms of macroscopic conservation equations for mass and energy. This approach leads to a so-called lumped parameter formulation (also called "0D" formulation).

In order to do so, the boundary ∂Ω l (t) is assumed to be partitioned under the form ∂Ω l (t) = i ∂Ω l,i Γ ls (t) and the following integral quantities are defined:

m l (t) = ρ l V l (t) (11)
Tl (t) = 1 V l (t) Ω l (t)
T l ( r, t)dV ( 12)

ṁls (t) = ρ s Γ ls v ls ( r, t) • n ls dS (13) φi l (t) = - 1 S i l ∂Ω l,i λ l ∇T l ( r, t) • ndS (14) φls l (t) = - 1 S ls (t) Γ ls (t) λ l ∇T l ( r, t) • n ls dS ( 15
)
where S ls (resp. S i l ) is the area of Γ ls (resp. ∂Ω l,i ).

125

In this way, using Reynolds transport and Gauss theorem, the integral formulation of Eqs. 1 and 4 is obtained as:

dm l (t) dt = ṁls (t) (16)
C p,l m l (t) d Tl dt (t) + ṁls (t) Tl (t) -T fus. = - i φi l (t)S i l (17)
φls l (t)S ls (t)

The same integration process can be carried out for the solid domain to obtain an integral formulation as:

dm s (t) dt = -ṁls (t) (18)
C p,s m s (t) d Ts dt (t) + ṁls (t) T fus. -Ts (t) = - i φi s (t)S i s (19) + φls s (t)S ls (t)
with similar notations as for the liquid phase and the conductive heat flux defined as:

φls s (t) = - 1 S ls (t) Γ ls (t) λ s ∇T s ( r, t) • n ls dS (20)
Finally, the interface Eq. 10 is written as:

ṁls (t) = S ls (t) ∆H fus. φls l (t) -φls s (t) (21)
Obviously, the practical use of such 0D equations requires associated closure relations or Neumann boundary conditions associated with the average heat fluxes φi l (t) and φi s (t). In the remainder, we will consider that adequate closure laws (i.e., in the form of Nusselt-Rayleigh correlations for natural convection) are available for the liquid phase 0D formulation associated with such a natural convection problem and we will focus our attention on the treatment of the solid phase.

1D slab solid approximation

In this study, the problem is further simplified under the hypothesis of a 1D slab fusion front as depicted in Figure 1. Note that the liquid shape is not specified except for its interface with the solid. The solid thickness e s (t) is given by e s = m s S ls ρ s . The x-axis origin is chosen on the external boundary of the liquid, the coordinate of the solid-liquid interface is given by e l (t) while the external boundary of the solid is at x s (t) = e l (t) + e s (t) In this case, the set of equations for the solid Eqs. 7, 8 and 9 can be simplified as:

ρ s C p,s ∂T s ∂t (x, t) = λ s ∂ 2 T s ∂x 2 (x, t) for x ∈ [e l (t), x s (t)] (22) conditions in T s (x, t) and/or -λ s ∂T s ∂x (x, t) for x = x s (t) (23) T s (x, t) = T fus. for x = e l (t) (24)
In the same way, the lumped solid conservation equations Eq. 20 is:

C p,s m s d Ts dt + ṁls T fus. -Ts = -φbc s + φls s S (25)
with S = S ls = S bc s and the average heat fluxes are:

φls s (t) = -λ s ∂T s ∂x (e l (t), t) (26) φbc s (t) = -λ s ∂T s ∂x (x s (t), t) (27)
In the remainder, the system of equations to be solved will consist in the set of ordinary differential equations (ODE) composed of:

• the mass conservation equations for both liquid and solid phases Eqs. 17

and 19 coupled by the mass flow rate associated with the plane fusion front equation Eq. 21;

• the lumped liquid energy equation Eq. 18 with appropriate boundary conditions and/or closure laws on ∂Ω l (t);

• the lumped solid energy equation Eq. 25 with heat fluxes φls s (t) and φbc s (t) that are given by model-dependent equations discussed in the next section.

Note that these equations are valid only when both liquid and solid phases are present but, actually, the modelling presented here (and numerical results) also encompasses the cases when the system is completely liquid or solid. For the sake of conciseness, the equation modifications associated with these system states are not given as they are straightforward. The only additional difficulty is related to the transitions from one state to the other (appearance or disappearance of one phase) and its treatment in the numerical time integration; some explanations are given in Appendix A.

Front propagation models

This section presents the different models pertaining to the closure of the ODE system of equations with solid boundary heat fluxes φls s (t) and φbc s (t).

Reference model

In order to compare different approximate models, a reference model has been constructed. A reference solution is obtained by solving explictly the 1D heat equation of Eq. 22 (with boundary conditions given by Eqs. 23 and 24) in order to calculate the average heat fluxes φls s (t) and φbc s (t). More precisely, the numerical scheme to treat the moving interface proceeds from time t n to time t n+1 = t n + ∆t as follows:

• Eq. 22 is integrated from t n to t n+1 based on an Euler implicit scheme i.e.

ρ s C p,s Ts (x, t n+1 ) -T s (x, t n ) ∆t = λ s ∂ 2 Ts ∂x 2 (x, t n+1 ) for x ∈ [e l (t n ), x s (t n )] (28)
conditions in Ts (x, t n+1 ) and/orλ s ∂T s ∂x (x, t n+1 ) for x = x s (t n ) (29)

Ts (x, t n+1 ) = T fus. for x = e l (t n ) (30) and a mesh-based discretization of this Laplacian equation.

• the ODE system is integrated from t n to t n+1 considering that ∀t ∈

[t n , t n+1 ], φls s (t) = -λ s ∂T s ∂x (e l (t n ), t n+1 ) (31) φbc s (t) = -λ s ∂T s ∂x (x s (t n ), t n+1 ) (32)
• from the previous integration, one obtains e l (t n+1 ) that is used to update the mesh associated with the solid slab and calculate T s (x, t n+1 ) by projection of Ts (x, t n+1 ) onto this new mesh over [e l (t n+1 ), x s (t n+1 )] and renormalization to ensure that

x s (t n+1 ) e l (t n+1 )
T s (x, t n+1 )dx = e s (t n+1 ) Ts (t n+1 ).

In practice, Eqs. 28, 29 and 30 are discretized using linear finite elements over a uniform mesh. Following [START_REF] Carey | Derivative calculation from finite element solutions[END_REF], special care was simply taken in the evaluation of the heat fluxes at the boundary from this finite element solution in order to obtain superconvergent estimates for essential boundary conditions. 185

Approximate models

In this paper, the previous reference model is used to assess the behavior of different approximations of φls s (t) and φbc s (t) that are used in integral models. In any case, these models give analytical closures for the ODE system.

Quadratic profile based model

A simple adequate choice, completely defined by these three constraints, is to use a second-order polynomial as it is solution of the stationary 1D slab heat conduction equation with an internal source. Such an approximation is used in the MAAP source term code [START_REF] Fauske | MAAP4 -modular accident analysis program for LWR power plants[END_REF] or the PROCOR software platform [START_REF] Tellier | Phenomenological analyses of corium propagation in LWRs: the PROCOR software platform[END_REF] in the framework of light water reactor severe accident analysis. In this case, Ts (x, t) is written as:

Ts (x, t) = a x s -x e s 2 + b x s -x e s + c (34)
with

c = T s (x s , t) (35) b = 6 Ts (t) -2T fus. -4c (36) a = 3 Ts (t) - b 2 -c (37)
The boundary heat fluxes are then given by φls

s (t) = λ s 2a + b e s = -
that are equivalent to Eqs. 38 and 39.

Parametric model

In a second class approximation, when the steady-state of the system is known, inequalities relating the different heat fluxes can be obtained and may be used to construct a parametric approximation. For instance, let us consider a system in liquid-solid state under fusion for which it is known that, at final steady-state, the fusion is complete. Then, assuming that the fusion front monotonously propagates with a monotous evolution of the solid average temperature, while Ts < T fus. , the two following inequalities hold during all the fusion transient: • if φls s = φls l , all the power goes to the solid heating until it reaches Ts = T fus. ; afterwards, the fusion front propagates with φls s = φbc s (i.e. Ts is constant);

dm l dt ≥ 0 (
• if φls s = φls l Sts+ φbc s 1+Sts
, Ts is constant while the fusion front propagates. Note that this approximation is handy for analytical calculations and is often used in "order of magnitude" evaluations such as in [START_REF] Seiler | A phenomenological analysis of melt progression in the lower head of a pressurized water reactor[END_REF].

Based on these relations, a parametric approximation of φls s may be used under the form: (

φls s =      min φls l , ω φls l + (1 -ω)
) 46
where ω is a dimensionless parameter ∈ [0, 1]. For Neumann or Robin boundary conditions applied on the solid surface at x s , this approximation of φls s can be supplemented by an approximation of φbc s or T s (x s , t) in order to close the ODE system of equations for specific configurations. For the results presented in Section 4, a linear profile approximation was used for T s (x s , t) i.e. T s (x s , t) = 2 Ts (t) -T fus. . Such a parametric model may be useful in Monte-Carlo based sensitivity analyses considering that the treatment of the ω parameter as a random variable is a sound choice for capturing the overall model uncertainty. It was developed in the PROCOR platform for this purpose [START_REF] Tellier | Phenomenological analyses of corium propagation in LWRs: the PROCOR software platform[END_REF].

Steady-state conduction assumption

As a last approximation, it is a common practice in "order of magnitude" evaluations such as in [START_REF] Seiler | A phenomenological analysis of melt progression in the lower head of a pressurized water reactor[END_REF] to assume that the temperature profile in the solid is, at any time, the steady-state one i.e. a linear profile in such a way that φls

s (t) = φbc s (t) = -λ s T s (x s , t) -T fus. e s (47)
In this approximation, the solution of the lumped solid energy equation Eq. 25 is

no longer calculated as it is not needed. Note that with such an approximation, this equation would reduce to

C p,s m s d Ts dt + ṁls T fus. -Ts = 0 (48)
In the general case, the inconsistency between this energy balance and the steady-state linear profile approximation gives rather poor results (this discrepancy will be illustrated by numerical results in the next section). It is only when the Stefan number is sufficiently small [START_REF] Crepeau | Approximate solutions to the Stefan problem with internal heat generation[END_REF] (i.e. the second term in Eq. 48 is small) that this approximation is of interest.

Numerical results

In this section, the different approximations detailed in the previous sections will be compared with the reference model in order to discuss their shortcomings. This discussion on the models validity is not exhaustive and should not be taken as a general conclusion on the model performances. However, for the sake of genericity, the cases of study are selected with respect to a nondimensionalization of the equations presented hereafter.

Nondimensionalization of the equations

This nondimensionalization is limited to the front propagation (Eq. 21) and solid-related equations (Eqs. 22 and 25) as the treatment of the conduction in the solid in such a fusion problem is the focus of this work.

The following notations are first introduced:

• scaling quantities: time t ⋆ , length x ⋆ , temperature T ⋆ , heat flux φ ⋆ ;

• dimensionless variables:

τ = t t⋆ , χ = x x⋆ , χ s = x s x⋆ , ǫ s = es x⋆ , ǫ l = e l x⋆ , ν = t⋆ x⋆ ṁls ρ s S , θ = T s -T fus. T ⋆ , ψ = φ φ ⋆ ; • dimensionless numbers: Fourier number F o = λ s t⋆ ρ s C p,s x 2 ⋆ , Biot number Bi = φ⋆x⋆ λ s T ⋆ , Stefan number St = C p,s T ⋆ ∆H fus. .
In this way, the front equation Eq. 21 is written as:

ν = F oBiSt ψls l -ψls s (49)
The heat equation and associated boundary conditions Eqs. 22, 23 and 24 take the form:

∂θ ∂τ = F o ∂ 2 θ ∂χ 2 for χ ∈ [ǫ l , χ s ] (50)
conditions in θ(χ, τ ) and/or -

1 Bi ∂θ ∂χ (χ, τ ) for χ = χ s (51) θ(χ, τ ) = 0 for χ = ǫ l (t) (52)
Finally, the lumped solid energy conservation equation Eq. 25 is rewritten as: As mentioned earlier, this nondimensionalization is incomplete as the coupling with the liquid equations (through ψls l ) is not treated but it is still useful in a qualitative manner. Indeed, one can notice that the dimensionless front velocity ν is directly proportional to the product BiSt while the dimensionless heat fluxes are inversely proportional to Bi; as a consequence, in the next section, the two following sensitivity analyses are considered: varying Bi while maintaining BiSt constant and varying St while maintaining Bi constant.

ǫ s d θ dτ -ν θ = F oBi ψls s -ψbc s (

Cases of study

The numerical results presented here are all based on the cylindrical geometry depicted in Figure 2 with a small aspect ratio (height over diameter). The system is heated from below with an imposed power φin l S; Rayleigh-Benard convection in the liquid controls the power partitioning between the lateral boundary (heat flux φlat l and surface S lat l ) and the solid/liquid interface. For these convective heat exchanges at the liquid boundaries, Nusselt-Rayleigh correlations are used: the Churchill & Chu and Globe & Dropkin correlations for the lateral and axial heat transfers respectively (see for instance [START_REF] Tellier | Transient stratification modelling of a corium pool in a LWR vessel lower head[END_REF] for more details on such closures for corium pool lumped models). with T fus the material liquidus temperature, h 0 a heat transfer coefficient representative of the heat exchange with water (nucleate boiling), T 0 the water saturation temperature and ǫ the emissivity of the solid surface. In order to consider the conditions on ∂Ω bc s in the reference model, the radiative condition has been linearized in such a way that Eq. 29 is written in any case as a Robin condition:

-λ s ∂T s ∂x (x, t n+1 ) = h T s (x, t n+1 ) -T 0 for x = x s (t n ) (57) with h =      h 0 for convective bc ǫσ (T s (x, t n ) + T 0 ) (T s (x, t n )) 2 + (T 0 ) 2 for radiative bc (58)
For the sake of simplicity, pure iron has been considered for the steel material.

For the physical properties of solid and liquid phases, temperature-independent values were obtained from the PROCOR code [START_REF] Tellier | Phenomenological analyses of corium propagation in LWRs: the PROCOR software platform[END_REF] data (based on TOLBIAC-ICB code [START_REF] Spindler | Simulation of MCCI with the TOLBIAC-ICB code based on the phase segregation model[END_REF] subroutines) evaluated at T fus = 1811K.

Then, the test case conditions are defined by the following:

• the cylinder radius is 1.4m and the total mass of steel is fixed so that when the system is completely solid, its height is e max s = 0.3m;

• at t = 0 -, the system is at steady-state for a given φin,0

l = 0.4MW/m 2
and convective cooling on top with h 0 = 10 4 W/m 2 /K, T 0 = 372.78K in such a way that the system is composed of a solid phase on top of a liquid phase;

• at t = 0 + , the boundary conditions are modified as follows: φin l = 3 × φin,0 l and convective heat transfer on the top boundary is replaced by radiative heat transfer.

The transient that follows leads to the complete fusion of the solid. In the context of severe accident analysis, the complete fusion of the solid is an important event because when the steel layer is completely liquid, the heat flux on its lateral boundary largely increases and becomes a major threat for the vessel integrity. This so-called "focusing effect" phenomenon (related to the high conductivity of this steel layer and its small aspect ratio) is of major concern in the evaluation of the in-vessel retention strategy where the aim is to contain the corium materials within the reactor pressure vessel (see, for instance, [START_REF] Rempe | Potential for invessel retention through ex-vessel flooding[END_REF]).

Based on Section 4.1, the dimensionless numbers Bi and St were evaluated with x ⋆ = e max s , T ⋆ = T fus -1000 (1000K being an estimate of Ts (t = 0)) and

φ ⋆ = φin,0 l .
From there, in addition to this nominal case, four additional tests have been defined by modifying the solid physical properties:

• in two of these tests (denoted γ Bi = 0.5 and γ Bi = 5.0), both λ s and ∆H fus.

have been modified in order to decrease or increase Bi (multiplicative factor γ Bi ) while maintaining BiSt constant. Note that the Bi value cannot be decreased much otherwise, the initial state (as given by the steady-state of the system for the conditions at t = 0 -) and steady-state correspond to a completely solid steel layer.

• in the two other tests (denoted γ St = 0.1 and γ St = 5.0), ∆H fus. has been modified in order to decrease or increase St (multiplicative factor γ St ) while maintaining Bi constant.

These test cases and associated Bi, St values are given in Table 1 along with the initial solid thickness e s (0) and the time t f at which the solid fusion is completed. The different models and associated parameters that are compared are summarized in Table 2. For the parametric model, while the ω parameter was varied from 0 to 1 with 0.1 increment, the results reported hereafter are limited to ω = 0, ω = 1 and the ω value giving the closest values to the R(10 For all models and all cases, the ODE system of equations is solved with an Euler explicit scheme with a prescribed timestep of δt = 0.1s. For the reference model, the timestep ∆t in the coupling scheme between the ODE system and the 1D conduction equation (see Section 3.1) is set to ∆t = 5.0s (except in the γ Bi = 0.5 case, ∆t = 2.5s). These values of δt and ∆t are sufficiently small in order to ensure a proper model comparison (for the sake of conciness, the convergence with respect to these time discretizations is not shown here).

Results and discussion

For all cases, taking R(10 -4 ) result as the reference, for any other model M, the comparison is presented in terms of:

• the absolute error on the fusion time ∆t f = t f (M)t f R(10 -4 ) ;

• the relative error in L 1 -norm for any quantity q ∈ φls s , φbc s , Ts , φls l , e s defined as:

ε q L1 = q (M) -q R(10 -4 ) 1 q (R(10 -4 )) 1 (59)
with

• 1 calculated over the time interval [0, max t f (M) , t f R(10 -4 ) ].
These results for the nominal case are given in Table 3 where it can be seen that R(10 -3 ) results are very close to R(10 -4 ) demonstrating the convergence of this reference solution in terms of the spatial mesh size.

This nominal case gives a general view on the performances of the different approximate models for this typical configuration of interest in light water reactor severe accident analysis. The discrepancies of the Q model are limited, especially for the most important quantities t f and e s . Regarding the parametric model P, while as expected the performances largely depend on the ω parameter, for the "optimal" ω value of 0.3, the agreement on t f and e s is comparable with the Q model but the error on the heat fluxes at the solid boundaries are larger.

Finally, the steady-state conduction approximation appears, as expected, as a very poor approximation with t f underestimated by about 50%.

Varying Bi (maintaining BiSt constant)

The results for γ Bi = 0.5 and γ Bi = 5.0 cases are presented in Tables 4 and5 respectively in order to discuss the effect of increasing Bi (while maintaining BiSt constant).

Regarding the reference model, the discrepancy between R(10 -3 ) and R(10 -4 ) increases when Bi is increased: indeed, when the conductivity decreases, the transient thermal gradients within the solid are more pronounced and localized in such a way that a finer mesh at the boundary is needed.

Then, when Bi is increased, the "optimal" value of ω for the parametric model P decreases (from 0.4 to 0.1 when going from γ Bi = 0.5 to γ Bi = 5.0), while the L model results are largely improved. For both models, this trends can be related to the modification in the power partitioning between the solid heating and the fusion front propagation when Bi is varied. This variation of the power partitioning can be clearly illustrated by considering the following normalized quantity ∈ [0, 1]:

Γ ls φ = φls l -φls s φls l -φbc s ( 60
)
that represents, for a given power input from the liquid at the interface and a given top boundary heat exchange, the fraction that is "used" for the front propagation. 1 -Γ ls φ is then the fraction that goes into the solid heat balance (right hand side of Eq. 25). This quantity is depicted in Figure 3 for both γ Bi = 0.5 and γ Bi = 5.0 cases as calculated by the R(10 -4 ) model (note that, at t = 0 + , due to the continuity of φls l and φls s , Γ ls φ is zero). It can be seen that, in average, Γ ls φ becomes greater during the fusion transient when Bi is increased. This is consistent with the performance trends for the parametric and steadystate conduction models; indeed, with the parametric model, Γ ls φ is maximum for ω = 0 while the steady-state conduction assumption gives Γ ls φ = 1. Finally, the Q model results are only slightly affected by Bi modification and there is no specific trends (see Tables 3, 4 and5). The results for γ St = 0.1 and γ St = 5.0 cases are presented in Tables 6 and7 respectively in order to discuss the effect of increasing St (while maintaining Bi constant).

First of all, one can notice that the effect of refining the spatial mesh size from 10 -3 to 10 -4 is very limited and almost insensitive to the St value.

Then, when St is increased, the "optimal" value of ω for the parametric model P increases (from 0.1 to 0.5 when going from γ St = 0.1 to γ St = 5.0), while the L model results are deteriorated. Alike previous cases where Bi was varied, this trend is related to the power partitioning between the solid heating and the fusion front propagation illustrated in Figure 4 in terms of Γ ls φ (t): when St increases, Γ ls φ (t) is decreased. In addition, when St becomes small (i.e. for long fusion transient), the sensitivity of the parametric model results to the ω value becomes more important as can be seen from the error amplitude between ω = 0 and ω = 1 cases. While ω can still be adjusted to get results close to the reference (in the γ St = 0.1 case, with ω = 0.08, the discrepancy on the fusion time is reduced to 39.7s), the constant power partitioning hypothesis of this model becomes inappropriate. This is further illustrated in Figure 5 where the Ts (t) evolution as calculated by R(10 -4 ) and P(0.1) is given. The evolution of Ts (t) can be separated around 1500s into two time regions exhibiting very different slopes; as a consequence, in such a case, an "optimal" ω value is bound to underestimate (resp. overestimate) the slope of Ts (t) at the beginning (resp. end) of the fusion transient.

Approximate models limitations

To complete this analysis, the shortcomings of the different approximate models are illustrated in more detail by additional results.

First, an important feature of the Q model can be observed at the beginning of the transient: because of the quadratic profile asumption fitted on the average temperature Ts and both boundary conditions, the discontinuity at t = 0 of the boundary condition on ∂Ω bc s leads to a non-physical discontinuity of the heat flux φls s at the solid-liquid interface as depicted in Figure 6 for the γ St = 5.0 case. This discontinuity can induce a spurious behavior of the front propagation at the beginning of the transient. Indeed, as shown on Figure 7 for the γ St = 5.0 case, for the first 16.6s, Q model predicts solidification at the interface instead of fusion: the Q model cannot ensure the correct direction of the front at all time.

Then, as explained in Section 3.2.3, an inherent feature of the L model is that the solid average temperature Ts decreases during a fusion transient. As shown in Figure 8 in the γ Bi = 0.5, this decreasing trend can be very pronounced and leads to negative temperature that completely invalidate the transient solid heat balance with such a model. Finally, as mentioned before, the parametric model performance is limited by the fixed value of ω that imposes a constant power partition at the interface.

In addition, it is also inherently limited by the hypothesis d Ts dt ≥ 0 that may not be always verified: as illustrated in Figure 9

Figure 1 :

 1 Figure 1: 1D front case and notations

  190A first approximate solution is constructed by replacing the direct discretization of Eq. 22 by a prescribed temperature profile Ts (x, t) that satisfies the boundary conditions Eqs. 23 and 24 and preserve the average temperature Ts (t) Ts (x, t)dx = Ts (t)

where

  St s = C p,s T fus. -Ts (t) ∆H fus. . N.B. similar relations can be obtained in the opposite case when it is assumed that a solidification front monotonously propagates with dm l dt ≤ 0 and d Tl dt ≤ 0. These relations give bounds corresponding to the two extreme cases of power partitioning at the liquid-solid interface between the front propagation and the solid heating:

φls l Sts+ φbc s 1+Sts while

 1+Sts Ts < T fus.φbc s when Ts = T fus.

  , τ ).Then, the dimensionless time is taken as the heat conduction characteristic time in such a way that F o = 1 and the nondimensionalized equations only depends on Bi and St numbers.

Figure 2 :

 2 Figure 2: Geometrical configuration for the different cases and associated notations

For

  the sake of clarity, further discussion of the results in terms of model comparison is splitted into two parts: a comparative analysis between γ Bi = 0.5 and γ Bi = 5.0 cases in Section 4.3.1, γ St = 0.1 and γ St = 5.0 cases in Section 4.3.2. Finally, Section 4.3.3 provides additional results obtained on these transients in order to clearly highlight the shortcomings of the different approximate models.

Figure 3 :

 3 Figure 3: Γ ls φ (t)γ Bi = 0.5 and γ Bi = 5.0 cases -R(10 -4 ) model

Finally

  , the Q model results are slightly deteriorated when St increases; the behavior of the Q model in the γ St = 5.0 case is further discussed in Section 4.3.3. ∆t f (s) εL 1

Figure 4 :

 4 Figure 4: Γ ls φ (t)γ St = 0.1 and γ St = 5.0 cases -R(10 -4 ) model.

Figure 5 :

 5 Figure 5: Ts (t) and Tl (t)γ St = 0.1 -R(10 -4 ) and P(0.1) models.

Figure 6

 6 Figure 6: φls l (t) and φls s (t)γ St = 5.0 -R(10 -4 ) and Q models.

Figure 7 :

 7 Figure 7: e l (t) and e s (t)γ St = 5.0 -R(10 -4 ) and Q models.

Figure 8 :

 8 Figure 8: Tl (t) and Ts (t)γ Bi = 0.5 -R(10 -4 ) and L models.

Figure 9 :

 9 Figure 9: Tl (t) and Ts (t)γ St = 5.0 -R(10 -4 ) and P(0.1) models.

  In this paper, the integral modelling of a two-phase Stefan problem with convection has been discussed from the point of view of solid boundary heat flux closure relations. The comparison of different approximate models has been carried out with respect to a "reference" model that combines the same integral formulation of the conservation equations as the approximate ones with a mesh-based solution of the 1D heat equation in the solid. Numerical results were discussed for fusion transients based on a typical configuration of interest for light water reactor severe accident analysis. A short parametric study based on nondimensionalization considerations in terms of Biot and Stefan numbers was presented in order to highlight the limitations of the different approximations.On the overall, the quadratic profile based model (equivalent to two-point Hermite approximations for the integrals that define the average temperatures and boundary heat fluxes) appears as a good approximation while the applicability of the parametric model and the steady state conduction assumption is shown

  More precisely, H 1,1 (T s ) and H 0,0 ∂T s ∂x (trapezoidal-rule) interpolation are used where, for any function f defined over [e l , x s ], H ν,µ (f ) denotes the Hermite interpolation based on values d i f dx

							λ s	6 Ts (t) -2T s (x s , t) -4T fus. e s	(38)
		φbc s (t) = λ s	b e s = λ s	6 Ts (t) -2T fus. -4T s (x s , t) e s	(39)
	Obviously if Ts (t) =	T fus. + T s (x s , t) 2	, a linear temperature profile is obtained
	i.e. a = 0 in such a way that the steady state of the system without an internal
	source is correctly represented.
	x s e l	∂T s ∂x	dx = T s (x s , t) -T fus. ≈	e s 2	∂T s ∂x	(e l , t) +	∂T s ∂x	(x s , t)

Actually, this approximation is equivalent to the one used in

[START_REF] An | Lumped parameter model for one-dimensional melting in a slab with volumetric heat generation[END_REF]

and based on two-point Hermite-type approximations

[START_REF] Phillips | Explicit forms for certain hermite approximations[END_REF]

for T s (x, t) and

∂T s ∂x (x, t) integrals. i (e l ) 0≤i≤ν and d i f dx i (x s ) 0≤i≤µ

. In this way, the integral approximations are given by:

Ts (t) ≈ 1 2 T fus. + T s (x s , t) + e s

12

∂T s ∂x (e l , t) -∂T s ∂x (x s , t) (40)

Table 1 :

 1 Test cases description

	case	Bi	St	e s (0) (m)	t f (s) †
	nominal	4.115	2.105	0.158	1371
	γ Bi = 0.5	2.057	4.209	0.267	1721
	γ Bi = 5.0 20.573	0.421	0.037	907
	γ St = 0.1	4.115	0.201	0.158	6053
	γ St = 5.0	4.115	10.523	0.158	974
	† as calculated with R(10 -4 ) model.

Table 2 :

 2 Notations for the different models considered in the numerical tests.

-4 ) model. notation model R(∆x) reference model with target mesh size ∆x ∈ 10 -4 , 10 -3 (in m) Q quadratic profile-based approximate model P(ω) parametric approximate model with parameter ω ∈ {0, . . . 1} L steady-state linear profile-based approximate model

Table 3 :

 3 Relative differences -nominal case -reference is R(10 -4 ).

Table 4 :

 4 Relative differencesγ Bi = 0.5 case -reference is R(10 -4 ).

		∆t f (s)	φls s	φbc s	εL 1 (%) Ts	φls l	es
	R(10 -3 )	-5.0	1.01	2.11	0.77	0.20	0.18
	Q	0.1	4.41	2.82	0.64	0.04	0.77
	P (0.0)	-163.1	48.32	86.55	38.65	6.52	18.96
	P (0.1)	-76.4	31.82	57.72	17.85	3.28	13.18
	P (1.0)	651.9	193.88	416.76	102.70	15.86	86.18
	L	-163.1	52.75	93.27	151.53	6.66	22.66

Table 5 :

 5 Relative differencesγ Bi = 5.0 case -reference is R(10 -4 ).

	1
	0.9
	0.8
	0.7
	0.6
	0.5
	0.4
	0.3
	0.2
	0.1
	0
	0	200	400	600	800	1 000	1 200	1 400	1 600	1 800
					time (s)

Table 6 :

 6 Relative differencesγ St = 0.1 case -reference is R(10 -4 ).

		∆t f (s)	φls s	φbc s	εL 1 (%) Ts	φls l	es
	R(10 -3 )	0.0	0.07	0.64	0.15	0.04	0.30
	Q	-5.1	2.99	15.99	3.55	1.23	5.59
	P (0.0)	-81.6	12.30	87.56	22.48	5.12	15.95
	P (0.5)	5.1	3.90	35.69	8.26	2.15	13.26
	P (1.0)	127.4	34.08	236.54	35.80	14.83	74.45
	L	-749.6	92.75	159.52	146.72	29.77	80.11

Table 7 :

 7 Relative differencesγ St = 5.0 case -reference is R(10 -4 ).

	1
	0.9
	0.8
	0.7
	0.6
	0.5
	0.4
	0.3
	0.2
	0.1
	0
	0	500	1 000	1 500	2 000	2 500	3 000
				time (s)

Acknowledgments

This work has been carried out within the framework of the PROCOR platform development funded by CEA, EDF and AREVA.

to be limited to certain fusion transients depending on the power partitioning at the liquid-solid interface between the front propagation and the solid heating.

Appendix A. ODE system numerical integration

The ODE systems are solved using the Apache Commons Math library [START_REF]Apache commons math[END_REF] in the JAVA language. A set of additional top-level classes have been added in the PROCOR platform (kernel package) in order to facilitate the description of such ODE systems in the context of lumped mass and energy conservation equations. In particular, based on the discrete event handling tools of this library (see Apache Commons Math ode package), a model can explicitly declare the different states of its underlying ODEs in such a way that the detection and treatment of the associated transitions can be taken care of properly during the integration process. The state transitions are described through "switching" functions (see Section 6.3 in [START_REF] Gupta | A review of recent developments in solving ODEs[END_REF] for instance) and, through a root-finding algorithm, these functions are used within the time integration process iin order to refine the time step (independently of the time integration scheme) in order to adequately capture the state transitions. Using this library and this formalism, many different integration methods can be used in a transparent way (in this study, an explicit Euler scheme was used). In addition to the three states "liquid only", "liquid-solid" and "solid only", two additional intermediate states have been introduced as transitional states in order to stabilize the front propagation and avoid numerical problems: when a phase appears (surface temperature reaches T fus ), in these intermediate states, the appearing phase is supposed to remain constant at T fus in such a way that only the mass conservation equation is integrated for this appearing phase. T s (0, t) = T fus.

T s (e l (t), t) = T l (e l (t), t) = T fus.

T s (e l (t), t) = T l (e l (t), t) = T fus. m a (resp. m d ) mass threshold associated to a phase appearance (resp. disappearance); P up out power exchanged by the system at its upper boundary;

P dwn in power received by the system at its lower boundary.